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:
lsqnonlin and ode

Subject: lsqnonlin and ode

From: Franz

Date: 4 Aug, 2011 14:47:23

Message: 1 of 8

Hi,

i have a system of 3 differential equations, which are solved by an ode45 solver.
Additionally I want to fit 9 parameters of these equations to experimental data, by using the "levenberg-marquardt"-algorithm of lsqnonlin. The experimental data is imported via xlsread. An error prompt isn't appearing, if I start the code.
Solving the equations and fitting the parameters lasts very long.

My question is: Did I make a mistake in programming or has onyone a proposal for improving the performance of the code?

This is my code:
FILE lsqparafit1.m

function lsqparafit1

clear global, clear all, clear clc
global c_v0 HPLC tspan

%read Excel file
HPLC=xlsread('Messwerte.xls','gewaschen','A82:D96');

c_v0(1) = HPLC (1,4); % Cellulose
c_v0(2) = HPLC (1,3); % Cellobiose
c_v0(3) = HPLC (1,2); % Glucose

%Initial Parameter values
KM=22;
kR1=20;
kR2=5;
kR3=265;
K1Ig=0.07;
K2Ig=0.03;
K3Ig=4;
K1Icb=0.02;
K2Icb=130.0;

par0 = [KM,kR1,kR2,kR3,K1Ig,K2Ig,K3Ig,K1Icb,K2Icb];

%Interval for the ode solver
tspan = HPLC(:,1)'; % Matrix transponieren

%choosing levenberg-marquardt
options = optimset('Algorithm','levenberg-marquardt');
%Supplying the basis values par0 to lsqnonlin
paropt = lsqnonlin(@minimiz,par0,[],[],options);
%solving ode with optimized parameters
[t,c_v] = ode45(@lsqparafit2,tspan,c_v0,' ',paropt);

function lsq=minimiz(par)
global c_v0 HPLC tspan

%solving ode with initial values par0
[t,c_v] = ode45(@lsqparafit2,tspan,c_v0,' ',par);

%addition residues
lsq = sum(abs(c_v(:,1)-HPLC(:,4)))+sum(abs(c_v(:,2)-HPLC(:,3)))+sum(abs(c_v(:,3)-HPLC(:,2)));

FILE lsqparafit2.m

function dc_v=lsqparafit2(t,c_v,par)

% Substratgehalt (Cellulose)
CF0=0.6;
SL0=100;
Cs0=CF0*SL0;
.........

%Parameters
KM= par(1);
kR1= par(2);
kR2= par(3);
kR3= par(4);
K1Ig= par(5);
K2Ig= par(6);
K3Ig= par(7);
K1Icb= par(8);%0.015;
K2Icb= par(9);%132.0;

.......

r1=(kR1*Arf*Ceb1*RS*c_v(1))/(1+(c_v(2)/K1Icb)+(c_v(3)/K1Ig)+(Cx/K1Ix));
r2=(kR2*Arf*Ceb1*RS*c_v(1))/(1+(c_v(2)/K2Icb)+(c_v(3)/K2Ig)+(Cx/K2Ix));
r3=(kR3*Arf*Cef2*c_v(2))/(KM*(1+(c_v(3)/K3Ig)+(Cx/K3Ix))+c_v(2));


dCs_dt=-r1-r2;
dCcb_dt=1.056*r1-r3;
dCg_dt=1.111*r2+1.053*r3;

dc_v=[dCs_dt;dCcb_dt;dCg_dt];

%%%%%%%%%%%%%%%%%%%%%%%
END of code


It would be nice if somebody could help me!

Thankyou:)

ciao

Subject: lsqnonlin and ode

From: Alan Weiss

Date: 4 Aug, 2011 14:58:28

Message: 2 of 8

On 8/4/2011 10:47 AM, Franz wrote:
> %addition residues
> lsq =
> sum(abs(c_v(:,1)-HPLC(:,4)))+sum(abs(c_v(:,2)-HPLC(:,3)))+sum(abs(c_v(:,3)-HPLC(:,2)));
>

This is a mistake in your code. Look at the documentation:
http://www.mathworks.com/help/toolbox/optim/ug/lsqnonlin.html
In the Description, it says "Rather than compute the value ||f(x)||^2,
lsqnonlin requires the user-defined function to compute the
vector-valued function f(x) = [f_1(x) f_2(x) ... f_n(x)].

By passing a scalar instead of a vector, you slow the solution tremendously.

You might also want to turn on iterative display (set the Display option
to 'iter').

Alan Weiss
MATLAB mathematical toolbox documentation

Subject: lsqnonlin and ode

From: Alan Weiss

