项目调度问题的一些matlab开发的工具箱

Matlab鲁棒控制工具箱（Robust Control Toolbox）

目 录
一、引言 2
1、工具箱函数简介 2
1）不确定元素 2
2）不确定矩阵和系统 2
3）不确定模型的控制 3
4）不确定模型间的互连 3
5）模型降阶 3
6）鲁棒性和最坏情况分析 4
7）参数依赖系统的鲁棒性分析 4
8）控制器综合 4
9）μ综合 5
10）采样系统 5
11）增益调度 5
12）频域响应模型 5
13）公用函数 6
14）LMI函数 6
15）LMI特征 6
16）LMI求解 6
17）结果验证 6
18）修改LMI对象 7
2、不确定性建模 7
3、最坏情况下的性能分析 9
4、MIMO系统的鲁棒控制器设计 11
5、模型降阶及逼近 13
6、作者简介 16
7、参考文献 17
代码1 ACC Benchmark Problem 18
代码2 NASA HiMAT 18
二、多变量回路成形设计 19
三、模型降阶 20
1、Hankel奇异值 20
2、加性误差方法和乘性误差方法 22
3、modreal函数 26
4、ncfmr函数 28
5、参考文献 30
四、鲁棒性分析 30
1、不确定性建模 30
2、鲁棒性分析 34
3、MIMO系统的鲁棒性分析 38
4、最坏情况下增益分析 43




















Matlab鲁棒控制工具箱（Robust Control Toolbox）

鲁棒控制工具箱提供了一系列的函数和工具以支持带有不确定元素的多输入多输出控制系统的设计。在该工具箱的帮助下，你可以建立带有不确定参数和动态特性的LTI模型，也可以分析MIMO系统的稳定性裕度和最坏情况下的性能。
该工具箱提供了一系列的控制器分析和综合函数，能够分析最坏情况下的性能及确定最坏情况下的参数值。利用模型降阶函数能够对复杂模型进行简化。同时提供了先进的鲁棒控制方法，如H2、H∞、LMI、μ分析等。
一、引言
1、工具箱函数简介
1）不确定元素
ucomplex 创建不确定复数参量
ucomplexm 创建不确定复数矩阵
udyn 创建未定义结构的不确定动态系统对象
ultidyn 创建线性时不变不确定性对象
ureal 创建不确定实数参量
2）不确定矩阵和系统
diag 对不确定矩阵和系统对角化
randatom 创建随机不确定性atom对象
randumat 创建随机不确定性矩阵
randuss 创建随机不确定性状态空间模型
ufrd 创建不确定性频域响应数据对象（ufrd），或者将其它模型转
换为ufrd对象
umat 创建不确定性矩阵
uss 定义不确定性状态空间模型，或者将LTI对象转换为不确定
性状态空间模型
rss 创建随机稳定连续时间状态空间模型
3）不确定模型的控制
actual2normalized 对于给定的atom对象，计算其与标准值间的归一化距离
gridureal 将ureal对象均匀网格化
isuncertain 判断是否为不确定性系统
lftdata 将不确定对象分解为固定的规范型和固定的不确定性部分
normalized2actual 将正规化坐标系

中的atom值转换为实际值
repmat 复制和命名不确定矩阵
simplify 简化不确定对象的表达式
squeeze 将umat对象去掉一维
usample 产生随机不确定对象
usubs 替换不确定对象中的不确定参数
uss/ssbal 将不确定系统中的状态/不确定性标量化
4）不确定模型间的互连
iconnect 创建空白的互连对象
icsignal 创建icsignal对象
imp2exp 将线性关系转换为输入-输出关系
stack 将不确定矩阵、模型或者数组压入堆中
sysic 互连确定与不确定矩阵或系统
5）模型降阶
balancmr 利用均方根法求降阶模型
bstmr 利用Schur法求随机模型的降阶模型
hankelmr 为降阶前的最小Hankel测度
hankelsv 计算稳定/不稳定系统、连续/离散系统的Hankel奇异值
imp2ss 将脉冲响应模型转换为近似状态空间模型
modreal 模态形式的实现
ncfmr 归一化的互质降阶模型
reduce 利用Hankel奇异值法降阶
schurmr 利用Schur法求降阶模型
slowfast 将状态空间模型按照快-慢原则分解
stabsep 将状态空间模型按照稳定/不稳定原则分解
6）鲁棒性和最坏情况分析
cpmargin 计算包含反馈控制器的闭环系统的稳定性裕度
gapmetric 计算两个系统间的gap、nugap距离上限
loopmargin 分析反馈回路
loosens 分析包含反馈控制器的闭环系统的灵敏度
mussv 计算所构造的奇异值（μ）的界限
mussvextract 从mussv创建的数据结构中提出muinfo对象
ncfmargin 计算反馈回路的归一化稳定性裕度
popov 执行Popov鲁棒性检验
robopt 创建可供robuststab和robustperf选择的对象
robustperf 计算不确定多变量系统的鲁棒性能裕度
robuststab 计算不确定多变量系统的鲁棒稳定性裕度
wcgain 计算不确定系统的最坏情况下的增益界限
wcgopt 创建可供wcgain、wcsens和wcmargin使用的对象
wcmargin 计算最坏情况下反馈系统的增益/幅值裕度
wcnorm 计算最坏情况下不确定矩阵的范数
wcsens 计算最坏情况下反馈回路的灵敏度和补灵敏度函数
7）参数依赖系统的鲁棒性分析
aff2pol 将仿射参数依赖系统（Parameter-Dependent Systems，P-系统）
转换为多胞模型
decay 多胞模型或者P-系统的二次衰减速率
ispsys 是否为P-系统
pdlstab 评估多胞模型或者P-系统的鲁棒稳定性
pdsimul P-系统沿参数变化轨迹的时域响应
polydec 计算多胞坐标
psinfo 多胞模型或者P-系统的设置参数
psys 设置线性的多胞或者参数依赖系统
pvec 设置不确定性向量或者时变参数
pvinfo 参数向量的参数
quadperf 计算多胞模型或者P-系统的二次型H∞指标
guadstab 评估多胞模型或者P-系统的二次型稳定性
8）控制器综合
augw 为加权混合灵敏度回路成形设计创

