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

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

%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)
%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

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

[IPSHITA]

there is an error in the text on the line global t c1meas c2meas c3meas c4meas c5meas

as it is showing that it is not inside the function

please solve the problem

[Torsten]

The code has errors, but the line you refer to shouldn't cause problems.

It's the script part of the code, so the line doesn't need to be inside a function.

I rearranged the order in which script and functions appear in the code above to make it "work".

Answer 1

[rickert]

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

Answer 2

[Star Strider]

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.