使用命令行识别线性模型- MATLAB & Simulink - 卡塔尔世界杯8强比赛直播

使用命令行识别线性模型

介绍

目标

从多输入/单输出(MISO)数据中估计和验证线性模型，以找到最能描述系统动力学的模型。

完成本教程后，您将能够使用命令行完成以下任务:

创建数据对象来表示数据。
图数据。
通过从输入和输出信号中去除偏移量来处理数据。
从数据中估计和验证线性模型。
模拟和预测模型输出。

请注意

本教程使用时域数据来演示如何估计线性模型。同样的工作流程也适用于拟合频域数据。

数据描述

本教程使用的是数据文件co2data.mat，其中包含两个实验的双输入和单输出(MISO)时域数据来自稳态，操作员从平衡值扰动。

在第一个实验中，操作员在两个输入端都引入了脉冲波。在第二个实验中，操作员在第一个输入中引入脉冲波，在第二个输入中引入阶跃信号。

准备数据

将数据加载到MATLAB工作空间

加载数据。

负载co2data

该命令将以下五个变量加载到MATLAB工作区中:

Input_exp1和Output_exp1分别是第一次实验的输入和输出数据。
Input_exp2和Output_exp2分别是第二次实验的输入和输出数据。
时间时间向量是从0到1000分钟，以相等的增量递增吗0.5分钟。

对于这两个实验，输入数据由两列值组成。第一列的值是化学物质消耗的速率(单位是千克每分钟)，第二列的值是电流(单位是安培)。输出数据是二氧化碳产生速率的单列(单位为毫克/分钟)。

绘制输入/输出数据

绘制两个实验的输入和输出数据。

情节(时间、Input_exp1时间,Output_exp1)传说(输入1的,“输入2”,“输出1”)图图(Time,Input_exp2,Time,Output_exp2)图例(输入1的,“输入2”,“输出1”)

第一个图显示了第一个实验，其中操作员对每个输入施加脉冲波，以扰动其稳态平衡。

第二个图显示了第二个实验，其中操作员对第一个输入施加脉冲波，对第二个输入施加步进信号。

从数据中去除均衡值

绘制数据，如绘制输入/输出数据，表明输入和输出具有非零平衡值。在本教程的这一部分中，您将从数据中减去均衡值。

你从每个信号中减去平均值的原因是，通常情况下，你建立了描述偏离物理平衡时的响应的线性模型。对于稳态数据，可以合理地假设信号的平均水平对应于这样的平衡。因此，你可以在零附近寻找模型，而不用在物理单位中建模绝对平衡水平。

放大图可以看到，操作员对输入施加扰动的最早时刻发生在稳态条件25分钟之后(或前50个样本之后)。因此，前50个样本的平均值代表了平衡条件。

使用以下命令去除两个实验中输入和输出的平衡值:

Input_exp1 = Input_exp1 -…(大小(Input_exp1, 1), 1) *意味着(Input_exp1 (1:50,:));Output_exp1 = Output_exp1 -…意思是(Output_exp1 (1:50,:));Input_exp2 = Input_exp2 -…(大小(Input_exp2, 1), 1) *意味着(Input_exp2 (1:50,:));Output_exp2 = Output_exp2 -…意思是(Output_exp2 (1:50,:));

使用对象表示数据进行系统识别

系统识别工具箱™数据对象，iddata和idfrd，将数据值和数据属性封装到单个实体中。可以使用“系统标识工具箱”命令方便地将这些数据对象作为单个实体进行操作。

在本教程的这一部分，您将创建两个职业iddata对象，两个实验各一个。你用实验1的数据进行模型估计，用实验2的数据进行模型验证。您使用两个独立的数据集，因为您使用一个数据集进行模型估计，另一个数据集进行模型验证。

请注意

当两个独立的数据集不可用时，可以将一个数据集分成两个部分，假设每个部分包含足够的信息来充分表示系统动态。

创建iddata对象

你一定已经把样本数据加载到MATLAB中了^®工作区，如加载数据到MATLAB工作区。

使用这些命令创建两个数据对象，泽和zv:

t = 0.5;%采样时间为0.5分钟泽= iddata (Output_exp1 Input_exp1, Ts);zv = iddata (Output_exp2 Input_exp2, Ts);

泽包含实验1和zv包含实验2的数据。Ts为采样时间。

的iddata构造函数需要三个实参用于时域数据，语法如下:

data_obj = iddata(输出、输入、sampling_interval);

的属性iddata对象,使用得到命令。例如，输入这个命令来获取估计数据的属性:

(泽)

ans = struct带字段:Domain: 'Time' Name: " OutputData: [2001x1 double] y: 'Same as OutputData' OutputName: {'y1'} OutputUnit: {"} InputData: [2001x2 double] u: 'Same as InputData' InputName: {2x1 cell} InputUnit: {2x1 cell} Period: [2x1 double] InterSample: {2x1 cell} Ts: 0.5000 Tstart: 0.5000 SamplingInstants: [2001x1 double] TimeUnit: 'seconds' ExperimentName: 'Exp1'注释:{}UserData: []

