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

1、利用Matlab实现H8控制Prof.Dr.-Ing.F.AllgwerInstituteforSystemsTheoryandAutomaticControlhttp:/www.ist.uni-stuttgart.de/education/courses/robust1引言控制器设计原理容易理解，难点在于编程。这里简单介绍Matlab里面几个相关函数的用法，希望能帮助你设计第一个控制器。Matlab提供了很多设计函数，与设计相关的几个重要的工具箱有：ControlSystemToolbox，mu-AnalysisandSynthesisToolbox(mu-tools)，RobustCont

2、rolToolbox(RCT)和LMIControlToolbox。Matlab7.0之后的版本中，LMI和mu-tools都包含在RCTv3.0.1中，Matlab7.0之前的版本中这些工具箱是独立的。本文中用到的函数都写在了一个m文件中(见附录)，也可以从网站下载。利用混合S/KS问题说明相关函数的用法。首先回顾这个问题。2回路成形传递函数混合S/KS问题可用图1来说明。wI<sKS从w到Z的闭环传递函数=可以表示为广义过程模型P(s)(见图2)可以表示为原文：ControlinMatlab.2011/12/8译假设上面几个状态空间变量具有如下的形式：于是可以得到P(s)的一个可能的

3、状态空间实现形式:比0-BSCBs-BSD0Aks00ks00A0B0-DSCDs-DSD0Cks00D曉00-CI-Dws和wks为调整参数。一种选择方法为:拆/'Z+TTr丄；=;=const,s+其中Avl为允许的最大稳态误差，为期望带宽，M为灵敏度峰值(一般情况下A=0.01,M=2)。从控制器设计的方面来说，“t的倒数为回路成形期望灵敏度的上限，'二影响控制器的输出u。有些情况下希望对补灵敏度函数丁=GI'(I+GI')1进行回路成形设计(在图1中增加一个输出3=丁)。一种选择方法为:Wt=s+/M乂4s+此函数与'涵数成山=",轴对

4、称，见图3所示。图中参数设置为A=0.01(=-40dB),M=2(=6dB),5=1iad/tie<?.-4-JD2410101D101DFrequencyrad/secFigure3:InverseofweightingfunctionsWsandWt3子系统的实现在Matlab中有几种方式得到G,WS和W。例如ControlSystemToolbox提SKS供的ss,tf和zpk等函数。Mu-tools也提供了诸如pck,nd2sys,zp2sys等函数,也可以用mksys和tree等方法。需要注意的是,Mu-tools提供了一种与ControlSystemToolbox不一样的表

5、达方式：系统矩阵(systemmatrix)。ControlSystemToolbox里面可以写成W：；IA,E八,D),对mu-tools则不适用。Mu-tools：AJBJCJD=ssdata(Gcst):Gmu=pck(AJBC,D).A,BjC,D=unpck(Gmu):Gcst=ss(A3B,C?D).4广义系统P的实现广义系统P也有多种产生方式。下面我列举了五种：(1) 直接写出传递函数矩阵，见2节。我倾向于利用mu-tools。如果后面需要转化为状态空间模型，需要利用minreal之类的函数得到最小化实现。重要的函数有：sbs(side-by-side),abv(above),m

6、mult(multiply),minv(inverse)。(2) 写出状态空间系数矩阵A,B,C,D，然后利用3节中的指令：P=pck(A,B,C,D)。(3) 利用sysic函数(systeminterconnect)o先在一个m文件中利用mu-tools设置了子系统，然后利用该函数将子系统互连起来。(4) 利用sconnect函数，这是LMI-tools提供的函数。子系统、输入、输出以参数的形式传递，sconnect返回互连的系统(5) 利用RCTv3.0.1提供的iconnect函数，该函数功能与sysic类似这些方法中，我倾向于sysic和iconnect函数，因为它们使用灵活，能够搭

7、建方法1和方法2难于搭建的复杂系统。一般情况下，得到均衡实现形式以避免数值计算问题是一个很好的想法，均衡实现形式可以利用ControlSystemToolbox提供的balreal实现。5控制器设计用、.设计问题是寻找一控制器K,使之稳定系统G,并使下列的范数最小：|旧（只K）|g=|WkSsKS有很多种得到Hw控制器的方法，例如hinfsyn,hinfric,hinflmi，这些函数将P作为输入，并以系统矩阵（mu-tools）方式表达。RCTv3.0.1提供了mixsyn函数，耳将二作为输入（H'F为补灵敏度函数权重），并不需要事先得到广义系统模型P。这些方法的主要区别在于是否用到

