利用matlab實(shí)現(xiàn)H-infinity魯棒控制

1、利用Matlab實(shí)現(xiàn)H8控制Prof.Dr.-Ing.F.AllgwerInstituteforSystemsTheoryandAutomaticControlhttp:/www.ist.uni-stuttgart.de/education/courses/robust1引言控制器設(shè)計(jì)原理容易理解，難點(diǎn)在于編程。這里簡單介紹Matlab里面幾個(gè)相關(guān)函數(shù)的用法，希望能幫助你設(shè)計(jì)第一個(gè)控制器。Matlab提供了很多設(shè)計(jì)函數(shù)，與設(shè)計(jì)相關(guān)的幾個(gè)重要的工具箱有：ControlSystemToolbox，mu-AnalysisandSynthesisToolbox(mu-tools)，RobustCont

2、rolToolbox(RCT)和LMIControlToolbox。Matlab7.0之后的版本中，LMI和mu-tools都包含在RCTv3.0.1中，Matlab7.0之前的版本中這些工具箱是獨(dú)立的。本文中用到的函數(shù)都寫在了一個(gè)m文件中(見附錄)，也可以從網(wǎng)站下載。利用混合S/KS問題說明相關(guān)函數(shù)的用法。首先回顧這個(gè)問題。2回路成形傳遞函數(shù)混合S/KS問題可用圖1來說明。wI<sKS從w到Z的閉環(huán)傳遞函數(shù)=可以表示為廣義過程模型P(s)(見圖2)可以表示為原文：ControlinMatlab.2011/12/8譯假設(shè)上面幾個(gè)狀態(tài)空間變量具有如下的形式：于是可以得到P(s)的一個(gè)可能的

3、狀態(tài)空間實(shí)現(xiàn)形式:比0-BSCBs-BSD0Aks00ks00A0B0-DSCDs-DSD0Cks00D曉00-CI-Dws和wks為調(diào)整參數(shù)。一種選擇方法為:拆/'Z+TTr丄；=;=const,s+其中Avl為允許的最大穩(wěn)態(tài)誤差，為期望帶寬，M為靈敏度峰值(一般情況下A=0.01,M=2)。從控制器設(shè)計(jì)的方面來說，“t的倒數(shù)為回路成形期望靈敏度的上限，'二影響控制器的輸出u。有些情況下希望對(duì)補(bǔ)靈敏度函數(shù)丁=GI'(I+GI')1進(jìn)行回路成形設(shè)計(jì)(在圖1中增加一個(gè)輸出3=丁)。一種選擇方法為:Wt=s+/M乂4s+此函數(shù)與'涵數(shù)成山=",軸對(duì)

4、稱，見圖3所示。圖中參數(shù)設(shè)置為A=0.01(=-40dB),M=2(=6dB),5=1iad/tie<?.-4-JD2410101D101DFrequencyrad/secFigure3:InverseofweightingfunctionsWsandWt3子系統(tǒng)的實(shí)現(xiàn)在Matlab中有幾種方式得到G,WS和W。例如ControlSystemToolbox提SKS供的ss,tf和zpk等函數(shù)。Mu-tools也提供了諸如pck,nd2sys,zp2sys等函數(shù),也可以用mksys和tree等方法。需要注意的是,Mu-tools提供了一種與ControlSystemToolbox不一樣的表

5、達(dá)方式：系統(tǒng)矩陣(systemmatrix)。ControlSystemToolbox里面可以寫成W：；IA,E八,D),對(duì)mu-tools則不適用。Mu-tools：AJBJCJD=ssdata(Gcst):Gmu=pck(AJBC,D).A,BjC,D=unpck(Gmu):Gcst=ss(A3B,C?D).4廣義系統(tǒng)P的實(shí)現(xiàn)廣義系統(tǒng)P也有多種產(chǎn)生方式。下面我列舉了五種：(1) 直接寫出傳遞函數(shù)矩陣，見2節(jié)。我傾向于利用mu-tools。如果后面需要轉(zhuǎn)化為狀態(tài)空間模型，需要利用minreal之類的函數(shù)得到最小化實(shí)現(xiàn)。重要的函數(shù)有：sbs(side-by-side),abv(above),m