要了解更多关于此数据对象属性的信息，请参见iddata参考页面。

要修改数据属性，可以使用点表示法或集命令。例如，要分配通道名和标记图轴的单元，在MATLAB命令窗口中键入以下语法:

%设置时间单位为分钟泽。TimeUnit =“最小值”;设置输入通道的名称泽。InputName = {“ConsumptionRate”,“当前”};%设置输入变量的单位泽。InputUnit = {公斤/分钟的,“一个”};%设置输出通道的名称泽。OutputName =“生产”;设置输出通道单位泽。OutputUnit =毫克/分钟的;%设置验证数据属性zv。TimeUnit =“最小值”;zv。InputName = {“ConsumptionRate”,“当前”};zv。InputUnit = {公斤/分钟的,“一个”};zv。OutputName =“生产”;zv。OutputUnit =毫克/分钟的;

你可以验证InputName的属性泽被更改，或“索引”到此属性:

ze.inputname;

提示

属性名称，例如InputUnit，都不区分大小写。还可以缩写以。开头的属性名输入或输出用u为输入和y为输出在属性名中。例如,OutputUnit相当于yunit。

在数据对象中绘制数据

你可以画出iddata对象的使用情节命令。

绘制估计数据。

情节(泽)

底部的坐标轴显示输入ConsumptionRate和当前的，顶部的坐标轴显示输出ProductionRate。

在一个新的图形窗口中绘制验证数据。

数字%打开一个新的MATLAB Figure窗口情节(zv)绘制验证数据

放大图，可以看到实验过程放大了第一个输入(ConsumptionRate)的2倍，并放大第二个输入(当前的)的10倍。

选择数据的一个子集

在开始之前，创建一个新的数据集，其中只包含原始估计和验证数据集的前1000个样本，以加快计算速度。

Ze1 =泽(1:1000);Zv1 = zv (1:1000);

有关索引到的更多信息iddata对象，见相应参考页。

估计脉冲响应模型

为什么估计阶跃和频率响应模型?

频率响应和阶跃响应为非参数可以帮助你理解系统动态特性的模型。这些模型不是用带有可调参数的紧凑数学公式来表示的。相反，它们由数据表组成。

在本教程的这一部分中，您将使用数据集估计这些模型泽。你一定已经创建了。泽，如中所述创建iddata对象。

来自这些模型的响应图显示了系统的以下特征:

从第一个输入到输出的响应可能是一个二阶函数。
从第二个输入到输出的响应可能是一阶函数或过阻尼函数。

估计频率响应

系统识别工具箱产品为估计频率响应提供了三个函数:

etfe使用傅立叶分析计算经验传递函数。
水疗中心估计传递函数使用谱分析的固定频率分辨率。
spafdr让你指定一个可变的频率分辨率来估计频率响应。

使用水疗中心来估计频率响应。

通用电气= spa(泽);

将频率响应绘制为波德图。

波德(Ge)

振幅在0.54 rad/min的频率处达到峰值，这表明第一个输入-输出组合可能存在共振行为(复极)-ConsumptionRate来生产。

在这两个图中，相位滚动迅速，这表明两个输入/输出组合都存在时间延迟。

估计经验步骤响应

为了从数据中估计阶跃响应，首先从数据中估计一个非参数脉冲响应模型(FIR滤波器)，然后绘制其阶跃响应图。

%模型估计Mimp =冲动(Ze1、60);%阶跃响应步骤(Mimp)

第一个输入/输出组合的阶跃响应提示超调，这表明物理系统中存在欠阻尼模式(复极)。

从第二个输入到输出的阶跃响应显示没有超调，这表明要么是一阶响应，要么是具有实极点的高阶响应(过阻尼响应)。

阶跃响应图表明系统中存在非零延迟，这与您在创建的波德图中得到的快速相位滚离一致估计经验步骤响应。

多输入系统的延迟估计

为什么估计延迟?

为了识别参数黑箱模型，必须指定输入/输出延迟作为模型顺序的一部分。

如果你从实验中不知道你的系统的输入/输出延迟，你可以使用系统识别工具箱软件来估计延迟。

使用ARX模型结构估算延迟

在单输入系统的情况下，你可以在脉冲-响应图上读到延迟。然而，在多输入系统的情况下，例如本教程中的一个，你可能无法分辨哪个输入导致了输出的初始变化，你可以使用延迟命令。

该命令通过估计一个具有一定时延范围的低阶离散时间ARX模型来估计动态系统中的时延，然后选择与最佳拟合相对应的时延。

ARX模型结构是最简单的黑箱参数结构之一。在离散时间中，ARX结构是一个形式如下的差分方程:

$\begin{array}{l} y (t) + {一个}_{1} y (t - 1) + \dots + {一个}_{n 一个} y (t - n_{一个}) = \\ b_{1} u (t - n_{k}) + \dots + b_{n b} u (t - n_{k} - n_{b} + 1) + e (t) \end{array}$

