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

1、西安交通大学自动化系利用Matlab实现廿控制Prof. Dr.-Ing.F.AllgwerInstitute for Systems Theory and Automatic Control http:/www.ist.uni-stuttgart.de/education/courses/robust 1引言Ha控制器设计原理容易理解，难点在于编程。这里简单介绍 Matlab里面几 个相关函数的用法，希望能帮助你设计第一个HF空制器。Matlab提供了很多Ho计函数，与Ho计相关的几个重要的工具箱有： Control System Toolbox, mu-Analysis and Synth

2、esis Toolbox(mu-tools), Robust Control Toolbox (RCT)和 LMI Control Toolbox。Matlab7.0 之后的版本中,LMI 和mu-tools都包含在RCT v3.0.1中，Matlab 7.0之前的版本中这些工具箱是独立 I本文中用到的函数都写在了一个 m文件中(见附录)，也可以从网站下载。利用混合S/KS问题说明HW目关函数的用法。首先回顾这个问题。2回路成形传递函数混合S/KS问题可用图1来说明。Figure 1: Mixed S/KS problem.从w到z的闭环传递函数丁 = E(尸,茨)可以表示为WsSVksKSF

3、igure 2: General control problem.广义过程模型P(s)(见图2)可以表示为原文：Hoo Control in Matlab. 2011/12/8 译假设上面几个状态空间变量具有如下的形式:Bks于是可以得到P(s)的一个可能的状态空间实现形式:As 00000cks0-DSC0ADC0-cBs -BSD 100瓦0IksBDk3-DWS和Wks为调整参数。一种选择方法为：s/M +v S ; H const.s + 3。A其中A1为允许的最大稳态误差，0为期望带宽，M为灵敏度峰值(一般情况 下A=0.01, M=2)o从控制器设计的方面来说， Hs的倒数为回路成

