利用matlab实现H-infinity鲁棒控制

西安交通大学自动化系
利用Matlab实现H∞控制
Prof. Dr.-Ing.F.Allgwer
Institute for Systems Theory and Automatic Control
http://www.ist.uni-stuttgart.de/education/courses/robust
1 引言
H∞控制器设计原理容易理解，难点在于编程。

这里简单介绍Matlab里面几个相关函数的用法，希望能帮助你设计第一个H∞控制器。

Matlab提供了很多H∞设计函数，与H∞设计相关的几个重要的工具箱有：Control System Toolbox，mu-Analysis and Synthesis Toolbox(mu-tools)，Robust Control Toolbox（RCT）和LMI Control Toolbox。

Matlab7.0之后的版本中，LMI 和mu-tools都包含在RCT v3.0.1中，Matlab 7.0之前的版本中这些工具箱是独立的。

本文中用到的函数都写在了一个m文件中（见附录），也可以从网站下载。

利用混合S/KS问题说明H∞相关函数的用法。

首先回顾这个问题。

2 回路成形传递函数
混合S/KS问题可用图1来说明。

从w到z的闭环传递函数可以表示为
广义过程模型P(s)（见图2）可以表示为
假设上面几个状态空间变量具有如下的形式：
于是可以得到P(s)的一个可能的状态空间实现形式：
W S和W KS为调整参数。

一种选择方法为：
其中A<1为允许的最大稳态误差，为期望带宽，M为灵敏度峰值（一般情况下A=0.01，M=2）。

从控制器设计的方面来说，的倒数为回路成形期望灵敏度的上限，影响控制器的输出u。

有些情况下希望对补灵敏度函数进行回路成形设计（在图1中增加一个输出）。

一种选择方法为：
此函数与函数成轴对称，见图3所示。

图中参数设置为A=0.01（=-40dB），M=2（=6dB），。

3 子系统的实现
在Matlab中有几种方式得到G，W S和W KS。

例如Control System Toolbox提供的ss，tf和zpk等函数。

Mu-tools也提供了诸如pck，nd2sys，zp2sys等函数，也可以用mksys和tree等方法。

需要注意的是，Mu-tools提供了一种与Control System Toolbox不一样的表达方式：系统矩阵（system matrix）。

Control System
Toolbox里面可以写成，对mu-tools则不适用。

Mu-tools：
4 广义系统P的实现
广义系统P也有多种产生方式。

下面我列举了五种：
（1）直接写出传递函数矩阵，见2节。

我倾向于利用mu-tools。

如果后面需要转化为状态空间模型，需要利用minreal之类的函数得到最小化实现。

重要的函数有：sbs（side-by-side），abv（above），mmult（multiply），minv（inverse）。

（2）写出状态空间系数矩阵A,B,C,D，然后利用3节中的指令：P=pck （A,B,C,D）。

（3）利用sysic函数（system interconnect）。

先在一个m文件中利用mu-tools 设置了子系统，然后利用该函数将子系统互连起来。

（4）利用sconnect函数，这是LMI-tools提供的函数。

子系统、输入、输出以参数的形式传递，sconnect返回互连的系统。

（5）利用RCT v3.0.1提供的iconnect函数，该函数功能与sysic类似。

这些方法中，我倾向于sysic和iconnect函数，因为它们使用灵活，能够搭建方法1和方法2难于搭建的复杂系统。

一般情况下，得到均衡实现形式以避免数值计算问题是一个很好的想法，均衡实现形式可以利用Control System Toolbox
提供的balreal实现。

5 控制器设计
设计问题是寻找一控制器K，使之稳定系统G，并使下列的H∞范数最小：
有很多种得到H∞控制器的方法，例如hinfsyn，hinfric，hinflmi，这些函数将P作为输入，并以系统矩阵（mu-tools）方式表达。

RCT v3.0.1提供了mixsyn
函数，将作为输入（为补灵敏度函数权重），并不需要事
先得到广义系统模型P。

这些方法的主要区别在于是否用到Riccati方程和γ迭代或者线性矩阵不等式来求解最优化问题。

LMI方法不需要求解Riccati时设定的假设条件。

另外还有些指令：ncfsyn和loopsyn（对开环传递函数L=GK进行H∞回路成形设计），hinfmix和msfsyn（多目标）。

具体请查阅手册。

6 结果分析
当控制器设计好后，需要对结果进行分析，这时可以利用Control System Toolbox提供的函数，例如利用lsim，step（阶跃响应），bode（伯德图），sigma （奇异值），freqresp（频域响应）等函数对传递函数S，KS，T，K，GK进行分析。

Mu-tools提供的函数有：trsp（时域响应），frsp（频域响应），vsvd（奇异值），vplot。

7 结论
从前面可见，有很多种进行H∞控制器设计的方法。

为了避免混淆，我建议尽可能使用mu-tools和RCT。

如果你知道Control System Toolbox存在某个函数，那么mu-tools也很可能存在具有相同功能的函数，只不过名称稍有些不同。