y (t)表示时间的输出t,u (t)表示时间的输入t,n_一个为极点数，n_b的数量。b参数(等于0的个数加1)，n_k输入前的样本数量是否影响系统的输出(称为延迟或死时间的模型)，和e (t)是白噪声干扰。

延迟假设n_一个=n_b=2噪音e是白的还是无关紧要的，而估计呢n_k。

要估计这个系统中的延迟，输入:

延迟(泽)

Ans = 5 10

这个结果包括两个数字，因为有两个输入:第一个输入的估计延迟为5数据样本，第二次输入的估计延迟为10数据样本。因为实验的样本时间是0.5Min，这对应于a2.5-min delay在第一个输入影响输出之前，和一个5.0第二次输入影响输出之前的-min延迟。

使用可选方法估计延迟

估计系统时延有两种可选方法:

绘制输入输出数据的时间图，并读取输入的第一个变化与输出的第一个变化之间的时间差。这种方法仅适用于单输入/单输出系统;在多输入系统的情况下，你可能无法判断是哪个输入导致了输出的初始变化。
用1标准偏差置信区域绘制数据的脉冲响应。你可以用脉冲响应第一次出现在置信区域之外的时间来估计时间延迟。

使用ARX模型结构估计模型订单

为什么要估计模型顺序?

模型秩序定义模型复杂性的一个或多个整数。一般来说，模型的顺序与极点的数量、零的数量和响应延迟(输出响应输入之前的样本数量表示的时间)有关。模型顺序的具体含义取决于模型结构。

要计算参数黑盒模型，必须提供模型顺序作为输入。如果你不知道你的系统的阶数，你可以估计它。

完成本节的步骤后，你会得到以下结果:

对于第一个输入/输出组合:n_一个= 2,n_b=2，和延迟n_k= 5。
对于第二个输入/输出组合:n_一个= 1,n_b=1，和延迟n_k= 10。

稍后，您将通过指定模型顺序值来探索不同的模型结构，这些模型顺序值是围绕这些初始估计的轻微变化。

用于估计模型顺序的命令

在本教程的这一部分中，使用struc,arxstruc,selstruc对一定范围的模型订单和延迟估计和比较简单多项式(ARX)模型，并根据模型的质量选择最佳订单。

下面的列表描述了使用每个命令的结果:

struc为的指定范围创建模型顺序组合的矩阵n_一个,n_b,n_k值。
arxstruc的输出。struc，系统地为每个模型订单估计一个ARX模型，并将模型输出与测量输出进行比较。arxstruc返回损失函数对于每个模型，它是预测误差的平方和的规范化和。
selstruc的输出。arxstruc，打开ARX模型结构选择窗口，如图所示，帮助您选择模型订单。
你用这个图来选择最适合的模型。
- 横轴为参数总数-n_一个+n_b。
- 纵轴，叫无法解释的输出方差(%),输出的部分不是用模型来解释ARX模型预测误差参数的数量显示在水平轴上。
  的预测误差为验证数据输出与模型领先一步预测输出之间差异的平方和。
- n_k是延迟。
三个矩形在图中以绿色、蓝色和红色突出显示。每种颜色表示一种最佳拟合准则，如下所示:
- Red -最佳拟合将验证数据输出与模型输出之间的差的平方和最小化。这个矩形表示整体的最佳拟合。
- Green -最佳适合最小化Rissanen MDL标准。
- Blue - Best fit将赤池AIC准则最小化。
在本教程中，无法解释的输出方差(%)对于从4到20的参数的组合数量，值保持大致不变。这样的恒定表明模型性能在更高的阶数下并没有改善。因此，低阶模型可能同样能很好地拟合数据。

请注意

当你使用相同的数据集进行估计和验证时，使用MDL和AIC标准来选择模型顺序。这些准则弥补了由于使用过多参数而导致的过拟合。有关这些标准的更多信息，请参见selstruc参考页面。

第一个输入输出组合的模型顺序

在本教程中，系统有两个输入和一个输出，你可以独立估计每个输入/输出组合的模型顺序。你可以同时估计两个输入的延迟，也可以一次估计一个输入的延迟。

对于第一个输入/输出组合，尝试以下顺序组合是有意义的:

n_一个=2:5
n_b=1:5
n_k=5

这是因为你估计的非参数模型估计脉冲响应模型证明第一个输入/输出组合的响应可能具有二阶响应。同样,在多输入系统的延迟估计，这个输入/输出组合的延迟估计为5。

估计第一个输入/输出组合的模型顺序:

使用struc创建一个可能的模型订单矩阵。
```
NN1 = struc (2:5 1:5 5);
```
使用selstruc计算输入订单的ARX模型的损失函数NN1。
```
selstruc (arxstruc(泽(:,:1),zv (:,: 1), NN1))
```
请注意

泽(:,:1)选择数据中的第一个输入。

这个命令打开交互式ARX Model Structure Selection窗口。

请注意

