Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Parameter estimation

Subject: Parameter estimation

From: Yo

Date: 14 Mar, 2012 19:35:28

Message: 1 of 6

Hi there.

Well, I'm trying to use fminsearch to find the minimum of an objective function.
The problem here is that the objective function is a set of differential equation that have to be solved numerically.... and some experimental values of them.

Let's say I have the following nonlinear system (and represents some physical phenomena)

dx1/dt = a*exp( b/x2) -x1;
dx2/dt = b*sin(x1*x2);

and some experimetal dada stored in the vectors xep1 and xexp2

I want to find the minimum of the square of the residuals

R = (xep-x)'(xep-x)

where
xexp = [xexp1 ; xexp2];
x = [x1 ; x2];

Any clue how to do this???

Subject: Parameter estimation

From: Roger Stafford

Date: 15 Mar, 2012 02:54:33

Message: 2 of 6

"Yo" wrote in message <jjqru0$apl$1@newscl01ah.mathworks.com>...
> Well, I'm trying to use fminsearch to find the minimum of an objective function.
> The problem here is that the objective function is a set of differential equation that have to be solved numerically.... and some experimental values of them.
> .......
- - - - - - - - - -
  If, as you seem to indicate, your question is to be taken with respect to an arbitrary set of differential equations, I seriously doubt if there exists any shortcut to finding the best parameters. You would have to write an ambitious objective function which is able to receive a set of parameters, then use them to call on an ode solver to evaluate your set of differential equation with those parameters from beginning to end, and finally return with a computed square residual as its objective value to be minimized. This objective function would be called by 'fminsearch', or whatever minimization function you select, as many times as it takes to find its way to a valid minimum. I hope you have a lot of time to devote to such a project. I would expect it to be a very time-consuming process.

  You could probably help things along time-wise if the accuracy is initially set in the ode solver and the optimizing function at a low accuracy level to get rough solutions to begin with. Then later you could use these results on a second run to make fairly close initial estimates for a second run with better accuracy. I very much doubt if the functions in the Optimization Toolbox were designed with a possible combined operation with an ode solver in mind.

Roger Stafford

Subject: Parameter estimation

From: Marc

Date: 17 Mar, 2012 04:40:27

Message: 3 of 6

"Yo" wrote in message <jjqru0$apl$1@newscl01ah.mathworks.com>...
> Hi there.
>
> Well, I'm trying to use fminsearch to find the minimum of an objective function.
> The problem here is that the objective function is a set of differential equation that have to be solved numerically.... and some experimental values of them.
>
> Let's say I have the following nonlinear system (and represents some physical phenomena)
>
> dx1/dt = a*exp( b/x2) -x1;
> dx2/dt = b*sin(x1*x2);
>
> and some experimetal dada stored in the vectors xep1 and xexp2
>
> I want to find the minimum of the square of the residuals
>
> R = (xep-x)'(xep-x)
>
> where
> xexp = [xexp1 ; xexp2];
> x = [x1 ; x2];
>

I do this is on a regular basis only I am fitting 20-300 parameters.

My odes are really stiff, so I use ode15s. I am also measuring the input and output only, so I am not mapping ode15s results if that makes sense.

rTime = [0 RTIME];
[rt,Ufinal] = ode15s(@KM, rTime, U);

Hence, I measure U (initial) and Ufinal(end,:)

So my actual function looks something like this

function [yFit] = KMSmall2a(xIn,Parm)

temp =xIn(1);
pressure = xIn(2);...
U = xIn(5:21);
% bunch of code
a1 = Parm(1);
a2 = Parm(2);
.....
aX = Parm(x);
% bunch more code
rTime = [0 RTIME];
[rt,Ufinal] = ode15s(@KM, rTime, U);

yOut = sum(Ufinal(end, 1:17)).*Ufinal(end,1:17)./100;

yFit = zeros(15,1);

yFit(1) = yOut(1);
yFit(2) = yOut(2);
yFit(3) = yOut(3);
yFit(4) = yOut(4);
yFit(5) = yOut(5); ....

yFit(15) = yOut(15);
%% here is my function with ODEs
function Rate = KM(rT,U)
% some code
% some code
end
end

So I use a nested function to create a function that takes my measured x's and my parameters to spit on calculated y's.

Now, fminsearch or the optimization toolbox fmincon, etc. etc. want a function and an initial guess at your parameters.

Here is where it becomes an art.

I set up another function to accept a guess at the parameter set to then spit out a calculated f (funcDiff) here. Within this function is where I put my x and y data that I measured. In another file I create the data in .mat files which then this function uses to "bring in" when called.
 
I use something like this for lsqnonlin
%----------------------------------------------------------------------
function funcDiff = ModelFitFunc(Parm)

%% Load data
 
Data = load 'xxxxx.mat';

m = size(Data);
funcDiff = zeros(m(1),1);

for i = 1:m(1)
    xIn = lData(i,1:21)';
    yMeas = Data(i,22:36)';

    yCalc = KMSmall2a(xIn,Parm);

    funcDiff(i) = sum((yCalc - yMeas).^2);
end

end

For fminsearch, it requires you to difine the function output (ie. a single value not a vector).

So you can see I create my funcOut here for fminsearch to minimize..

I use this for fminsearch and some other solvers
%---------------------------------------------------------------------------
function funcOut = ModelFitFuncSingle(Parm)
%% Load data
 
Data = load 'xxxxx.mat';