如果你清楚自己的目的，但不知道函数的名字，那么建议你浏览帮助手册中的函数索引。

希望此处简短的介绍能够帮助你进行H∞控制器的设计。

好运！
附录
%下面的代码展示了如何在Matlab中进行H-infinity控制器的设计。

此处举的例子与混合
%S/KS问题有些不同。

此处用到的模型和权重函数见Skogestad和Postlethwaite，1996，
%ed.1,p.60.权重函数并不是“最优”的。

%
%大部分函数来自mu-tools，一些来自lmi-tools。

mu-tools和lmi-tools均包含在RCT
%v3.0.1中。

%-Jorgen Johnsen 14.12.06
%-------------------------------------------------------------------------
%建立子系统
%------------------------------------------------------------------------- %Plant:G=200/((10s+1)(0.05s+1)^2)
%方法1：直接方法，利用mu-tools
G=nd2sys(1,conv([10,1],conv([0.05,1],[0.05 1])),200);
%方法2：control system toolbox
s=tf('s');
Gcst=200/((10*s+1)*(0.05*s+1)^2);
[a,b,c,d]=ssdata(balreal(Gcst));
G=pck(a,b,c,d);
%权重：Ws=(s/M+w0)/(s+w0*A),Wks=1
M=1.5;w0=10;A=1.e-4;
Ws=nd2sys([1/M w0],[1 w0*A]);
Wks=1;
%------------------------------------------------------------------------- %建立广义系统P
%------------------------------------------------------------------------- %方法0：直接方法
%/z1\ /Ws -Ws*G\ /r\
%|z2| =|0 Wks | | |
%\ v/ \I -G / \u/
%传递函数表达方法
Z1=sbs(Ws,mmult(-1,Ws,G));
Z2=sbs(0,Wks);
V=sbs(1,mmult(-1,G));
P0=abv(Z1,Z2,V);
%通常情况下P0并不是最小实现，所以需要降阶
[a,b,c,d]=unpck(P0);
[ab,bb,cb,db]=ssdata(balreal(minreal(ss(a,b,c,d))));
P0=pck(ab,bb,cb,db); %此时得到变量为System类型
%------------------------------------------------------------------------- %建立广义系统P
%------------------------------------------------------------------------- %方法1：直接方法
%/z1\ /W1 -W1*G\ /r\
%|z2| =|0 W2 | | |
%\ v/ \I -G / \u/
%子系统的ss实现
[A,B,C,D]=unpck(G);
[A1,B1,C1,D1]=unpck(Ws);
[A2,B2,C2,D2]=unpck(Wks);
%计算不同子系统的输入、输出变量的个数
n1=size(A1,1);[q1,p1]=size(D1);
n2=size(A2,1);[q2,p2]=size(D2);
n=size(A,1);[q,p]=size(D);%原文此处为[p,q]=size(D);
%全系统的ss实现
Ap=[A1 zeros(n1,n2) -B1*C;
zeros(n2,n1) A2 zeros(n2,n);
zeros(n,n1) zeros(n,n2) A];
Bp=[B1 -B1*D;
zeros(n2,p) B2;
zeros(n,p) B];
Cp=[C1 zeros(q1,n2) -D1*C;
zeros(q2,n1) C2 zeros(q2,n);
zeros(q,n1) zeros(q,n2) -C];
Dp=[D1 -D1*D;
zeros(q2,p) D2;
eye(p) -D];
%得到均衡实现形式，以减少可能产生的计算问题
[Apb,Bpb,Cpb,Dpb]=ssdata(balreal(ss(Ap,Bp,Cp,Dp)));
P1=pck(Apb,Bpb,Cpb,Dpb); %P1与P0数值相近，但符号有差别%-------------------------------------------------------------------------
%建立广义系统P
%-------------------------------------------------------------------------
%方法2：利用sysic函数
systemnames='G Ws Wks';
inputvar='[r(1);u(1)]';%所有输入均为标量，r(2)为两维信号outputvar='[Ws;Wks;r-G]';
input_to_G='[u]';
input_to_Ws='[r-G]';
input_to_Wks='[u]';
sysoutname='P2';
cleanupsysic='yes';
sysic
%-------------------------------------------------------------------------
%建立广义系统P
%-------------------------------------------------------------------------
%方法3：利用sconnect函数
inputs='r(1);u(1)';
outputs='Ws;Wks;e=r-G';
K_in=[];%无控制器
G_in='G:u';
Ws_in='Ws:e';
Wks_in='Wks:u';
[P3,r]=sconnect(inputs,outputs,K_in,G_in,G,Ws_in,Ws,Wks_in,Wks);
%-------------------------------------------------------------------------
%建立广义系统P
%-------------------------------------------------------------------------
%方法4：利用iconnect函数
%注意1：不再使用mu-tools System表达形式
%注意2：iconnet函数仅适用于Robust Control Toolbox v3.0.1及以上版本r=icsignal(1);
u=icsignal(1);
ws=icsignal(1);
wks=icsignal(1);
e=icsignal(1);
y=icsignal(1);
M=iconnect;
M.Input=[r;u];
M.Output=[ws;wks;e];
M.Equation{1}=equate(e,r-y);
M.Equation{2}=equate(y,ss(A,B,C,D)*u);
M.Equation{3}=equate(ws,ss(A1,B1,C1,D1)*e);
M.Equation{4}=equate(wks,ss(A2,B2,C2,D2)*u);
[ab,bb,cb,db]=ssdata(balreal(M.System));
P4=pck(ab,bb,cb,db);
%-------------------------------------------------------------------------
%控制器的设计
%-------------------------------------------------------------------------
%下面使用的方法均基于广义系统P的System矩阵表达形式
%选择你喜欢的控制器设计方法和广义系统P
%选择广义系统P
P=P1; %P=P0;P=P1;P=P2;P=P3;P=P4;
%然后设置一些参数（测量变量个数，输入变量个数，gamma限制）nmeas=1;nu=1;gmn=0.5;gmx=20;tol=0.001;
%控制器设计：选择你喜欢的设计方法，不需要的注释掉
[K,CL,gopt]=hinfsyn(P,nmeas,nu,gmn,gmx,tol);
[gopt,K]=hinflmi(P,[nmeas,nu],0,tol); CL=starp(P,K,nmeas,nu);
[gopt,K]=hinfric(P,[nmeas,nu],gmn,gmx);CL=starp(P,K,nmeas,nu);
%利用RCT v3.0.1进行控制器的设计
%通常不需要系统的传递函数表达式，但更需要ss类型
[a,b,c,d]=unpck(G); Gcst=ss(a,b,c,d);
[a,b,c,d]=unpck(Ws); Wscst=ss(a,b,c,d);
[a,b,c,d]=unpck(Wks); Wkscst=ss(a,b,c,d);
[K,CL,gopt]=mixsyn(Gcst,Wscst,Wkscst,[]);
[a,b,c,d]=ssdata(balreal(K)); K=pck(a,b,c,d);
[a,b,c,d]=ssdata(balreal(CL)); CL=pck(a,b,c,d);
%-------------------------------------------------------------------------
%分析结果
%-------------------------------------------------------------------------
%注意：此处使用mu-tools函数。

