我写的代码如下:function ODEfunctiongroup_boke_wo
% 动力学ODE方程模型的参数估计
% 此例数据只有t,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10
% k1-k12 为待拟合参数
% The variables y here are y(1)=x1, y(2)=x2, y(3)=x3,y(4)=x4,y(5)=x5,y(6)=x6,y(7)=x7,y(8)=x8,y(9)=x9,y(10)=x10.
%
clear all;
clc;
global t0
k0 = [5 5 5 5 5 5 5 5 5 5 5 5]; % 参数初值
lb = [0 0 0 0 0 0 0 0 0 0 0 0]; % 参数下限
ub = [+inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf]; % 参数上限
x0 = [2276.45 20.11 0.029 0.0010 0.19 0.23 0.12 0.045 119.24 0.089]; %初始状态 %动力学数据:
% t x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
ExpData = ...
[0 2276.45 20.11 0.03 0.00 0.19 0.23 0.12 0.04 119.24 0.09
4 2003.47 67.60 0.22 0.00 0.57 0.19 0.18 0.29 746.86 0.16
8 2032.79 77.54 0.29 0.00 0.93 0.15 0.24 0.50 1101.89 0.23
12 2149.09 56.71 0.21 0.00 0.94 0.15 0.24 0.47 1176.67 0.35
16 2011.55 42.71 0.15 0.01 0.89 0.14 0.26 0.47 995.61 0.33
20 1839.29 48.31 0.15 0.01 0.79 0.13 0.27 0.56 864.48 0.52
24 2003.07 34.93 0.19 0.02 0.88 0.11 0.28 0.79 737.48 0.60
28 1923.94 37.12 0.16 0.02 1.07 0.13 0.30 0.96 670.42 0.60
32 1790.84 60.29 0.19 0.03 1.34 0.14 0.52 1.54 672.09 1.16
36 1909.44 45.17 0.19 0.03 1.25 0.13 0.48 0.90 710.46 1.29
40 1700.47 45.07 0.21 0.03 1.06 0.14 0.33 0.44 672.98 0.78
];
t0 = ExpData(:,1);
yexp = ExpData(:,2:11); % yexp: 对应实验数据[x1 x2 x3 x4 x5 x6 x7 x8 x9 x10]
%% 使用函数fmincon()进行参数估计
% [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,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('\tk9 = %.4f\n',k(9))
% fprintf('\tk10 = %.4f\n',k(10))
% fprintf('\tk11 = %.4f\n',k(11))
% fprintf('\tk12 = %.4f\n',k(12))
% fprintf('The sum of the squares is: %.1e\n\n',fval)
% k_fmincon = k;
% 这一步通常被省略,通过反复迭代初始值得到最优解,加上后可以降低对初始值的依赖。
% 以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
% 需要指出,这种方法并非在所有场合均有效,但有时确实可以改善求解效果。
%% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
Output(k,ci,resnorm)
%% 以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
% k0 = k_fmincon;
% [k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
% lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);
% ci = nlparci(k,residual,jacobian);
% fprintf('\n\n以fmincon()的结果为初值,使用函数lsqnonlin()估计得到的参数值为:\n')
% Output(k,ci,resnorm)
%% ----------------------------------------------------
tspan = t0'; % 指定点时间,也可增加步数,只需多调用一次ode45函数
[t x] = ode45(@KineticEqs,tspan,x0,[],k);
y(:,1:10) = x(:,1:10);
figure(1)
plot(tspan,y(:,1),'b',ExpData(:,1),yexp(:,1),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y1(x1)');
figure(2)
plot(tspan,y(:,2),'b',ExpData(:,1),yexp(:,2),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y2(x2)');
figure(3)
plot(tspan,y(:,3),'b',ExpData(:,1),yexp(:,3),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y3(x3)');
figure(4)
plot(tspan,y(:,4),'b',ExpData(:,1),yexp(:,4),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y4(x4)');
figure(5)
plot(tspan,y(:,5),'b',ExpData(:,1),yexp(:,5),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y5(x5)');
figure(6)
plot(tspan,y(:,6),'b',ExpData(:,1),yexp(:,6),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y6(x6)');
figure(7)
plot(tspan,y(:,7),'b',ExpData(:,1),yexp(:,7),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y7(x7)');
figure(8)
plot(tspan,y(:,8),'b',ExpData(:,1),yexp(:,8),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y8(x8)');
figure(9)
plot(tspan,y(:,9),'b',ExpData(:,1),yexp(:,9),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y9(x9)');
figure(4)
plot(tspan,y(:,10),'b',ExpData(:,1),yexp(:,10),'or'),legend('计算值','实验值','Location','best')
xlabel('时间');ylabel('y10(x10)');
end
%% ------------------------------------------------------------------
function Output(k,ci,resnorm)
% Output.m
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))
fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))
fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))
fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))
fprintf('\tk9 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9))
fprintf('\tk10 = %.4f ± %.4f\n',k(10),ci(10,2)-k(10))
fprintf('\tk11 = %.4f ± %.4f\n',k(11),ci(11,2)-k(11))
fprintf('\tk12 = %.4f ± %.4f\n',k(12),ci(12,2)-k(12))
fprintf('The sum of the squares is: %.1e\n\n',resnorm)
end
%% ------------------------------------------------------------------
% function f = ObjFunc4Fmincon(k,x0,yexp)
% global t0
% tspan = t0';
% [t x] = ode45(@KineticEqs,tspan,x0,[],k);
% y(:,1:10) = x(:,1:10); % 对应实验数据 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
% f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2) ...
% + sum((y(:,3)-yexp(:,3)).^2) + sum((y(:,4)-yexp(:,4)).^2)+sum((y(:,5)-yexp(:,5)).^2)+sum((y(:,6)-yexp(:,6)).^2)
% +sum((y(:,7)-yexp(:,7)).^2)+sum((y(:,8)-yexp(:,8)).^2)+sum((y(:,9)-yexp(:,9)).^2)+sum((y(:,10)-yexp(:,10)).^2); % 计算平方和,供fmincon调用
% end
%% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)
global t0
tspan = t0';
[t x] = ode45(@KineticEqs,tspan,x0,[],k);
y(:,1:10) = x(:,1:10); % 对应实验数据 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f5 = y(:,5) - yexp(:,5);
f6 = y(:,6) - yexp(:,6);
f7 = y(:,7) - yexp(:,7);
f8 = y(:,8) - yexp(:,8);
f9 = y(:,9) - yexp(:,9);
f10 = y(:,10) - yexp(:,10);
f = [f1; f2; f3; f4;f5;f6;f7;f8;f9;f10]; % 理论值与实验值的差值,残差
% figure(5)
% subplot(2,2,1)
% plot(f1);
% subplot(2,2,2)
% plot(f2);
% subplot(2,2,3)
% plot(f3);
% subplot(2,2,4)
% plot(f4);
end
%% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k) % 微分方程组
dxdt = ...
[(-k(1)*x(1));
(k(1)*x(1)-k(2)*x(3)-k(5)*x(3)-k(7)*x(3)-k(9)*x(3));
(k(2)*x(3)-k(3)*x(4)+k(13)*x(5));
(k(7)*x(3)-k(8)*x(5)-k(12)*x(5)-k(13)*x(5));
(k(9)*x(3)-k(10)*x(6)-k(11)*x(6)-k(14)*x(6));
(k(5)*x(3)-k(6)*x(7)+k(14)*x(6));
(k(8)*x(5)+k(11)*x(6));
(k(3)*x(4));
(k(12)*x(5)+k(10)*x(6));
(k(6)*x(7))
];
end
% code end
错误提示:索引超出数组元素的数目(12)。
出错 ODEfunctiongroup_boke_wo>KineticEqs (line 159)
(k(2)*x(3)-k(3)*x(4)+k(13)*x(5));
出错 odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
出错 ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
出错 ODEfunctiongroup_boke_wo>ObjFunc4LNL (line 131)
[t x] = ode45(@KineticEqs,tspan,x0,[],k);
出错 lsqnonlin (line 215)
initVals.F = feval(funfcn{3},xCurrent,varargin{:});
出错 ODEfunctiongroup_boke_wo (line 54)
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);
原因: Failure in initial objective function evaluation. LSQNONLIN cannot continue.