6、mult(multiply),minv(inverse)。(2) 寫出狀態(tài)空間系數(shù)矩陣A,B,C,D，然后利用3節(jié)中的指令：P=pck(A,B,C,D)。(3) 利用sysic函數(shù)(systeminterconnect)o先在一個(gè)m文件中利用mu-tools設(shè)置了子系統(tǒng)，然后利用該函數(shù)將子系統(tǒng)互連起來。(4) 利用sconnect函數(shù)，這是LMI-tools提供的函數(shù)。子系統(tǒng)、輸入、輸出以參數(shù)的形式傳遞，sconnect返回互連的系統(tǒng)(5) 利用RCTv3.0.1提供的iconnect函數(shù)，該函數(shù)功能與sysic類似這些方法中，我傾向于sysic和iconnect函數(shù)，因?yàn)樗鼈兪褂渺`活，能夠搭

7、建方法1和方法2難于搭建的復(fù)雜系統(tǒng)。一般情況下，得到均衡實(shí)現(xiàn)形式以避免數(shù)值計(jì)算問題是一個(gè)很好的想法，均衡實(shí)現(xiàn)形式可以利用ControlSystemToolbox提供的balreal實(shí)現(xiàn)。5控制器設(shè)計(jì)用、.設(shè)計(jì)問題是尋找一控制器K,使之穩(wěn)定系統(tǒng)G,并使下列的范數(shù)最?。簗舊（只K）|g=|WkSsKS有很多種得到Hw控制器的方法，例如hinfsyn,hinfric,hinflmi，這些函數(shù)將P作為輸入，并以系統(tǒng)矩陣（mu-tools）方式表達(dá)。RCTv3.0.1提供了mixsyn函數(shù)，耳將二作為輸入（H'F為補(bǔ)靈敏度函數(shù)權(quán)重），并不需要事先得到廣義系統(tǒng)模型P。這些方法的主要區(qū)別在于是否用到

8、Riccati方程和Y迭代或者線性矩陣不等式來求解最優(yōu)化問題。LMI方法不需要求解Riccati時(shí)設(shè)定的假設(shè)條件。另外還有些指令：ncfsyn和loopsyn（對(duì)開環(huán)傳遞函數(shù)L=GK進(jìn)行Hw回路成形設(shè)計(jì)），hinfmix和msfsyn（多目標(biāo)）。具體請(qǐng)查閱手冊(cè)。6結(jié)果分析當(dāng)控制器設(shè)計(jì)好后，需要對(duì)結(jié)果進(jìn)行分析，這時(shí)可以利用ControlSystemToolbox提供的函數(shù)，例如利用lsim，step（階躍響應(yīng)）,bode（伯德圖）,sigma（奇異值），freqresp（頻域響應(yīng)）等函數(shù)對(duì)傳遞函數(shù)S，KS，T，K，GK進(jìn)行分析。Mu-tools提供的函數(shù)有：trsp（時(shí)域響應(yīng)）,frsp（頻域響

9、應(yīng)）,vsvd（奇異值）,vplot。7結(jié)論從前面可見，有很多種進(jìn)行Hw控制器設(shè)計(jì)的方法。為了避免混淆，我建議盡可能使用mu-tools和RCT。如果你知道ControlSystemToolbox存在某個(gè)函數(shù),那么mu-tools也很可能存在具有相同功能的函數(shù)，只不過名稱稍有些不同。如果你清楚自己的目的，但不知道函數(shù)的名字，那么建議你瀏覽幫助手冊(cè)中的函數(shù)索引。希望此處簡短的介紹能夠幫助你進(jìn)行Hw控制器的設(shè)計(jì)。好運(yùn)！附錄%下面的代碼展示了如何在Matlab中進(jìn)行H-infinity控制器的設(shè)計(jì)。此處舉的例子與混合%S/KS問題有些不同。此處用到的模型和權(quán)重函數(shù)見Skogestad和Postlet

10、hwaite,1996，%ed.1,p.6O.權(quán)重函數(shù)并不是“最優(yōu)”的。%大部分函數(shù)來自mu-tools,一些來自lmi-tools。mu-tools和lmi-tools均包含在RCT%v3.0.1中。%-JorgenJohnsen14.12.06%建立子系統(tǒng)%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);%權(quán)重：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;%建立廣義系統(tǒng)P%方法0：直接方法%/z1/Ws-Ws*G/r%|z2|=|0Wks|%v/I-G/u/%傳遞函數(shù)表達(dá)方法Z1=sbs(Ws,mmult(-1,Ws,G);Z2=sbs(0,Wks);V=sbs(1,mmult(-1,G);P0=abv(Z1,Z2,V);%通常情況下P0并不是最小實(shí)現(xiàn)，所以需要降階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);%此時(shí)得到變量為System類型%建立廣義系統(tǒng)P%方法1：直接方法%/z1/W1-W1*G/r%|z2|=|0W2|%v/I-G/u/%子系統(tǒng)的ss實(shí)現(xiàn)A,B,C,D=unpck(G);A1,B1,C1,D1=unpck(Ws);A2,B2,C2,D2=unpck(Wks);%計(jì)算不同子系統(tǒng)的輸入、輸出變量的個(gè)數(shù)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);%全系統(tǒng)的ss實(shí)現(xiàn)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;%得到均衡實(shí)現(xiàn)形式，以減少可能產(chǎn)生的計(jì)算問題Apb,Bpb,Cpb,Dpb=ssdata(bal

14、real(ss(Ap,Bp,Cp,Dp);Pl=pck(Apb,Bpb,Cpb,Dpb);%Pl與P0數(shù)值相近，但符號(hào)有差別%建立廣義系統(tǒng)P%方法2：利用sysic函數(shù)systemnames='GWsWks'inputvar='r(1);u(1)'%所有輸入均為標(biāo)量，r(2)為兩維信號(hào)outputvar='Ws;Wks;r-G'input_to_G='u'input_to_Ws='r-G'input_to_Wks='u'sysoutname='P2'cleanupsysic='

15、;yes'sysic%建立廣義系統(tǒng)P%方法3：利用sconnect函數(shù)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);%建立廣義系統(tǒng)P%方法4：利用iconnect函數(shù)%注意1：不再使用mu-toolsSystem表達(dá)形式%注意2：iconnet函數(shù)僅適用于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);%控制器的設(shè)計(jì)%下面使用的方法均基于廣義系統(tǒng)P的System矩陣表達(dá)形式%選擇你喜歡的控制器設(shè)計(jì)方法和廣義系統(tǒng)P%選擇廣義系統(tǒng)PP=P1;%P=P0;P=P1;P=P2;P=P3;P=P4;%然后設(shè)置一些參數(shù)(測量變量個(gè)數(shù)，輸入變量個(gè)數(shù)，gamma限制)nmeas=1;nu=1;gmn=0.5;gmx=20;tol=0.001;%控制器設(shè)計(jì)：選擇你喜歡的設(shè)計(jì)方法，不需要的注釋掉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進(jìn)行控制器的設(shè)計(jì)%通常不需要系統(tǒng)的傳遞函數(shù)表達(dá)式，但更需要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);%分析結(jié)果%注意：此處使用mu-tools函數(shù)。也可以使用controltoolbox指令，例如用series和%feedback進(jìn)行子系統(tǒng)連接，sigma或freqresp,svd和bode畫奇異值，step和lsim分析%時(shí)域響應(yīng)%畫(Weighted)閉環(huán)系統(tǒng)的奇異值w=logspace(-4,6,50);CLw=vsvd(frsp(CL,w);figure(1)vplot('liv,m',CLw)title(

20、9;singularvaluesofweightedclosedloopsystem')%得到傳遞函數(shù)矩陣type,out,in,n=minfo(G);I=eye(out);S=minv(madd(I,mmult(G,K);%靈敏度函數(shù)T=msub(I,S);%補(bǔ)靈敏度函數(shù)KS=mmult(K,S);%G的輸入GK=mmult(G,K);%回路傳遞函數(shù)%頻域的奇異值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')%是否滿足設(shè)計(jì)要求？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')%最后，階躍響應(yīng)reference=1;tfinal=1;step=0.01;y=trsp(T,reference,tfinal,step);u=trsp(KS,reference,tfinal,step);figure(6)