MATLAB计算动力学参数时出现误差大的问题

作者 ygqsgc: 来源: 小木虫 300 6 举报帖子

昨天修改了某个地方导致今天求解参数误差极大，但我看了半天也没找到原因，也忘记昨天修改的什么地方了
程序如下，图1是昨天的，图2 是今天的，差距有些大，希望大神能够解惑

CODE:

function Kinetics2

clear all; clc

format long

k0=[0 0 0 0 0 0 0 0];

tspan=[0 2 4 6 8 10 15 20 25 30 40 50 60 75 90 105 120];             %参数初值

lb=[0 0 0 0 0 0 0 0];  ub = [1 1 1 1 1 1 1 1]*1e1;  %lb、ub：参数下限和上限

Y0=[0 0];

Kinetics1=[0        0        0

2        0.101406296        0.025750304

4        0.148139381        0.077208075

6        0.171971874        0.097848253

8        0.200101464        0.131543106

10        0.223088343        0.149258466

15        0.268349285        0.206242074

20        0.319194899        0.249534539

25        0.352139308        0.280899049

30        0.373849138        0.309100951

40        0.4145939        0.348550324

50        0.442947836        0.381294016

60        0.460584916        0.395704095

75        0.483226647        0.421979897

90        0.499672966        0.430587898

105        0.50291736        0.452012783

120        0.518449036        0.452107899];

yexp=Kinetics1(1:17,2:3);            

%使用函数fmincon()进行参数估计

[k,fval,flag]=fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],Y0,yexp);

fprintf('\n 使用函数fmincon()估计得到的参数值为:\n')

fprintf('\tk1=%.4f\n',k(1)), fprintf('\tk2=%.4f\n',k(2)), fprintf('\tk3=%.4f\n',k(3)),fprintf('\tk4=%.4f\n',k(4)), fprintf('\tk5=%.4f\n',k(5)),fprintf('\tk6=%.4f\n',k(6)),fprintf('\tk7=%.4f\n',k(7)),fprintf('\tk8=%.4f\n',k(8));

fprintf(' The sum of the squares is:%.1e\n\n', fval);

k_fmincon=k;

ts=0:1:max(tspan);

[ts ys] = ode45(@func,ts,Y0,[],k);

yy = [Kinetics1(:,2) Kinetics1(:,3)];

plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo',ts,ys(:,2),'r',tspan,yy(:,2),'ro'),

legend('Y1的计算值','Y1的实验值','Y2的计算值','Y2的实验值','Location','best');

function f=ObjFunc4Fmincon(k,Y0,yexp)

tspan=[0 2 4 6 8 10 15 20 25 30 40 50 60 75 90 105 120];

[t,Y]=ode45(@func,tspan,Y0,[],k);

y(:,1)=Y(:,1);y(:,2:3)=Y(:,1:2);

f=sum((y(:,2)-yexp(:,1)).^2)+sum((y(:,3)-yexp(:,2)).^2)

function dYdt=func(t,Y,k)

k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);k7=k(7);k8=k(8);

f1=(k1.*k2.*k3.*k4.*(1-Y(1)).*6.25.*(1-Y(2)).*3.125-k5.*k6.*k7.*k8.*Y(1).*6.25.*Y(2).*3.125)./(((k1.*k2.*k7+k1.*k2.*k4).*(1-Y(1)).*6.25+(k5.*k3.*k4+k2.*k3.*k4).*(1-Y(2)).*3.125+(k5.*k6.*k7+k1.*k6.*k7).*Y(1).*6.25+(k2.*k6.*k8+k5.*k7.*k8).*Y(2).*3.125+(k1.*k3.*k4+k1.*k2.*k3).*(1-Y(1)).*6.25.*(1-Y(2)).*3.125+(k1.*k6.*k7+k5.*k6.*k4).*(1-Y(1)).*6.25.*Y(1).*6.25+(k2.*k4.*k8+k5.*k3.*k8).*(1-Y(2)).*3.125+(k6.*k7.*k8+k5.*k6.*k8).*Y(1).*6.25.*Y(2).*3.125).*6.125);

f2=(k1.*k2.*k3.*k4.*(1-Y(1)).*6.25.*(1-Y(2)).*3.125-k5.*k6.*k7.*k8.*Y(1).*6.25.*Y(2).*3.125)./(((k1.*k2.*k7+k1.*k2.*k4).*(1-Y(1)).*6.25+(k5.*k3.*k4+k2.*k3.*k4).*(1-Y(2)).*3.125+(k5.*k6.*k7+k1.*k6.*k7).*Y(1).*6.25+(k2.*k6.*k8+k5.*k7.*k8).*Y(2).*3.125+(k1.*k3.*k4+k1.*k2.*k3).*(1-Y(1)).*6.25.*(1-Y(2)).*3.125+(k1.*k6.*k7+k5.*k6.*k4).*(1-Y(1)).*6.25.*Y(1).*6.25+(k2.*k4.*k8+k5.*k3.*k8).*(1-Y(2)).*3.125+(k6.*k7.*k8+k5.*k6.*k8).*Y(1).*6.25.*Y(2).*3.125).*3.125);

dYdt=[f1;f2];

今日热帖