建增广系统模型
h2hinfsyn 极点配置约束下的混合H2/H∞设计
h2syn 针对LTI模型设计H2控制器
hinfsyn 针对LTI模型设计H∞最优控制器
loopsyn H∞回路成形控制器设计
ltrsyn LQG回路传递函数恢复控制器设计
mkfilter 创建Bessel，Butterworth，Chebyshev，RC滤波器
mixsyn H∞混合灵敏度控制器设计
ncfsyn H∞归一化互质因子控制器设计
sigma 画出LTI反馈回路的奇异值
makeweight H∞混合灵敏度中的权重（mixsyn，augw）
9）μ综合
cmsclsyn 常值矩阵μ综合
dksyn 利用μ综合D-K迭代法进行鲁棒控制器设计
dkitopt 创建供dksyn使用的对象
drawmag 鼠标交互操作函数
fitfrd 针对状态空间模型，合理化D标度的频域响应
fitmagfrd 针对稳定、最小相位的状态空间模型，合理化标量幅值数据
genphase 针对单输入单输出、实有理、最小相位传递函数，合理化的
幅值响应数据
10）采样系统
sdhinfnorm 计算采样系统的L2范数
sdhinfsyn 采样系统的H∞控制器设计
sdlsim 带有反馈回路的采样系统的时域响应
11）增益调度
hinfgs 增益调度H∞控制器设计
12）频域响应模型
frd/loglog frd对象的log-log度量
frd/semilogx frd对象的半对数（semilog）度量
frd/rcond frd对象的互逆条件判断
frd/schur frd对象的Schur分解
frd/svd frd对象的奇异值分解
13）公用函数
bilin 多变量频域双线性变换（s域、z域）
dmplot 对增益和相位裕度进行说明
mktito 将LTI系统分解为两输入-两输出系统
sectf 状态空间模型的双线性变换
skewdec 创建反对称矩阵
symdec 创建对称矩阵
14）LMI函数
getlmis LMI系统的内部描述
lmiedit 用matlab语言编辑或显示LMI系统
lmiterm 增加新项给现有的LMI对象
lmivar 在现有的LMI系统中设定矩阵变量
newlmi 增加新的LMI对象给现有的LMI系统
setlmis 初始化LMI系统
15）LMI特征
dec2mat 从决策变量向量中提取出矩阵值
decinfo 描述矩阵变量和决策变量中的联系
decnbr LMI系统中的决策变量个数
lmiinfo 现有LMI对象的信息
lminbr LMI系统中的LMI对象的个数
mat2dec 从矩阵中提取决策变量向量
matnbr LMI系统中矩阵变量的个数
16）LMI求解
defcx 为mincx对象设定cTx对象
feasp 求解给定的LMI系统
gevp LMI约束下的广义特征值
mincx LMI约束下最小化线性对象
17）结果验证
evallmi 针对给定的决策变量评估LMI
showlmi 待评估LMI对象的左边项和右边项
18）修改LMI对象
dellmi 从现有LMI系统中移除LMI对象
delmvar 从带求解的LMI问题中移除矩阵变量
setmvar 举例矩阵变量和待评估的LMI项
2、不确定性建模
鲁棒控制的核心是对不确定LTI系统的建模。系统的不确定性来

源于系统参数不能精确获得，或者变化范围很大，比如系统零极点位置未知、增益未知。也可能存在结构未知。
利用鲁棒控制工具箱可以利用matlab对象构造不确定LTI系统。
例1：ACC Bechmark Problem
考虑如图1所示的Bechmark Problem，两辆卡车由弹簧连接，不考虑摩擦力。

图1 ACC Benchmark Problem
图1的结构框图见图2所示，图2中的一些传递函数为

其中，参数m1,m2,k具有20%的不确定性，即
m1 = 1±0.2，m2 = 1±0.2，k = 1±0.2

图2 "ACC Benchmark" Two-Cart System Block Diagram y1 = P(s) u1
最里面的虚线框传递函数矩阵为

其中，输入为u1, u2, 输出为y1, y2。
利用matlab对不确定系统P建模如下：
% Create the uncertain real parameters m1, m2, & k
m1 = ureal('m1',1,'percent',20);
m2 = ureal('m2',1,'percent',20);
k = ureal('k',1,'percent',20);

s = zpk('s'); % create the Laplace variable s
G1 = ss(1/s^2)/m1; % Cart 1
G2 = ss(1/s^2)/m2; % Cart 2

% Now build F and P
F = [0;G1]*[1 -1]+[1;-1]*[0,G2];
P = lft(F,k) % close the loop with the spring k
系统P为SISO、不确定、状态空间模型，具有4个状态和3个不确定参数，m1,m2和k，利用下面的命令得到标称模型：
zpk(P.nominal)
得到返回结果：
Zero/pole/gain:
1
--------------
s^2 (s^2 + 2)
如果不确定系统P与一LTI控制器C连接，

连接形式见图3所示。

图3 控制器连接
得到闭环控制系统y1 = T(s) u1 。针对三个不确定参数m1,m2和k进行五次蒙特卡洛仿真，观察0~0.1s内的控制结果，程序如下：
% Create the uncertain real parameters m1, m2, & k
m1 = ureal('m1',1,'percent',20);
m2 = ureal('m2',1,'percent',20);
k = ureal('k',1,'percent',20);

s = zpk('s'); % create the Laplace variable s
G1 = ss(1/s^2)/m1; % Cart 1
G2 = ss(1/s^2)/m2; % Cart 2

% Now build F and P
F = [0;G1]*[1 -1]+[1;-1]*[0,G2];
P = lft(F,k) % close the loop with the spring k

zpk(P.nominal)
C=100*ss((s+1)/(.001*s+1))^3 % LTI controller
T=feedback(P*C,1); % closed-loop uncertain system
step(usample(T,5),.1);
结果如图4所示：

图4 Monte Carlo Sampling of Uncertain System's Step Response
该例完整代码见代码1。
3、最坏情况下的性能分析
针对不确定参数可能的变化范围，所设计的控制系统应该都能够满足稳定性和性能指标的要求。利用usample指令能够针对不确定参数进行蒙特卡洛仿真，但不可能将所有情况都仿真到。为了得到最坏情况下的性能，利用蒙特卡洛仿真次数会非常大。
鲁棒控制提供了强有力的鲁棒性能分析工具，能够直接计算最坏情况下的性能上下限，而不需要进行随机采样仿真。点击最坏情况下鲁棒性分析指令。
参考例1，闭环系统为
T=feedback(P*C,1); % closed-loop uncertain system
不确定状态空间模型T有三个不确定性参数：m1,m2和k，每个