也可以使用control toolbox指令，例如用series 和
%feedback进行子系统连接，sigma或freqresp，svd和bode画奇异值，step和lsim 分析
%时域响应
%画(Weighted)闭环系统的奇异值
w=logspace(-4,6,50);
CLw=vsvd(frsp(CL,w));
figure(1)
vplot('liv,m',CLw)
title('singular values of weighted closed loop system')
%得到传递函数矩阵
[type,out,in,n]=minfo(G);
I=eye(out);
S=minv(madd(I,mmult(G,K))); %灵敏度函数
T=msub(I,S); %补灵敏度函数
KS=mmult(K,S); %G的输入
GK=mmult(G,K); %回路传递函数
%频域的奇异值
Sw=vsvd(frsp(S,w));
Tw=vsvd(frsp(T,w));
Kw=vsvd(frsp(K,w));
KSw=vsvd(frsp(KS,w));
GKw=vsvd(frsp(GK,w));
%画奇异值
%注意：如果愿意的话，可以将vplot用plot代替，单位是dB。

figure(2)
vplot('liv,lm',Sw,'-',Tw,'--',GKw,'-.');
title('\sigma(S(jw))(solid),\sigma(T(jw))(dashed) and \sigma(GK(jw))(dashdot),') xlabel('Frequency [rad/sec]');ylabel('Amplitude')
figure(3)
vplot('liv,lm',Kw)
title('\sigma(K(jw))')
xlabel('Frequency [rad/sec]');ylabel('Amplitude')
%是否满足设计要求？
Sd=minv(Ws); Sdw=vsvd(frsp(Sd,w)); %期望的灵敏度
KSd=minv(Wks); KSdw=vsvd(frsp(KSd,w)); %期望的输出
figure(4)
vplot('liv,lm',Sw,'-',Sdw,'--');
title('\sigma(S(jw))(solid) and \sigma(Ws^{-1}(jw))(dashed),')
xlabel('Frequency [rad/sec]');ylabel('Amplitude')
figure(5)
vplot('liv,lm',KSw,'-',KSdw,'--');
title('\sigma(KS(jw))(solid) and \sigma(Wks^{-1}(jw))(dashed),')
xlabel('Frequency [rad/sec]');ylabel('Amplitude')
%最后，阶跃响应
reference=1;tfinal=1;step=0.01;
y=trsp(T,reference,tfinal,step);
u=trsp(KS,reference,tfinal,step);
figure(6)
subplot(2,1,1)
vplot('iv,d',y)
title('Step response');ylabel('y')
subplot(2,1,2)
vplot('iv,d',u)
ylabel('y');xlabel('time')。