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

15 ビュー (過去 30 日間)
vinutha paravastu
vinutha paravastu 2019 年 8 月 18 日
編集済み: Torsten 2024 年 3 月 1 日 12:04
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
  2 件のコメント
IPSHITA
IPSHITA 2024 年 3 月 1 日 8:49
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
Torsten 2024 年 3 月 1 日 9:35
編集済み: Torsten 2024 年 3 月 1 日 12:04
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".

サインインしてコメントする。

回答 (2 件)

rickert
rickert 2019 年 8 月 18 日
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

Star Strider
Star Strider 2019 年 8 月 18 日
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.

カテゴリ

Help Center および File ExchangePhysics についてさらに検索

タグ

タグが未入力です。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by