4、形期望灵敏度的上限，”屋影响控制器的输出Uo有些情况下希望对补灵敏度函数丁 = GI(I + GA3T进行回路成形设计(在图1中增加一个输出3 = W)。一种选择方法为:TT7 + 30/AIWt =As + 3。此函数与函数成2 = 30.轴对称，见图3所示。图中参数设置为A=0.01(=-40dB), M=2 (=6dB),1 ri/sec,oia-2D2101010Frequency rad.7secFigure 3: Inverse of weighting functions 11 s and 0。3子系统的实现在Matlab中有几种方式得到 G, WS和 Wks。例如Control

5、 System Toolbox提 供的ss, tf和zpk等函数。Mu-tools也提供了诸如pck, nd2sys, zp2sys等函数， 也可以用mksys和tree等方法。需要注意的是，Mu-tools提供了一种与Control System Toolbox不一样的表达方式：系统矩阵(system matrix)。Control System Toolbox里面可以写成Gcst =SS (A jB , C jD)对mu-tools则不适用。Mlltnnk A,B,C,D=ssdata(Gcst); Gmu=pck(A,B,C,D),mutooisA,BfC,D=unpck(Gmu); G

6、cst = ss(AfB,C,D).4广义系统P的实现广义系统P也有多种产生方式。下面我列举了五种：(1)直接写出传递函数矩阵，见2节。我倾向于利用mu-tools。如果后面需 要转化为状态空间模型，需要利用 minreal之类的函数得到最小化实现。重要的 函数有：sbs (side-by-side), abv (above), mmult (multiply ), minv (inverse) 0(2)写出状态空间系数矩阵A,B,C,D,然后利用 3节中的指令：P=pck(A,B,C,D)。(3)利用 sysic函数(system interconnect。先在一个 m文件中利用 mu-to

7、ols 设置了子系统，然后利用该函数将子系统互连起来。(4)利用sconnect函数，这是LMI-tools提供的函数。子系统、输入、输出 以参数的形式传递，sconnect返回互连的系统。(5)利用RCT v3.0.1提供的iconnect函数，该函数功能与sysic类似。这些方法中，我倾向于sysic和iconnect函数，因为它们使用灵活，能够搭 建方法1和方法2难于搭建的复杂系统。一般情况下，得到均衡实现形式以避免 数值计算问题是一个很好的想法，均衡实现形式可以利用 Control System Toolbox提供的balreal实现5控制器设计兄X S/KS设计问题是寻找一控制器 K

8、,使之稳定系统G,并使下列的H8范数最小:|瓦（产,女川8 = /Us有很多种得到H0t空制器的方法，例如 hinfsyn, hinfric , hinflmi ,这些函数将P作为输入，并以系统矩阵（mu-tools）方式表达。RCT v3.0.1提供了 mixsyn函数，将G，RS，口 AW作为输入（为补灵敏度函数权重），并不需要事 先得到广义系统模型P。这些方法的主要区别在于是否用到 Riccati方程和丫迭 代或者线性矩阵不等式来求解最优化问题。 LMI方法不需要求解Riccati时设定 的假设条件。另外还有些指令：ncfsyn和loopsyn （对开环彳专递函数L=GK进行H3回路成

9、形设计），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。如果你

10、知道Control System Toolbox存在某个函数， 那么mu-tools也很可能存在具有相同功能的函数，只不过名称稍有些不同。如果 你清楚自己的目的，但不知道函数的名字，那么建议你浏览帮助手册中的函数索希望此处简短的介绍能够帮助你进行 H问空制器的设计。好运！ I%下面的代码展示了如何在Matlab中进行H-infinity控制器的设计。此处举的例 子与混合%S/KS问题有些不同。此处用到的模型和权重函数见 Skogestadffi Postlethwaite, 1996,%ed.1,p.60权重函数并不是“最优”的。 %大部分函数来自 mu-tools, 一些来自lmi-tool

11、s。mu-tools和lmi-t001s均包含在 RCT%v3.0.1 中。%-Jorgen Johnsen 14.12.06% %建立子系统%Plant:G=200/(10s+1)(0.05s+1)A2)%方法1：直接方法，利用 mu-toolsG=nd2sys(1,conv(10,1,conv(0.05,1,0.05 1),200);%方法 2： control system toolboxs=tf(s);Gcst=200/(10*s+1)*(0.05*s+1)A2);a,b,c,d=ssdata(balreal(Gcst);G=pck(a,b,c,d);%权重：Ws=(s/M+w0)/(

12、s+w0*A),Wks=1M=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

13、);P0=pck(ab,bb,cb,db); %t匕时得到变量为 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);%全

14、系统的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(A

15、p,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:禾1J用sconnect函数inputs=r(1);u(1);outputs=Ws;Wks;e=r-G

16、;K_in=;%无控制器G_in=G:u;Ws_in=Ws:e;Wks_in=Wks:u;P3,r=sconnect(inputs,outputs,K_in,G_in,GWs_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)

17、;M=iconnect;M.Input=r;u;M.Output=ws;wks;e;M.Equation1=equate(e,r-y);M.Equation2=equate(y,ss(A,B,C,D)*u);M.Equation3=equate(ws,ss(A1,B1,C1,D1)*e);M.Equation4=equate(wks,ss(A2,B2,C2,D2)*u);ab,bb,cb,db=ssdata(balreal(M.System);P4=pck(ab,bb,cb,db);%控制器的设计%下面使用的方法均基于广义系统 P的System矩阵表达形式%选择你喜欢的控制器设计方法和广义系统

18、 P%选择广义系统PP=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);%利用RC

19、T 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函数。也

20、可以使用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

21、); % 灵敏度函数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代替，单位是dBofigure(2)vplot(liv,lm,Sw,-,Tw,-,GKw,-.);title(sigma(S(jw)(solid),sigma(T(jw)(dashed) a

22、nd 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(WsA-1(jw)(dashed),)xlabel(Frequency rad/sec);ylabel(Amplitude) figure(5)vplot(liv,lm,KSw,-,KSdw,-);title(sigma(KS(jw)(solid) and sigma(WksA-1(jw)(dashed),) xlabel(Frequency rad/sec);ylabel(Amplitude)%最后，阶跃