双罐系统

目的是对实验室规模的双罐系统的下罐液位进行建模，如图1所示。

图1:双坦克系统的示意图。

输入输出数据

我们通过加载可用的输入输出数据开始建模工作，这些数据使用下面的IDNLGREY模型结构进行模拟，并向输出添加噪声。twotankdata。mat文件包含一个包含3000个输入输出样本的数据集，使用0.2秒的采样率(Ts)生成。输入u(t)是施加在泵上的电压[V]，它产生了流入上部油箱的流量。在这个上罐底部有一个相当小的洞，流出的水进入下罐，两个罐系统的输出y(t)就是下罐的液位[m]。我们创建一个IDDATA对象z来保存坦克数据。为了记账和记录，我们还指定了通道名称和单元。此步骤是可选的。

负载(fullfile (matlabroot“工具箱”，“识别”，“iddemos”，“数据”，“twotankdata”));Z = iddata(y, u, 0.2，“名字”，“两个坦克”）;集(z,“InputName”，“泵电压”，“InputUnit”，“V”，.．.“OutputName”，“降低水箱水位”，“OutputUnit”，“米”，.．.“Tstart”0,“TimeUnit”，“年代”）;

用于估计的输入输出数据显示在一个图窗口中。

图(“名字”[z。的名字':输入输出数据']);情节(z);

图2:两罐系统的输入输出数据。

双坦克系统建模

下一步是指定一个描述两个坦克系统的模型结构。为此，设x1(t)和x2(t)分别表示上水箱和下水箱的水位。对于每个水箱，基本物理(质量平衡)表明，水量的变化取决于流入和流出的差(i = 1,2):

d/dt (Ai*xi(t)) = Qini(t) - Qouti(t)

其中Ai [m^2]为槽i的截面积，Qini(t)和Qouti(t) [m^3/s]为槽i在t时刻的流入和流出量。

对于上槽，假定流入与施加在泵上的电压成正比，即Qin1(t) = k*u(t)。由于上部水箱的出水口较小，可以应用伯努利定律，即流出量与水位的平方根成正比，更准确地说为:

Qout1(t) = a1*√(2*g*x1(t))

其中a1为出口孔截面积，g为重力常数。对于下槽，流入等于上槽流出量，即Qin2(t) = Qout1(t)，流出量由伯努利定律给出:

Qout2(t) = a2*√(2*g*x2(t))

其中a2为出口孔的截面积。

综合这些事实，可以得出以下状态空间结构:

d / dt x1 (t) = 1 / A1 * (k * u (t) - A1 * sqrt (2 * g * x1 (t))) d / dt x2 (t) = 1 / A2 * (A1 * 12 + (2 * g * x1 (t)) - A2 * 12 + (2 * g * x2 (t))) y (t) = x2 (t)

两罐C MEX模型文件

这些方程接下来被放入一个带有6个参数(或常数)的C mex文件，A1, k, A1, g, A2和A2。C MEX-文件通常比用MATLAB语言编写的相应文件要复杂一些，但是C MEX建模通常在执行速度方面具有明显的优势，特别是对于更复杂的模型。提供了一个C mex模板文件(见下文)来帮助用户构造代码。对于大多数应用程序，定义输出的数量并在模板中输入描述dx和y的代码行就足够了。一个IDNLGREY C mex文件应该总是被构造成返回两个输出:

Dx:状态空间方程的右手边y:输出方程的右手边

并且它应该接受3+Npo(+1)个输入参数，如下所示:

T:当前时间x: T时刻的状态向量([]对于静态模型)u: T时刻的输入向量([]对于时间序列模型)p1, p2，…， pNpo:单个参数(可以是实标量、列向量或二维矩阵);这里Npo是参数对象的数量，对于带有标量参数的模型，它与参数Np FileArgument的数量相一致:模型文件的可选输入

在我们的两个坦克系统中有6个标量参数，因此C MEX建模文件的输入参数的数量应该是3+Npo = 3+6 = 9。后面的第10:th参数可以省略，因为在这个应用程序中没有使用可选的FileArgument。

编写C MEX建模文件通常分为四个步骤:

1.包含c库和输出数量的定义。2.写出计算状态方程右边的函数compute_dx。3.写出计算输出方程的右边的函数compute_y。4.编写主接口函数，其中包括基本的错误检查功能、创建和处理输入和输出参数的代码，以及对compute_dx和compute_y的调用。

让我们查看双坦克系统的C MEX源文件(除了一些注释)，并在此基础上更详细地讨论这四个项目。

图3:cmex两罐系统的源代码。

1.通常包含两个c库mex.h和math.h，以提供对许多与mex相关的和数学函数的访问。每个建模文件使用标准的C-define声明输出的数量:

/*包含库。*/ #include "mex.h" #include "math.h"

/*在这里指定输出的数量。*/ #定义NY

2 - 3。接下来，我们在文件中找到更新状态的函数compute_dx和输出compute_y。这两个函数都包含参数列表，在位置1处计算输出(dx或y)，之后依次是计算状态方程的右侧和输出方程所需的所有变量和参数。

这些函数的第一步是解压缩将在后续方程中使用的模型参数。任何有效的变量名(输入实参列表中使用的变量名除外)都可以用于提供各个参数的物理意义名称。

和C语言一样，数组的第一个元素存储在0号位置。因此，C中的dx[0]对应MATLAB®中的dx(或在它是标量的情况下只是dx)，输入u[0]对应u(或u(1))，参数A1[0]对应A1，以此类推。

两个坦克模型文件涉及平方根计算。这是通过包含数学C库math.h实现的。数学库实现了最常见的三角函数(sin, cos, tan, asin, acos, atan等)，指数函数(exp)和对数函数(log, log10)，平方根函数(sqrt)和函数幂函数(pow)，以及绝对值计算(fabs)。每当使用任何math.h函数时，必须包含math.h库;否则可以省略。有关C数学库的更多细节，请参见“非线性灰盒模型识别教程:创建IDNLGREY模型文件”。

4.主界面函数应该几乎总是具有相同的内容，对于大多数应用程序来说，不需要进行任何修改。原则上，唯一可能考虑更改的部分是调用compute_dx和compute_y的地方。对于静态系统，可以省略对compute_dx的调用。在其他情况下，可能希望只传递状态和输出方程中引用的变量和参数。例如，在双罐系统的输出方程中，只使用了一种状态，我们可以很好地将输入参数列表缩短为:

Void compute_y(double *y, double *x)

并调用compute_y as:

compute_y (y、x);

还可以扩展compute_dx和compute_y的输入参数列表，以包括在接口函数中推断出的进一步变量，如状态的数量和参数的数量。

一旦模型源文件完成，它必须被编译，这可以从MATLAB命令提示符使用mex命令完成;参见“help mex”。(这里省略了这一步。)

在开发特定于模型的C MEX文件时，通过复制IDNLGREY C MEX模板文件开始工作通常是有用的。该模板包含基本源代码以及关于如何为特定应用程序定制代码的详细说明。在MATLAB命令提示符中输入以下命令可以显示模板文件的位置。

fullfile(matlabroot， 'toolbox'， 'ident'， ' nlient '， 'IDNLGREY_MODEL_TEMPLATE.c')

关于IDNLGREY C MEX模型文件的更多细节，请参见“创建IDNLGREY模型文件”示例。

创建两个坦克IDNLGREY模型对象

下一步是创建一个IDNLGREY对象来描述两个坦克系统。为了方便起见，我们还设置了一些关于输入和输出(名称和单位)的记账信息。

文件名=“twotanks_c”；描述模型结构的文件。Order = [1 1 2];模型订单[ny nu nx]。参数= {0.5;0.0035;0.019;.．.9.81;0.25;0.016};初始参数。InitialStates = [0;0.1);初始状态的初始值。Ts = 0;时间连续系统。nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts，.．.“名字”，“两个坦克”）;集(nlgr,“InputName”，“泵电压”，“InputUnit”，“V”，.．.“OutputName”，“降低水箱水位”，“OutputUnit”，“米”，.．..．.“TimeUnit”，“年代”）;

我们继续通过命令SETINIT和SETPAR添加关于状态的名称、单位和模型参数的信息。此外，状态x1(t)和x2(t)都是不能为负的槽位，因此我们还通过“最小值”属性指定x1(0)和x2(0) >= 0。事实上，我们也知道所有的模型参数都应该是严格正的。因此，我们将所有参数的“Minimum”属性设置为一个小的正数值(eps(0))。这些设置意味着约束估计将在接下来的估计步骤中执行(也就是说，估计的模型将是一个所有输入的约束都被遵守的模型)。

NLGR = setinit(NLGR，“名字”, {“水箱上水位”“降低水箱水位”});NLGR = setinit(NLGR，“单位”, {“米”“米”});NLGR = setinit(NLGR，“最低”， {0 0});%积极的水平!NLGR = setpar(NLGR，“名字”, {“上槽区”.．.“泵恒”.．.“上槽出口区域”.．.“引力常数”.．.“下油罐区”.．.“下油罐出口区域”});NLGR = setpar(NLGR，“单位”, {“m ^ 2”' m ^ 3 / (s * V) '“m ^ 2”“m / s ^ 2)”“m ^ 2”“m ^ 2”});NLGR = setpar(NLGR，“最低”, num2cell (eps (0) * (6,1)));所有参数> 0!

两个罐的截面积(A1和A2)可以相当准确地确定。因此，我们将这些和g视为常量，并验证通过命令GETPAR为所有6个参数正确设置了“Fixed”字段。总而言之，这意味着3个模型参数将被估计。

nlgr.Parameters(1)。固定= true;nlgr.Parameters(4)。固定= true;nlgr.Parameters(5)。固定= true;getpar (nlgr“固定”）

ans x1 = 6单元阵列{[1]}{[0]}{[0]}{[1]}{[1]}{[0]}

初始双坦克模型的性能

在估计自由参数k, a1和a2之前，我们用初始参数值模拟系统。我们使用默认的微分方程求解器(带有自适应步长调整的龙格-库塔45求解器)，并将绝对和相对误差容差设置为相当小的值(分别为1e-6和1e-5)。注意，当使用两个输入参数调用COMPARE命令时，因为default将估计所有初始状态，而不管是否有初始状态被定义为“Fixed”。为了只估计空闲的初始状态，调用COMPARE与第三和第四个输入参数，如下所示:COMPARE (z, nlgr， 'init'， 'm');由于坦克模型的初始状态默认为“固定”，因此此命令不会执行初始状态估计。

nlgr.SimulationOptions.AbsTol = 1e-6;nlgr.SimulationOptions.RelTol = 1e-5;比较(z, nlgr);

图4:初始两坦克模型的真实输出与模拟输出的比较。

模拟和真实的输出显示在一个图窗口，可以看到，适合不是那么令人印象深刻。

参数估计

为了提高拟合性，接下来用NLGREYEST对3个自由参数进行估计。(因为默认情况下，所有初始状态的'Fixed'字段都是假的，所以在对估计器的调用中不会对初始状态进行估计。)

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

估计双坦克模型的性能

为了研究估计模型的性能，对它进行了模拟(这里重新估计了初始状态)。

比较(z, nlgr);

图5:两种坦克模型的真实输出与模拟输出的比较。

计算结果与模拟结果吻合较好。然而，剩下的一个问题是，是否可以使用一个更简单的线性模型结构来精确地描述两个坦克系统。为了回答这个问题，让我们试着将数据与一些标准线性模型结构相匹配，然后使用COMPARE来查看这些模型如何很好地捕捉坦克的动力学。

Nk = delayest(z);Arx22 = arx(z， [2 2 nk]);二阶线性ARX模型。Arx33 = arx(z， [3 3 nk]);三阶线性ARX模型。Arx44 = arx(z， [4 4 nk]);四阶线性ARX模型。Oe22 = oe(z， [2 2 nk]);二阶线性OE模型。Oe33 = oe(z， [3 3 nk]);三阶线性OE模型。Oe44 = oe(z， [4 4 nk]);四阶线性OE模型。Sslin = ssest(z);%状态空间模型(自动确定顺序)比较(z, nlgr“b”arx22,“m -”arx33,”男:“arx44,“m——”，.．.oe22,“g -”oe33,“旅客:”oe44,“g——”sslin,的r -）;

图6:两种坦克模型的真实输出和模拟输出的比较。

对比图清楚地表明，线性模型不能反映两个坦克系统的所有动力学。另一方面，估计的非线性IDNLGREY模型显示出与真实输出的良好拟合。此外，IDNLGREY模型参数也很好地符合用于生成真实输出的参数。在接下来的显示计算中，我们使用GETPVEC命令，该命令返回从包含IDNLGREY对象的模型参数的结构数组创建的参数向量。

disp (' True估计参数向量'）;Ptrue = [0.5;0.005;0.02;9.81;0.25;0.015);流(' %1.4f %1.4f\n'[ptrue ';getpvec (nlgr) '));

True估计参数向量0.5000 0.5000 0.0050 0.0049 0.0200 0.0200 9.8100 9.8100 0.2500 0.2500 0.0150 0.0147

使用PE得到的预测误差很小，看起来非常像随机噪声。

图;体育(z, nlgr);

图7:得到了估计的IDNLGREY双坦克模型的预测误差。

让我们也研究一下，如果输入电压逐步从5增加到6、7、8、9和10 V会发生什么。我们通过从5伏特的固定偏移量开始调用不同的指定步长振幅的STEP来做到这一点。创建的专用选项集简化了步骤响应配置stepDataOptions：

图(“名字”, [nlgr。的名字':步骤响应']);T = (-20:0.1:80)';Opt = stepDataOptions(“InputOffset”5,“StepAmplitude”6);步骤(nlgr t“b”、选择);持有在step振幅= 7;步骤(nlgr t‘g’、选择);step振幅= 8;步骤(nlgr t“r”、选择);step振幅= 9;步骤(nlgr t“米”、选择);step振幅= 10;步骤(nlgr t“k”、选择);网格在；传奇('5 -> 6v '，'5 -> 7v '，'5 -> 8 v '，'5 -> 9v '，'5 -> 10v '，.．.“位置”，“最佳”）;

图8:得到了估计的IDNLGREY双坦克模型的阶跃响应。

最后使用PRESENT命令，我们得到关于估计模型的摘要信息:

礼物(nlgr);

nlgr =由'twotanks_c' (MEX-file)定义的连续时间非线性灰盒模型:dx/dt = F(t, u(t)， x(t)， p1，…， p6) y(t) = H(t, u(t)， x(t)， p1，…，p6) + e(t) with 1 input(s), 2 state(s), 1 output(s), and 3 free parameter(s) (out of 6). Inputs: u(1) Pump voltage(t) [V] States: Initial value x(1) Upper tank water level(t) [m] xinit@exp1 0 (fixed) in [0, Inf] x(2) Lower tank water level(t) [m] xinit@exp1 0.1 (fixed) in [0, Inf] Outputs: y(1) Lower tank water level(t) [m] Parameters: Value Standard Deviation p1 Upper tank area [m^2] 0.5 0 (fixed) in ]0, Inf] p2 Pump constant [m^3/(s*V)] 0.00488584 0.0259032 (estimated) in ]0, Inf] p3 Upper tank outlet area [m^2] 0.0199719 0.0064682 (estimated) in ]0, Inf] p4 Gravity constant [m/(s^2)] 9.81 0 (fixed) in ]0, Inf] p5 Lower tank area [m^2] 0.25 0 (fixed) in ]0, Inf] p6 Lower tank outlet area [m^2] 0.0146546 0.0776058 (estimated) in ]0, Inf] Name: Two tanks 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 "Two tanks". Fit to estimation data: 97.35% FPE: 2.419e-05, MSE: 2.414e-05 More information in model's "Report" property. Model Properties

结论

在这个例子中，我们展示了:

1.如何使用C mex文件进行IDNLGREY建模;提供了一个相当简单的示例，其中非线性状态空间建模显示出良好的潜力

额外的信息

有关使用系统识别工具箱™识别动态系统的更多信息，请访问系统识别工具箱产品信息页面。