使用最小平方的非线性 ODE 优化答案

【问题标题】：nonlinear ODEs optimization with leastsq使用最小平方的非线性 ODE 优化
【发布时间】：2021-11-30 06:28:55
【问题描述】：

[更新] 我正在研究非线性 ODE 系统优化并将其拟合到实验数据。我有一个由 5 个模型 ODE 组成的系统，必须通过 17 个参数进行优化。我的方法是计算求解的 ODE 和实验数据之间的差异 - 函数差异，然后使用最小平方求解器最小化差异并找到最佳参数，如下代码：

//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);

现在的结果是：

运行此脚本会导致这样的错误：

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”——这些都没有解决问题。在 ode 中使用“修复”方法出现此错误：

ode: rksimp exit with state 3.

请指教如何解决这个问题？

附：文件“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

【问题讨论】：

尝试强制“僵硬”的方法。除此之外，leastsq 很有可能使用参数的非物理值进行调用。在 leastsq 调用中添加约束。
在 leastsq 中添加了约束并尝试了“stiff”方法 - 没有结果，但给出了不同的错误：``` lsode-- at t (=r1) with step h (=r2), Corrector不与 abs(h) = hmin 收敛，其中 r1 为：0.1366396046954D+01 和 r2：0.6917767912662D-16 反复收敛失败（可能是错误的雅可比提供或雅可比类型或公差的错误选择）```
你能用新代码更新问题吗，包括约束和实际数据分配到Dat ？
请在上面找到更新
我在代码中做了修改（这样就可以直接在Scilab中执行脚本了。你应该在residual函数中显示参数的值，以了解ode求解器是否失败，因为奇怪的值。它可以帮助增加更严格的界限。

标签： optimization ode least-squares scilab

【解决方案1】：

在 Scilab 中，leastsq（基于 optim）非常差，并且不具有全局收敛性，不像 ipopt 可以作为原子模块使用。像这样安装它：

--> atomsInstall sci_ipopt

我已通过以下方式修改了您的脚本

为这种生物动力学保留 ode 的经典用法，即“僵硬”（使用 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 的有限差分梯度非常不准确。

希望对你有所帮助...

【讨论】：

非常感谢。花了点时间看懂了代码，认识了 ipopt 工具箱（顺便，很高兴认识这个模块的作者）。我已经稍微简化了我的模型 ODE，结果令人满意。再次感谢！
太棒了！仅供参考，它是什么样的模型？
使用 Monod 动力学的 Diauxic 细菌生长模型。