通过辨识对复杂系统进行线性逼近
本示例演示了如何通过线性模型辨识，获得一个复非线性系统的线性逼近。该方法基于选择一个能激励系统的输入信号。通过将线性模型拟合到非线性模型对所选输入信号的仿真响应，可得到线性逼近。
本示例使用了 Simulink®, Control System Toolbox™ 和 Simulink Control Design™。
简介
在许多情况下，线性模型是通过在某些局部条件下对更复杂的非线性系统进行简化而得到的。例如，飞机动态特性的高保真模型可以通过一个详细的 Simulink 模型来描述。为了加快此类系统的仿真速度、研究其在工作点附近的局部行为或设计补偿器，一种常用的方法是对其进行线性化处理。如果我们围绕某个工作点对原始模型进行解析线性化，一般而言，所得模型的阶数会与原始模型的状态数相同或接近。对于该阶数在分析或控制系统设计中应适用的输入类型而言，该阶数可能过高。因此，我们可以考虑一种替代方法，该方法的核心是从系统仿真中收集投入产出数据，并利用这些数据推导出恰好具有正确阶数的线性模型。
F14 模型的解析线性化
以 F14 模型为例。这虽然已经是一个线性模型，但其中包含导数模块和可能影响其输出性质的扰动源。我们可以按照以下方式，将其在单个输入端口与两个输出端口之间进行“线性化”：
open_system('idF14Model') IO = getlinio('idF14Model') syslin = linearize('idF14Model',IO)
3x1 vector of Linearization IOs:
--------------------------
1. Linearization input perturbation located at the following signal:
- Block: idF14Model/Pilot
- Port: 1
- Signal Name: Stick Input
2. Linearization output measurement located at the following signal:
- Block: idF14Model/Gain5
- Port: 1
- Signal Name: Angle of Attack
3. Linearization output measurement located at the following signal:
- Block: idF14Model/Pilot G force
- Port: 1
- Signal Name: Pilot G force
syslin =
A =
Transfer Fcn Derivative Transfer Fcn Derivative1
Transfer Fcn -0.6385 0 689.4 0
Derivative 1 -1e+05 0 0
Transfer Fcn -0.00592 0 -0.6571 0
Derivative1 0 0 1 -1e+05
Actuator Mod 0 0 1.424 0
Alpha-sensor 0.001451 0 0 0
Stick Prefil 0 0 0 0
Pitch Rate L 0 0 1 0
Proportional 0 0 -0.8156 0
Actuator Mod Alpha-sensor Stick Prefil Pitch Rate L
Transfer Fcn -1280 0 0 0
Derivative 0 0 0 0
Transfer Fcn -137.7 0 0 0
Derivative1 0 0 0 0
Actuator Mod -20 2.986 -39.32 -1.67
Alpha-sensor 0 -2.526 0 0
Stick Prefil 0 0 -22.52 0
Pitch Rate L 0 0 0 -4.144
Proportional 0 -1.71 22.52 0.9567
Proportional
Transfer Fcn 0
Derivative 0
Transfer Fcn 0
Derivative1 0
Actuator Mod -3.864
Alpha-sensor 0
Stick Prefil 0
Pitch Rate L 0
Proportional 0
B =
Stick Input
Transfer Fcn 0
Derivative 0
Transfer Fcn 0
Derivative1 0
Actuator Mod 0
Alpha-sensor 0
Stick Prefil 1
Pitch Rate L 0
Proportional 0
C =
Transfer Fcn Derivative Transfer Fcn Derivative1
Angle of Att 0.001451 0 0 0
Pilot G forc -3106 3.106e+08 7.083e+04 -7.081e+09
Actuator Mod Alpha-sensor Stick Prefil Pitch Rate L
Angle of Att 0 0 0 0
Pilot G forc 0 0 0 0
Proportional
Angle of Att 0
Pilot G forc 0
D =
Stick Input
Angle of Att 0
Pilot G forc 0
Continuous-time state-space model.
syslin 是一个具有 2 个输出、1 个输入和 9 个状态的模型。这是因为在原始系统中，从“摇杆输入”输入端到两个输出端的线性化路径共有 9 种状态。我们可以使用 operpoint 来验证这一点：
operpoint('idF14Model')
ans = Operating point for the Model idF14Model. (Time-Varying Components Evaluated at time t=0) States: ---------- x _ (1.) idF14Model/Actuator Model 0 (2.) idF14Model/Aircraft Dynamics Model/Transfer Fcn.1 0 (3.) idF14Model/Aircraft Dynamics Model/Transfer Fcn.2 0 (4.) idF14Model/Controller/Alpha-sensor Low-pass Filter 0 (5.) idF14Model/Controller/Pitch Rate Lead Filter 0 (6.) idF14Model/Controller/Proportional plus integral compensator 0 (7.) idF14Model/Controller/Stick Prefilter 0 (8.) idF14Model/Dryden Wind Gust Models/Q-gust model 0 (9.) idF14Model/Dryden Wind Gust Models/W-gust model 0 0 Inputs: None ----------
在保持所选方波 ("Stick input") 输入响应精度的同时，能否减少该命令的执行次数？
准备身份辨识数据
对模型进行仿真，并在 0:30 秒的时间范围内记录输入方波 (u) 以及输出 "Angle of attack" (y1) 和 "Pilot G force" (y2)。这些数据经插值处理后，转换为间隔均匀的时间向量（采样时间为 0.0444 秒），并存储在 "idF14SimData.mat" 文件中。
load idF14SimData Z = iddata([y1, y2],u,Ts,'Tstart',0); Z.InputName = 'Stick input'; Z.InputUnit = 'rad/s'; Z.OutputName = {'Angle of attack', 'Pilot G force'}; Z.OutputUnit = {'rad', 'g'}; t = Z.SamplingInstants; subplot(311) plot(t,Z.y(:,1)), ylabel('Angle of attack (rad)') title('Logged Input-Output Data') subplot(312) plot(t,Z.y(:,2)), ylabel('Pilot G force (g)') subplot(313) plot(t,Z.u), ylabel('Stick input (rad/s)') axis([0 30 -1.2 1.2]) xlabel('Time (seconds)')
状态空间模型的估计
使用 ssest 命令估计 2 到 4 阶的状态空间模型。我们将估计配置为使用“仿真”重点，并选择不估计模型的扰动分量。
opt = ssestOptions('Focus','simulation'); syslin2 = ssest(Z, 2, 'DisturbanceModel', 'none', opt); syslin3 = ssest(Z, 3, 'DisturbanceModel', 'none', opt); syslin4 = ssest(Z, 4, 'DisturbanceModel', 'none', opt);
将线性化模型 syslin 与三个辨识模型的拟合效果与数据进行比较。请注意，syslin 属于 SS 模型，而 syslin2、syslin3 和 syslin4 属于 IDSS 模型。
syslin.InputName = Z.InputName;
syslin.OutputName = Z.OutputName; % reconcile names to facilitate comparison
clf
compare(Z, syslin, syslin2, syslin3, syslin4)
图中显示，三阶模型 (syslin3) 作为飞机在默认 (t=0) 运行条件下的动力学线性逼近，表现相当不错。它与数据的吻合度略好于通过解析线性化 (syslin) 得到的结果。如果原始模型 idF14Model 是线性的，那么为什么其线性化结果 syslin 无法与数据实现 100% 的拟合？这有两个原因：
测量到的输出值会受到阵风的影响，这意味着记录下来的输出值并非仅仅是 Stick 输入函数。有一些干扰因素正在影响它。
"Pilot G force" 模块使用了导数模块，其线性化取决于时间常数 "c" 的值。"c" 应取较小的数值（我们采用 1e-5），但不能为零。c 的非零值会在线性化过程中引入逼近误差。
让我们来看一下模型 syslin3 的参数，该模型似乎很好地捕捉到了响应情况：
syslin3
syslin3 =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3
x1 -1.006 2.029 0.5842
x2 -8.284 -19.39 5.611
x3 2.784 12.63 -6.956
B =
Stick input
x1 0.2614
x2 -5.512
x3 3.606
C =
x1 x2 x3
Angle of att -8.841 0.5347 1.402
Pilot G forc -86.42 15.85 66.12
D =
Stick input
Angle of att 0
Pilot G forc 0
K =
Angle of att Pilot G forc
x1 0 0
x2 0 0
x3 0 0
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: none
Number of free coefficients: 18
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using SSEST on time domain data "Z".
Fit to estimation data: [98.4;97.02]%
FPE: 2.367e-05, MSE: 0.1103
通过简化与估计实现模型简化
我们还可以采用以下方法：降低线性化模型 syslin 的阶数，并优化降阶模型的参数，使其能最好地拟合数据 Z。为了确定一个合适的降阶值，我们使用 hsvd：
[S, BalData] = hsvd(syslin); clf; bar(S)
柱状图显示，对于第 4 态及更高态，奇异值都相当小。因此，三阶可能是在简化过程中最优的选择。
sysr = balred(syslin,3,BalData)
opt2 = bodeoptions; opt2.PhaseMatching = 'on';
clf; bodeplot(sysr,syslin,opt2)
sysr =
A =
x1 x2 x3
x1 -2.854 -7.61 -54.04
x2 -0.9714 2.341 9.123
x3 0.6979 -7.203 -24.08
B =
Stick input
x1 -137.7
x2 -869
x3 -506.7
C =
x1 x2 x3
Angle of att -0.0005063 -0.0008826 -0.001016
Pilot G forc -0.005926 -0.04692 -0.1646
D =
Stick input
Angle of att -0.03784
Pilot G forc -1.617
Continuous-time state-space model.
波德图显示，在 10 rad/s 以内，模型再现性良好。如 compare 图所示，sysr 能够模拟出与原始 9 态模型同样良好的响应：
compare(Z, sysr, syslin)
让我们对 sysr 的参数进行调整，以改善其对数据的拟合效果。在此估计中，我们选择了“莱文伯格-马夸特”搜索法，并将允许的最大迭代次数设为 10。这些选择是经过一番试错后做出的。我们还启用了估计进度显示功能。
opt.Display = 'on'; opt.SearchMethod = 'lm'; opt.SearchOptions.MaxIterations = 10; sysr2 = ssest(Z, sysr, opt) compare(Z, sysr2)
sysr2 =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3
x1 -4.048 -7.681 -54.01
x2 -0.4844 1.549 8.895
x3 -0.2398 -6.777 -25.78
B =
Stick input
x1 -137.7
x2 -869
x3 -506.7
C =
x1 x2 x3
Angle of att -0.0003361 -0.0004964 0.00215
Pilot G forc -0.01191 -0.03599 -0.1434
D =
Stick input
Angle of att 0.003022
Pilot G forc 0.6438
K =
Angle of att Pilot G forc
x1 0 0
x2 0 0
x3 0 0
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: yes
Disturbance component: none
Number of free coefficients: 20
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using SSEST on time domain data "Z".
Fit to estimation data: [98.78;97.03]%
FPE: 1.434e-05, MSE: 0.1097
优化后的模型 sysr2 与 F14 模型的响应吻合得相当好（第一个输出结果的吻合度约为 99%，第二个输出结果的吻合度约为 97%）。
结论
我们提出了一种替代的解析线性化方法，用于获得复杂系统的线性逼近。这些结果是针对特定输入信号得出的，严格来说，仅适用于该输入信号。为了提高研究结果对各种输入特征的适用性，我们可以使用不同类型的输入进行多次仿真。然后，我们可以将生成的数据集合并为一个多试验数据集（请参阅 iddata/merge），并将其用于估计。为了方便起见，在这个示例中，我们使用了一个复数系统，但该系统是线性的。这种方法的真正优势将在非线性系统中显现出来。
我们还提出了一种方法，可以在降低线性系统阶数的同时，确保简化后的模型能够忠实再现原始 Simulink 模型的仿真响应。简化模型 sysr 的作用是为估计模型 sysr2 提供一个初始估计值。该方法还突显了一个事实，即任何线性系统（包括不同类别的线性系统）均可作为估计的初始模型。
bdclose('idF14Model')