Rissanen MDL和Akaike AIC准则产生等价的结果，都在图上用蓝色矩形表示。

红色矩形代表最适合的整体，这发生在n_一个= 2,n_b= 3,n_k= 5。红色和蓝色矩形之间的高度差是不显著的。因此，可以选择对应最低模型阶数和最简单模型的参数组合。
点击蓝色矩形，然后点击选择要选择那个订单组合:
n_一个=2
n_b=2
n_k=5
要继续，在MATLAB命令窗口中按任意键。

第二次投入产出组合的模型顺序

对于第二个输入/输出组合，可以尝试以下顺序组合:

n_一个=1:3
n_b=1:3
n_k=10

这是因为你估计的非参数模型估计脉冲响应模型表明第二次输入/输出组合的响应可能具有一阶响应。同样,在多输入系统的延迟估计，这个输入/输出组合的延迟估计为10。

估计第二个输入/输出组合的模型顺序:

使用struc创建一个可能的模型订单矩阵。
```
NN2 = struc (1:3, 1:3, 10);
```
使用selstruc计算输入订单的ARX模型的损失函数NN2。
```
selstruc (arxstruc(泽(:,:2),zv (:,: 2), NN2))
```
这个命令打开交互式ARX Model Structure Selection窗口。

请注意

Akaike AIC和整体最佳拟合标准产生了等价的结果。两者都由图上相同的红色矩形表示。

所有矩形之间的高度差是不显著的，所有这些模型顺序都会导致相似的模型性能。因此，可以选择对应最低模型阶数和最简单模型的参数组合。
点击最左边的黄色矩形，然后点击选择选择最低的顺序:n_一个= 1,n_b= 1,n_k= 10。
要继续，在MATLAB命令窗口中按任意键。

估计转移函数

指定传递函数的结构

在本教程的这一部分，你估计一个连续时间传递函数。您必须已经准备好您的数据，如准备数据。你可以使用以下估计模型顺序的结果来指定模型的顺序:

对于第一个输入/输出组合，使用:
- 两极，对应na = 2在ARX模型。
- 延迟5，对应nk = 5ARX模型中的样本(或2.5分钟)。
对于第二个输入/输出组合，使用:
- 一极，对应na = 1在ARX模型中
- 延迟为10，对应于nk = 10ARX模型中的样本(或5分钟)。

你可以估计这些阶的传递函数使用特遣部队命令。，您也可以选择查看进度报告显示选项在的创建的选项集中tfestOptions命令。

选择= tfestOptions (“显示”,“上”);

将模型指令和延迟收集成要传递给的变量特遣部队。

Np = [2 1];ioDelay = [2.5 5];

估计传递函数。

mtf =特遣部队(Ze1、np、[],ioDelay,选择);

查看模型的系数。

mtf

mtf =从输入“ConsumptionRate”到输出“Production”:-1.382 s + 0.0008305 exp(-2.5*s) * ------------------------- s^2 + 1.014 s + 5.444e-12从输入“Current”到输出“Production”:5.93 exp(-5*s) * ---------- s + 0.5858连续时间识别传递函数。参数化:极点数:[2 1]零数:[1 0]自由系数数:6使用“tfdata”，“getpvec”，“getcov”来表示参数及其不确定性。状态:使用时域数据“Ze1”上的TFEST进行估计。拟合估计数据:78.92% FPE: 14.22, MSE: 13.97

模型的显示显示超过85%的拟合估计数据。

验证模型

在本教程的这一部分中，您将创建一个图来比较实际输出和模型输出使用比较命令。

比较(Zv1 mtf)

比较显示约60%适合。

残留分析

使用渣油命令来评估残差的特征。

渣油(Zv1 mtf)

残差显示出高度的自相关性。这并不出乎意料，因为模型mtf没有任何组件来单独描述噪声的性质。要对测量的输入输出动力学和扰动动力学进行建模，您需要使用一个包含描述噪声成分的元素的模型结构。您可以使用bj,党卫军和过程命令，它们分别创建多项式、状态空间和过程结构的模型。其中，这些结构包含捕获噪声行为的元素。

评估过程模型

指定过程模型的结构

在本教程的这一部分中，您估计了一个低阶、连续时间的传递函数(过程模型)。系统识别工具箱产品支持最多三个极点(其中可能包含欠阻尼极点)、一个零、一个延迟元素和一个积分器的连续时间模型。

您必须已经准备好您的数据，如准备数据。

你可以使用以下估计模型顺序的结果来指定模型的顺序:

对于第一个输入/输出组合，使用:
- 两极，对应n_一个= ARX模型中的2。
- 延迟5，对应n_k= ARX模型中的5个样本(或2.5分钟)。
对于第二个输入/输出组合，使用:
- 一极，对应n_一个= ARX模型中的1。
- 延迟为10，对应于n_k= ARX模型中的10个样本(或5分钟)。

请注意

因为离散时间ARX模型估计的零个数和它的连续时间对应对象之间没有关系，所以你没有对零个数的估计。在本教程中，您可以为第一个输入/输出组合指定一个零，而为第二个输出组合指定不为零。

