主要内容

三种生态种群系统:MATLAB和C MEX-File时间序列建模

这个例子展示了如何创建非线性灰盒时间序列模型。时间序列模型是不使用任何测量输入的模型。研究了三种理想的生态系统,其中两种物种:

  • 竞争同样的食物,或者:

  • 处于捕食者-被捕食者的情境中

该实例展示了基于MATLAB和C mex文件的建模。

生态人口系统

在所有三种被调查的种群系统中,我们对两种物种的种群如何随时间变化很感兴趣。为了对此进行建模,设x1(t)和x2(t)表示在时间t时各自物种的个体数量。设l1和l2分别表示与x1(t)和x2(t)相关的出生率,两者均假设为随时间变化的常数。物种的死亡率取决于食物的可获得性,如果有捕食者存在,则取决于被吃掉的风险。通常(而且通常),物种i(i = 1或2)的死亡率可以写成ui(x1(t), x2(t)),其中ui(.)是某个适当的函数。在实践中,这意味着ui(x1(t), x2(t))*xi(t)种i的动物每一个时间单位死亡。这些表述的净效果可以总结为状态空间类型的模型结构:(时间序列):

d

——x1(t) = l1*x1(t) - u1(x1(t), x2(t))*x1(t)

dt

d

——x2(t) = l2*x2(t) - u2(x1(t), x2(t))*x2(t)

dt

在这里,选择这两种状态作为输出是很自然的,即我们令y1(t) = x1(t)和y2(t) = x2(t)。

. 1。两个物种争夺同样的食物

如果两个物种争夺同样的食物,那么控制食物供应的是该物种的总体数量,进而控制它们的死亡率。一种简单而常见的方法是假设死亡率可以写成:

Ui (x1(t), x2(t)) = gi + di*(x1(t) + x2(t))

对于两个物种(i = 1或2),其中gi和di为未知参数。这就形成了状态空间结构:

d

——x1(t) = (l1-g1)*x1(t) - d1*(x1(t)+x2(t))*x1(t)

dt

d

——x2(t) = (l2-g2)*x2(t) - d2*(x1(t)+x2(t))*x2(t)

dt

这种结构的一个直接问题是l1、g1、l2和g2不能分别被识别。我们只能希望p1 = l1-g1和p3 = l2-g2。同样令p2 = d1, p4 = d2,得到重新参数化的模型结构:

d

——x1(t) = p1*x1(t) - p2*(x1(t)+x2(t))*x1(t)

dt

d

——x2(t) = p3*x2(t) - p4*(x1(t)+x2(t))*x2(t)

dt

Y1 (t) = x1(t)

Y2 (t) = x2(t)

在第一个人口示例中,我们求助于MATLAB文件建模。然后将上面的方程输入到MATLAB文件preys_m中。M,与以下内容。

函数[dx, y] = preys_m(t, x, u, p1, p2, p3, p4, varargin)

两个物种争夺同一种食物。

%输出方程。

Y = [x(1);...猎物种类

x(2)……猎物种类2。

];

状态方程。

Dx = [p1*x(1)-p2*(x(1)+x(2))*x(1);...猎物种类

p3 * x (2) p4 * x (x (1) + (2)) * x(2)……猎物种类2。

];

MATLAB文件以及初始参数向量、适当的初始状态和一些管理信息将作为输入输入到IDNLGREY对象构造函数。注意,参数和初始状态的初始值都被指定为结构数组,其中分别包含Npo(参数对象的数量=参数的数量,如果所有参数都是标量)和Nx(状态的数量)元素。通过这些结构数组,可以完全将非默认属性值(“名称”、“单位”、“值”、“最小值”、“最大值”和“固定值”)分配给每个参数和初始状态。在这里,我们将每个初始状态的“Minimum”值赋值为零(总体为正!),并指定在默认情况下估计两个初始状态。

