A双槽系统

目标是模拟一个实验室规模的两罐系统的下罐的液位，示意图如图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]为t时刻储罐i的流入量和流出量。

对于上部罐，假设流入量与施加到泵上的电压成正比，即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模型文件

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

Dx:状态空间方程(s)的右边(s) y:输出方程(s)的右边(s)

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

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

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

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

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

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

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

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

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

/*指定输出的数量。*/ #define NY 1

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

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

和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库;否则可以省略。参见“非线性灰盒模型识别教程:创建IDNLGREY模型文件”了解C数学库的更多细节。

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

Void compute_y(double *y, double *x)

调用compute_y为:

compute_y (y、x);

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

一旦模型源文件完成，必须对其进行编译，这可以在MATLAB命令提示符中使用mex命令完成;看到“帮助墨西哥人”。(此处省略此步骤。)

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

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

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

创建两个坦克IDNLGREY模型对象

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

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

我们继续通过命令setit和SETPAR添加关于状态的名称和单元以及模型参数的信息。此外，状态x1(t)和x2(t)都是不能为负的坦克级别，因此我们也通过“Minimum”属性指定x1(0)和x2(0) >= 0。事实上，我们也知道所有的模型参数应该是严格正的。因此，我们将所有参数的“最小值”属性设置为一个小的正值(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命令时，默认情况下将估计所有初始状态，而不管任何初始状态是否被定义为“Fixed”。为了只估计自由初始状态(s)，用第三个和第四个输入参数调用COMPARE，如下所示:由于坦克模型的初始状态默认都是“固定的”，这个命令不会执行初始状态估计。

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

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

模拟的和真实的输出显示在绘图窗口中，可以看到，匹配不是那么令人印象深刻。

参数估计

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

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

估计两坦克模型的性能

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

比较(z, nlgr);

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

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

nk =延迟(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 = ss (z);%状态空间模型(顺序自动确定)比较(z, nlgr“b”arx22,“m -”arx33,”男:“arx44,“m——”，．..oe22,“g -”oe33,“旅客:”oe44,“g——”sslin,的r -）;

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

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

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

估计参数矢量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 V增加到6 V、7 V、8 V、9 V和10 V，会发生什么情况。我们通过调用从固定偏移量5伏特开始的不同指定阶跃振幅的STEP来实现这一点。由。创建的专用选项集简化了阶跃响应配置stepDataOptions：

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

图8:得到了估计的IDNLGREY两罐模型的阶跃响应。

最后使用PRESENT命令，得到估计模型的概要信息:

礼物(nlgr);

nlgr =由'twotanks_c'定义的连续时间非线性灰箱模型(MEX-file): dx/dt = F(t, u(t)， x(t)， p1，…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.

结论

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

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

额外的信息

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