参数波动范围为1±20%。利用下列指令检验闭环系统T是否对这三个参数的所有组合具有鲁棒稳定性：
[StabilityMargin,Udestab,REPORT] = robuststab(T);
REPORT
得到下面的报告结果：
Uncertain System is robustly stable to modeled uncertainty.
-- It can tolerate up to 308% of the modeled uncertainty.
-- A destabilizing combination of 500% of the modeled uncertainty exists,
causing an instability at 0.56 rad/s.
-- Sensitivity with respect to uncertain element ...
'k' is 22%. Increasing 'k' by 25% leads to a 6% decrease in the margin.
'm1' is 57%. Increasing 'm1' by 25% leads to a 14% decrease in the margin.
'm2' is 61%. Increasing 'm2' by 25% leads to a 15% decrease in the margin.
报告说明了对于20%的波动，系统是鲁棒稳定的。Udestab指令返回使系统稳定性降低5倍的不确定参数的组合：
Udestab =
k: 1.1364e-010
m1: 1.1364e-010
m2: 1.9992
通过上面的仿真得到，在闭环系统稳定的前提下，不确定参数波动范围能够达到308%~500%，远远超过之前预期的20%。在20%的严格约束下，由于参数变动对闭环系统性能的影响是怎样的？下面的代码计算最坏情况下T的峰值增益，以及峰值增益时的频率和参数值：
[PeakGain,Uwc] = wcgain(T);
Twc=usubs(T,Uwc); % Worst case closed-loop system T
Trand=usample(T,4); % 4 random samples of uncertain system T
bodemag(Twc,'r',Trand,'b-.',{.5,50}); % Do bode plot
legend('T_{wc} - worst-case',...
'T_{rand} - random samples',3);


图5 Uncertain System Closed-Loop Bode Plots
4、MIMO系统的鲁棒控制器设计
利用鲁棒控制工具箱提供的指令函数，能够对MIMO LTI模型进行鲁棒控制器的设计。
例2 利用loopsyn设计控制器
loopsyn是一种功能强大但设计过程简洁的控制器设计函数。对于给定的LTI模型，事先指定期望的开环系统的频域响应波形，loopsyn能够计算一稳定的控制器来最佳逼近指定的波形。
以文献[8]中的两输入两输出的NASA HiMAT飞机模型为例，见图6所示。

图6 Aircraft Configuration and Vertical Plane Geometry
输入变量为 和 ，输出变量为α和θ，6个状态变量分别为

该模型的状态空间方程用matlab语言表示为：
% NASA HiMAT model G(s)
ag =[ -2.2567e-02 -3.6617e+01 -1.8897e+01 -3.2090e+01 3.2509e+00 -7.6257e-01;
9.2572e-05 -1.8997e+00 9.8312e-01 -7.2562e-04 -1.7080e-01 -4.9652e-03;
1.2338e-02 1.1720e+01 -2.6316e+00 8.7582e-04 -3.1604e+01 2.2396e+01;
0 0 1.0000e+00 0 0 0;
0 0 0 0 -3.0000e+01 0;
0 0 0 0 0 -3.0000e+01];
bg = [ 0 0;
0 0;

0 0;
0 0;
30 0;
0 30];
cg = [ 0 1 0 0 0 0;
0 0 0 1 0 0];
dg = [ 0 0;
0 0];
G=ss(ag,bg,cg,dg);
指定系统的带宽为10rad/s，于是期望传递函数为Gd(s)=10/s，然后loopsyn(G,Gd)得到回路成形控制器G，以满足期望的回路波形Gd。
s=zpk('s'); w0=10; Gd=w0/(s+.001);
[K,CL,GAM]=loopsyn(G,Gd); % Design a loop-shaping controller K

% Plot the results
sigma(G*K,'r',Gd,'k-.',Gd/GAM,'k:',Gd*GAM,'k:',{.1,30})
figure ;T=feedback(G*K,eye(2));
sigma(T,ss(GAM),'k:',{.1,30});grid

图7 MIMO Robust Loop Shaping with loopsyn(G,Gd)
得到的 说明了得到的最优波形与期望波形的逼近程度，它是闭环传递函数T=feedback(G*K,eye(2))的峰值增益。此例中，γ=1.6023=4db。所设计的回路波形与期望波形的差别在4db内。
此例代码见代码2。
5、模型降阶及逼近
有时候复杂的模型并不需要精确的控制。一些优化方法（例如基于H2、H∞、μ综合等）得到的控制器阶数至少等于控制对象状态变量的数量。鲁棒控制工具箱提供了模型降阶工具，以便得到原复杂模型和控制器的低阶近似模型。点击见鲁棒控制工具箱提供的模型降阶工具函数。
最重要的模型降阶方法是最小化加性、乘性、规范化的互质因子模型误差。这三种方法可以用reduce指令来实现。
以上节中的NASA HiMAT飞机模型为例，该模型含有8个状态，得到的最优回路成形控制器有16个状态。利用模型降阶，在不影响稳定性及闭环系统性的情况下去掉一些状态变量。
针对上节例子，在命令行中输入：
hankelsv(K,'ncf','log');
得到对数坐标系下的NCF Hankel奇异值，见图8所示。

图8 Hankel Singular Values of Coprime Factorization of K
理论证明，可以去掉含有远小于ncfmargin(G,K)的NCF Hankel奇异值的控制器状态变量而不会造成不稳定。
计算ncfmargin(G,K)，并加入到所画的Hankel奇异值图中：
hankelsv(K,'ncf','log');v=axis;
hold on; plot(v(1:2), ncfmargin(G,K)*[1 1],'--'); hold off

图9 HiMAT控制器K的16个NCF Hankel奇异值中有4个小于ncfmargin(G,K)
在此例中，能够去掉控制器K中16个状态变量中的4个，计算降阶后的包含12个状态变量控制器的指令为：
K1=reduce(K,11,'errortype','ncf');
结果用以下指令画出：
sigma(G*K1,'b',G*K,'r--',{.1,30});

图10 HiMAT with 11-State Controller K1 vs. Original 16-State Controller K
6、作者简介
Professor Andy Packard is with the Faculty of Mechanical Engineering at the University of California, Berkeley. His research interests include robustness issues in control analysis and design, linear algebra and numerical algorithms in control problems, applications of system theory to aerospace problems, flight control, and control of fluid flow.
Professor Gary Balas is with the Facu

lty of Aerospace Engineering & Mechanics at the University of Minnesota and is president of MUSYN Inc. His research interests include aerospace control systems, both experimental and theoretical.
Dr. Michael Safonov is with the Faculty of Electrical Engineering at the University of Southern California. His research interests include control and decision theory.
Dr. Richard Chiang is employed by Boeing Satellite Systems, El Segundo, CA. He is a Boeing Technical Fellow and has been working in the aerospace industry over 25 years. In his career, Richard has designed 3 flight control laws, 12 spacecraft attitude control laws, and 3 large space structure vibration controllers, using modern robust control theory and the tools he built in this toolbox. His research interests include robust control theory, model reduction, and in-flight system identification. Working in industry instead of academia, Richard serves a unique role in our team, bridging the gap between theory and reality.
The linear matrix inequality (LMI) portion of Robust Control Toolbox? software was developed by these two authors:
Dr. Pascal Gahinet is employed by The MathWorks. His research interests include robust control theory, linear matrix inequalities, numerical linear algebra, and numerical software for control.
Professor Arkadi Nemirovski is with the Faculty of Industrial Engineering and Management at Technion, Haifa, Israel. His research interests include convex optimization, complexity theory, and nonparametric statistics.