Date: 4 Aug, 2011 14:59:57

Message: 3 of 8

On 8/4/2011 10:58 AM, Alan Weiss wrote:
> On 8/4/2011 10:47 AM, Franz wrote:
>> %addition residues
>> lsq =
>> sum(abs(c_v(:,1)-HPLC(:,4)))+sum(abs(c_v(:,2)-HPLC(:,3)))+sum(abs(c_v(:,3)-HPLC(:,2)));
>>
>>
>
> This is a mistake in your code. Look at the documentation:
> http://www.mathworks.com/help/toolbox/optim/ug/lsqnonlin.html
> In the Description, it says "Rather than compute the value ||f(x)||^2,
> lsqnonlin requires the user-defined function to compute the
> vector-valued function f(x) = [f_1(x) f_2(x) ... f_n(x)].
>
> By passing a scalar instead of a vector, you slow the solution
> tremendously.
>
> You might also want to turn on iterative display (set the Display option
> to 'iter').
>
> Alan Weiss
> MATLAB mathematical toolbox documentation

Oh, and I forgot to say that you should probably set the DiffMinChange
option to a larger value than the default; try 1e-4.

Alan Weiss
MATLAB mathematical toolbox documentation

Subject: lsqnonlin and ode

From: Franz

Date: 5 Aug, 2011 14:58:10

Message: 4 of 8

Thank you very much for your help.
I have already changed the options of optimset and odeset.
But I still don't know (even after reading the manual of lsqnonlin) how I have to pass the parameters to lsqnonlin.
I have two matrices, one with measured data , and one with data from the ode45 -solver. And I want to minimize the difference of these matrices, by passing the difference to lsqnonlin and thereby fit the initial parameters.
I´m afraid I just began to use Matlab, but this software alreaady helped me solving several differential equations and it would be great if I could also solve this parameter fitting problem.
Perhaps you have some kind of example or a tip for me again.

Thank you!

Franz




Alan Weiss <aweiss@mathworks.com> wrote in message <j1ec5d$qgo$2@newscl01ah.mathworks.com>...
> On 8/4/2011 10:58 AM, Alan Weiss wrote:
> > On 8/4/2011 10:47 AM, Franz wrote:
> >> %addition residues
> >> lsq =
> >> sum(abs(c_v(:,1)-HPLC(:,4)))+sum(abs(c_v(:,2)-HPLC(:,3)))+sum(abs(c_v(:,3)-HPLC(:,2)));
> >>
> >>
> >
> > This is a mistake in your code. Look at the documentation:
> > http://www.mathworks.com/help/toolbox/optim/ug/lsqnonlin.html
> > In the Description, it says "Rather than compute the value ||f(x)||^2,
> > lsqnonlin requires the user-defined function to compute the
> > vector-valued function f(x) = [f_1(x) f_2(x) ... f_n(x)].
> >
> > By passing a scalar instead of a vector, you slow the solution
> > tremendously.
> >
> > You might also want to turn on iterative display (set the Display option
> > to 'iter').
> >
> > Alan Weiss
> > MATLAB mathematical toolbox documentation
>
> Oh, and I forgot to say that you should probably set the DiffMinChange
> option to a larger value than the default; try 1e-4.
>
> Alan Weiss
> MATLAB mathematical toolbox documentation

Subject: lsqnonlin and ode

From: Alan Weiss

Date: 5 Aug, 2011 17:42:40

Message: 5 of 8

On 8/5/2011 10:58 AM, Franz wrote:
> Thank you very much for your help.
> I have already changed the options of optimset and odeset.
> But I still don't know (even after reading the manual of lsqnonlin) how
> I have to pass the parameters to lsqnonlin.
> I have two matrices, one with measured data , and one with data from the
> ode45 -solver. And I want to minimize the difference of these matrices,
> by passing the difference to lsqnonlin and thereby fit the initial
> parameters.
> I´m afraid I just began to use Matlab, but this software alreaady helped
> me solving several differential equations and it would be great if I
> could also solve this parameter fitting problem.
> Perhaps you have some kind of example or a tip for me again.
>
> Thank you!
>
> Franz

The objective function for lsqnonlin is (measured - ode45solution).
lsqnonlin will attempt to minimize the sum of squares of these differences.

To pass parameters from one function to another (i.e., between lsqnonlin
and ode45), see
http://www.mathworks.com/help/toolbox/optim/ug/brhkghv-7.html

Good luck,

Alan Weiss
MATLAB mathematical toolbox documentation

Subject: lsqnonlin and ode (2)

From: Frans

Date: 8 Dec, 2011 19:50:08

Message: 6 of 8

