"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 = (xepx)'(xepx)
>
> where
> xexp = [xexp1 ; xexp2];
> x = [x1 ; x2];
>
I do this is on a regular basis only I am fitting 20300 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','levenbergmarquardt');
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.30.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.