使用idproc命令创建两个模型结构，每个模型对应一个输入/输出组合:

midproc0 = idproc ({“P2ZUD”,“P1D”},“TimeUnit”,“分钟”);

细胞数组{“P2ZUD”、“P1D”}指定每个输入/输出组合的模型结构:

“P2ZUD”表示具有两个极点的传递函数(P2)，一个零(Z)，欠阻尼(复共轭)极点(U)和延迟(D)。
“P1D”表示具有一极点的传递函数(P1)和延迟(D)。

该示例还使用了TimeUnit参数来指定所使用的时间单位。

查看模型结构和参数值

查看两个得到的模型。

midproc0

midproc0 =流程模型与2输入:y = G11 (s) u1 + G12 (s) u2从输入1输出1:1 + Tz * s G11 (s) = Kp  * ---------------------- * exp (Td * s) 1 + 2 *ζ* Tw * s + (Tw * s) ^ 2 Kp =南Tw Td =南Tz = =南ζ=南南从输入2输出1:Kp G12 (s ) = ---------- * exp (Td * s) 1 + Tp1 * s Kp Td = =南Tp1 =南南参数化:{“P2DUZ”}{‘P1D}很多免费的系数:8使用“getpvec”、“getcov”参数及其不确定性。状态:由直接构建或转换创建。不估计。

参数值设置为南因为它们还没有估算出来。

指定时间延迟的初始猜测

将模型对象的延时属性设置为2.5最小值和5为每个输入/输出对的最小值作为初始猜测。此外，设置延迟的上限，因为有良好的初始猜测。

midproc0.Structure .Td(1,1)。值= 2.5;midproc0.Structure .Td(1、2)。值= 5;midproc0.Structure .Td(1,1)。最大= 3;midproc0.Structure .Td(1、2)。最大= 7;

请注意

在设置延迟约束时，必须以实际时间单位(在本例中为分钟)来指定延迟，而不是样本的数量。

利用过程估计模型参数

过程是一个迭代过程模型的估计量，这意味着它使用迭代非线性最小二乘算法来最小化成本函数。的成本函数为误差平方和的加权和。

根据它的论点，过程估计不同的黑盒多项式模型。您可以使用过程，例如，估计线性连续时间传递函数、状态空间或多项式模型结构的参数。使用过程，您必须提供一个具有未知参数和估计数据作为输入参数的模型结构。

对于本教程的这一部分，您必须已经定义了模型结构，如中所述指定过程模型的结构。使用midproc0作为模型结构和Ze1作为估计数据:

midproc = proc (Ze1 midproc0);礼物(midproc)

midproc =流程模型与2输入:y = G11 (s) u1 + G12 (s) u2从输入输出“生产”:“ConsumptionRate”1 + Tz * s G11 (s) = Kp  * ---------------------- * exp (Td * s) 1 + 2 *ζ* Tw * s + (Tw * s) ^ 2 Kp = -1.1807 + / - 0.29986太瓦= 1.6437 + / - 714.6ζ= 16.036 + / - 6958.9 Td = 2.426 + / - 64.276 Tz = -109.19 + / - 63.733从输入输出“生产”“当前”:Kp G12 (s ) = ---------- * exp (Td * s) 1 + Tp1 * s Kp = 10.264 + / - 0.048404 Tp1 = 2.049 + / - 0.054901 Td = 4.9175 + / - 0.034374参数化:{'P2DUZ'} {'P1D'}自由系数数:8使用"getpvec"， "getcov"表示参数及其不确定性。状态:终止条件:达到最大迭代次数..迭代次数:20，函数求值次数:279使用PROCEST对时域数据“Ze1”进行估计。拟合估计数据:86.2% FPE: 6.081, MSE: 5.984模型“Report”属性中的更多信息。

与离散时间多项式模型不同，连续时间模型可以让你估计延迟。在这种情况下，估计的延迟值与您指定的初始猜测不同2.5和5,分别。的参数估计值中存在较大的不确定性G_1里面(s)表示来自输入的动态1对输出的贡献没有很好地捕获。

要了解更多关于估计模型的信息，请参见流程模型。

验证模型

在本教程的这一部分中，您创建了一个比较实际输出和模型输出的图。

比较(Zv1 midproc)

图显示，模型输出与测量输出合理一致:模型与验证数据之间有60%的一致性。

使用渣油进行残差分析。

渣油(Zv1 midproc)

第二次输入和残差之间的相互关系是显著的。自相关图显示置信区域以外的值，表明残差是相关的。

将迭代参数估计算法改为Levenberg-Marquardt。

选择= procestOptions;Opt.SearchMethod =“lm”;midproc1 = proc (Ze1 midproc,选择);

调整算法属性或指定初始参数猜测并重新运行估计可能会改善模拟结果。添加噪声模型可能会改善预测结果，但不一定会改善仿真结果。

用噪声模型估计一个过程模型