My problem is similar: i'm struggling with passing on variables through functions.
Somehow my script with my ODE's does not sent the result back to the other sccript.
Working with nargin and nargout did not help.

Adjusting the previous post with al my constants&symbols yielded a very long calculation time, which is why I post my own code.

Case:
I have three ODE's with two parameters. I want to fit the result of the last ODE to my measurement. I tried to do it in three files:

file 1: Fitterfrans.m

clear global
close all
clear all
clc

%Initial Parameter values
V10=100*10^-6;
V20= 60*10^-6;

par0 = [V10 V20];

%choosing levenberg-marquardt
options = optimset('DiffMinChange',1e-4,'display','iter');

%Supplying the basis values par0 to lsqnonlin
paropt = lsqnonlin(@tanksinseriesaan,par0,[],[],options)

end of scipt 1
File 2: tanksinseriesaan.m

function rss = tanksinseriesaan(par,v,V1,V2,V3)

% Load measurement data
    load('D:\My Documents\Artikel 2 RTD\work up\3 mm - 06-Dec-2011.mat')

    for n=5
            a=Et12(:,n);
            a(a == 0)=[NaN];
            Eexp=a; % this gives the y data that we want to fit
            b=theta(:,n);
            b(b == 0)=[NaN];
            tauexp=b; % this gives the x data that we want to fit
    end

 % Load constants
     v=15*10^-6; % Volumetric flow rate
     Vtot=167.0176*10^-6; % Reactor volume at 1 mm discspacing

% ODE
    %Interval for the ode solver
    tspan = tauexp(:,1); % equal to tspan of the fitted measurement
   
    %define parameter to vary
    V1=par(1); % volume 1
    V2=par(2); % volume 1
    V3=Vtot-V1-V2; % total = volume 1 + 2 +3
    
    % ODE
    [Cmod] = ode45(@tanksinseries,tspan,[0 0 0],v,V1,V2,V3);

        % convert Concentration of tank 3 into E(t) diagram

            taavg=trapz(t,(1-Cmod)); % Residence time from integral of 1-F!
            taavg=taavg(:,3); % Residence time after third tank
            Cdiff=diff(Cmod(:,3)); % delta y
            tdiff=diff(t(:,1)); % delta x
            diffinal=Cdiff./tdiff; % dydx = afgeleide = E(t)
            Edim=taavg.*diffinal;
       
rss=Edim(:,1)-Eexp(:,1);

end of scipt 2
file 3: tanksinseries.m

function diff=tanksinseries(v,V1,V2,V3)

diff = zeros(3,1);

end of script 3

% the three ODE's:
   diff(1)=v/V1*1-v/V1*C(1);
   diff(2)=v/V2*C(1)-v/V2*C(2);
   diff(3)=v/V3*C(2)-v/V3*C(3);

%
By applying the pause command I can see that eveything goes well, untill

Subject: lsqnonlin and ode (2.1)

From: Frans

Date: 8 Dec, 2011 20:14:09

Message: 7 of 8

Apologize for my previous post.
I have the same problem as before: combining ODE with LSQNONLIN results in long waiting time.
Opposite as posted befor my script is working now, but it is slow.

Script 1: naame fitterfrans.m

clear global
close all
clear all
clc

%Initial Parameter values
V10=100*10^-6;
V20= 60*10^-6;
par0 = [V10 V20];

%choosing levenberg-marquardt
options = optimset('DiffMinChange',1e-4);

%Supplying the basis values par0 to lsqnonlin
paropt = lsqnonlin(@tanksinseriesaan,par0,[],[],options)

end of script 1.
script 2. name: tanksinseriesaan.m

function rss = tanksinseriesaan(par)

% Load measurement data
    load('D:\My Documents\Artikel 2 RTD\work up\3 mm - 06-Dec-2011.mat')

    for n=5
            a=Et12(:,n);
            a(a == 0)=[NaN];
            Eexp=a;
            b=theta(:,n);
            b(b == 0)=[NaN];
            tauexp=b;
    end

 % Load constants
     v=15*10^-6; % Volumetric flow rate
     Vtot=167.0176*10^-6; % Reactor volume at 1 mm discspacing

% ODE
    %Interval for the ode solver
    tspan = tauexp(:,1); % equal to tspan of measurement
    options =odeset('RelTol',1e-5,'AbsTol',1e-5);
   
    %define parameter to vary
    V1=par(1); % volume 1
    V2=par(2); % volume 1
    V3=Vtot-V1-V2; % total = volume 1 + 2 +3
    
    C0=[0 0 0];
    % ODE
    [t,diff] = ode45(@tanksinseries,tspan,C0,options,v,V1,V2,V3);
    
        % convert Concentration of tank 3 into E(t) diagram
            taavg=trapz(t,(1-Cmod)); % Residence time from integral of 1-F!
            taavg=taavg(:,3); % Residence time after third tank
            Cdiff=diff(Cmod(:,3)); % delta y
            tdiff=diff(t(:,1)); % delta x
            diffinal=Cdiff./tdiff; % dydx = afgeleide = E(t)
            Edim=taavg.*diffinal;
       