7、参考文献
[1] Boyd, S.P., El Ghaoui, L., Feron, E., and Balakrishnan, V., Linear Matrix Inequalities in Systems and Control Theory, Philadelphia, PA, SIAM, 1994.
[2] Dorato, P. (editor), Robust Control, New York, IEEE Press, 1987.
[3] Dorato, P., and Yedavalli, R.K. (editors), Recent Advances in Robust Control, New York, IEEE Press, 1990.
[4] Doyle, J.C., and Stein, G., "Multivariable Feedback Design: Concepts for a Classical/Modern Synthesis," IEEE Trans. on Automat. Contr., 1981, AC-26(1), pp. 4-16.
[5] El Ghaoui, L., and Niculescu, S., Recent Advances in LMI Theory for Control, Philadelphia, PA, SIAM, 2000.
[6] Lehtomaki, N.A., Sandell, Jr., N.R., and Athans, M., "Robustness Results in Linear-Quadratic Gaussian Based Multivariable Control Designs," IEEE Trans. on Automat. Contr., Vol. AC-26, No. 1, Feb. 1981, pp. 75-92.
[7] Safonov, M.G., Stability and Robustness of Multivariable Feedback Systems, Cambridge, MA, MIT Press, 1980.
[8] Safonov, M.G., Laub, A.J., and Hartmann, G., "Feedback Properties of Multivariable Systems: The Role and Use of Return Difference Matrix," IEEE Trans. of Automat. Contr., 1981, AC-26(1), pp. 47-65.
[9] Safonov, M.G., Chiang, R.Y., and Flashner, H., "H∞ Control Synthesis for a Large Space Structure," Proc. of American Contr. Conf., Atlanta, GA, June 15-17, 1988.
[10] Safonov, M.G., and Chiang, R.Y., "CACSD Using the State-Space L∞ Theory -- A Design Example," IEEE Trans. on

Automatic Control, 1988, AC-33(5), pp. 477-479.
[11] Sanchez-Pena, R.S., and Sznaier, M., Robust Systems Theory and Applications, New York, Wiley, 1998.
[12] Skogestad, S., and Postlethwaite, I., Multivariable Feedback Control, New York, Wiley, 1996.
[13] Wie, B., and Bernstein, D.S., "A Benchmark Problem for Robust Controller Design," Proc. American Control Conf., San Diego, CA, May 23-25, 1990; also Boston, MA, June 26-28, 1991.
[14] Zhou, K., Doyle, J.C., and Glover, K., Robust and Optimal Control, Englewood Cliffs, NJ, Prentice Hall, 1996.

代码1 ACC Benchmark Problem
% Create the uncertain real parameters m1, m2, & k
m1 = ureal('m1',1,'percent',20);
m2 = ureal('m2',1,'percent',20);
k = ureal('k',1,'percent',20);

s = zpk('s'); % create the Laplace variable s
G1 = ss(1/s^2)/m1; % Cart 1
G2 = ss(1/s^2)/m2; % Cart 2

% Now build F and P
F = [0;G1]*[1 -1]+[1;-1]*[0,G2];
P = lft(F,k) % close the loop with the spring k

zpk(P.nominal)
C=100*ss((s+1)/(.001*s+1))^3 % LTI controller
T=feedback(P*C,1); % closed-loop uncertain system
step(usample(T,5),.1);

[StabilityMargin,Udestab,REPORT] = robuststab(T);
REPORT
代码2 NASA HiMAT
% NASA HiMAT model G(s)
ag =[ -2.2567e-02 -3.6617e+01 -1.8897e+01 -3.2090e+01 3.2509e+00 -7.6257e-01;
9.2572e-05 -1.8997e+00 9.8312e-01 -7.2562e-04 -1.7080e-01 -4.9652e-03;
1.2338e-02 1.1720e+01 -2.6316e+00 8.7582e-04 -3.1604e+01 2.2396e+01;
0 0 1.0000e+00 0 0 0;
0 0 0 0 -3.0000e+01 0;
0 0 0 0 0 -3.0000e+01];
bg = [ 0 0;
0 0;
0 0;
0 0;
30 0;
0 30];
cg = [ 0 1 0 0 0 0;
0 0 0 1 0 0];
dg = [ 0 0;
0 0];
G=ss(ag,bg,cg,dg);

s=zpk('s'); w0=10; Gd=w0/(s+.001);
[K,CL,GAM]=loopsyn(G,Gd); % Design a loop-shaping controller K

% Plot the results
sigma(G*K,'r',Gd,'k-.',Gd/GAM,'k:',Gd*GAM,'k:',{.1,30})
figure ;T=feedback(G*K,eye(2));
sigma(T,ss(GAM),'k:',{.1,30});grid

hankelsv(K,'ncf','log');v=axis;
hold on; plot(v(1:2), ncfmargin(G,K)*[1 1],'--'); hold off
figure
K1=reduce(K,12,'errortype','ncf');
sigma(G*K1,'b',G*K,'r--',{.1,30});
二、多变量回路成形设计
点击见鲁棒控制工具箱多变量回路成形设计函数。
当系统不确定性不是很大时，可以设计高增益、高性能的反馈控制器。高回路增益尤其是幅值大于1时，能够减小系统不确定性的影响，降低整个系统对噪声的敏感性。但是，如果模型的不确定性很大，甚至不能得知增益的正负符号，那就不能采用大反馈增益设计方法，以免造成系统不稳定。因此，模型的不确定性是是反馈设计问题一个基本的限制因素。
1）乘性不确定性
对于模型G的逼近模型G0，

模型的G0的乘性不确定性 定义为：

或者，等价地有

模型的不确定性来源于多个方面。可能存在未建模的时延、寄生电容、执行器时间常数不能准确获得、机械系统中的高频振动形式等，都是模型不确定性的来源。这些类型的不确定性低频下很小，高频则会显著增大。
在单输入单输出系统中， 是个关键阈值，高于这个值就不能获得模型足够的信息来设计控制器。当模型不确定性变为200%时，模型不能提供实际系统的相位信息，这就意味着此时最稳妥的稳定系统的方法是使反馈增益小于1。
2）范数和奇异值
MIMO系统的传递函数为矩阵，增益通常用奇异值、H2范数、H∞范数来表示。
H2和H∞范数的物理含义为：H2范数表示模型G的脉冲响应的能量，H∞范数表示模型G在所有频域、所有输入下的峰值增益。
奇异值的定义：对于r阶矩阵 ，定义 为 所有特征值的非负均方根值，并排序： 。如果r
最大的奇异值 有时候定义为 。
如果A是n×n矩阵，第n个奇异值（即最小的奇异值）定义为 。
3）奇异值的性质
奇异值一些有用的性质如下：

