输入-输出数据集

在本例中，数据变量存储在twotankdata.mat文件将被使用。它包含使用0.2秒采样时间生成的一个双罐系统的3000个输入输出数据样本。输入u(t)是施加在泵上的电压[V]，它产生流入上油箱。在这个上罐底部有一个相当小的孔，流出物进入下罐，两个罐系统的输出y(t)就是下罐的液位[m]。我们创建一个IDDATA对象z封装加载的数据:

该数据集也用于非线性灰盒示例“双罐系统:时间连续SISO系统的C MEX-File建模”，其中假设了描述系统的物理方程。这个例子是关于黑盒模型的，因此没有假设实际物理定律的知识。

负载twotankdataZ = iddata(y, u, 0.2，“名字”,“双坦克系统”）;

分割数据和绘图

将这个数据集分成3个相同大小的子集，每个子集包含1000个样本。

我们将用z1，并对其进行评估z2而且z3．z1类似于z2，但不是z3．因此在某种意义上，z2可以帮助我们评估模型的插值能力，而z3能帮我们评估他们的推断能力。

Z1 = z(1:1000);Z2 = z(1001:2000);Z3 = z(2001:3000);情节(z1、z2 z3)传说(“估计”,“验证1”,《验证2》）

回归量选择

估计非线性黑盒模型的第一步是选择回归集。回归量是基于时滞输入/输出变量的简单公式，最简单的情况是变量被一小组连续值滞后。例如，如果“u”是一个输入变量的名称，“y”是一个输出变量的名称，那么一个示例回归集可以是{y(t-1)， y(t-2)， u(t)， u(t-1)， u(t-2)}，其中“t”表示时间变量。另一个涉及多项式项的例子可能是{abs(y(t-2))， y(t-2)*u(t-1)， u(t-4)^2}。延迟变量中更复杂、任意的公式也是可能的。

在系统标识工具箱™中，不需要显式地创建回归器集;只需指定它们的形式(如多项式顺序)、贡献变量及其滞后，就可以隐式生成一个大型集合。规范的对象linearRegressor,polynomialRegressor而且customRegressor都用于此目的。例如R = linearRegressor ({' y ', ' u '},{[1 2 4],[2] 10})意味着回归量集{y y (t - 1), (2), y(第四节),u (2), u (t-10)}。类似地，R =多项式回归({'y'}，{[1 2]}，3)表示变量'y'中滞后1和2的三阶多项式回归量，即集合{y(t-1)^3, y(t-2)^3}。还有更多的配置选择，例如只使用绝对值，和/或允许在回归公式中混合滞后。

使用顺序矩阵的线性回归规范

当模型仅使用基于连续滞后的线性回归器时，通常更容易使用顺序矩阵来指定它们，而不是使用linearRegressor对象。序矩阵[na nb nk]定义了回归公式中使用的过去输出、过去输入和输入延迟的数量。这与估计线性ARX模型时使用的思想是相同的(参见ARX, IDPOLY)。例如，NN =[2 3 4]意味着输出变量使用滞后(1,2)，而输入变量使用滞后(4,5,6)，导致回归集{y(t-1)， y(t-2)， u(t-4)， u(t-5)， u(t-6)}。

通常模型顺序是通过试错来选择的。有时，首先尝试线性ARX模型是有用的，以便获得关于最佳顺序的提示。因此，让我们首先尝试确定线性ARX模型的最佳阶数，使用函数arxstruc而且selstruc．

V = arxstruc(z1,z2,struc(1:5, 1:5,1:5));%用赤池信息准则(AIC)选择最佳顺序nn = selstruc(V，“另类投资会议”）

Nn = 5 13 3

AIC准则选择nn = [na nb nk] =[5 1 3]，这意味着，在所选的ARX模型结构中，y(t)由6个回归量y(t-1)，y(t-1)，…，y(t-5)和u(t-3)预测。我们将尝试在估计非线性模型时使用这些值。

非线性ARX (IDNLARX)模型

在IDNLARX模型中，模型输出是回归函数的非线性函数，如y(t) = Fcn(y(t-1)，y(t-1)，…，y(t-5)， u(t-3))。

IDNLARX模型通常使用以下语法进行估计:

M = NLARX(数据，订单，Fcn)。

它类似于线性ARX命令，附加参数“Fcn”指定非线性函数的类型。非线性ARX模型通过使用两阶段转换产生输出:将训练数据(输入-输出信号)映射到一组回归函数。从数值上讲，回归集是一个双矩阵R，每个回归数有一列数据。2.将回归函数映射到单个输出y = Fcn®，其中Fcn()是一个标量非线性函数。

为了估计一个IDNLARX模型，在选择了模型阶数之后，我们需要选择使用的非线性映射函数Fcn()的类型。让我们首先尝试一个小波网络(idWaveletNetwork)函数。

mw1 = nlarx(z1，[5 1 3]， idWaveletNetwork);

小波网络使用小波单元的组合来将回归量映射到模型的输出。该模型mw1是一个@idnalrx对象。它将非线性映射函数(此处为idWaveletNetwork对象)存储在其OutputFcn财产。使用这些单位的数量可以提前指定，或者我们可以让估计算法自动确定一个最优值。属性NonlinearFcn。NumberOfUnits存储模型中选择的单元数量。

NLFcn = mw1.OutputFcn;NLFcn.NonlinearFcn.NumberOfUnits

Ans = 3

通过比较模型的响应来检查结果mw1到数据集中的测量输出z1．

比较(z1, mw1);

使用模型属性

第一个模型不太令人满意，因此让我们尝试修改模型的一些属性。与其让估计算法自动选择idWaveletNetwork函数中的单元数，不如显式地指定这个数字。

mw2 = nlarx(z1，[5 1 3]， idWaveletNetwork(8));

默认InputName和OutputName已被使用。

对象。InputName对象。OutputName

Ans = 1×1单元格数组{'u1'} Ans = 1×1单元格数组{'y1'}

IDNLARX模型的回归函数可以使用GETREG命令查看。回归量是由对象。InputName,对象。OutputName还有时间延迟。

getreg(对象)

ans = 6×1单元阵列{y1 (t - 1)的}{‘y1 (2)} {y1(条t - 3)的}{y1(第四节)的}{‘y1 (t-5)} {u1(条t - 3)的}

注意顺序矩阵等价于使用specification对象创建的线性回归集，如下所示:

R =线性回归([z1.OutputName;z1。InputName] {1:5, 3});mw2_a = nlarx(z1, R, idWaveletNetwork(8));getreg (mw2_a)

ans = 6×1单元阵列{y1 (t - 1)的}{‘y1 (2)} {y1(条t - 3)的}{y1(第四节)的}{‘y1 (t-5)} {u1(条t - 3)的}

mw2_a本质上是相同的模型对象．规范对象的使用允许更大的灵活性;当您不想使用连续延迟，或者输出变量的最小延迟大于1时，可以使用它。

将回归函数分配给映射函数

IDNLARX模型的输出是其回归函数的静态函数。通常这个映射函数有一个线性块，一个非线性块，外加一个偏置项(输出偏置)。模型输出是这些块的输出和偏置的总和。为了降低模型复杂度，只选择一个回归量子集进入非线性块和线性块。映射函数的线性和非线性组件的回归函数的分配是使用模型的RegressorUsage属性执行的。

RegUseTable = mw2。RegressorUsage

RegUseTable = 6×2 table y1:LinearFcn y1:NonlinearFcn ____________ _______________ y1(t-1)真真y1(t-2)真真y1(t-3)真真y1(t-4)真真y1(t-5)真真u1(t-3)真真真

RegUseTable是一个表，它显示了哪些回归器被用作小波网络的线性和非线性组件的输入。每一行对应一个回归量。假设我们只想使用非线性分量中输入变量u1的回归函数。可以通过以下方式进行配置。

RegNames = regustable . properties . rownames;I = contains(RegNames，‘u1’）;RegUseTable。（“日元:NonlinearFcn”)(~I) = false;

作为一个例子，考虑一个模型，除了使用的线性回归函数外，还使用了2阶多项式回归函数对象．首先，生成一个多项式回归规范集。

P =多项式回归({“日元”,‘u1’0:3}}, {1:2, 2)

P = Order 2回归变量y1, u1 Order: 2变量:{'y1' 'u1'} lag: {[1 2] [0 1 2 3]} UseAbsolute: [0 0 0] AllowVariableMix: 0 AllowLagMix: 0 TimeVariable: 't'

现在添加规格P到模型的Regressors属性。

Mw3 = mw2;%从mw2结构开始进行新的估计mw3。Regressors = [mw3.Regressors;P];将多项式集添加到模型结构中

最后，通过最小化提前1步的预测误差，更新模型参数以拟合数据。

Mw3 = nlarx(z1, Mw3)

mw3 = 非线性ARX模型，1输出1输入输入:u1输出:y1回归量:1。变量y1 u1 2中的线性回归。变量y1, u1中的二阶回归函数输出函数:8单元小波网络采样时间:0.2秒状态:在时域数据“两罐系统”上使用NLARX估计。拟合估计数据:96.76%(预测焦点)FPE: 3.499e-05, MSE: 3.338e-05

评估估计模型

不同的模型，包括线性模型和非线性模型，都可以在同一个COMPARE命令中进行比较。

Mlin = arx(z1，[5 1 3]);%线性ARX模型％对象,比较(z1, mlin mw3)

可以通过在验证数据集上运行COMPARE命令来执行模型验证。首先，通过检查它们的模拟响应与测量输出信号的匹配程度来验证模型z2．

对象,比较(z2 mlin mw3)

类似地，针对数据集z3验证估计的模型。

对象,比较(z3 mlin mw3)

还可以对模型残差进行分析，以验证模型质量。我们期望残差是白色的(除滞后0外与自身不相关)并且与输入信号不相关。残差是超前一步的预测误差。通过使用RESID命令，可以用无关紧要值的边界来显示残差。例如，在验证数据集上z2时，三种已识别模型的残差如下所示。

对象,渣油(z2 mlin mw3)传奇显示

残差都与输入信号不相关z2(z2.InputData)。然而，它们对模型表现出一些自相关性mlin而且mw3．该模型对象在两项(自相关和相互关联)测试中表现更好。

函数图可以用来查看各种IDNLARX模型的非线性响应。该图本质上显示了非线性映射函数Fcn()的特征。

mw3情节(对象)绘制非线性响应作为所选回归量的函数

基于Sigmoid网络的非线性ARX模型

该模型可以使用其他非线性函数来代替小波网络。尝试idSigmoidNetwork函数，它使用sigmoid单位的和(加上一个线性项和一个偏移项)将回归函数映射到输出。将网络配置为使用8个单位的sigmoid(按不同数量进行扩展和转换)。

ms1 = nlarx(z1，[5 1 3]， idSigmoidNetwork(8));比较(z1, ms1)

现在在数据集z2上计算新模型。

比较(z2, ms1);

在数据集z3上。

比较(ms1 z3);

IDNLARX模型中的其他非线性映射函数

idCustomNetwork类似于idSigmoidNetwork，但是用户可以在MATLAB文件中定义自己的单位函数，而不是sigmoid单位函数。本例使用gaussunit函数。

类型gaussunit.m

函数(f, g) = gaussunit (x) % gaussunit单位函数用于idCustomNetwork(例子)% % (f, g) = gaussunit (x) % % x:单位函数变量% f:单位函数值% g: df / dx %:单位活跃的范围(g (x)显著不为零的区间[-])% %单位函数必须“矢量”:对于一个向量或矩阵x,输出% f和g参数必须具有相同的大小x,计算中的元素。作者:张清华f = exp(-x.*x);如果nargout>1 g = - 2*x .* f;A = 0.2;结束

对于基于定制网的模型，只使用第三、第五和第六个回归量作为定制网函数非线性成分的输入。

Mc1_ini = idnlarx(“日元”,‘u1’， [5 1 3]， idCustomNetwork(@gaussunit));使用= false(6,1);使用([3 5 6])= true;mc1_ini。RegressorUsage{:，2} =使用;MC1 = nlarx(z1,mc1_ini);

树配分函数是另一个可以用于非线性ARX模型的映射函数。使用idTreePartition函数估计一个模型:

mt1 = nlarx(z1， [5 1 3]， idTreePartition);

比较模型的响应哪而且mt1数据z1．

比较(z1, MC1, mt1)

使用回归算法从统计和机器学习工具箱

如果您可以访问统计和机器学习工具箱™，您还可以使用高斯过程(GP)函数和袋装或增强树集成来创建您的非线性ARX模型。GPs由idGaussianProcess对象表示，而回归树集合由idTreeEnsemble对象表示。GPs默认使用平方指数核。在创建对象时可以指定其他内核类型。在本例中，我们使用'Matern32'内核。

如果存在(“fitrgp”,“文件”(2)警告=警告(“关闭”,“统计数据:套索:MaxIterReached”）;mp1 = nlarx(z1， [5 1 3]， idGaussianProcess(“Matern32”）)创建200个回归树的增强集合En = idTreeEnsemble;En.EstimationOptions.FitMethod =“lsboost-resampled”；En.EstimationOptions.NumLearningCycles = 200;En.EstimationOptions.ResampleFraction = 1;mp2 = nlarx(z1， [0 70 3]， En)比较(z1, mp1, mp2)%比较模型拟合估计数据警告(警告)结束

mp1 =  1输出1输入的非线性ARX模型输入:u1输出:y1回归:变量y1, u1中的线性回归输出函数:使用Matern32核的高斯过程函数。采样时间:0.2秒状态:在时域数据上使用NLARX估计“两罐系统”。拟合估计数据:97.1%(预测焦点)FPE: 2.738e-05, MSE: 2.675e-05 mp2 =  1输出1输入的非线性ARX模型输入:u1输出:y1回归函数:变量u1中的线性回归函数输出函数:袋装回归树集合采样时间:0.2秒状态:在时域数据上使用NLARX估计“双坦克系统”。拟合估计数据:91.14% MSE: 0.0002503

使用深度学习工具箱中的网络对象

如果您有深度学习工具箱™可用，还可以使用神经网络来定义IDNLARX模型的映射函数Fcn()。这个网络必须是前馈网络(没有循环组件)。

在这里，我们将创建一个输入数量未知的单输出网络(由网络对象中的输入大小为零表示)。在估计过程中，输入的数量被自动设置为回归量的数量。

如果存在(“feedforwardnet”,“文件”) = = 2%仅在“深度学习工具箱”可用时运行Ff = feedforwardnet([5 20]);%使用前馈网络将回归量映射到输出ff.layers{2}。transferFcn =“radbas”；ff.trainParam.epochs = 50;创建一个包裹前馈网络的神经网络映射函数。OutputFcn = idFeedforwardNetwork(ff);mn1 = nlarx(z1，[5 1 3]， OutputFcn);估计网络参数%比较(z1, mn1)%比较适合估计数据结束

Hammerstein-Wiener (IDNLHW)模型

Hammerstein-Wiener模型最多由3个块组成:一个线性传递函数块前面是一个非线性静态块和/或后面是另一个非线性静态块。如果第一个非线性静态块不存在，则称为Wiener模型，如果第二个非线性静态块不存在，则称为Hammerstein模型。

IDNLHW模型通常使用以下语法进行估计:

M = NLHW(数据，订单，输入非线性，输出非线性)。

其中Orders = [nb bf nk]为线性传递函数的阶数和时延。input非线性和output非线性在线性块的输入和输出端为两个非线性块指定了非线性函数。M是一个@idnlhw模型。线性输出误差(OE)模型对应于平凡(单位增益)非线性的情况，即如果IDNLHW模型的输入和输出非线性函数都是单位增益，则模型结构等效于线性OE模型。

分段线性非线性函数的Hammerstein-Wiener模型

idpiecewislinear(分段线性)非线性函数在一般IDNLHW模型中是有用的。

注意，在Orders = [nb nf nk]中，nf指定线性传递函数的极点数，与IDNLARX模型的“na”阶有一定关系。“nb”与数字if为零有关。

mhw1 = nlhw(z1， [1 5 3]， idpiecewislinear, idpiecewislinear);

可以像以前一样使用COMPARE来计算结果。

比较(z1, mhw1)

模型属性可以通过PLOT命令可视化。

单击方框图，选择输入非线性(UNL)、线性块或输出非线性(YNL)视图。

对于线性块，可以在阶跃响应、波德图、脉冲响应和极点零映射中选择图的类型。

情节(mhw1)

这两个分段线性函数的断点可以被检验。这是一个两行矩阵，第一行对应idpiecewislinear函数的输入，第二行对应idpiecewislinear函数的输出。

mhw1.InputNonlinearity.BreakPoints

ans =列1至7 2.6795 3.4262 4.1729 4.9195 5.6324 6.3599 7.0874 -2.0650 -1.6316 -0.2396 1.4602 2.7554 3.6889 4.2632列8至10 7.8148 8.3543 9.1260 4.5074 4.4653 6.2634

含饱和和死区非线性的Hammerstein-Wiener模型

饱和函数和死区函数是物理学启发的非线性函数。它们可用于描述执行器/传感器饱和场景和干摩擦效应。

mhw2 = nlhw(z1， [1 5 3]， idDeadZone, idSaturation);比较(z1, mhw2)

非线性的缺失由IDUNITGAIN对象或空的“[]”值表示。单位增益只是一个输入信号到输出的馈电，没有任何修改。

mhw3 = nlhw(z1， [1 5 3]， idDeadZone, idUnitGain);%无输出非线性mhw4 = nlhw(z1， [1 5 3]， []，idSaturation);%无输入非线性

估计后可以检验死区和饱和函数的极限值。

mhw2.InputNonlinearity。ZeroInterval mhw2.OutputNonlinearity.LinearInterval

Ans = -0.9982 4.7565 Ans = 0.2601 0.5848

在估计模型时，可以指定死区零区间(ZeroInterval)或饱和函数线性区间(LinearInterval)的初始猜测。

mhw5 = nlhw(z1， [1 5 3]， idDeadZone([0.1 0.2])， idSaturation([-1 1])));

Hammerstein-Wiener模型-指定属性

估计算法选项可以通过NLHWOPTIONS命令指定。

opt = nlhwOptions();opt.SearchMethod =“玲娜”；opt.SearchOptions.MaxIterations = 7;mhw6 = nlhw(z1， [1 5 3]， idDeadZone, idSaturation, opt);

一个已经估计的模型可以通过更多的估计迭代来改进。

根据估计数据进行评估z1,mhw7应该有比mhw6．

opt.SearchOptions.MaxIterations = 30;Mhw7 = nlhw(z1, mhw6, opt);比较(z1, mhw6, mhw7)

Hammerstein-Wiener模型-使用其他非线性函数

非线性函数分段线性(idpiecewislinear)、饱和(idSaturation)和死区(idDeadZone)主要设计用于IDNLHW模型。它们只能估计单变量非线性。除了IDNLARX模型外，其他更通用的非线性函数SigmoidNetwork (idSigmoidNetwork)、CustomNetwork (idCustomNetwork)和WaveletNetwork (idWaveletNetwork)也可以用于IDNLHW模型。然而，由于IDNLHW模型的估计需要可微非线性函数，因此不能使用不可微函数Tree-partition (idTreePartition)、多层神经网络(idFeedforwardNetwork)、高斯过程函数(idGaussianProcess)和回归树集成(idTreeEnsemble)。