MATLAB Answers

0

parameter estimation for a set of differential equations using fmincon function in matlab

Asked by vinutha paravastu on 18 Aug 2019
Latest activity Answered by Star Strider
on 18 Aug 2019
hello all i am new to matlab pls help in rectifying my errors in matlab code
basically i am developing kinetic model for a reaction and i need to estimate parameters based on the experimental values and i used fmincon function for optimizing parameters and used minimum sum of least squares as objective function
i used the following codes
1.kinetics1 function code
function dcdt = kinetics1(p,t,c)
%Parameters
k1=p(1);
k2=p(2);
k3=p(3);
k4=p(4);
k5=p(5);
k6=p(6);
k7=p(7);
k8=p(8);
k9=p(9);
k10=p(10);
k11=p(11);
k12=p(12);
k13=p(13);
k14=p(14);
% mass balances
dcdt(1) = (k1+k2+k3+k4)*c(1);
dcdt(2) = (k1*c(1))-(k5+k6+k7)*c(2);
dcdt(3) = (k2*c(1))-(k8+k9+k11+k12)*c(3);
dcdt(4) = (k3*c(1)+k6*c(2)+k11*c(3)+k8*c(3))-(k13+k14+k10)*c(4);
dcdt(5)=(k4*c(1)+k7*c(2)+k9*c(3)+k12*c(3)+k13*c(4)+k14*c(4)+k10*c(4));
end
2.simulate function code
function c = simulate(p)
global t c1meas c2meas c3meas c4meas c5meas
c=zeros(length(t),5);
c(1,1)=c1meas(1);
c(1,2)=c2meas(2);
c(1,3)=c3meas(3);
c(1,4)=c4meas(4);
c(1,5)=c5meas(5);
c0=[c(1,1),c(1,2),c(1,3),c(1,4),c(1,5)];
for i=length(t)-1
ts=[t(i),t(i+1)];
[t,c]= ode45(@(t,c)kinetics1(p,t,c),ts,c0);
c(i+1,1)=c0(1);
c(i+1,2)=c0(2);
c(i+1,3)=c0(3);
c(i+1,4)=c0(4);
c(i+1,5)=c0(5);
end
end
3.Objective function code
%define objective
function obj=objective(p)
global c1meas c2meas c3meas c4meas c5meas
%simulate model
c=simulate(p);
obj = sum(((c(:,1)-c1meas)./c1meas).^2)
+ sum(((c(:,2)-c2meas)./c2meas).^2)+sum(((c(:,2)-c3meas)./c3meas).^2)+sum(((c(:,2)-c4meas)./c4meas).^2)+sum(((c(:,2)-c5meas)./c5meas).^2);
end
4. my final script vinutha3
global t c1meas c2meas c3meas c4meas c5meas
t= [0.88 0.96 0.98 1.04 1.05];
c1meas=[0.211 0.066 0.17 0.455 0.088];
c2meas =[0.666 0.165 0.083 0.047 0.009];
c3meas=[0.302 0.093 0.155 0.341 0.094]
c4meas=[0.237 0.084 0.177 0.404 0.082]
c5meas=[0.686 0.16 0.072 0.041 0.008]
%parameters initial guess
k1=1;
k2=1;
k3=1;
k4=1;
k5=1;
k6=1;
k7=1;
k8=1;
k9=1;
k10=1;
k11=1;
k12=1;
k13=1;
k14=1;
p0=[k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14];
% show initial objective
disp(['Initial SSE Objective: ' num2str(objective(p0))])
% optimize parameters
% no linear constraints
A = [];
b = [];
Aeq = [];
beq = [];
nlcon = [];
lb=[];
ub=[];
% options = optimoptions(@fmincon,'Algorithm','interior-point');
p = fmincon(@objective,p0,A,b,Aeq,beq,lb,ub,nlcon); %,options);
% show final objective
disp(['Final SSE Objective: ' num2str(objective(p))])
% optimized parameter values
k1=p(1);
k2=p(2);
k3=p(3);
k4=p(4);
k5=p(5);
k6=p(6);
k7=p(7);
k8=p(8);
k9=p(9);
k10=p(10);
k11=p(11);
k12=p(12);
k13=p(13);
k14=p(14);
disp(['k1: ' num2str(k1)])
disp(['k2: ' num2str(k2)])
disp(['k3: ' num2str(k3)])
disp(['k4: ' num2str(k4)])
disp(['k5: ' num2str(k5)])
disp(['k6: ' num2str(k6)])
disp(['k7: ' num2str(k7)])
disp(['k8: ' num2str(k8)])
disp(['k9: ' num2str(k9)])
disp(['k10: ' num2str(k10)])
disp(['k11: ' num2str(k11)])
disp(['k12: ' num2str(k12)])
disp(['k13: ' num2str(k13)])
disp(['k14: ' num2str(k14)])
% calculate model with updated parameters
ci = simulate(p0);
cp = simulate(p)
i am getting following errors
Error using odearguments (line 93)
@(T,C)KINETICS1(P,T,C) must return a column vector.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in simulate (line 12)
[t,c]= ode45(@(t,c)kinetics1(p,t,c),ts,c0);
Error in objective (line 6)
c=simulate(p);
Error in vinutha3 (line 25)
disp(['Initial SSE Objective: ' num2str(objective(p0))])
plssssssss help meeeeeee

  0 Comments

