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

# 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 wrote in message ... > 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 " 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.