文件名=“preys_m”;描述模型结构的文件。Order = [2 0 2];模型订单[ny nu nx]。参数= struct(“名字”, {“生存因素,物种1”“死亡因素,物种1”...“生存因素,物种2”“死亡因素,物种2”},...“单位”, {1 /年的1 /年的1 /年的1 /年的},...“价值”, {1.8 0.8 1.2 0.8},...“最低”, {-Inf -Inf -Inf},...“最大”, {Inf Inf Inf Inf},...“固定”,{假假假假假});估计所有4个参数。InitialStates = struct(“名字”, {“种群,物种1”“种群,物种2”},...“单位”, {“大小(千)”“大小(千)”},...“价值”, {0.2 1.8},...“最低”, {0 0},...“最大”, {Inf Inf},...“固定”, {false false});估计两个初始状态。Ts = 0;时间连续系统。nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts,...“名字”“两个物种争夺同一种食物”...“OutputName”, {“种群,物种1”“种群,物种2”},...“OutputUnit”, {“大小(千)”“大小(千)”},...“TimeUnit”“年”);

PRESENT命令可以用来查看初始模型的信息:

礼物(nlgr);
nlgr = 'preys_m'定义的连续时间非线性灰盒模型(MATLAB文件):dx/dt = F(t, x(t), p1,…, p4) y(t) = H(t, x(t), p1,…,p4) + e(t) with 0 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 4). States: Initial value x(1) Population, species 1(t) [Size (in t..] xinit@exp1 0.2 (estimated) in [0, Inf] x(2) Population, species 2(t) [Size (in t..] xinit@exp1 1.8 (estimated) in [0, Inf] Outputs: y(1) Population, species 1(t) [Size (in thousands)] y(2) Population, species 2(t) [Size (in thousands)] Parameters: Value p1 Survival factor, species 1 [1/year] 1.8 (estimated) in [-Inf, Inf] p2 Death factor, species 1 [1/year] 0.8 (estimated) in [-Inf, Inf] p3 Survival factor, species 2 [1/year] 1.2 (estimated) in [-Inf, Inf] p4 Death factor, species 2 [1/year] 0.8 (estimated) in [-Inf, Inf] Name: Two species competing for the same food Status: Created by direct construction or transformation. Not estimated. More information in model's "Report" property. Model Properties

由信用证。输入输出数据

接下来,我们加载(模拟的,但噪声已被破坏)数据,并创建一个IDDATA对象,描述两个物种争夺同一种食物的特定情况。该数据集包含201个数据样本,涵盖20年的进化。

负载(fullfile (matlabroot“工具箱”“识别”“iddemos”“数据”“preydata”));Z = iddata(y, [], 0.1,“名字”“两个物种争夺同一种食物”);集(z,“OutputName”, {“种群,物种1”“种群,物种2”},...“Tstart”0,“TimeUnit”“年”);

出具。初始两种模型的性能

对初始模型的模拟清楚地表明,它不能处理真实的人口动态。见图示。对于时间序列类型的IDNLGREY模型,请注意模型输出是由初始状态决定的。

比较(z, nlgr, 1);

图中包含2个轴对象。坐标轴对象1包含2个line类型的对象。这些对象代表验证数据(种群,物种1),两个物种竞争同一食物:67.36%。坐标轴对象2包含两个line类型的对象。这些对象代表验证数据(种群,物种2),两个物种竞争同一食物:67.83%。

图1:初始两种模型的真实输出与模拟输出的比较。

各。参数估计

为了克服初始模型相当差的性能,我们继续使用NLGREYEST对4个未知参数和2个初始状态进行估计。使用NLGREYESTOPTIONS指定估计选项;在本例中,“Display”被设置为“on”,这意味着估计进度信息显示在进度窗口中。你可以使用NLGREYESTOPTIONS来指定基本的算法属性,如“GradientOptions”,“SearchMethod”,“MaxIterations”,“Tolerance”,“Display”。

opt = nlgreyestOptions;opt.Display =“上”;opt.SearchOptions.MaxIterations = 50;NLGR = nlgreyest(z, NLGR, opt);

本。估计两物种模型的性能

参数和初始状态的估计值与用于生成真实输出数据的估计值非常一致:

disp (' True估计参数向量');
估计参数向量
Ptrue = [2;1;1;1);流(' %6.3f %6.3f\n'[ptrue ';getpvec (nlgr) '));
2.000 2.004 1.000 1.002 1.000 1.018 1.000 1.010
disp (' ');

               
disp (“真实的估计初始状态”);
估计初始状态
X0true = [0.1;2);流(' %6.3f %6.3f\n'[x0true ';cell2mat (getinit (nlgr,“价值”))));
0.100 0.101 2.000 1.989

为了进一步评估模型的质量(并说明与初始模型相比的改进),我们还模拟了估计的模型。在一个图窗口中将模拟输出与真实输出进行比较。可以看出,估计的模型是很好的。

比较(z, nlgr, 1);

图中包含2个轴对象。坐标轴对象1包含2个line类型的对象。这些对象代表验证数据(种群,物种1),两个物种竞争同一食物:98.42%。坐标轴对象2包含两个line类型的对象。这些对象代表验证数据(种群,物种2),两个物种竞争同一食物:97.92%。

图2:两种估计模型的真实输出与模拟输出的比较。

PRESENT提供了关于估计模型的进一步信息,例如关于参数不确定性和其他估计相关量,如损失函数和Akaike的FPE(最终预测误差)度量。

礼物(nlgr);
nlgr = 'preys_m'定义的连续时间非线性灰盒模型(MATLAB文件):dx/dt = F(t, x(t), p1,…, p4) y(t) = H(t, x(t), p1,…,p4) + e(t) with 0 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 4). States: Initial value x(1) Population, species 1(t) [Size (in t..] xinit@exp1 0.100729 (estimated) in [0, Inf] x(2) Population, species 2(t) [Size (in t..] xinit@exp1 1.98855 (estimated) in [0, Inf] Outputs: y(1) Population, species 1(t) [Size (in thousands)] y(2) Population, species 2(t) [Size (in thousands)] Parameters: ValueStandard Deviation p1 Survival factor, species 1 [1/year] 2.00429 0.00971109 (estimated) in [-Inf, Inf] p2 Death factor, species 1 [1/year] 1.00235 0.00501783 (estimated) in [-Inf, Inf] p3 Survival factor, species 2 [1/year] 1.01779 0.0229598 (estimated) in [-Inf, Inf] p4 Death factor, species 2 [1/year] 1.0102 0.0163506 (estimated) in [-Inf, Inf] Name: Two species competing for the same food Status: Termination condition: Near (local) minimum, (norm(g) < tol).. Number of iterations: 5, Number of function evaluations: 6 Estimated using Solver: ode45; Search: lsqnonlin on time domain data "Two species competing for the same food". Fit to estimation data: [98.42;97.92]% FPE: 7.747e-09, MSE: 0.0001743 More information in model's "Report" property. Model Properties

责任。一个经典的捕食者-猎物系统

假设第一个物种生活在第二个物种上。物种1的食物可得性与x2(t)(物种2的个体数)成正比,这意味着物种1的死亡率随着x2(t)的增加而降低。这个事实可以用简单的表达式来表达:

U1 (x1(t), x2(t)) = g1 - a1*x2(t)

其中g1和a1是未知参数。同样,当第一个物种的个体数量增加时,物种2的死亡率也会增加,例如:

U2 (x1(t), x2(t)) = g2 + a2*x1(t)

g2和a2是两个未知参数。利用上面假设的线性出生率,可以得到状态空间结构:

d

——x1(t) = (l1-g1)*x1(t) + a1*x2(t)*x1(t)

dt

d

——x2(t) = (l2-g2)*x2(t) - a2*x1(t)*x2(t)

dt

正如在前面的总体示例中一样,在这里也不可能唯一地识别这六个单独的参数。通过与上例相同的重参数化,即p1 = l1-g1, p2 = a1, p3 = l2-g2, p4 = a2,可得到如下模型结构:

d

——x1(t) = p1*x1(t) + p2*x2(t)*x1(t)

dt

d

——x2(t) = p3*x2(t) - p4*x1(t)*x2(t)

dt

Y1 (t) = x1(t)

Y2 (t) = x2(t)

从估计的角度来看,这更合适。

这次我们将此信息输入到名为predprey1_c.c的C mex文件中。模型文件的结构是标准的IDNLGREY C mex -文件(参见标题为“创建IDNLGREY模型文件”或idnlgreydemo2.m的示例),带有状态和输出更新函数compute_dx和compute_y,如下所示。

Void compute_dx(double *dx, double t, double *x, double *p,

常量mxArray *auxvar

/*检索模型参数。* /

双*p1, *p2, *p3, *p4;

P1 = p[0];/*生存因素,捕食者。* /

P2 = p[1];/*死亡因素,捕食者。* /

P3 = p[2];/*生存因素,猎物。* /

P4 = p[3];/*死亡因素,猎物。* /

/* x[0]:捕食者物种。* /

/* x[1]:猎物种类。* /

Dx [0] = p1[0]*x[0]+p2[0]*x[1]*x[0];

Dx [1] = p3[0]*x[1]-p4[0]*x[0]*x[1];

/*输出方程。* /

Void compute_y(双*y,双* t,双*x,双**p,

常量mxArray *auxvar

/* y[0]:捕食者物种。* /

/* y[1]:猎物种类。* /

Y [0] = x[0];

Y [1] = x[1];

因为模型是时间序列类型的,所以compute_dx和compute_y都不包括u在输入参数列表中。事实上,predprey1_c.c的主接口函数甚至没有声明u,尽管空的u([])总是通过IDNLGREY方法传递给predprey1_c。

编译后的C mex文件,以及初始参数向量、适当的初始状态和一些管理信息,将作为输入参数提供给IDNLGREY对象构造函数:

文件名=“predprey1_c”;描述模型结构的文件。Order = [2 0 2];模型订单[ny nu nx]。参数= struct(“名字”, {“生存因素,捕食者”“死亡因素,捕食者”...“生存因素,猎物”“死亡因素,猎物”},...“单位”, {1 /年的1 /年的1 /年的1 /年的},...“价值”, {-1.1 0.9 1.1 0.9},...“最低”, {-Inf -Inf -Inf},...“最大”, {Inf Inf Inf Inf},...“固定”,{假假假假假});估计所有4个参数。InitialStates = struct(“名字”, {“人口,捕食者”“人口,猎物”},...“单位”, {“大小(千)”“大小(千)”},...“价值”, {1.8 1.8},...“最低”, {0 0},...“最大”, {Inf Inf},...“固定”, {false false});估计两个初始状态。Ts = 0;时间连续系统。nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts,...“名字”“经典的1个捕食者- 1个猎物系统”...“OutputName”, {“人口,捕食者”“人口,猎物”},...“OutputUnit”, {“大小(千)”“大小(千)”},...“TimeUnit”“年”);

接下来通过PRESENT命令查看捕食者和猎物模型。

礼物(nlgr);
nlgr =由'predprey1_c' (MEX-file)定义的连续时间非线性灰盒模型:dx/dt = F(t, x(t), p1,…, p4) y(t) = H(t, x(t), p1,…,p4) + e(t) with 0 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 4). States: Initial value x(1) Population, predators(t) [Size (in t..] xinit@exp1 1.8 (estimated) in [0, Inf] x(2) Population, preys(t) [Size (in t..] xinit@exp1 1.8 (estimated) in [0, Inf] Outputs: y(1) Population, predators(t) [Size (in thousands)] y(2) Population, preys(t) [Size (in thousands)] Parameters: Value p1 Survival factor, predators [1/year] -1.1 (estimated) in [-Inf, Inf] p2 Death factor, predators [1/year] 0.9 (estimated) in [-Inf, Inf] p3 Survival factor, preys [1/year] 1.1 (estimated) in [-Inf, Inf] p4 Death factor, preys [1/year] 0.9 (estimated) in [-Inf, Inf] Name: Classical 1 predator - 1 prey system Status: Created by direct construction or transformation. Not estimated. More information in model's "Report" property. Model Properties

B.2。输入输出数据

我们的下一步是加载(模拟的,尽管噪声已经损坏)数据并创建一个IDDATA对象来描述这个特定的捕食者-猎物情况。该数据集还包含了涵盖20年进化的201个数据样本。

负载(fullfile (matlabroot“工具箱”“识别”“iddemos”“数据”“predprey1data”));Z = iddata(y, [], 0.1,“名字”“经典的1个捕食者- 1个猎物系统”);集(z,“OutputName”, {“人口,捕食者”“人口,猎物”},...“Tstart”0,“TimeUnit”“年”);

B.3。初始经典捕食者-被捕食者模型的表现

对初始模型的模拟表明,该模型不能准确地处理真实的种群动态。查看图窗口。

比较(z, nlgr, 1);

图中包含2个轴对象。坐标轴对象1包含2个line类型的对象。这些对象代表验证数据(种群,捕食者),经典1捕食者- 1被捕食者系统:4.842%。坐标轴对象2包含两个line类型的对象。这些对象表示验证数据(种群,猎物),经典1捕食者- 1猎物系统:4.523%。

图3:原始经典捕食者-被捕食者模型的真实输出与模拟输出的比较。

B.4。参数估计

为了提高初始模型的性能,我们继续使用NLGREYEST对4个未知参数和2个初始状态进行估计。

nlgr = nlgreyest(z, nlgr, nlgreyestOptions)“显示”“上”));

B.5。估计的经典捕食者-被捕食者模型的性能

参数和初始状态的估计值与生成真实输出数据时的估计值非常接近:

disp (' True估计参数向量');
估计参数向量
Ptrue = [-1;1;1;1);流(' %6.3f %6.3f\n'[ptrue ';getpvec (nlgr) '));
-1.000 -1.000 1.000 1.000 1.000 1.000 1.000 0.999
disp (' ');

               
disp (“真实的估计初始状态”);
估计初始状态
X0true = [2;2);流(' %6.3f %6.3f\n'[x0true ';cell2mat (getinit (nlgr,“价值”))));
2.000 2.002 2.000 2.000

为了进一步评估模型的质量(并说明与初始模型相比的改进),我们还模拟了估计的模型。在一个图窗口中将模拟输出与真实输出进行比较。再次可以看出,估计的模型是相当好的。

比较(z, nlgr, 1);

图中包含2个轴对象。坐标轴对象1包含2个line类型的对象。这些对象代表验证数据(种群,捕食者),经典1捕食者- 1猎物系统:98.08%。坐标轴对象2包含两个line类型的对象。这些对象代表验证数据(种群,猎物),经典的1捕食者- 1猎物系统:97.97%。

图4:经典捕食者-被捕食者模型的真实输出与模拟输出的比较。

正如预期的那样,PE返回的预测误差很小,具有随机性质。

体育(z, nlgr);

图中包含2个轴对象。axis对象1包含一个类型为line的对象。这些对象代表z(种群,捕食者),经典的1捕食者- 1被捕食者系统。Axes对象2包含一个类型为line的对象。这些对象代表z(种群,猎物),经典的1捕食者- 1猎物系统。

图5:用估计的IDNLGREY经典捕食者-食饵模型得到的预测误差。

C.1。具有猎物拥挤的捕食者-猎物系统

最后的种群研究也致力于1捕食者和1猎物系统,不同之处在于我们在这里引入了一个术语-p5*x2(t)^2,表示由于拥挤导致的猎物生长迟缓。从上一个例子中重新参数化的模型结构,然后用这个拥挤项进行补充:

d

——x1(t) = p1*x1(t) + p2*x2(t)*x1(t)

dt

d

——x2(t) = p3*x2(t) - p4*x1(t)*x2(t) - p5*x2(t)^2

dt

Y1 (t) = x1(t)

Y2 (t) = x2(t)

这些方程的解释基本上与上面相同,除了在没有捕食者的情况下,猎物数量的增长将受到抑制。

新的建模情况反映在一个名为predprey2_c.c的C mex文件中,它几乎与predprey1_c.c相同。除了与模型参数数量相关的更改(5而不是4),主要的区别在于状态更新函数compute_dx:

/*状态方程。* /

Void compute_dx(double *dx, double t, double *x, double *p,

常量mxArray *auxvar

/*检索模型参数。* /

双*p1, *p2, *p3, *p4, *p5;

P1 = p[0];/*生存因素,捕食者。* /

P2 = p[1];/*死亡因素,捕食者。* /

P3 = p[2];/*生存因素,猎物。* /

P4 = p[3];/*死亡因素,猎物。* /

P5 = p[4];/*拥挤因素,猎物。* /

/* x[0]:捕食者物种。* /

/* x[1]:猎物种类。* /

Dx [0] = p1[0]*x[0]+p2[0]*x[1]*x[0];

dx [1] = p3 [0] * x [1] p4 [0] * x [0] * [1] p5[0] *战俘(x [1], 2);

注意,添加的延迟期计算为p5[0]*pow(x[1],2)。函数pow(。,.) is defined in the C library math.h, which must be included at the top of the model file through the statement #include "math.h" (this is not necessary in predprey1_c.c as it only holds standard C mathematics).

编译后的C mex文件,以及初始参数向量、适当的初始状态和一些管理信息,将作为输入参数提供给IDNLGREY对象构造函数:

文件名=“predprey2_c”;描述模型结构的文件。Order = [2 0 2];模型订单[ny nu nx]。参数= struct(“名字”, {“生存因素,捕食者”“死亡因素,捕食者”...“生存因素,猎物”“死亡因素,猎物”...“拥挤因素,猎物”},...“单位”, {1 /年的1 /年的1 /年的1 /年的1 /年的},...“价值”, {-1.1 0.9 1.1 0.9 0.2},...“最低”, {-Inf -Inf -Inf -Inf},...“最大”, {Inf Inf Inf Inf Inf},...“固定”,{假假假假假假假});估计所有5个参数。InitialStates = struct(“名字”, {“人口,捕食者”“人口,猎物”},...“单位”, {“大小(千)”“大小(千)”},...“价值”, {1.8 1.8},...“最低”, {0 0},...“最大”, {Inf Inf},...“固定”, {false false});估计两个初始状态。Ts = 0;时间连续系统。nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts,...“名字”“1个捕食者- 1个猎物系统显示拥挤”...“OutputName”, {“人口,捕食者”“人口,猎物”},...“OutputUnit”, {“大小(千)”“大小(千)”},...“TimeUnit”“年”);

通过输入模型对象的名称(nlgr),命令窗口中显示模型的基本信息。请注意,与前面一样,present(nlgr)提供了更全面的模型摘要。

nlgr
nlgr =由'predprey2_c' (MEX-file)定义的连续时间非线性灰盒模型:dx/dt = F(t, x(t), p1,…, p5) y(t) = H(t, x(t), p1,…,p5) + e(t) with 0 input(s), 2 state(s), 2 output(s), and 5 free parameter(s) (out of 5). Name: 1 predator - 1 prey system exhibiting crowding Status: Created by direct construction or transformation. Not estimated.

C.2。输入输出数据

接下来,我们加载(模拟的,尽管噪声受到了破坏)数据并创建一个IDDATA对象来描述这种拥挤类型的捕食者-被捕食者情况。该数据集包含201个数据样本,涵盖20年的进化。

负载(fullfile (matlabroot“工具箱”“识别”“iddemos”“数据”“predprey2data”));Z = iddata(y, [], 0.1,“名字”“1个捕食者- 1个猎物系统显示拥挤”);集(z,“OutputName”, {“人口,捕食者”“人口,猎物”},...“Tstart”0,“TimeUnit”“年”);

C.3。具有猎物拥挤的初始捕食者-被捕食者模型的表现

对初始模型的模拟清楚地表明,它不能应付真实的人口动态。见图。

CLF compare(z, nlgr, 1);

图中包含2个轴对象。坐标轴对象1包含2个line类型的对象。这些对象表示验证数据(种群,捕食者),1捕食者- 1捕食者系统显示拥挤:62.34%。坐标轴对象2包含两个line类型的对象。这些对象表示验证数据(种群,猎物),1个捕食者- 1个猎物系统显示拥挤:46.38%。

图6:初始捕食者-被捕食者模型的真实输出与模拟输出的比较。

C.4。参数估计

为了提高初始模型的性能,我们对5个未知参数和2个初始状态进行了估计。

nlgr = nlgreyest(z, nlgr, nlgreyestOptions)“显示”“上”));

C.5。具有猎物拥挤的捕食者-食饵估计模型的性能

参数和初始状态的估计值再次非常接近于用于生成真实输出数据的值:

disp (' True估计参数向量');
估计参数向量
Ptrue = [-1;1;1;1;0.1);流(' %6.3f %6.3f\n'[ptrue ';getpvec (nlgr) '));
-1.000 -1.000 1.001 1.000 1.002 1.000 1.002 0.100 0.101
disp (' ');

               
disp (“真实的估计初始状态”);
估计初始状态
X0true = [2;2);流(' %6.3f %6.3f\n'[x0true ';cell2mat (getinit (nlgr,“价值”))));
2.000 2.003 2.000 2.002

为了进一步评估模型的质量(并说明与初始模型相比的改进),我们还模拟了估计的模型。在一个图窗口中将模拟输出与真实输出进行比较。再次可以看出,估计的模型是相当好的。

比较(z, nlgr, 1);

图中包含2个轴对象。坐标轴对象1包含2个line类型的对象。这些对象表示验证数据(种群,捕食者),1捕食者- 1捕食者系统显示拥挤:97.53%。坐标轴对象2包含两个line类型的对象。这些对象表示验证数据(种群,猎物),1个捕食者- 1个猎物系统显示拥挤:97.36%。

图7:捕食者-食饵模型在被捕食者拥挤情况下的真实输出与模拟输出的比较。

我们通过PRESENT返回的模型信息来总结第三个总体示例。

礼物(nlgr);
nlgr =由'predprey2_c' (MEX-file)定义的连续时间非线性灰盒模型:dx/dt = F(t, x(t), p1,…, p5) y(t) = H(t, x(t), p1,…,p5) + e(t) with 0 input(s), 2 state(s), 2 output(s), and 5 free parameter(s) (out of 5). States: Initial value x(1) Population, predators(t) [Size (in t..] xinit@exp1 2.00281 (estimated) in [0, Inf] x(2) Population, preys(t) [Size (in t..] xinit@exp1 2.00224 (estimated) in [0, Inf] Outputs: y(1) Population, predators(t) [Size (in thousands)] y(2) Population, preys(t) [Size (in thousands)] Parameters: Value Standard Deviation p1 Survival factor, predators [1/year] -0.999914 0.00280581 (estimated) in [-Inf, Inf] p2 Death factor, predators [1/year] 1.00058 0.00276684 (estimated) in [-Inf, Inf] p3 Survival factor, preys [1/year] 1.0019 0.00272154 (estimated) in [-Inf, Inf] p4 Death factor, preys [1/year] 1.00224 0.00268423 (estimated) in [-Inf, Inf] p5 Crowding factor, preys [1/year] 0.101331 0.0005023 (estimated) in [-Inf, Inf] Name: 1 predator - 1 prey system exhibiting crowding Status: Termination condition: Change in cost was less than the specified tolerance.. Number of iterations: 8, Number of function evaluations: 9 Estimated using Solver: ode45; Search: lsqnonlin on time domain data "1 predator - 1 prey system exhibiting crowding". Fit to estimation data: [97.53;97.36]% FPE: 4.327e-08, MSE: 0.0004023 More information in model's "Report" property. Model Properties

结论

本例展示了如何基于MATLAB和MEX模型文件进行IDNLGREY时间序列建模。

Baidu
map