我写的代码如下：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.