m = size(Data);
funcDiff = zeros(m(1),1);

for i = 1:m(1)
    xIn = lData(i,1:21)';
    yMeas = Data(i,22:36)';

    yCalc = KMSmall2a(xIn,Parm);

    funcDiff(i) = sum((yCalc - yMeas).^2);
end

funcOut = sum(funcDiff.^2);
end

Here is a how some of my code looks to call the above... I left in some hints showing how easy it can be to then use other functions from the optimization toolbox....

Parm = ParmOpt;
%options = optimset('Algorithm','levenberg-marquardt');
tic

%[ParmOpt,resnorm] = lsqnonlin(@ModelFitFunc,Parm,[],[], options);
%ParmOpt = fminunc(@ModelFitFuncSingle,Parm)
ParmOpt = fminsearch(@ModelFitFuncSingle,Parm)

fOutOld = ModelFitFunc(Parm)
fOut = ModelFitFunc(ParmOpt)
toc

What I am not showing you is how I optimize my function by weighting y values. Hopefully you can figure that out.

For ode equations, lsqnonlin and fminunc are faster but I run into issues with the numerical jacobians being singular and ending up with poorer fits at times and early exit from the program.

fminsearch uses a nelder mead simplex algorithm and takes anywhere from minutes to hours.

You should be able to estimate how long this will take by tic toc your function out.

If it takes roughly 0.01 sec to run my ode function and I have 30 data points, it takes 0.3-0.4s for ONE call from fminsearch to my function. So, patience is a virtue.

Hope this helps you think about how best to set up your problem.

Subject: Parameter estimation

From: John Hedengren

Date: 20 Mar, 2012 05:00:18

Message: 4 of 6

Another option for solving parameter estimation is to simultaneously solve the objective function and model equations. Some benchmark testing shows that the simultaneous approach can solve much larger systems than the "shooting method" as used with ode15s or other forward-stepping integrators.

Below is some additional information on using the APM MATLAB package for solving large-scale differential and algebraic equation systems in parameter estimation or nonlinear control.

http://www.mathworks.com/matlabcentral/fileexchange/32068-optimization-nonlinear-control-and-estimation-toolbox

There is one example in this package of a pendulum motion that is used to estimate the mass and length of the pendulum - it may be a good starting point for the development of your application. We have bi-weekly webinars that give tutorials on getting starting:

http://apmonitor.com/wiki/index.php/Main/ApplicationWebinars

-John

"Yo" wrote in message <jjqru0$apl$1@newscl01ah.mathworks.com>...
> Hi there.
>
> Well, I'm trying to use fminsearch to find the minimum of an objective function.
> The problem here is that the objective function is a set of differential equation that have to be solved numerically.... and some experimental values of them.
>
> Let's say I have the following nonlinear system (and represents some physical phenomena)
>
> dx1/dt = a*exp( b/x2) -x1;
> dx2/dt = b*sin(x1*x2);
>
> and some experimetal dada stored in the vectors xep1 and xexp2
>
> I want to find the minimum of the square of the residuals
>
> R = (xep-x)'(xep-x)
>
> where
> xexp = [xexp1 ; xexp2];
> x = [x1 ; x2];
>
> Any clue how to do this???

Subject: Parameter estimation

From: Yo

Date: 30 Mar, 2012 22:00:28

Message: 5 of 6

Thanks, now I've tried different approaches and is not that hard. The key is to define an error function to minimize and that's all.

Thank you all!

Subject: Parameter estimation

From: Marc

Date: 5 Apr, 2012 05:22:25

Message: 6 of 6

"John Hedengren" <john@hedengren.net> wrote in message <jk92t2$bqn$1@newscl01ah.mathworks.com>...
> Another option for solving parameter estimation is to simultaneously solve the objective function and model equations. Some benchmark testing shows that the simultaneous approach can solve much larger systems than the "shooting method" as used with ode15s or other forward-stepping integrators.
>
> Below is some additional information on using the APM MATLAB package for solving large-scale differential and algebraic equation systems in parameter estimation or nonlinear control.
>
> http://www.mathworks.com/matlabcentral/fileexchange/32068-optimization-nonlinear-control-and-estimation-toolbox
>
> There is one example in this package of a pendulum motion that is used to estimate the mass and length of the pendulum - it may be a good starting point for the development of your application. We have bi-weekly webinars that give tutorials on getting starting:
>
> http://apmonitor.com/wiki/index.php/Main/ApplicationWebinars
>
> -John
>
> "Yo" wrote in message <jjqru0$apl$1@newscl01ah.mathworks.com>...
> > Hi there.
> >
> > Well, I'm trying to use fminsearch to find the minimum of an objective function.
> > The problem here is that the objective function is a set of differential equation that have to be solved numerically.... and some experimental values of them.
> >
> > Let's say I have the following nonlinear system (and represents some physical phenomena)
> >
> > dx1/dt = a*exp( b/x2) -x1;
> > dx2/dt = b*sin(x1*x2);
> >
> > and some experimetal dada stored in the vectors xep1 and xexp2
> >
> > I want to find the minimum of the square of the residuals
> >
> > R = (xep-x)'(xep-x)
> >
> > where
> > xexp = [xexp1 ; xexp2];
> > x = [x1 ; x2];
> >
> > Any clue how to do this???

John

I would be really interested in how you solve both the model equations and the objective function simultaneously.

I pretty much showed you mine, now you show me yours.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us