本教程的这一部分展示了如何估计一个过程模型，并在估计中包括一个噪声模型。包含噪声模型通常会改善模型预测结果，但不一定会改善仿真结果。

使用以下命令指定一阶ARMA噪声:

选择= procestOptions;Opt.DisturbanceModel =“ARMA1”;midproc2 = proc (Ze1 midproc0,选择);

你可以输入“距离”而不是“DisturbanceModel”。属性名称不区分大小写，只需要包含名称中唯一标识属性的部分。

将得到的模型与实测数据进行比较。

比较(Zv1 midproc2)

图显示，模型输出与验证数据输出保持了合理的一致性。

执行残留分析。

渣油(Zv1 midproc2)

残差图显示自相关值在置信范围内。因此，添加一个噪声模型会产生不相关的残差。这表明了一个更精确的模型。

估计黑箱多项式模型

用于估计多项式模型的模型顺序

在本教程的这一部分中，您将估计几种不同类型的黑箱、输入输出多项式模型。

您必须已经准备好您的数据，如准备数据。

你可以用以下先前估计模型阶数的结果来指定多项式模型的阶数:

对于第一个输入/输出组合，使用:
- 两极，对应n_一个= ARX模型中的2。
- 一个零，对应n_b= ARX模型中的2。
- 延迟5，对应n_k= ARX模型中的5个样本(或2.5分钟)。
对于第二个输入/输出组合，使用:
- 一极，对应n_一个= ARX模型中的1。
- 无零，对应n_b= ARX模型中的1。
- 延迟为10，对应于n_k= ARX模型中的10个样本(或5分钟)。

估计一个线性ARX模型

关于ARX模型。对于单输入/单输出系统(SISO)， ARX模型结构为:

$\begin{array}{l} y (t) + {一个}_{1} y (t - 1) + \dots + {一个}_{n 一个} y (t - n_{一个}) = \\ b_{1} u (t - n_{k}) + \dots + b_{n b} u (t - n_{k} - n_{b} + 1) + e (t) \end{array}$

y (t)表示时间的输出t,u (t)表示时间的输入t,n_一个为极点数，n_b是0的个数加1，n_k输入前的样本数量是否影响系统的输出(称为延迟或死时间的模型)，和e (t)是白噪声干扰。

ARX模型结构不区分单个输入/输出路径的极点:将ARX方程除以一个包含两极的facebook就说明了这一点一个出现在两个输入的分母中。因此，可以设置n_一个到每个输入/输出对的极点之和，等于3.在这种情况下。

系统标识工具箱产品估计参数 ${一个}_{1} \dots {一个}_{n}$ 和 $b_{1} \dots b_{n}$ 使用您指定的数据和模型顺序。

使用ARX估计ARX模型。使用arx用快速、非迭代的方法计算多项式系数arx:

马克思= arx (Ze1,“na”3,“注”(2 - 1),“朝鲜”10 [5]);礼物(马克思)%显示模型参数%带有不确定性信息

marx =离散时间ARX模型:A(z)y(t) = B(z)u(t) + e(t) A(z) = 1 - 1.027 (+/- 0.02907) z^-1 + 0.1678 (+/- 0.02583) z^-3 B1(z) = 1.86 (+/- 0.189) z^-5 - 1.608 (+/- 0.1888) z^-6 B2(z) = 1.612 (+/- 0.07392) z^-10采样时间:0.5分钟参数化:多项式阶:na=3 nb=[2 1] nk=[5 10]自由系数数:6使用"polydata"、"getpvec"、"getcov"表示参数及其不确定性。状态:在时域数据“Ze1”上使用ARX进行估计。拟合估计数据:90.7%(预测焦点)FPE: 2.768, MSE: 2.719模型“Report”属性中的更多信息。

MATLAB估计多项式一个,B1,B2。每个模型参数的不确定度计算为1个标准差，并显示在每个参数值旁边的括号中。

或者，你也可以使用下面的简写语法，将模型顺序指定为单个向量:

marx = arx(Ze1，[3 2 1 5 10]);

访问模型数据。你估计的模型，马克思，是离散时间idpoly对象。要获取此模型对象的属性，可以使用得到功能:

(马克思)

A: [1 -1.0267 0.1678 0.0129] B: {[0 0 0 0 1.8599 -1.6084] [0 0 0 0 0 0 0 0 0 0 0 0 1.6118]} C: 1 D: 1 F: {[1] [1]} IntegrateNoise: 0变量:'z^-1' IODelay:[0 0]结构:[1x1 pmodel。polynomial] NoiseVariance: 2.7436 InputDelay: [2x1 double] OutputDelay: 0 InputName: {2x1 cell} InputUnit: {2x1 cell} InputGroup: [1x1 struct] OutputName: {'Production'} OutputUnit: {'mg/min'} OutputGroup: [1x1 struct]注释:[0x1 string] UserData: [] Name: " Ts: 0.5000 TimeUnit: 'minutes' SamplingGrid: [1x1 struct]报告:[1x1 idresults.arx]

