【MIMO稳定裕度】基于数据驱动的多输入多输出系统稳定裕度分析
最近一直在忙着写论文，只能说要写一篇高水平论文确实不容易，要一直反复来回修改调整，要求比较高，所以没太有时间和精力写博客，这两天结束了初稿，又正好是假期，出来冒个泡。
本次分享的主题是：多输入多输出系统的稳定裕度求解方法。
0. 问题描述
对于多输入多输出系统，经典的Bode图分析稳定裕度不再适用，目前市面上对于多输入多输出系统用的比较广泛的办法有：
- 回差阵奇异值法；
- μ分析法；
其中回差阵奇异值法相对比较保守，但是两种方法都是基于模型的，也就是必须得有建立好的状态空间模型才能求解稳定裕度。然而，最近却了解到一种不基于模型，基于实验数据来解算稳定裕度的方法，并且求解过程简单，感觉在实际中更加有效。以下首先介绍基于模型的SM求解算法，然后介绍基于数据驱动的SM求解方法，最后再补充介绍一下μ分析法求解SM法。
1. 基于模型的MIMO稳定裕度分析
对于多输入多输出闭环系统，在被控对象输入端注入一个复对角乘性扰动 I m + Δ ( s ) \mathbf{I}_m+\mathbf{\Delta }\left( s \right) Im+Δ(s) 来定义多变量增益和时延裕度（GM and DM）。
当求解增益裕度时，
Δ ( s ) = d i a g ( r 1 , r 2 , . . . , r m ) , ∣ r i ∣ ⩽ r max , i = 1 , 2 , . . . m \mathbf{\Delta }\left( s \right) =\mathrm{diag}\left( r_1, r_2, ... , r_m \right) , \left| r_i \right|\leqslant r_{\max}, i=1,2,...m Δ(s)=diag(r1,r2,...,rm),∣ri∣⩽rmax,i=1,2,...m当求解时延裕度时， Δ ( s ) = d i a g ( e − s τ 1 − 1 , e − s τ 2 − 1 , . . . , e − s τ m − 1 ) , 0 ⩽ τ i ⩽ τ max , i = 1 , 2 , . . . m \mathbf{\Delta }\left( s \right) =\mathrm{diag}\left( e^{-s\tau _1-1}, e^{-s\tau _2-1}, ... , e^{-s\tau _m-1} \right) , 0\leqslant \tau _i\leqslant \tau _{\max}, i=1,2,...m Δ(s)=diag(e−sτ1−1,e−sτ2−1,...,e−sτm−1),0⩽τi⩽τmax,i=1,2,...m对上图利用环路变换和小增益定理进行转换可得 
R i ( s ) = ( I m + ( H ( s ) K G ( s ) ) − 1 ) − 1 \mathbf{R}_{\mathrm{i}}\left( s \right) =\left( \mathbf{I}_m+\left( \mathbf{H}\left( s \right) \mathbf{KG}\left( s \right) \right) ^{-1} \right) ^{-1} Ri(s)=(Im+(H(s)KG(s))−1)−1此时的增益裕度和时延裕度的最大值为 r max = ( 1 − ε 1 ) ∥ R i ( s ) ∥ − 1 τ max = ( 1 − ε 2 ) ∥ s R i ( s ) ∥ − 1 r_{\max}=\left( 1-\varepsilon _1 \right) \left\| \mathbf{R}_i\left( s \right) \right\| ^{-1} \\ \tau _{\max}=\left( 1-\varepsilon _2 \right) \left\| s\mathbf{R}_i\left( s \right) \right\| ^{-1} rmax=(1−ε1)∥Ri(s)∥−1τmax=(1−ε2)∥sRi(s)∥−1对于 H ∝ H_{\propto} H∝范数 ∥ R i ( s ) ∥ \left\| R_i\left( s \right) \right\| ∥Ri(s)∥求通过sigma函数求解最大奇异值得到，然而可转换成对应的增益裕度和时延裕度。 d b 2 m a g ( − r max + 1 ) ⩽ G M ⩽ d b 2 m a g ( r max + 1 ) 0 ⩽ D M ⩽ τ max db2mag\left( -r_{\max}+1 \right) \leqslant GM\leqslant db2mag\left( r_{\max}+1 \right) \\ 0\leqslant DM\leqslant \tau _{\max} db2mag(−rmax+1)⩽GM⩽db2mag(rmax+1)0⩽DM⩽τmax在验证上述结果的时候，直接在反馈回路上加上增益环节和时延环节即可，-r*[1 0;0 1]+[1 0;0 1]，经过验证，该裕度取值合适，较不保守。
2. 基于数据驱动的MIMO稳定裕度分析
以上方法是基于有确切模型的时候求解稳定裕度，然而对于实际情况来说，大部分时候是不知道具体模型的，此时使用上述方法得到的稳定裕度就会有较大的偏差。
对上一节的控制回路进行调整，可得到基于扫频数据的稳定裕度求解方法。
上图表示的意思是，通过在z点注入扫频信号，在w点获取对应的输出数据，此时借助频率响应综合识别软件包（CIFER） 的帮助即可得到对应的 R i ( s ) \boldsymbol{R}_i\left( s \right) Ri(s)，具体的求解步骤如下：
最终得到奇异值曲线如下所示，对最大奇异值根据上节的公式进行转换，可得到对应的增益裕度为6.17dB，时延裕度为0.1618s。经过验证，数值合理，较不保守。
上述方法的优势之处为：
- 简单；
- 无模型基于实验数据获得；
- 在不中断环路的情况下进行在线评估SM；
- 对于含时滞环节的不确定性能更好的处理；
劣势在于：虽然扫频试验得到的数据驱动SMs具有一定的保守性，但在所得到的数据驱动SMs的扰动下，系统仍能保持稳定。
3. μ分析法求解稳定裕度
顺带补充一下μ分析法求解稳定裕度的方法，同样在输入端引入乘性不确定性：
将不确定性从系统中分离出来，并转换为如下形式：
其中 P = − ( I + K G ) − 1 K G \boldsymbol{P}=-\left( \boldsymbol{I}+\boldsymbol{KG} \right) ^{-1}\boldsymbol{KG} P=−(I+KG)−1KG 这样，原系统转化为适合鲁棒稳定性分析的标准形式。如果知道了回路开环传递函数KG ，就可以通过该式计算出P(s)的频率特性，进而利用Matlab中的 μ-tool得到系统的结构奇异曲线。并找出系统的μ峰值。
上图μ峰值为3.75，转换成经典的增益裕度和相角裕度分别为2.0492dB和15.2904°。经过验证，此结果合理，系统在此范围内仍能稳定。
4. 总结
1、分别阐述了基于理论模型和基于实验数据的稳定裕度计算方法；
2、可根据实际情况合理选择计算方法；在理论上可以采用基于模型的计算方法（小增益 L 2 L_2 L2方法和μ分析法），在实践上可以采用基于扫频数据驱动的测量方法。