这两个性质非常重要，因为它们决定了矩阵A的最大和最小奇异值，从而决定了当输入沿所有方向时的最大和最小增益。
对于稳定的连续时间LTI系统G(s)，H2范数和H∞范数用频域奇异值来定义。
H2范数：

H∞范数：

（此小节内容matlab帮助文档不全）
三、模型降阶
1、Hankel奇异值
针对设计复杂系统的鲁棒控制器，模型降阶具有下面优点：
1、根据待设计控制系统的指标，简化模型；
2、在保留了系统绝大部分动态特性的情况下，采用降阶模型能够加速仿真过程；
3、如果利用LQG、H∞设计方法，得到的控制器阶数会非常高，而且可能不是很必要。对控制律的阶数进行降阶有时能够在不影响控制性能的情况下降低控制律的复杂性。
鲁棒控制工具箱模型降阶函数提供了两种模型降阶方法：
（1）加性误差方法。在某种误差准则下，降阶模型有加性误差限制。
（2）乘性误差方法。在某种误差准则下，降阶模型有乘性误差限制。
误差用峰值增益穿越频率（H∞范数）来度量，误差的界限是被忽略的Hankel奇异值的函数。
在控制理论中，特征值决定了系统的稳定性，Hankel奇异值则决定了每个状态变量的能量。保留大能量的状态变量则意味着保留了系统大部分的动态特性，例如稳定性、频域和时域响应特性等。该工具箱中用到的模型降阶方法都是基于Hankel奇异值。
对于一个稳定的状态空间系统(A,B,C,D)，它的Hankel奇异值定义如下：

其中，P和Q分别表示能控和能观矩阵，并且满足

例如，输入以下指令：
rand('state',1234); randn('state',5678);
G = rss(30,4,3);
hankelsv(G)
得到图11所示的Hankel奇异值。

图11 Hankel奇异值
图11说明了大部分的能量储存在状态1~15中。后面将会看到如何利用模型降阶工具保留15个状态。
2、加性误差方法和乘性误差方法
鲁棒控制工具箱将系统不确定性分为加性不确定性和乘性不确定性。于是模型降阶误差分为两种：加性误差和乘性误差。
设得到的原系统G的降阶模型为Gred，加性误差上限为：

其中， 表示原系统G的第i个Hankel奇异值。
乘性误差上限为：