您可以使用点表示法访问由这些属性存储的信息。例如，你可以通过找到的根来计算模型的离散极点一个多项式。

marx_poles =根(marx.a)

马克思极点= 0.7953 0.2877 -0.0564

在这种情况下，您访问一个多项式使用marx.a。

该模型马克思使用三个离散极点描述系统动力学。

提示

你也可以使用极直接计算模型的极点。

学习更多的知识。要了解关于估计多项式模型的更多信息，请参见输入-输出多项式模型。

有关访问模型数据的更多信息，请参见数据提取。

估计状态空间模型

对状态空间模型。一般状态空间模型结构为:

$\begin{array}{l} \frac{d x}{d t} = 一个 x (t) + B u (t) + K e (t) \\ y (t) = C x (t) + D u (t) + e (t) \end{array}$

y (t)表示时间的输出t,u (t)表示时间的输入t,x (t)状态向量在时间t,e (t)是白噪声干扰。

必须指定一个整数作为模型顺序(状态向量的维度)，以估计状态空间模型。默认情况下，延迟等于1。

系统识别工具箱产品估计状态空间矩阵一个,B,C,D,K使用指定的模型订单和数据。

状态空间模型结构是快速估计的好选择，因为它只包含两个参数:n是杆数(的大小一个矩阵),nk是延迟。

使用n4sid估计状态空间模型。使用n4sid命令指定模型订单的范围，并评估几个状态空间模型的性能(订单2到8):

mn4sid = n4sid (Ze1 2:8,“InputDelay”9 [4]);

该命令使用快速、非迭代(子空间)方法，并打开以下图。使用这个图来决定哪些状态对输入/输出行为提供了显著的相对贡献，哪些状态提供了最小的贡献。

纵轴是对每个状态对模型的输入/输出行为贡献多少的相对度量(协方差矩阵奇异值的对数)。横轴对应模型顺序n。这块地建议n = 3，用红色矩形表示。

的选择的顺序订单字段显示推荐的模型订单，3.在这种情况下，默认为。的方法可以更改顺序选择选择的顺序下拉列表。的值应用选择的顺序字段并通过单击关闭订单选择窗口应用。

默认情况下,n4sid使用状态空间形式的自由参数化。的值来估计标准形式SSParameterization财产“规范”。，还可以指定输入到输出延迟(在示例中)“InputDelay”财产。

mCanonical = n4sid (Ze1 3“SSParameterization”,“规范”,“InputDelay”9 [4]);礼物(mCanonical);%显示模型属性

mCanonical =离散时间识别状态空间模型:x (t + Ts) = x (t) + B K u (t) + e (t) y (t) = C x (t) + D u (t) + e (t) = (x1, x2) x3 x1 1 0 1 0 x2 0 0 x3 0.0737 + / - 0.05919 - -0.6093 + / - 0.1626 - 1.446 + / - 0.1287 B = ConsumptionR当前x1 1.844 + / - 0.175 - 0.5633 + / - 0.122 x2 1.063 + / - 0.1673 2.308 + / - 0.1222 x3 0.2779 + / - 0.09615 - 1.878 + / - 0.1058 C = (x1, x2) x3生产1 0 0 D = ConsumptionR当前生产0 0 K = 0.8674 + / - 0.03169 x1 x2 0.6849 + / - 0.04145 x3 0.5105 + / - 0.04352输入延迟(采样周期):4 9样时间:0.5分钟参数化:带指数的CANONICAL形式:3。馈通:无扰动分量:估计自由系数数量:12使用“idssdata”，“getpvec”，“getcov”表示参数及其不确定性。状态:使用时域数据“Ze1”上的N4SID进行估计。拟合估计数据:91.39%(预测焦点)FPE: 2.402, MSE: 2.331模型“Report”属性中的更多信息。

请注意

mn4sid和mCanonical离散时间模型。要估计一个连续时间模型，请设置“t”财产0在估计命令中，或使用党卫军命令:

mCT1 = n4sid(Ze1, 3，“t”0,“InputDelay”， [2.5 5]) mCT2 = ssest(Ze1, 3，“InputDelay”(2.5 - 5))

学习更多的知识。要了解关于估计状态空间模型的更多信息，请参见状态空间模型。

Box-Jenkins模型的估计

关于Box-Jenkins模型。一般Box-Jenkins (BJ)结构为:

$y (t) = \sum_{我 = 1}^{n u} \frac{B_{我} (问)}{F_{我} (问)} u_{我} (t - n k_{我}) + \frac{C (问)}{D (问)} e (t)$

要估计一个BJ模型，需要指定参数n_b,n_f,n_c,n_d,n_k。

ARX模型结构没有区分单个输入/输出路径的极点，而BJ模型在将扰动的极点和零与系统动力学的极点和零分开建模方面提供了更大的灵活性。

使用polyest估计一个BJ模型。您可以使用聚来估计BJ模型。聚是一种迭代方法，具有以下通用语法:

Polyest (data，[na nb nc nd nf nk]);

要估计BJ模型，输入:

na = 0;Nb = [2 1];数控= 1;nd = 1;Nf = [2 1];Nk = [5 10];mbj = polyest(Ze1，[na nb nc nd nf nk]);

这个命令指定nf = 2,nb = 2,nk = 5对于第一个输入，和nf = nb = 1和nk = 10对于第二个输入。

显示模型信息。

礼物(mbj)

mbj =离散BJ模型:y (t) = (B (z) / F (z)] u (t) + (C (z) / D (z)] e (t) B1 (z) = 1.823 (+ / - 0.1792) z ^ 5 - 1.315 (+ / - 0.2367) z ^ 6 B2 (z) = 1.791 (+ / - 0.06431) z ^ -10 C (z) = 1 + 0.1068 (+ / - 0.04009) z ^ 1 D (z) = 1 - 0.7452 (+ / - 0.02694) z ^ 1 F1 (z) = 1 - 1.321 (+ / - 0.06936) z ^ 1 + 0.5911 (+ / - 0.05514) z ^ 2 F2 (z) = 1 - 0.8314 (+ / - 0.006441) z ^ 1样品时间:0.5分钟参数化:多项式订单:nb =[2 1]数控= 1 nd = 1 nf = [2 1] nk = 10[5]很多免费的系数:8使用“polydata”、“getpvec”、“getcov”参数及其不确定性。状态:终止条件:接近(局部)最小值，(norm(g) < tol)..迭代次数:6，函数求值次数:13使用POLYEST对时域数据“Ze1”进行估计。拟合估计数据:90.75%(预测焦点)FPE: 2.733, MSE: 2.689模型“Report”属性中的更多信息。

每个模型参数的不确定度计算为1个标准差，并显示在每个参数值旁边的括号中。

的多项式C和D分别给出噪声模型的分子和分母。

提示

或者，你也可以使用下面的简写语法，将顺序指定为单个向量:

mbj = bj(Ze1，[2 1 1 1 2 1 5 10]);

bj是的一个版本聚具体估计了BJ模型结构。

学习更多的知识。要了解更多关于识别输入输出多项式模型(如BJ)的信息，请参见输入-输出多项式模型。

比较模型输出和测量输出

将ARX、状态空间和Box-Jenkins模型的输出与测量输出进行比较。

比较(Zv1,马克思,mn4sid mbj)

比较将验证数据集中的测量输出与模型的模拟输出进行对比。来自验证数据集的输入数据作为模型的输入。

对ARX、状态空间和Box-Jenkins模型执行残差分析。

渣油(Zv1,马克思,mn4sid mbj)

这三种模型都能很好地模拟输出，并且具有不相关的残差。因此，选择ARX模型，因为它是三个输入输出多项式模型中最简单的，并充分捕捉了过程动态。

模拟和预测模型输出

模拟模型输出

在本教程的这一部分中，您将模拟模型输出。您必须已经创建了连续时间模型midproc2，如中所述评估过程模型。

模拟模型输出需要以下信息:

向模型输入值
模拟的初始条件(也称为初始状态)

例如，以下命令使用iddata和idinput命令来构造一个输入数据集，并使用sim卡模拟模型输出:

%为模拟创建输入U = iddata([]，idinput([200 2])，“t”, 0.5,“TimeUnit”,“最小值”);%模拟设置初始条件为零的响应ysim_1 = sim (midproc2 U);

为了最大化模型的模拟响应与相同输入的测量输出之间的拟合，可以计算测量数据对应的初始条件。通过使用可以获得最佳拟合初始条件findstates在估计模型的状态空间版本上。以下命令估计初始状态x0从数据集中Zv1:

计算初始状态所需的模型的%状态空间版本midproc2_ss = ids (midproc2);x0 = findstates (midproc2_ss Zv1);

接下来，使用从数据中估计的初始状态来模拟模型。

%的模拟输入Usim = Zv1 ([],::);选择= simOptions (“InitialCondition”, x0);ysim_2 = sim (midproc2_ss、Usim选择);

在一个图上比较模拟输出和测量输出。

图([ysim_2的阴谋。y, Zv1.y])传说({模型输出的,“测量”})包含(“时间”), ylabel (“输出”)

预测未来产出

许多控制设计应用程序要求您使用过去的输入/输出数据来预测动态系统的未来输出。

例如,使用预测提前5步预测模型响应:

预测(midproc2 Ze1 5)

将预测的输出值与测量的输出值进行比较。的第三个参数比较指定提前五步预测。当不指定第三个参数时，比较假设一个无限的预测视界，并转而模拟模型输出。

比较(Ze1 midproc2 5)

使用体育计算预测误差犯错的预测输出之间midproc2以及测量的输出。然后，绘制误差谱光谱命令。

(错)= pe (midproc2 Zv1);谱(spa(呃,[],logspace (2200)))

预测误差用1标准偏差置信区间绘制。由于扰动的高频性质，在高频处误差更大。