rss=[Edim(:,1)-Eexp(:,1)];

end of script 2.
script 3: name tanksinseries.m

function diff=tanksinseries(t,diff,v,V1,V2,V3)
diff = zeros(3,1);
diff(1)=v/V1*1-v/V1*diff(1);
diff(2)=v/V2*diff(1)-v/V2*diff(2);
diff(3)=v/V3*diff(2)-v/V3*diff(3);

Subject: lsqnonlin and ode (2.1)

From: Torsten

Date: 9 Dec, 2011 07:36:08

Message: 8 of 8

On 8 Dez., 21:14, "Frans " <f.vissc...@tue.nl> wrote:
> Apologize for my previous post.
> I have the same problem as before: combining ODE with LSQNONLIN results in long waiting time.
> Opposite as posted befor my script is working now, but it is slow.
>
> Script 1: naame fitterfrans.m
>
> clear global
> close all
> clear all
> clc
>
> %Initial Parameter values
> V10=100*10^-6;
> V20= 60*10^-6;
> par0 = [V10 V20];
>
> %choosing levenberg-marquardt
> options = optimset('DiffMinChange',1e-4);
>
> %Supplying the basis values par0 to lsqnonlin
> paropt = lsqnonlin(@tanksinseriesaan,par0,[],[],options)
>
> end of script 1.
> script 2. name: tanksinseriesaan.m
>
> function rss = tanksinseriesaan(par)
>
> %   Load measurement data
>     load('D:\My Documents\Artikel 2 RTD\work up\3  mm - 06-Dec-2011.mat')
>
>     for n=5
>             a=Et12(:,n);
>             a(a == 0)=[NaN];
>             Eexp=a;
>             b=theta(:,n);
>             b(b == 0)=[NaN];
>             tauexp=b;
>     end
>
>  % Load constants
>      v=15*10^-6;             % Volumetric flow rate
>      Vtot=167.0176*10^-6;    % Reactor volume at 1 mm discspacing
>
> % ODE
>     %Interval for the ode solver
>     tspan = tauexp(:,1);                            % equal to tspan of measurement
>     options =odeset('RelTol',1e-5,'AbsTol',1e-5);
>
>     %define parameter to vary
>     V1=par(1);      % volume 1
>     V2=par(2);      % volume 1
>     V3=Vtot-V1-V2;  % total = volume 1 + 2 +3
>
>     C0=[0 0 0];
>     % ODE
>     [t,diff] = ode45(@tanksinseries,tspan,C0,options,v,V1,V2,V3);
>
>         % convert Concentration of tank 3 into E(t) diagram
>             taavg=trapz(t,(1-Cmod)); % Residence time from integral of 1-F!
>             taavg=taavg(:,3);        % Residence time after third tank
>             Cdiff=diff(Cmod(:,3));   % delta y
>             tdiff=diff(t(:,1));      % delta x
>             diffinal=Cdiff./tdiff;   % dydx = afgeleide = E(t)
>             Edim=taavg.*diffinal;
>
> rss=[Edim(:,1)-Eexp(:,1)];
>
> end of script 2.
> script 3: name tanksinseries.m
>
> function diff=tanksinseries(t,diff,v,V1,V2,V3)
> diff = zeros(3,1);
> diff(1)=v/V1*1-v/V1*diff(1);
> diff(2)=v/V2*diff(1)-v/V2*diff(2);
> diff(3)=v/V3*diff(2)-v/V3*diff(3);

Two remarks:
1. Load the measurement data once - not every time lsqnonlin calls
'tanksinseriesaan'.
2. I guess your system of differential equations is set up
incorrectly ; I think it should read
 > function diff=tanksinseries(t,y,v,V1,V2,V3)
> diff = zeros(3,1);
> diff(1)=v/V1*1-v/V1*y(1);
> diff(2)=v/V2*y(1)-v/V2*y(2);
> diff(3)=v/V3*y(2)-v/V3*y(3);
If this is the case, use dsolve before you start your calculation in
order to solve it analytically,
insert the solution 'by hand' in 'tanksinseriesaan' and use it to
calculate the required residuals
for lsqnonlin.

Best wishes
Torsten.

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