How to fit the exp data with model data using some other optimisation function
显示 更早的评论
global k x0
%Auranchytrium sp. temperature=250 oC
time=[0 5 10 15 20 25 30 35 40 45 50 55 60];
xs=[1 0.635 0.27 0.205 0.14 0.115 0.09 0.07 0.05 0.0525 0.055 0.0575 0.06]; % wt% solid
x1=[0.3215 0.20425 0.087 0.066 0.045 0.037 0.029 0.0225 0.016 0.01675 0.0175 0.01825 0.019]; % wt% protein
x2=[0.0622 0.0396 0.017 0.01285 0.0087 0.00715 0.0056 0.00435 0.0031 0.00325 0.0034 0.00355 0.0037]; % wt % carbohydrates
x3=[0.613 0.3895 0.166 0.126 0.086 0.0705 0.055 0.043 0.031 0.0325 0.034 0.0355 0.037]; % wt % lipids
x4=[0 0.165 0.33 0.375 0.42 0.435 0.45 0.46 0.47 0.4525 0.455 0.4575 0.46];% wt% aqueous phase
x5=[0 0.2 0.40 0.425 0.45 0.45 0.45 0.455 0.46 0.46 0.46 0.46 0.46]; % wt% bio crude
x6=[0 0 0 0 0 0.005 0.01 0.015 0.02 0.0225 0.025 0.0275 0.03]; % wt% gas
%xO= Initial guess
x=[x1' x2' x3' x4' x5' x6'];
x0=[0.3215 0.0622 0.6163 0 0 0 ];
p0=[0.085 0.1 0.06 0.06 0.045 0.08 0.085 0.08 0.0003 0.002];
pL=[0 0 0 0 0 0 0 0 0 0];
pU=[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5];
options = optimoptions(@lsqcurvefit,'Algorithm','levenberg-marquardt','MaxIter',10000,'TolX',1e-12);
[p,resnorm,res,EXITFLAG,OUTPUT,LAMBDA,jocob]=lsqcurvefit(@fun72500,p0,time,x,pL,pU);
ci=nlparci(p,res,jocob);
k1p=p(1);k1c=p(2);k1l=p(3);k2p=p(4);k2c=p(5);k2l=p(6);k3=p(7);k4=p(8);k5=p(9);k6=p(10);
k=[k1p k1c k1l k2p k2c k2l k3 k4 k5 k6 ];
fprintf('\n\n The estimated parameters values are;')
fprintf('\n k1p=%.4f,k1c=%.4f,k1l=%.4f,k2p=%.4f,k2c=%.4f,k2l=%.4f, k3=%.4f,k4=%.4f, k5=%.7f, k6=%.7f,',k1p,k1c,k1l,k2p,k2c,k2l,k3, k4, k5, k6)
y_m=fun72500(k,time);
figure(1)
plot(time,x(:,1),'*',time,y_m(:,1),'b',time,x(:,2),'*',time,y_m(:,2),'r',time,x(:,3),'*',time,y_m(:,3),'g',time,x(:,4),'*',time,y_m(:,4),'y',time,x(:,5),'*',time,y_m(:,5),'m',time,x(:,6),'*',time,y_m(:,6),'c');
legend({'Exp data','protein','Exp data','carbohydrates','Exp data','lipids','Exp data','aqueous phase ','Exp data','bio crude','Exp data','Gas'},'Location','northwest','NumColumns',3);
xlabel('Time(min)')
ylabel('wt.%')
function y_m=fun72500(p1,time)
global k x0
k=p1;
time;
[t,y_m]=ode23s(@model72500,time,x0);
end
function dydt=model72500(t,x)
global k x0
k1p=k(1);k1c=k(2);k1l=k(3);k2p=k(4);k2c=k(5);k2l=k(6);k3=k(7);k4=k(8);k5=k(9);k6=k(10);
y1=x(1);y2=x(2);y3=x(3);y4=x(4); y5=x(5); y6=x(6); % y1:protein; y2:carb; y3: lipid; y4: AP; y5: BC, y6: gas
dxpdt= -(k1p+k2p)*y1;
dxcdt= -(k1c+k2c)*y2;
dxldt= -(k1l+k2l)*y3;
dxaqdt=k1p*y1+k1c*y2+k1l*y3+k3*y5-(k4+k5)*y4;
dxbcdt=k2p*y1+k2c*y2+k2l*y3+k4*y4-(k3+k6)*y5;
dxgdt=(k5*y4+k6*y5);
dydt=[dxpdt;dxcdt;dxldt;dxbcdt;dxaqdt;dxgdt];
end
Any better idea to rearrange the code like giving values in function OR using function like fmicon(code for the same(because getting error while implementing)
4 个评论
Mohammad Nidal
2021-5-9
Walter Roberson
2021-5-9
Start by getting rid of the global variables.
Mohammad Nidal
2021-5-9
Walter Roberson
2021-5-9
I could, but I will not. The examples in the documentation are pretty clear, and you should easily be able to learn from them how to create appropriate anonymous functions.
回答(0 个)
类别
在 帮助中心 和 File Exchange 中查找有关 Physics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!