使用 leastsq 进行非线性 ODE 优化
nonlinear ODEs optimization with leastsq
[已更新] 我正在研究非线性 ODE 系统优化并将其拟合到实验数据。我有一个包含 5 个模型 ODE 的系统，必须通过 17 个参数对其进行优化。我的方法是计算求解的 ODE 和实验数据之间的差异 - 函数差异，然后使用 leastsq 求解器最小化差异并找到最佳参数，如下代码：
//RHSs of ODEs to be fitted:
function dx=model3(t,x,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H)
X=x(1);
S=x(2);
A=x(3);
DO=x(4);
V=x(5);`
qs=((q_Smax*S/(S+Ks))*Kia/(Kia+A));
qsof=(p_Amax*qs/(qs+Kap));
qsox=(qs-qsof)*DO/(DO+Ko);
qsa=(q_Amax*A/(A+Ksa))*(Kis/(qs+Kis));
pa=qsof*Yas;
qa=pa-qsa;
qo=(qsox-qm)*Yos+qsa*Yoa;
u=(qsox-qm)*Yem+qsof*Yxsof+qsa*Yxa;
dx(1)=u*X-F*X/V;
dx(2)=(F*(Sf-S)/V)-qs*X;
dx(3)=qsa*X-(F*A/V);
dx(4)=200*(100-DO)-qo*X*H;
dx(5)=F;
endfunction
//experimental data:
//Dat=fscanfMat('dane_exper_III_etap.txt');
Dat = [
0 30 1.4 24.1 99 6884.754
1 35 0.2 23.2 89 6959.754
2 40 0.1 21.6 80 7034.754
3 52 0.1 19.5 67 7109.754
4 61 0.1 18.7 70 7184.754
5 66 0.1 16.4 79 7259.754
6 71 0.1 15 94 7334.754
7 74 0 14.3 100 7409.754
8 76 0 13.8 100 7484.754
9 78 0 13.4 100 7559.754
9.5 79 0 13.2 100 7597.254
10 79 0 13.5 100 7634.754]
t=Dat(:,1);
x_exp(:,1)=Dat(:,2);
x_exp(:,2)=Dat(:,3);
x_exp(:,3)=Dat(:,4);
x_exp(:,4)=Dat(:,5);
x_exp(:,5)=Dat(:,6);
global MYDATA;
MYDATA.t=t;
MYDATA.x_exp=x_exp;
MYDATA.funeval=0;
//calculating differences between calculated values and experimental data:
function f=Differences(k)
global MYDATA
t=MYDATA.t;
x_exp=MYDATA.x_exp;
Kap=k(1); //g/L
Ksa=k(2); //g/L
Ko=k(3); //g/L
Ks=k(4); //g/L
Kia=k(5); //g/L
Kis=k(6); //g/L
p_Amax=k(7); //g/(g*h)
q_Amax=k(8); //g/(g*h)
qm=k(9);
q_Smax=k(10);
Yas=k(11); //g/g
Yoa=k(12);
Yxa=k(13);
Yem=k(14);
Yos=k(15);
Yxsof=k(16);
H=k(17);
x0=x_exp(1,:);
t0=0;
F=75;
Sf=500;
%ODEOPTIONS=[1,0,0,%inf,0,2,10000,12,5,0,-1,-1]
x_calc=ode('rk',x0',t0,t,list(model3,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H));
diffmat=x_calc'-x_exp;
//column vector of differences (concatenates 4 columns of the difference matrix)
f=diffmat(:);
MYDATA.funeval=MYDATA.funeval+1;
endfunction
// Initial guess
Kap=0.3; //g/L
Ksa=0.05; //g/L
Ko=0.1; //g/L
Ks=0.5; //g/L
Kia=0.5; //g/L
Kis=0.05; //g/L
p_Amax=0.4; //g/(g*h)
q_Amax=0.8; //g/(g*h)
qm=0.2;
q_Smax=0.6;
Yas=0.5; //g/g
Yoa=0.5;
Yxa=0.5;
Yem=0.5;
Yos=1.5;
Yxsof=0.22;
H=1000;
y0=[Kap;Ksa;Ko;Ks;Kia;Kis;p_Amax;q_Amax;qm;q_Smax;Yas;Yoa;Yxa;Yem;Yos;Yxsof;H];
yinf=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,100];
ysup=[%inf,%inf,%inf,%inf,%inf,%inf,3,3,3,3,3,3,3,3,3,3,10000];
[fopt,xopt,gopt]=leastsq(Differences,'b',yinf,ysup,y0);
现在的结果是：
0.2994018
0.0508325
0.0999987
0.4994088
0.5081272
0.
0.4004560
0.7050746
0.2774195
0.6068328
0.5
0.4926150
0.4053860
0.5255006
1.5018725
0.2193901
1000.0000
33591.642
运行 这个脚本导致这样的错误：
lsoda-- caution... t (=r1) and h (=r2) are
such that t + h = t at next step
(h = pas). integration continues
where r1 is : 0.5658105345269D+01 and r2 : 0.1884898700920D-17
lsoda-- previous message precedent given i1 times
will no more be repeated
where i1 is : 10
lsoda-- at t (=r1), mxstep (=i1) steps
needed before reaching tout
where i1 is : 500000
where r1 is : 0.5658105345270D+01
Excessive work done on this call (perhaps wrong jacobian type).
at line 27 of function Differences
我知道问题出在 ODE 求解步骤上。因此，我尝试更改 mxstep，同时将解决方法类型也改为 'adams'、'rk' 和 'stiff' - none，这解决了问题。在 ode 中使用 'fix' 方法我得到这个错误：
ode: rksimp exit with state 3.
请问如何解决？
P.S。文件中的实验数据'dane_exper_III_etap.txt':
0 30 1.4 24.1 99 6884.754
1 35 0.2 23.2 89 6959.754
2 40 0.1 21.6 80 7034.754
3 52 0.1 19.5 67 7109.754
4 61 0.1 18.7 70 7184.754
5 66 0.1 16.4 79 7259.754
6 71 0.1 15 94 7334.754
7 74 0 14.3 100 7409.754
8 76 0 13.8 100 7484.754
9 78 0 13.4 100 7559.754
9.5 79 0 13.2 100 7597.254
10 79 0 13.5 100 7634.754
在 Scilab 中 leastsq（基于 optim）非常差并且没有全局收敛特性，这与作为原子模块可用的 ipopt 不同。像这样安装：
--> atomsInstall sci_ipopt
我已经按照以下方式修改了你的脚本
- 保持这种生物动力学的经典用法，即“stiff”（使用 BDF 方法）。您使用的 Runge-Kutta 非常差，因为它是一种仅适用于温和 ODE 的显式方法。
- 使用
ipopt代替leastsq
- 在残差计算周围使用 try/catch/end 块以捕获对 ode 求解器的失败调用。
- 对残差使用一些权重。你应该玩玩它以提高合身性
- 使用严格的正下限而不是 0，因为某些参数的值非常低会使 ode 求解器失败。
- 添加一个绘图回调，它还会保存参数的当前值，以防您使用 ctrl-C 停止优化。
//RHSs of ODEs to be fitted:
function dx=model3(t,x,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H)
X=x(1);
S=x(2);
A=x(3);
DO=x(4);
V=x(5);
qs=((q_Smax*S/(S+Ks))*Kia/(Kia+A));
qsof=(p_Amax*qs/(qs+Kap));
qsox=(qs-qsof)*DO/(DO+Ko);
qsa=(q_Amax*A/(A+Ksa))*(Kis/(qs+Kis));
pa=qsof*Yas;
qa=pa-qsa;
qo=(qsox-qm)*Yos+qsa*Yoa;
u=(qsox-qm)*Yem+qsof*Yxsof+qsa*Yxa;
dx(1)=u*X-F*X/V;
dx(2)=(F*(Sf-S)/V)-qs*X;
dx(3)=qsa*X-(F*A/V);
dx(4)=200*(100-DO)-qo*X*H;
dx(5)=F;
endfunction
//calculating differences between calculated values and experimental data:
function [f,x_calc]=Differences(k, t, x_exp)
Kap=k(1); //g/L
Ksa=k(2); //g/L
Ko=k(3); //g/L
Ks=k(4); //g/L
Kia=k(5); //g/L
Kis=k(6); //g/L
p_Amax=k(7); //g/(g*h)
q_Amax=k(8); //g/(g*h)
qm=k(9);
q_Smax=k(10);
Yas=k(11); //g/g
Yoa=k(12);
Yxa=k(13);
Yem=k(14);
Yos=k(15);
Yxsof=k(16);
H=k(17);
x0=x_exp(1,:);
t0=0;
F=75;
Sf=500;
[x_calc]=ode("stiff",x0',t0,t,list(model3,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H));
diffmat=(x_calc'-x_exp)*residual_weight;
//column vector of differences (concatenates 4 columns of the difference matrix)
f=diffmat(:);
MYDATA.funeval=MYDATA.funeval+1;
endfunction
function [f,g]=normdiff2(k,new_k,t,x_exp)
try
res = Differences(k,t,x_exp)
if argn(1) == 2
JacRes = numderivative(list(Differences,t,x_exp),k)
g = 2*JacRes'*res;
end
f = sum(res.*res)
catch
f=%inf;
g=%inf*ones(k);
end
endfunction
function out=callback(param)
global MYDATA
if isfield(param,"x")
k = param.x;
MYDATA.k = k;
[f,x_calc]=Differences(k,t,x_exp)
plot_weight = diag(1./max(x_exp,'r'));
drawlater
clf
plot(t,x_exp*plot_weight,'-o')
plot(t,x_calc'*plot_weight,'-x')
legend X S A DO X
drawnow
end
out = %t;
endfunction
//experimental data:
//Dat=fscanfMat('dane_exper_III_etap.txt');
Dat = [
0 30 1.4 24.1 99 6884.754
1 35 0.2 23.2 89 6959.754
2 40 0.1 21.6 80 7034.754
3 52 0.1 19.5 67 7109.754
4 61 0.1 18.7 70 7184.754
5 66 0.1 16.4 79 7259.754
6 71 0.1 15 94 7334.754
7 74 0 14.3 100 7409.754
8 76 0 13.8 100 7484.754
9 78 0 13.4 100 7559.754
9.5 79 0 13.2 100 7597.254
10 79 0 13.5 100 7634.754]
t=Dat(:,1);
x_exp(:,1)=Dat(:,2);
x_exp(:,2)=Dat(:,3);
x_exp(:,3)=Dat(:,4);
x_exp(:,4)=Dat(:,5);
x_exp(:,5)=Dat(:,6);
global MYDATA;
MYDATA.funeval=0;
// Initial guess
Kap=0.3; //g/L
Ksa=0.05; //g/L
Ko=0.1; //g/L
Ks=0.5; //g/L
Kia=0.5; //g/L
Kis=0.05; //g/L
p_Amax=0.4; //g/(g*h)
q_Amax=0.8; //g/(g*h)
qm=0.2;
q_Smax=0.6;
Yas=0.5; //g/g
Yoa=0.5;
Yxa=0.5;
Yem=0.5;
Yos=1.5;
Yxsof=0.22;
H=100;
k0 = [Kap;Ksa;Ko;Ks;Kia;Kis;p_Amax;q_Amax;qm;q_Smax;Yas;Yoa;Yxa;Yem;Yos;Yxsof;H];
residual_weight = diag(1./[79,1.4, 24.1, 100, 7634.754]);
BIG = 1000;
SMALL = 1e-3;
problem = struct();
problem.f = list(normdiff2,t,x_exp);
problem.x0 = k0;
problem.x_lower = [SMALL*ones(16,1);100];
problem.x_upper = [BIG,BIG,BIG,BIG,BIG,BIG,3,3,3,3,3,3,3,3,3,3,10000]';
problem.int_cb = callback;
problem.params = struct("max_iter",200);
//
k = ipopt(problem)
这是 100 次迭代后的结果图（您可以更改 ipopt 选项中的值）。但是，不要指望正常终止
- 几乎可以确定你的参数集是无法识别的
- ODE 的有限差分梯度非常不准确。
希望对您有所帮助...
[已更新] 我正在研究非线性 ODE 系统优化并将其拟合到实验数据。我有一个包含 5 个模型 ODE 的系统，必须通过 17 个参数对其进行优化。我的方法是计算求解的 ODE 和实验数据之间的差异 - 函数差异，然后使用 leastsq 求解器最小化差异并找到最佳参数，如下代码：
//RHSs of ODEs to be fitted:
function dx=model3(t,x,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H)
X=x(1);
S=x(2);
A=x(3);
DO=x(4);
V=x(5);`
qs=((q_Smax*S/(S+Ks))*Kia/(Kia+A));
qsof=(p_Amax*qs/(qs+Kap));
qsox=(qs-qsof)*DO/(DO+Ko);
qsa=(q_Amax*A/(A+Ksa))*(Kis/(qs+Kis));
pa=qsof*Yas;
qa=pa-qsa;
qo=(qsox-qm)*Yos+qsa*Yoa;
u=(qsox-qm)*Yem+qsof*Yxsof+qsa*Yxa;
dx(1)=u*X-F*X/V;
dx(2)=(F*(Sf-S)/V)-qs*X;
dx(3)=qsa*X-(F*A/V);
dx(4)=200*(100-DO)-qo*X*H;
dx(5)=F;
endfunction
//experimental data:
//Dat=fscanfMat('dane_exper_III_etap.txt');
Dat = [
0 30 1.4 24.1 99 6884.754
1 35 0.2 23.2 89 6959.754
2 40 0.1 21.6 80 7034.754
3 52 0.1 19.5 67 7109.754
4 61 0.1 18.7 70 7184.754
5 66 0.1 16.4 79 7259.754
6 71 0.1 15 94 7334.754
7 74 0 14.3 100 7409.754
8 76 0 13.8 100 7484.754
9 78 0 13.4 100 7559.754
9.5 79 0 13.2 100 7597.254
10 79 0 13.5 100 7634.754]
t=Dat(:,1);
x_exp(:,1)=Dat(:,2);
x_exp(:,2)=Dat(:,3);
x_exp(:,3)=Dat(:,4);
x_exp(:,4)=Dat(:,5);
x_exp(:,5)=Dat(:,6);
global MYDATA;
MYDATA.t=t;
MYDATA.x_exp=x_exp;
MYDATA.funeval=0;
//calculating differences between calculated values and experimental data:
function f=Differences(k)
global MYDATA
t=MYDATA.t;
x_exp=MYDATA.x_exp;
Kap=k(1); //g/L
Ksa=k(2); //g/L
Ko=k(3); //g/L
Ks=k(4); //g/L
Kia=k(5); //g/L
Kis=k(6); //g/L
p_Amax=k(7); //g/(g*h)
q_Amax=k(8); //g/(g*h)
qm=k(9);
q_Smax=k(10);
Yas=k(11); //g/g
Yoa=k(12);
Yxa=k(13);
Yem=k(14);
Yos=k(15);
Yxsof=k(16);
H=k(17);
x0=x_exp(1,:);
t0=0;
F=75;
Sf=500;
%ODEOPTIONS=[1,0,0,%inf,0,2,10000,12,5,0,-1,-1]
x_calc=ode('rk',x0',t0,t,list(model3,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H));
diffmat=x_calc'-x_exp;
//column vector of differences (concatenates 4 columns of the difference matrix)
f=diffmat(:);
MYDATA.funeval=MYDATA.funeval+1;
endfunction
// Initial guess
Kap=0.3; //g/L
Ksa=0.05; //g/L
Ko=0.1; //g/L
Ks=0.5; //g/L
Kia=0.5; //g/L
Kis=0.05; //g/L
p_Amax=0.4; //g/(g*h)
q_Amax=0.8; //g/(g*h)
qm=0.2;
q_Smax=0.6;
Yas=0.5; //g/g
Yoa=0.5;
Yxa=0.5;
Yem=0.5;
Yos=1.5;
Yxsof=0.22;
H=1000;
y0=[Kap;Ksa;Ko;Ks;Kia;Kis;p_Amax;q_Amax;qm;q_Smax;Yas;Yoa;Yxa;Yem;Yos;Yxsof;H];
yinf=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,100];
ysup=[%inf,%inf,%inf,%inf,%inf,%inf,3,3,3,3,3,3,3,3,3,3,10000];
[fopt,xopt,gopt]=leastsq(Differences,'b',yinf,ysup,y0);
现在的结果是：
0.2994018
0.0508325
0.0999987
0.4994088
0.5081272
0.
0.4004560
0.7050746
0.2774195
0.6068328
0.5
0.4926150
0.4053860
0.5255006
1.5018725
0.2193901
1000.0000
33591.642
运行 这个脚本导致这样的错误：
lsoda-- caution... t (=r1) and h (=r2) are
such that t + h = t at next step
(h = pas). integration continues
where r1 is : 0.5658105345269D+01 and r2 : 0.1884898700920D-17
lsoda-- previous message precedent given i1 times
will no more be repeated
where i1 is : 10
lsoda-- at t (=r1), mxstep (=i1) steps
needed before reaching tout
where i1 is : 500000
where r1 is : 0.5658105345270D+01
Excessive work done on this call (perhaps wrong jacobian type).
at line 27 of function Differences
我知道问题出在 ODE 求解步骤上。因此，我尝试更改 mxstep，同时将解决方法类型也改为 'adams'、'rk' 和 'stiff' - none，这解决了问题。在 ode 中使用 'fix' 方法我得到这个错误：
ode: rksimp exit with state 3.
请问如何解决？
P.S。文件中的实验数据'dane_exper_III_etap.txt':
0 30 1.4 24.1 99 6884.754
1 35 0.2 23.2 89 6959.754
2 40 0.1 21.6 80 7034.754
3 52 0.1 19.5 67 7109.754
4 61 0.1 18.7 70 7184.754
5 66 0.1 16.4 79 7259.754
6 71 0.1 15 94 7334.754
7 74 0 14.3 100 7409.754
8 76 0 13.8 100 7484.754
9 78 0 13.4 100 7559.754
9.5 79 0 13.2 100 7597.254
10 79 0 13.5 100 7634.754
在 Scilab 中 leastsq（基于 optim）非常差并且没有全局收敛特性，这与作为原子模块可用的 ipopt 不同。像这样安装：
--> atomsInstall sci_ipopt
我已经按照以下方式修改了你的脚本
- 保持这种生物动力学的经典用法，即“stiff”（使用 BDF 方法）。您使用的 Runge-Kutta 非常差，因为它是一种仅适用于温和 ODE 的显式方法。
- 使用
ipopt代替leastsq - 在残差计算周围使用 try/catch/end 块以捕获对 ode 求解器的失败调用。
- 对残差使用一些权重。你应该玩玩它以提高合身性
- 使用严格的正下限而不是 0，因为某些参数的值非常低会使 ode 求解器失败。
- 添加一个绘图回调，它还会保存参数的当前值，以防您使用 ctrl-C 停止优化。
//RHSs of ODEs to be fitted:
function dx=model3(t,x,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H)
X=x(1);
S=x(2);
A=x(3);
DO=x(4);
V=x(5);
qs=((q_Smax*S/(S+Ks))*Kia/(Kia+A));
qsof=(p_Amax*qs/(qs+Kap));
qsox=(qs-qsof)*DO/(DO+Ko);
qsa=(q_Amax*A/(A+Ksa))*(Kis/(qs+Kis));
pa=qsof*Yas;
qa=pa-qsa;
qo=(qsox-qm)*Yos+qsa*Yoa;
u=(qsox-qm)*Yem+qsof*Yxsof+qsa*Yxa;
dx(1)=u*X-F*X/V;
dx(2)=(F*(Sf-S)/V)-qs*X;
dx(3)=qsa*X-(F*A/V);
dx(4)=200*(100-DO)-qo*X*H;
dx(5)=F;
endfunction
//calculating differences between calculated values and experimental data:
function [f,x_calc]=Differences(k, t, x_exp)
Kap=k(1); //g/L
Ksa=k(2); //g/L
Ko=k(3); //g/L
Ks=k(4); //g/L
Kia=k(5); //g/L
Kis=k(6); //g/L
p_Amax=k(7); //g/(g*h)
q_Amax=k(8); //g/(g*h)
qm=k(9);
q_Smax=k(10);
Yas=k(11); //g/g
Yoa=k(12);
Yxa=k(13);
Yem=k(14);
Yos=k(15);
Yxsof=k(16);
H=k(17);
x0=x_exp(1,:);
t0=0;
F=75;
Sf=500;
[x_calc]=ode("stiff",x0',t0,t,list(model3,Kap,Ksa,Ko,Ks,Kia,Kis,p_Amax,q_Amax,qm,q_Smax,Yas,Yoa,Yxa,Yem,Yos,Yxsof,H));
diffmat=(x_calc'-x_exp)*residual_weight;
//column vector of differences (concatenates 4 columns of the difference matrix)
f=diffmat(:);
MYDATA.funeval=MYDATA.funeval+1;
endfunction
function [f,g]=normdiff2(k,new_k,t,x_exp)
try
res = Differences(k,t,x_exp)
if argn(1) == 2
JacRes = numderivative(list(Differences,t,x_exp),k)
g = 2*JacRes'*res;
end
f = sum(res.*res)
catch
f=%inf;
g=%inf*ones(k);
end
endfunction
function out=callback(param)
global MYDATA
if isfield(param,"x")
k = param.x;
MYDATA.k = k;
[f,x_calc]=Differences(k,t,x_exp)
plot_weight = diag(1./max(x_exp,'r'));
drawlater
clf
plot(t,x_exp*plot_weight,'-o')
plot(t,x_calc'*plot_weight,'-x')
legend X S A DO X
drawnow
end
out = %t;
endfunction
//experimental data:
//Dat=fscanfMat('dane_exper_III_etap.txt');
Dat = [
0 30 1.4 24.1 99 6884.754
1 35 0.2 23.2 89 6959.754
2 40 0.1 21.6 80 7034.754
3 52 0.1 19.5 67 7109.754
4 61 0.1 18.7 70 7184.754
5 66 0.1 16.4 79 7259.754
6 71 0.1 15 94 7334.754
7 74 0 14.3 100 7409.754
8 76 0 13.8 100 7484.754
9 78 0 13.4 100 7559.754
9.5 79 0 13.2 100 7597.254
10 79 0 13.5 100 7634.754]
t=Dat(:,1);
x_exp(:,1)=Dat(:,2);
x_exp(:,2)=Dat(:,3);
x_exp(:,3)=Dat(:,4);
x_exp(:,4)=Dat(:,5);
x_exp(:,5)=Dat(:,6);
global MYDATA;
MYDATA.funeval=0;
// Initial guess
Kap=0.3; //g/L
Ksa=0.05; //g/L
Ko=0.1; //g/L
Ks=0.5; //g/L
Kia=0.5; //g/L
Kis=0.05; //g/L
p_Amax=0.4; //g/(g*h)
q_Amax=0.8; //g/(g*h)
qm=0.2;
q_Smax=0.6;
Yas=0.5; //g/g
Yoa=0.5;
Yxa=0.5;
Yem=0.5;
Yos=1.5;
Yxsof=0.22;
H=100;
k0 = [Kap;Ksa;Ko;Ks;Kia;Kis;p_Amax;q_Amax;qm;q_Smax;Yas;Yoa;Yxa;Yem;Yos;Yxsof;H];
residual_weight = diag(1./[79,1.4, 24.1, 100, 7634.754]);
BIG = 1000;
SMALL = 1e-3;
problem = struct();
problem.f = list(normdiff2,t,x_exp);
problem.x0 = k0;
problem.x_lower = [SMALL*ones(16,1);100];
problem.x_upper = [BIG,BIG,BIG,BIG,BIG,BIG,3,3,3,3,3,3,3,3,3,3,10000]';
problem.int_cb = callback;
problem.params = struct("max_iter",200);
//
k = ipopt(problem)
这是 100 次迭代后的结果图（您可以更改 ipopt 选项中的值）。但是，不要指望正常终止
- 几乎可以确定你的参数集是无法识别的
- ODE 的有限差分梯度非常不准确。
希望对您有所帮助...