[MATLAB数学相关] 关于带有微分方程组的参数最优化估计问题

[星竹石]

最近刚刚接触matlab， 打算用matlab的lsqnonlin函数去解决以下问题：
某一反应体系可以用如下常微分方程组描述
dc1/dt=k1*sqrt(c3)
dc2/dt=k2*c3
dc3/dt=-k1*sqrt(c3)-k2*c3-k3*c3

现c1,c2和c3随时间的变化已经通过实验得到，想拟和出合适的k1, k2和k3.

我自己编了一个小程序，但总是被提示“too many output arguments”, 不知道是怎么回事， 希望大家能帮助一下。
我的程序如下

function KineticsEst5

clear all
clc

k0 = [0.5 0.5 0.5]; % 参数初值
lb = [0 0 0]; % 参数下限
ub = [+inf +inf +inf]; % 参数上限
c0 = [212.650 474.378 208070];
KineticsData2;
yexp = ExpData(:,2:4); % yexp: 实验数据[c1 c2 c3]

% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],c0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
Output

% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,c0,yexp)
tspan = [0 1 2 3 5 7 9 12 14 16 19 21 23 26];
[t c] = ode45(@KineticEqs,tspan,c0,[],k);
y(:,1:3) = c(:,1:3);
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f = [f1; f2; f3];

% ------------------------------------------------------------------
function dcdt = KineticEqs(t,c,k)
dcdt = ...
[ ( k(1)*sqrt(c(3)) )
( k(2)*c(3) )
( -k(1)*sqrt(c(3))-k(2)*c(3)-k(3)*c(3) )
];

程序所调用的实验数据如下

% 动力学数据:
% t c1 c2 c3
ExpData = ...
[
0.000 212.650 474.378 208070.000
1.000 2787.164 1745.013 164290.000
2.000 5163.403 2734.322 128720.000
3.000 7536.437 3207.740 114670.000
5.000 10579.022 4986.831 31950.000
7.000 11776.221 6099.661 25140.000
9.000 12455.458 6839.463 12610.000
12.000 12850.537 7157.148 11690.000
14.000 12981.767 7436.120 7910.000
16.000 12982.335 7517.885 9520.000
19.000 13079.799 7527.819 8650.000
21.000 13520.587 7586.632 8200.000
23.000 13565.417 7526.917 7630.000
26.000 13673.061 7477.547 8980.000
]