Sign in to comment.

Tags

No tags entered yet.

2 Answers

Answer by rickert
on 18 Aug 2019

Just look at your error traceback!
Error using odearguments (line 93)
@(T,C)KINETICS1(P,T,C) must return a column vector.
But your kinetics1 function returns dcdt as a row vector. Either add a transpose, or initialize dcdt as e.g. NaNs:
function dcdt = kinetics1(p,t,c)
%Parameters
k1=p(1);
...
% mass balances
dcdt = NaN(5,1);
dcdt(1) = (k1+k2+k3+k4)*c(1);
dcdt(2) = (k1*c(1))-(k5+k6+k7)*c(2);
dcdt(3) = (k2*c(1))-(k8+k9+k11+k12)*c(3);
dcdt(4) = (k3*c(1)+k6*c(2)+k11*c(3)+k8*c(3))-(k13+k14+k10)*c(4);
dcdt(5)=(k4*c(1)+k7*c(2)+k9*c(3)+k12*c(3)+k13*c(4)+k14*c(4)+k10*c(4));
end

  0 Comments

Sign in to comment.


Answer by Star Strider
on 18 Aug 2019

You will likely find Monod kinetics and curve fitting helpful. (There are similar posts as well.)
One way to force ‘dcdt’ to become a column vector is to use a zeros call first:
function dcdt = kinetics1(p,t,c)
%Parameters
k1=p(1);
k2=p(2);
k3=p(3);
k4=p(4);
k5=p(5);
k6=p(6);
k7=p(7);
k8=p(8);
k9=p(9);
k10=p(10);
k11=p(11);
k12=p(12);
k13=p(13);
k14=p(14);
% mass balances
dcdt = zeros(5,1); % Define As Column Vector
dcdt(1) = (k1+k2+k3+k4)*c(1);
dcdt(2) = (k1*c(1))-(k5+k6+k7)*c(2);
dcdt(3) = (k2*c(1))-(k8+k9+k11+k12)*c(3);
dcdt(4) = (k3*c(1)+k6*c(2)+k11*c(3)+k8*c(3))-(k13+k14+k10)*c(4);
dcdt(5)=(k4*c(1)+k7*c(2)+k9*c(3)+k12*c(3)+k13*c(4)+k14*c(4)+k10*c(4));
end
You are doing parameter estimation, so fmincon is likely not the optimal function to use.
Also it is best to avoid global variables. See the documentation section on Passing Extra Parameters for the correct way to do ihat.

  0 Comments

Sign in to comment.