parameter estimation for a set of differential equations using fmincon function in matlab
15 ビュー (過去 30 日間)
古いコメントを表示
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
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
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
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
0 件のコメント
Star Strider
2019 年 8 月 18 日
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 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!