Parameter Estimation for a complex phyiscal system described by a set of differential equations

Hello.

I have a matLab program performing a time dependent numerical simulation for a complex phyiscal system, described by a set of differential equations.

let say there is 4 input parameters to be determined, and the output is voltage vesus time.

I want to do parameter estimation(determining their value automatically) for that 4 parameters,

by fitting the simulation result to the experiment result.

I imagine bayesian optimization could do the job, but I have no idea how to implement it to my matLab program.

Is there any product, module or application of matLab can perform this? Or any hints for the implmentation.

Any help would be appreciated. Thank you

[Star Strider]

I am not certain what you want to do. This search may provide some solutions.

[Sam Wong]

Let me make my question more clear.

I have a function doing high cost simulation,

function [timeS, voltageS] = startSim(a,b,c,d)

...

end

and, I have an experiment data in array [timeE, voltageE], the target of fitting.

My target is to do parameter fitting to automatically find the value of (a,b,c,d) that [timeS, voltageS] = [timeE, voltageE] in terms of mean squared error.

Here is the process I imagine

①Do simulation by calling function [timeS, voltageS] = startSim(a,b,c,d)

↓

②interpolate [timeS, voltageS] to the same time value [timeE, voltageE] , in order to calculate the mean squared error.

"Sum(voltageE - voltageS) for all timeE"

↓

③Call bayesopt for Bayesian Optimization to minimize "Sum(voltageE - voltageS) for all time", and obtain the cooresponding parameter (a,b,c,d)

However, I do not know how to implement ③BayesOpt to ①,②.

To learn using BayesOpt, I now simplify my case into fitting sin curves (of two parameters a,b) following the above 3 steps,

but now it says I didnt define the parameter well. The code is as below.

%exp-data

x0 = 0:1:6*pi;

y0 = 2*sin(x0/3);

%simulation parameter to be fitted, expected result would be a=2, b=0.33333

a = optimizableVariable('a',[-1,1],'Type','real');

b = optimizableVariable('b',[0,1],'Type','real');

%Do simulation calling function startSim

x=0:0.5:6*pi;

y = startSim(a,b)

% Bayesian Optimization

obj=@fun;

results = bayesopt(obj,[a,b])

%simulation program of sin curve with parameter a,b

function y = startSim(a,b)

x=0:0.5:6*pi;

y=a*sin(x*b);

end

%Definition of the object function

function obj = fun(x,y,x0,y0)

yi=interp1(x,y,x0,"spline");

obj = sum((yi-y0).^2);

end

[Star Strider]

The bayesopt function wants a single vector of parameters, so ‘fun’ has to present only that vector. I rarely use bayesopt, however the code likely needs to be something like this (given what I understand of it which is not much at this point) —

% Bayesian Optimization
obj = @(y)fun(x,y,x0,y0);
results = bayesopt(obj,[a,b])

since it would appear that the optimisation is with respect to ‘y’. I am certain that you understand where I am going with this, so choose the appropriate parameter to optimise in order to get the result you want.

.

[Sam Wong]

Thank you! with the Tips on https://www.mathworks.com/help/stats/variables-for-a-bayesian-optimization.html ,

now I can fit the sin curves. But in my real application.

My parameters are stored in structure as fields. e.g. param{1}.a

After I change all the a to param{1}.a following the format, the program cant run with error

Error using tabular/throwNDSubscriptError (line 503)
Subscripting into a table using one subscript (as in t(i)) is not supported. Specify a row subscript and a variable
subscript, as in t(rows,vars). To select variables, use t(:,i) or for one variable t.(i). To select rows, use t(i,:).

I think I need to change input from structure to array. Any help? Here is my code.

param{1}= Parameters;

%exp-data
x0 = 0:0.2:6*pi;
y0 = 2*sin(x0/3)+param{1}.c;
x=0:0.5:6*pi;
obj=@(param)fun(param,x,x0,y0);
results = bayesopt(obj,[param{1}.a,param{1}.b],'IsObjectiveDeterministic',true,'MaxObjectiveEvaluations',50);

function obj = fun(param,x,x0,y0)
   y = startSim(param);
   yi=interp1(x,y,x0,"spline");
   obj = sum((yi-y0).^2);
end    
function y = startSim(param)
     x=0:0.5:6*pi;
     y=param{1}.a*sin(x*param{1}.b)+param{1}.c;
end
function param= Parameters(~)
    param.a = optimizableVariable('aa',[-1,3],'Type','real');
    param.b = optimizableVariable('bb',[-1,1],'Type','real');
    param.c = 1.3;
end