8、Riccati方程和Y迭代或者线性矩阵不等式来求解最优化问题。LMI方法不需要求解Riccati时设定的假设条件。另外还有些指令：ncfsyn和loopsyn（对开环传递函数L=GK进行Hw回路成形设计），hinfmix和msfsyn（多目标）。具体请查阅手册。6结果分析当控制器设计好后，需要对结果进行分析，这时可以利用ControlSystemToolbox提供的函数，例如利用lsim，step（阶跃响应）,bode（伯德图）,sigma（奇异值），freqresp（频域响应）等函数对传递函数S，KS，T，K，GK进行分析。Mu-tools提供的函数有：trsp（时域响应）,frsp（频域响

9、应）,vsvd（奇异值）,vplot。7结论从前面可见，有很多种进行Hw控制器设计的方法。为了避免混淆，我建议尽可能使用mu-tools和RCT。如果你知道ControlSystemToolbox存在某个函数,那么mu-tools也很可能存在具有相同功能的函数，只不过名称稍有些不同。如果你清楚自己的目的，但不知道函数的名字，那么建议你浏览帮助手册中的函数索引。希望此处简短的介绍能够帮助你进行Hw控制器的设计。好运！附录%下面的代码展示了如何在Matlab中进行H-infinity控制器的设计。此处举的例子与混合%S/KS问题有些不同。此处用到的模型和权重函数见Skogestad和Postlet

10、hwaite,1996，%ed.1,p.6O.权重函数并不是“最优”的。%大部分函数来自mu-tools,一些来自lmi-tools。mu-tools和lmi-tools均包含在RCT%v3.0.1中。%-JorgenJohnsen14.12.06%建立子系统%Plant:G=200/(10s+l)(0.05s+l)A2)%方法1：直接方法，利用mu-toolsG=nd2sys(1,conv(10,1,conv(0.05,1,0.051),200);%方法2：controlsystemtoolboxs=tf('s');Gcst=200/(10*s+1)*(0.05*s+1)A2

11、);a,b,c,d=ssdata(balreal(Gcst);G=pck(a,b,c,d);%权重：Ws=(s/M+wO)/(s+wO*A),Wks=1M=1.5;w0=10;A=1.e-4;Ws=nd2sys(1/Mw0,1w0*A);Wks=1;%建立广义系统P%方法0：直接方法%/z1/Ws-Ws*G/r%|z2|=|0Wks|%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)

12、;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|=|0W2|%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,l);q,p=s

13、ize(D);%原文此处为p,q=size(D);%全系统的ss实现Ap=Alzeros(nl,n2)-Bl*C;zeros(n2,nl)A2zeros(n2,n);zeros(n,nl)zeros(n,n2)A;Bp=Bl-Bl*D;zeros(n2,p)B2;zeros(n,p)B;Cp=Clzeros(ql,n2)-Dl*C;zeros(q2,nl)C2zeros(q2,n);zeros(q,nl)zeros(q,n2)-C;Dp=Dl-Dl*D;zeros(q2,p)D2;eye(p)-D;%得到均衡实现形式，以减少可能产生的计算问题Apb,Bpb,Cpb,Dpb=ssdata(bal

14、real(ss(Ap,Bp,Cp,Dp);Pl=pck(Apb,Bpb,Cpb,Dpb);%Pl与P0数值相近，但符号有差别%建立广义系统P%方法2：利用sysic函数systemnames='GWsWks'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='

15、;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-toolsSystem表达形式%注意2：iconnet函数仅适用于Ro

16、bustControlToolboxv3.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.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

17、);ab,bb,cb,db=ssdata(balreal(M.System);P4=pck(ab,bb,cb,db);%控制器的设计%下面使用的方法均基于广义系统P的System矩阵表达形式%选择你喜欢的控制器设计方法和广义系统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=hi

18、nflmi(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);%利用RCTv3.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

19、,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函数。也可以使用controltoolbox指令，例如用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(

20、9;singularvaluesofweightedclosedloopsystem')%得到传递函数矩阵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);%画奇异值%注意：如果愿意的话，可以

21、将vplot用plot代替，单位是dB。figure(2)vplot('liv,lm',Sw,'-',Tw,'-',GKw,'-.');title('sigma(S(jw)(solid),sigma(T(jw)(dashed)andsigma(GK(jw)(dashdot),')xlabel('Frequencyrad/sec');ylabel('Amplitude')figure(3)vplot('liv,lm',Kw)title('sigma(K(jw)&

22、#39;)xlabel('Frequencyrad/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)andsigma(WsA-1(jw)(dashed),')xlabel('Freq

23、uencyrad/sec');ylabel('Amplitude')figure(5)vplot('liv,lm',KSw,'-',KSdw,'-');title('sigma(KS(jw)(solid)andsigma(WksA-1(jw)(dashed),')xlabel('Frequencyrad/sec');ylabel('Amplitude')%最后，阶跃响应reference=1;tfinal=1;step=0.01;y=trsp(T,reference,tfinal,step);u=trsp(KS,reference,tfinal,step);figure(6)