System Identification Toolbox 中的数据和模型对象
本示例演示了如何管理 System Identification Toolbox™ 中可用的数据和模型对象。系统辨识旨在根据数据构建模型。一个数据集表征为以下几项信息：输入和输出信号、采样时间、变量名称和单位等。同样，估计模型中包含各种信息 - 估计参数、其协方差矩阵、模型结构等等。
这意味着，将数据和模型相关的信息整合到对象中是合适且可取的。System Identification Toolbox 包含多个此类对象，本示例介绍了这些对象的基本特性。
IDDATA 对象
首先创建一些数据：
u = sign(randn(200,2)); % 2 inputs y = randn(200,1); % 1 output ts = 0.1; % The sample time
要将输入和输出收集到一个对象中，请执行以下操作：
z = iddata(y,u,ts);
只需输入数据名称，即可显示相关信息：
z
z =
Time domain data set with 200 samples.
Sample time: 0.1 seconds
Outputs Unit (if specified)
y1
Inputs Unit (if specified)
u1
u2
该数据由 plot 命令绘制为 iddata，与 plot(z) 中的情况相同。按任意键继续，并在各支线剧情之间切换。在此，我们将各通道分别绘制出来：
plot(z(:,1,1)) % Data subset with Input 1 and Output 1.
plot(z(:,1,2)) % Data subset with Input 2 and Output 1.
要获取输出和输入，请使用
u = z.u; % or, equivalently u = get(z,'u'); y = z.y; % or, equivalently y = get(z,'y');
要选择部分数据：
zp = z(48:79);
要选择第一个输出和第二个输入：
zs = z(:,1,2); % The ':' refers to all the data time points.
这些子选项可以组合使用：
plot(z(45:54,1,2)) % samples 45 to 54 of response from second input to the first output.
这些通道被赋予了 'y1'、'u2'等默认名称。可以通过以下方式将其更改为任意值：
set(z,'InputName',{'Voltage';'Current'},'OutputName','Speed');
等价地，也可以写成
z.inputn = {'Voltage';'Current'}; % Autofill is used for properties
z.outputn = 'Speed'; % Upper and lower cases are also ignored
对于记账和绘图，还可以设置单位：
z.InputUnit = {'Volt';'Ampere'};
z.OutputUnit = 'm/s';
z
z =
Time domain data set with 200 samples.
Sample time: 0.1 seconds
Outputs Unit (if specified)
Speed m/s
Inputs Unit (if specified)
Voltage Volt
Current Ampere
所有当前属性（与任何对象一样）都是通过 get 获取的：
get(z)
ans =
struct with fields:
Domain: 'Time'
Name: ''
OutputData: [200×1 double]
y: 'Same as OutputData'
OutputName: {'Speed'}
OutputUnit: {'m/s'}
InputData: [200×2 double]
u: 'Same as InputData'
InputName: {2×1 cell}
InputUnit: {2×1 cell}
Period: [2×1 double]
InterSample: {2×1 cell}
Ts: 0.1000
Tstart: 0.1000
SamplingInstants: [200×1 double]
TimeUnit: 'seconds'
ExperimentName: 'Exp1'
Notes: {}
UserData: []
除了前面讨论过的属性外，还有“周期” (Period) 属性，它表示输入的周期；如果周期为 inf，则表示非周期性输入：
z.Period
ans = Inf Inf
输入的采样间隔内行为可指定为 'zoh'（零阶保持，即分段常量）或 'foh'（一阶保持，即分段线性）。辨识程序利用这些信息来调整算法。
z.InterSample
ans =
2×1 cell array
{'zoh'}
{'zoh'}
可以通过“水平拼接”来添加通道（包括输入和输出），即 z = [z1 z2]：
z2 = iddata(rand(200,1),ones(200,1),0.1,'OutputName','New Output',... 'InputName','New Input'); z3 = [z,z2]
z3 =
Time domain data set with 200 samples.
Sample time: 0.1 seconds
Outputs Unit (if specified)
Speed m/s
New Output
Inputs Unit (if specified)
Voltage Volt
Current Ampere
New Input
让我们绘制 z3 的一些通道：
plot(z3(:,1,1)) % Data subset with Input 2 and Output 1.
plot(z3(:,2,3)) % Data subset with Input 2 and Output 3.
生成输入
命令 idinput 会生成典型的输入信号。
u = idinput([30 1 10],'sine'); % 10 periods of 30 samples u = iddata([],u,1,'Period',30) % Making the input an IDDATA object.
u =
Time domain data set with 300 samples.
Sample time: 1 seconds
Inputs Unit (if specified)
u1
将 SIM 应用于 iddata 输入后，将产生 iddata 输出。让我们使用 sim 来获取基于输入 u 的估计模型 m 的响应。我们还会根据模型的噪声动态特性，向模型响应中添加噪声。我们通过使用 "AddNoise" 仿真选项来实现这一点：
m = idpoly([1 -1.5 0.7],[0 1 0.5]); % This creates a model; see below. options = simOptions; options.AddNoise = true; y = sim(m,u,options) % simulated response produced as an iddata object
y =
Time domain data set with 300 samples.
Sample time: 1 seconds
Name: m
Outputs Unit (if specified)
y1
仿真输入 u 和输出 y 可以按以下方式合并为一个 iddata 对象：
z5 = [y u] % The output-input iddata.
z5 =
Time domain data set with 300 samples.
Sample time: 1 seconds
Name: m
Outputs Unit (if specified)
y1
Inputs Unit (if specified)
u1
有关 iddata 对象的更多信息，请参阅 help iddata。
线性模型对象
所有模型均以 MATLAB® 对象的形式提供。根据所用模型的类型不同，会有几种不同的对象，但这部分内容基本上是透明的。
load iddata1 m = armax(z1,[2 2 2 1]); % This creates an ARMAX model, delivered as an IDPOLY object
该模型的所有相关属性都被封装为一个对象（此处为 idpoly）。要显示它，只需输入其名称：
m
m =
Discrete-time ARMAX model: A(z)y(t) = B(z)u(t) + C(z)e(t)
A(z) = 1 - 1.531 z^-1 + 0.7293 z^-2
B(z) = 0.943 z^-1 + 0.5224 z^-2
C(z) = 1 - 1.059 z^-1 + 0.1968 z^-2
Sample time: 0.1 seconds
Parameterization:
Polynomial orders: na=2 nb=2 nc=2 nk=1
Number of free coefficients: 6
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARMAX on time domain data "z1".
Fit to estimation data: 76.38% (prediction focus)
FPE: 1.127, MSE: 1.082
Many of the model properties are directly accessible
m.a % The A-polynomial
ans =
1.0000 -1.5312 0.7293
通过 get 方法可以获取属性列表：
get(m)
A: [1 -1.5312 0.7293]
B: [0 0.9430 0.5224]
C: [1 -1.0587 0.1968]
D: 1
F: 1
IntegrateNoise: 0
Variable: 'z^-1'
IODelay: 0
Structure: [1×1 pmodel.polynomial]
NoiseVariance: 1.1045
InputDelay: 0
OutputDelay: 0
InputName: {'u1'}
InputUnit: {''}
InputGroup: [1×1 struct]
OutputName: {'y1'}
OutputUnit: {''}
OutputGroup: [1×1 struct]
Notes: [0×1 string]
UserData: []
Name: ''
Ts: 0.1000
TimeUnit: 'seconds'
SamplingGrid: [1×1 struct]
Report: [1×1 idresults.polyest]
使用 present 将参数协方差的估计值表示为各参数的 ±1 个标准差的不确定度值：
present(m)
m =
Discrete-time ARMAX model: A(z)y(t) = B(z)u(t) + C(z)e(t)
A(z) = 1 - 1.531 (+/- 0.01801) z^-1 + 0.7293 (+/- 0.01473) z^-2
B(z) = 0.943 (+/- 0.06074) z^-1 + 0.5224 (+/- 0.07818) z^-2
C(z) = 1 - 1.059 (+/- 0.06067) z^-1 + 0.1968 (+/- 0.05957) z^-2
Sample time: 0.1 seconds
Parameterization:
Polynomial orders: na=2 nb=2 nc=2 nk=1
Number of free coefficients: 6
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Termination condition: Near (local) minimum, (norm(g) < tol)..
Number of iterations: 3, Number of function evaluations: 7
Estimated using ARMAX on time domain data "z1".
Fit to estimation data: 76.38% (prediction focus)
FPE: 1.127, MSE: 1.082
More information in model's "Report" property.
使用 getpvec 获取所有模型参数（或仅自由参数）及其不确定度的扁平列表。使用 getcov 获取完整的协方差矩阵。
[par, dpar] = getpvec(m, 'free') CovFree = getcov(m,'value')
par =
-1.5312
0.7293
0.9430
0.5224
-1.0587
0.1968
dpar =
0.0180
0.0147
0.0607
0.0782
0.0607
0.0596
CovFree =
0.0003 -0.0003 0.0000 0.0007 0.0004 -0.0003
-0.0003 0.0002 -0.0000 -0.0004 -0.0003 0.0002
0.0000 -0.0000 0.0037 -0.0034 -0.0000 0.0001
0.0007 -0.0004 -0.0034 0.0061 0.0008 -0.0005
0.0004 -0.0003 -0.0000 0.0008 0.0037 -0.0032
-0.0003 0.0002 0.0001 -0.0005 -0.0032 0.0035
nf = 0、nd = 0 分别表示一般线性模型的阶数，而 ARMAX 模型是其特例。
该报告包含有关估计过程的信息：
m.Report m.Report.DataUsed % record of data used for estimation m.Report.Fit % quantitative measures of model quality m.Report.Termination % search termination conditions
ans =
Status: 'Estimated using ARMAX with prediction focus'
Method: 'ARMAX'
InitialCondition: 'zero'
Fit: [1×1 struct]
Parameters: [1×1 struct]
OptionsUsed: [1×1 idoptions.polyest]
RandState: [1×1 struct]
DataUsed: [1×1 struct]
Termination: [1×1 struct]
ans =
struct with fields:
Name: 'z1'
Type: 'Time domain data'
Length: 300
Ts: 0.1000
InterSample: 'zoh'
InputOffset: []
OutputOffset: []
ans =
struct with fields:
FitPercent: 76.3807
LossFcn: 1.0824
MSE: 1.0824
FPE: 1.1266
AIC: 887.1256
AICc: 887.4123
nAIC: 0.1192
BIC: 909.3483
ans =
struct with fields:
WhyStop: 'Near (local) minimum, (norm(g) < tol).'
Iterations: 3
FirstOrderOptimality: 7.2436
FcnCount: 7
UpdateNorm: 0.0067
LastImprovement: 0.0067
若要获取有关最小化的在线信息，请使用 'Display' 估计选项，其可选值有 'off'、'on' 和 'full'。这将启动一个进度查看器，显示有关模型估计进度的信息。
Opt = armaxOptions('Display','on'); m1 = armax(z1,[2 2 2 1],Opt);
线性模型的变体 - IDTF、IDPOLY、IDPROC、IDSS 和 IDGREY
线性模型有几种类型。上文是一个适用于多项式类型模型的 idpoly 版本的示例。通过使用相应的估计量（如 bj, oe, armax, arx 等），可以得到多项式型模型的不同变体，例如 Box-Jenkins 模型、输出误差模型、ARMAX 模型等。所有这些都以 idpoly 对象的形式呈现。
其他变体包括：idss（用于状态空间模型）；idgrey（用于用户定义的结构化状态空间模型）；idtf（用于传递函数模型）；以及 idproc（用于过程模型（增益+延迟+静态增益））。
用于评估模型的命令：bode, step, iopzmap, compare 等，均直接作用于模型对象，例如：
compare(z1,m1)
通过 idssdata、tfdata 和 zpkdata 可得到状态空间、传递函数以及零点/极点的变换：
[num,den] = tfdata(m1,'v')
num =
0 0.9430 0.5224
den =
1.0000 -1.5312 0.7293
'v' 表示 num 和 den 将作为向量返回，而不是作为元胞数组返回。元胞数组在处理多变量系统时非常有用。若要同时获取 num 和 den 值的 1 个标准差的不确定度，请使用：
[num, den, ~, dnum, dden] = tfdata(m1,'v')
num =
0 0.9430 0.5224
den =
1.0000 -1.5312 0.7293
dnum =
0 0.0607 0.0782
dden =
0 0.0180 0.0147
将辨识模型转换为 Control System Toolbox 的数值线性时不变 (LTI) 模型
这些对象还可直接与 Control System Toolbox™ 模型对象（如 tf、ss 和 zpk）建立连接，并且如果 Control System Toolbox 可用，则可将其转换为这些 LTI 对象。例如，tf 将一个 idpoly 对象转换为一个 tf 对象。
CSTBInstalled = exist('tf','class')==8; if CSTBInstalled % check if Control System Toolbox is installed tfm = tf(m1) % convert IDPOLY model m1 into a TF object end
tfm =
From input "u1" to output "y1":
0.943 z^-1 + 0.5224 z^-2
----------------------------
1 - 1.531 z^-1 + 0.7293 z^-2
Sample time: 0.1 seconds
Discrete-time transfer function.
将 IDLTI 模型转换为 Control Systems Toolbox 中的 LTI 模型时，噪声组件不会被保留。若要将噪声通道也作为 LTI 模型的常规输入，请使用 'augmented' 标志：
if CSTBInstalled tfm2 = tf(m1,'augmented') end
tfm2 =
From input "u1" to output "y1":
0.943 z^-1 + 0.5224 z^-2
----------------------------
1 - 1.531 z^-1 + 0.7293 z^-2
From input "v@y1" to output "y1":
1.051 - 1.113 z^-1 + 0.2069 z^-2
--------------------------------
1 - 1.531 z^-1 + 0.7293 z^-2
Input groups:
Name Channels
Measured 1
Noise 2
Sample time: 0.1 seconds
Discrete-time transfer function.
在模型 tfm2 中，噪声通道的名称为 v@y1。