其中， 表示原系统G相位矩阵的第i个Hankel奇异值。
1、加性误差方法
给定一LTI系统，下列的指令能够将系统降阶为任意期望的阶数。
rand('state',1234); randn('state',5678);
G = rss(30,4,3);
% balanced truncation to models with sizes 12:16
[g1,info1] = balancmr(G,12:16); % or use reduce
% Schur balanced truncation by specifying `MaxError'
[g2,info2] = schurmr(G,'MaxError',[1,0.8,0.5,0.2]);
sigma(G,'b-',g1,'r--',g2,'g-.')
图12为原系统G与降阶系统g1、g2的比较图。

图12 原系统G与降阶系统g1、g2的比较图
为了检验是否满足理论误差，输入：
norm(G-g1(:,:,1),'inf') % 2.0123
info1.ErrorBound(1) % 2.8529
或者用下面的指令画出降阶误差与理论上限值：

图13 降阶误差与理论上限值
3、乘性误差方法
大多数情况下，降阶函数bstmr能够将原模型与降阶模型在期望的穿越频率处的相对误差限制在一定的范围内，得到比加性方法更准确的降阶模型。这种特性当系统具有小阻尼极点时更为明显。
下列的指令对比了乘性误差模型与加性误差模型：
rand('state',1234); randn('state',5678); G = rss(30,1,1);
[gr,infor] = reduce(G,'algo','balance','order',7);
[gs,infos] = reduce(G,'algo','bst','order',7);
figure(1);bode(G,'b-',gr,'r--');
title('Additive Error Method')
figure(2);bode(G,'b-',gs,'r--');
title('Relative Error Method')

图14 原模型与加性误差模型

图15 原模型与乘性误差模型
从图14与图15中可以看到，乘性误差模型具有更好的相位逼近特性。
3、modreal函数
对于某些具有小阻尼零/极点的系统，bstmr函数得到更好的降阶模型，在期望的频率范围内乘性误差很小。加性误差方法，如balancmr，schurmr，或者hankelmr只关心总体的绝对峰值误差，忽略了小阻尼零/极点频率范围内的模型降阶误差。
在一些情况下，降阶模型需要保留jω轴上的极点，例如柔性结构的刚性动力学模型、积分型控制器。modreal函数提供了此类实现方法。
modreal函数将系统转化为它的模态形式，对角化系数矩阵A元素均为特征值，实数特征值为1×1形式，复数特征值为2×2形式，默认情况下所有的块按

照特征值幅值的大小升序排列，或者按照特征值的实部降序排列。指定jω轴上的极点将模型分为两个子系统：一个子系统仅包含jω轴上的动态特性，另一个则包含非jω轴上的动态特性。
rand('state',5678); randn('state',1234); G = rss(30,1,1);
[Gjw,G2] = modreal(G,1); % only one rigid body dynamics
G2.d = Gjw.d; Gjw.d = 0; % put DC gain of G into G2
subplot(211);sigma(Gjw);ylabel('Rigid Body')
subplot(212);sigma(G2);ylabel('Nonrigid Body')

图16 modreal函数将原系统转换为两个子系统
对G2可以进一步进行降阶。当G2降阶为Gred时，原系统最终的降阶模型为Gjw+Gred。
区别jω轴上的极点的方法已经集成在所有的模型降阶函数中（balancmr, schurmr, hankelmr, bstmr, hankelsv），使用者不必自己划分模型。
下面的指令将原来30个状态变量的模型降阶为8阶模型：
rand('state',5678); randn('state',1234); G = rss(30,1,1);
[gr,info] = reduce(G); % choose a size of 8 at prompt
bode(G,'b-',gr,'r--')
上面的指令并没有指定降阶模型的阶数，自动画出Hankel奇异值，见下图：

图17 Hankel奇异值

图18 伯德图
reduce函数默认的降阶算法为balancmr，它完成了将30个状态变量降阶为8个状态变量中绝大部分的工作，同时保留了刚性动态特性以便将来控制器的设计。
当给大型系统（状态变量个数>200）时，modreal是唯一的降阶方法。由于计算的复杂性，大部分基于Hankel奇异值的方法都不能得到很好的降阶模型。modreal函数将大型动态特性转换为模态形式，然后变为包含有50个左右状态变量的中间模型，然后可以利用基于Hankel奇异值的方法进一步对中间模型进行降阶，得到阶数更低的、适合控制器设计的降阶模型。
在航天工业中，典型的柔性航天器包含有240个状态变量，依次利用modreal和bstmr（或者其它加性误差函数）函数，最终得到7个变量的包含有刚性动态特性的三轴模型。这时候可以轻而易举地利用鲁棒控制设计方法来进行控制器的设计。
4、ncfmr函数
鲁棒控制工具箱提供了一种特殊的模型降阶函数：ncfmr，它利用给定模型的互质集来降阶，利用归一化的互质因子将包含有积分器的控制器转换成降阶模型。它不需要modreal函数进行前/后处理，但是不再保留积分器。
rand('state',5678); randn('state',1234);
K= rss(30,4,3); % The same model G used in the 1st example
[Kred,info2] = ncfmr(K);
sigma(K,Kred)
上面的指令并没有指定降阶模型的阶数，自动画出Hankel奇异值，见下图：

图19 Hankel奇异值
同时，命令行提示输入期望阶数，此例中输入15，得到下图：

图20
如果积分控制很重要，前面提到的除了ncfmr的方法都可以很好地保留原模型中的积分器。
5、参考文献
[15] Glover, K., "All

Optimal Hankel Norm Approximation of Linear Multivariable Systems, and Their L∞ - Error Bounds," Int. J. Control, Vol. 39, No. 6, 1984, pp. 1145-1193.
[16] Zhou, K., Doyle, J.C., and Glover, K., Robust and Optimal Control, Englewood Cliffs, NJ, Prentice Hall, 1996.
[17] Safonov, M.G., and Chiang, R.Y., "A Schur Method for Balanced Model Reduction," IEEE Trans. on Automat. Contr., Vol. 34, No. 7, July 1989, pp. 729-733.
[18] Safonov, M.G., Chiang, R.Y., and Limebeer, D.J.N., "Optimal Hankel Model Reduction for Nonminimal Systems," IEEE Trans. on Automat. Contr., Vol. 35, No. 4, April 1990, pp. 496-502.
[19] Safonov, M.G., and Chiang, R.Y., "Model Reduction for Robust Control: A Schur Relative Error Method," International J. of Adaptive Control and Signal Processing, Vol. 2, 1988, pp. 259-272.
[20] Obinata, G., and Anderson, B.D.O., Model Reduction for Control System Design, London, Springer-Verlag, 2001.
四、鲁棒性分析
1、不确定性建模
对于控制工程师来说，理解不确定性的影响是一件重要的事情。减少某些不确定性（初始条件、低频干扰等）而不显著引入其它形式干扰的影响（传感器噪声、建模不确定性）是反馈控制系统设计的主要工作。
闭环稳定性分析是目前处理处理初始条件的不确定性及随机小干扰的主要方式。
低频范围内的高增益反馈是处理未知漂移及输出干扰的主要方式。在这种情况下，在高频范围内需要用滤波器处理高频范围内的传感器噪声。
最后，当面临建模不确定性时，增益及相位裕度是保证灵敏度及稳定性的有效工具，其并不需要精确知道控制输入如何影响反馈变量。
鲁棒控制工具箱提供了内置函数能够轻松建立不确定性模型。主要的建模模块，称为不确定元素（uncertain elements）或者atoms，为不确定实参数和不确定LTI对象，能够在现有的过程模型中建立简单或者复杂的不确定性模型。
当建成不确定模型后，能够利用高层次的鲁棒分析工具来评估不确定性对闭环系统稳定性等方面造成的性能影响。
模型不确定性主要有两种：
（1）微分方程中参数的不确定性；
（2）频域的不确定性，有时用频域响应的绝对或相对不确定性来量化。
一个不确定参数有它的名字（与其它不确定参数相区别）和标称值，它的不确定性用以下途径来表述：
（1）与标称值存在加性偏离；
（2）在标称值某一范围内波动；
（3）在标称值某一比例内波动；
创建一个实参数，名称为bw，标称值为5，比例不确定性为10%：
bw = ureal('bw',5,'Percentage',10)
这样就创建了一个ureal对象。用get命令观察它的属性：
get(bw)
结果如下：
Name: 'bw'
NominalValue: 5
Mode: 'Percentage'
Range: [4.5000 5.5000]
PlusMinus: [-0.5000 0.5000]

Percentage: [-10 10]
AutoSimplify: 'basic'
注意到波动范围（Range属性）和加性偏离（PlusMinus属性）与比例（Percentage）属性相一致。
可以在状态空间模型或者传递函数模型中引入ureal对象属性的不确定实参数，得到不确定状态空间模型对象（uss）。例如，用不确定实参数bw对一阶系统建模，其带宽在4.5~5.5rad/s之间：
H = tf(1,[1/bw 1])
结果如下：
USS: 1 State, 1 Output, 1 Input, Continuous System
bw: real, nominal = 5, variability = [-10 10]%, 1 occurrence
所得到的结果H为一不确定系统，称为uss对象。H的标称值为一状态空间模型，其极点为-5：
pole(H.NominalValue)
结果：ans =-5
查看H的频域和阶跃响应：
bode(H,{1e-1 1e2});
step(H)


图21 不确定系统H的伯德图

图22 不确定系统H的阶跃响应
尽管H的带宽和时间常数存在波动，但在高频处均以20db/decade下降，不受bw的影响。利用ultidyn函数能够观察到更多的在高频处复杂的不确定性行为。
带宽是描述数学模型与实际对象之间差别的一种非正式的方法。经常听到“这个模型在8rad/s处性能很好”。例如低于5rad/s，高于30rad/s，所建立的模型并不能代表实际对象，在这期间则与实际对象很好地吻合。不确定LTI对象ultidyn能够用来描述这种情况。ultidyn对象用来描述未知的线性系统，只知道其在频域范围内幅值一致有界。当存在标称模型和频域滤波器时，ultidyn对象能够捕获到模型的不确定性。
假设当大于9rad/s时，系统H的性能偏离其标称一阶系统，例如低频处的相对误差为5%，在高频衰减处则增加至1000%。这种情况用ultidyn对象按照下述步骤来表述：
（1）用tf，ss，zpk等函数建立标称模型Gnom，Gnom也许已经存在有不确定参数。在此例中Gnom是H，其为一阶系统，并具有一个不确定时间常数。
（2）创建一个权重滤波器W，它的幅值表示每一频率处的相对不确定性。makeweight函数能够创建一阶权重，并能够指定低频、高频增益以及增益穿越频率。
（3）创建名为Delta的ultidyn的对象，其幅值限制为1。
不确定对象用下式来表示：G = Gnom*(1+W*Delta)。
如果W表示绝对（而不是相对）不确定性，则利用下式代替上式：G = Gnom + W*Delta。
用matlab语言来描述上述步骤：
Gnom = H;
W = makeweight(.05,9,10);
Delta = ultidyn('Delta',[1 1]);
G = Gnom*(1+W*Delta)
得到下面结果：
SS: 2 States, 1 Output, 1 Input, Continuous System
Delta: 1x1 LTI, max. gain = 1, 1 occurrence
bw: real, nominal = 5, variability = [-10 10]%, 1 occurrence
G也为不确定系统，其依赖Delta和W。利用bode函数观察G在频域[0.1 100]rad/s处的伯德图：
bode(G,{1e-1 1e2})

图23 系统G的伯德图
下一节中将会为G设计两个反馈控制器，并比较它们

的性能。
2、鲁棒性分析
下面为G设计反馈控制器。设计的目标：良好的稳态跟踪性能和抑制干扰的能力。G为一阶滞后模型，选择PI控制器。设期望的闭环系统阻尼系数为 ，自然频率为 ，KI和Kp（基于标称系统开环时间常数为0.2）的表达式为

按照以下步骤设计控制器：
（1）为了研究G的不确定性如何影响闭环系统的带宽，设计两组控制器，阻尼系数 均为0.707，自然频率 分别为3和7.5。
xi = 0.707;
wn = 3;
K1 = tf([(2*xi*wn/5-1) wn*wn/5],[1 0]);
wn = 7.5;
K2 = tf([(2*xi*wn/5-1) wn*wn/5],[1 0]);
注意到，K2在G不确定性最显著的区域具有与标称系统同样的带宽。后面将会看到，模型的变化将会带来闭环系统性能显著下降，这在意料之中。
（2）利用feedback函数形成闭环系统：
T1 = feedback(G*K1,1);
T2 = feedback(G*K2,1);
（3）观察闭环系统的阶跃响应：
tfinal = 3;
step(T1,'b',T2,'r',tfinal,20)

图24 闭环系统T1和T2的阶跃响应
从图中可以看到，T2的阶跃响应更快，这是由于K2设置了更高的闭环系统带宽，但模型不确定性造成的影响也更大。
利用robuststab函数检验对系统参数摄动的鲁棒性：
[stabmarg1,destabu1,report1] = robuststab(T1);
stabmarg1
结果如下：
ubound: 4.0241
lbound: 4.0241
destabfreq: 3.4959
输入：
[stabmarg2,destabu2,report2] = robuststab(T2);
stabmarg2
结果如下：
ubound: 1.2545
lbound: 1.2544
destabfreq: 10.5249
变量stabmarg给出了稳定频率区域的上下界。稳定区域大于1意味着对于所有的摄动值都是稳定的。如果小于1就意味着存在使系统不稳定的摄动值。变量report简要描述了鲁棒性分析结果：
report1
结果如下：
report1 =
Uncertain System is robustly stable to modeled uncertainty.
-- It can tolerate up to 402% of modeled uncertainty.
-- A destabilizing combination of 402% the modeled uncertainty
exists, causing an instability at 3.5 rad/s.
输入：
report2
结果如下：
report2 =
Uncertain System is robustly stable to modeled uncertainty.
-- It can tolerate up to 125% of modeled uncertainty.
-- A destabilizing combination of 125% the modeled uncertainty
exists, causing an instability at 10.5 rad/s.
此例中的两个系统均对不确定性稳定，但不确定性对其影响的程度不一样。为了描述不确定性对闭环系统性能的影响，利用wcgain函数计算最坏情况下闭环系统灵敏度函数S=1/(1+GK)的峰值增益。峰值增益与单位阶跃响应超调量有关。
计算灵敏度函数，并调用wcgain函数：
S1 = feedback(1,G*K1);
S2 = feedback(1,G*K2);
[maxgain1,wcu1] = wcgain(S1);
maxgain1
结果如下：
maxgain1 =
lbound: 1.8684
ubound: 1.9025
critfreq: 3.5152
输入

：
[maxgain2,wcu2] = wcgain(S2);
maxgain2
结果如下：
maxgain2 =
lbound: 4.6031
ubound: 4.6671
critfreq: 11.0231
变量maxgain给出了最坏情况下的灵敏度函数峰值增益的上限和下限值，以及最大增益处的频率值。变量wcu含有达到最坏情况的不确定变量的值。
利用usubs函数把不确定变量用最坏情况下的值来代替，并比较标称系统与最坏情况下的伯德图和阶跃响应：
bodemag(S1.NominalValue,'b',usubs(S1,wcu1),'b');
hold on
bodemag(S2.NominalValue,'r',usubs(S2,wcu2),'r');
hold off

图25 标称系统与最坏情况下系统的伯德图
tfinal = 3;
step(S1.NominalValue,'b',usubs(S1,wcu1),'b',tfinal)
hold on
step(S2.NominalValue,'r',usubs(S2,wcu2),'r',tfinal)
hold off

图26 标称系统与最坏情况下系统的阶跃响应图
从图中可以看到，虽然K2灵敏度优于K1，闭环系统的带宽过宽，已经到了不确定性主要集中的频域。此例中K2的最坏情况下的性能要逊于K1。
3、MIMO系统的鲁棒性分析
上一节中以一阶传递函数为例，研究了SISO不确定系统的鲁棒控制问题。类似地，也可以建立含有不确定系数矩阵的状态空间模型，利用前面介绍的工具对其进行分析。
考虑两输入、两输出、两状态的状态空间模型，其系数矩阵中含有不确定参数。首先创建一个不确定参数p。利用这个参数，构造不确定性A、C矩阵。B矩阵不包含不确定性：
p = ureal('p',10,'Percentage',10);
A = [0 p;-p 0];
B = eye(2);
C = [1 p;-p 1];
H = ss(A,B,C,[0 0;0 0]);
利用get函数观察系统H的属性：
get(H)
结果如下：
a: [2x2 umat]
b: [2x2 double]
c: [2x2 umat]
d: [2x2 double]
StateName: {2x1 cell}
Ts: 0
InputName: {2x1 cell}
OutputName: {2x1 cell}
InputGroup: [1x1 struct]
OutputGroup: [1x1 struct]
NominalValue: [2x2 ss]
Uncertainty: [1x1 atomlist]
Notes: {}
UserData: []
属性a，b，c，d和StateName与ss对象相同，属性InputName，OutputName，InputGroup和OutputGroup与系统对象（ss，zpk，tf，frd）相同。NomialValue为ss对象。
系统H不包含有执行器的动态特性。换句话说，执行器在所有频率处增益均为1。
然而，通道1的执行器在低频处具有10%的不确定性，在大于20rad/s的高频处建模不准确。通道2的执行器在低频处具有20%的不确定性，一直到45rad/s建模都很准确。利用ultidyn对象Delta1和Delta2，以及成形滤波器W1和W2，将此处描述的频域不确定性加入到模型H中：
W1 = makeweight(.1,20,50);
W2 = makeweight(.2,45,50);
Delta1 = ultidyn('Delta1',[1 1]);
Delta2 = ultidyn('Delta2',[1 1]);
G = H*blkdiag(1+W1*Delta1,1+W2*Delta2)
结果如下：
USS: 4 States, 2 Outputs, 2 Inputs, Continuous System
Delta1: 1x1 LTI

, max. gain = 1, 1 occurrence
Delta2: 1x1 LTI, max. gain = 1, 1 occurrence
p: real, nominal = 10, variability = [-10 10]%, 2 occurrences
注意，G是两输入、两输出的不确定系统，具有三个不确定元素：Delta1，Delta2和p。它具有四个状态变量，两个继承于系统H，两外两个来自成形滤波器W1和W2。
观察系统G在2s内的阶跃响应：
step(G,2)

图27 MIMO系统G的阶跃响应
观察系统G的伯德图，为了简洁起见，从振荡结束后开始作图：
bode(G,{13 100})

图28 MIMO系统G的伯德图
载入控制器，检验其是否为两输入、两输出：
load mimoKexample
size(K)
结果如下：
State-space model with 2 outputs, 2 inputs, and 9 states.
利用loopsens函数能够形成标准的过程模型/控制器反馈连接形式，包括输入、输出端的灵敏度和补灵敏度分析：
F = loopsens(G,K)
结果如下：
F =
Poles: [13x1 double]
Stable: 1
Si: [2x2 uss]
Ti: [2x2 uss]
Li: [2x2 uss]
So: [2x2 uss]
To: [2x2 uss]
Lo: [2x2 uss]
PSi: [2x2 uss]
CSo: [2x2 uss]
F.Ploes中包含有标称系统的极点。如果闭环标称系统是稳定的，F.Stable值为1。S表示灵敏度，T表示补灵敏度，L表示开环增益，下标i、o表示模型G的输入和输出。P表示模型，C表示控制器。
Ti的数学形式为

Lo为G*K，CSo数学形式为

F.PSi包含了输出对输入干扰的传递函数，用伯德图来查看：
bodemag(F.PSi,':',{1e-1 100})

图29 输出对输入干扰的传递函数的伯德图
可以用loopmargin函数查看某一时刻回路的增益和相位裕度等指标。它们针对标称系统，没有考虑系统G包含的不确定性。
考察每一输入的的稳定裕度、输出的稳定裕度及综合输入/输出稳定裕度：（Explore the simultaneous margins individually at the plant input, individually at the plant output, and simultaneously at both input and output）：
[I,DI,SimI,O,DO,SimO,Sim] = loopmargin(G,K);（此命令没有运行成功）
第三个输出变量为模型输入允许的即时增益和相位变化：
SimI
结果显示：
SimI =
GainMargin: [0.1180 8.4769]
PhaseMargin: [-76.5441 76.5441]
Frequency: 6.2287
结果说明，每个输入端增益可以大致在1/8~8之间波动，相位波动范围达到76°。
第六个输出变量为模型输出端的即时增益和相位裕度：
SimO
结果显示：
SimO =
GainMargin: [0.1193 8.3836]
PhaseMargin: [-76.3957 76.3957]
Frequency: 18.3522
注意到输出端即时波动范围与输入端数值近似，在多回路反馈控制中并不总是这种情况。
最后一个输出变量是综合输入/输出即时增益、相位稳定裕度。这种情况下同时考虑输入、输出变化，这时得到的稳定范围要比输入、输出端单独变

化要小：
Sim
结果显示：
Sim =
GainMargin: [0.5671 1.7635]
PhaseMargin: [-30.8882 30.8882]
Frequency: 18.3522
这些数据说明闭环系统具有鲁棒性，能够容忍输入、输出端同时进行大增益变化及30°的相位变化。
loopmargin函数能够分析标称、多回路系统的稳定裕度，不能用来分析含有ureal、ultidyn等对象的不确定系统。当分析复杂的不确定模型时，loopmargin函数计算得到的数据可能不是实际稳定裕度。此时可以利用robuststab函数计算不确定系统的稳定裕度。
例如，用robuststab函数计算含有不确定元素Delta1，Delta2，p的闭环控制系统的稳定裕度。利用前面得到的闭环控制系统：F = loopsens(G,K)。F中所有的元素，例如F.Si, F.To等，都具有相同的动态，因此具有相同的稳定特性：
[stabmarg,desgtabu,report] = robuststab(F.So);
stabmarg
结果如下：
stabmarg =
ubound: 2.2175
lbound: 2.2175
destabfreq: 13.7576
输入：report
结果如下：
report =
Uncertain System is robustly stable to modeled uncertainty.
-- It can tolerate up to 222% of the modeled uncertainty.
-- A destabilizing combination of 222% of the modeled uncertainty exists,
causing an instability at 13.8 rad/s.
-- Sensitivity with respect to uncertain element ...
'Delta1' is 51%. Increasing 'Delta1' by 25% leads to a 13% decrease in the margin.
'Delta2' is 55%. Increasing 'Delta2' by 25% leads to a 14% decrease in the margin.
'p' is 93%. Increasing 'p' by 25% leads to a 23% decrease in the margin.
分析结果也证实了loopmargin函数的分析。闭环系统对与Delta1，Delta2，p的变化具有相当的鲁棒性。该系统实际上能够容忍2倍于现在的不确定性而且能保持闭环稳定性。
4、最坏情况下增益分析
画出标准输出灵敏度函数的伯德图，结果显示在低频处所有通道具有良好的抑制干扰的能力：
bodemag(F.So.NominalValue,{1e-1 100})

图30 标准输出灵敏度函数的伯德图
利用norm函数计算频率响应矩阵最大奇异值的峰值：
[PeakNom,freq] = norm(F.So.NominalValue,'inf')
结果如下：
PeakNom =
1.1317
freq =
7.0483
可见，峰值为1.13，发生在频率36rad/s处。
当Delta1，Delta2和p变化时，最大输出灵敏度增益是多少？可以用wcgain函数来回答这个问题。
[maxgain,wcu] = wcgain(F.So);
maxgain
结果如下：
maxgain =
lbound: 2.1017
ubound: 2.1835
critfreq: 8.5546
分析结果表明，最坏情况下的增益介于2.1与2.2之间，峰值频率大概为8.5。
利用usubs函数将Delta1，Delta2，p替换成达到峰值增益2.1时的值。在输出补灵敏度函数中做此替换，并施加单位阶跃响应：
ste

p(F.To.NominalValue,usubs(F.To,wcu),5)

图31 补灵敏度函数的单位阶跃响应
结果显示，即使面临不确定参数最坏情况下的组合，输出灵敏度函数的响应并没有明显恶化。稳定时间由2s增加到4s，非对角耦合增加了2倍，但仍然非常小。