基于Scilab的Lotka Volterra模型参数估计

Question

我正在尝试使用Scilab对Lotka-Volterra模型进行参数估计（我是一个完全新手）。当我尝试运行脚本时，Scilab警告不连贯的减法。我想我的问题与this topic相同，但是那里的解决方案使用了Matlab函数。

这是我的脚本：

// 1. Create Lotka Volterra function

function [dY]=LotkaVolterra(t,X,c,n,m,e)
    IngestC = c * X(1) * X(2)
    GrowthP = n * X(1)
    MortC = m * X(2)
    dY(1) = GrowthP - IngestC
    dY(2) = IngestC * e - MortC
endfunction

// 2. Define the Nonlinear Least Squares functions

function f = Differences ( x ) 
    // Returns the difference between the simulated differential 
    // equation and the experimental data.
    c = x(1) ;n = x(2);m = x(3);e = x(4);y0 = y_exp(1,:);t0 = 0
    y_calc=ode(y0',t0,t,list(LotkaVolterra,c,n,m,e)) 
    diffmat = y_calc' - y_exp
    f = diffmat(:)
endfunction 

function val = L_Squares ( x ) 
    // Computes the sum of squares of the differences.
    f = Differences ( x ) 
    val = sum(f.^2)
endfunction 

// Experimental data 
t = [0:19]'; 
H=[20,20,20,12,28,58,75,75,88,61,75,88,69,32,13,21,30,2,153,148];
L=[30,45,49,40,21,8,6,5,10,20,33,34,30,21,14,8,4,4,14,38];
y_exp=[H',L'];


// compute the model cost function
function [f, g, ind] = modelCost (x, ind)
    f = L_Squares ( x )
    g = derivative ( L_Squares , x )
endfunction

// use of optim function with loops to avoid local minimum 
tic
i=0
fitminx=zeros(4,100);
fitminy=zeros(1,100);
for c=[0:0.1:1]
    for n=[0:0.1:1]
        for m=[0:0.1:1]
            for e=[0:0.1:1]
                i=i+1
                x0 = [c;n;m;e]
                [ fopt , xopt , gopt ] = optim ( modelCost , x0 )
                fitminx(:,i)=xopt;
                fitminy(:,i)=fopt;
            end
        end
    end 
end
[a,b]=min(fitminy)
fitminx(:,a)
toc

错误消息是：

lsoda--  at t (=r1), mxstep (=i1) steps   
needed before reaching tout
      where i1 is :        500                                                  
      where r1 is :   0.4145715729197D+01                                       
Attention : Le résultat est peut être inexact.

 !--error 9 
Soustraction incohérente.
at line       4 of function Differences called by :  
at line       2 of function L_Squares called by :  
at line      16 of function %R_ called by :  
at line      15 of function %deriv1_ called by :  
at line      58 of function derivative called by :  
at line       3 of function modelCost called by :  
                [ fopt , xopt , gopt ] = optim ( modelCost , x0 )

感谢您对我的问题的兴趣和时间（对不起我的英语）

Answer 1

我的答案here重复

问题在于求解器以某种方式达到了无法解决每个t上的颂歌并在某一点停止的点。因此，你的y_calc比y_exp小。

如果这对您来说不是问题，可以通过将diffmat函数的第6行上的Differences更改为

diffmat = y_calc' - y_exp(1:size(y_calc',1),:)

Answer 2

你的问题是在优化过程中c,n,m,e参数得到负值。只需在optim调用中添加约束，如下所示：

[fopt, xopt, gopt] = optim(modelCost, 'b', zeros(4,1), %inf*ones(4,1), x0)

Answer 3

问题在于

diffmat = y_calc' - y_exp

添加以下代码：

disp( "Y_calc dimensions:");
disp( size(y_calc'));

disp( "y_exp dimensions:");
disp( size(y_exp));    

实测：

Y_calc dimensions:    
91.    2.  

y_exp dimensions:    
20.    2.  

我对预期的行为和预期的矩阵大小一无所知，但这至少是错误的根本原因。