基于混合灵敏度的H∞鲁棒控制器的设计仿真

clcclose allK=31.31b0=50.88b1=57.81Ag=[0 1;-b0 -b1]Bg=[0;1]Cg=[K 0]Dg=[0];[num,den]=ss2tf(Ag,Bg,Cg,Dg) Gs=tf(num,den)[z,p,k]=tf2zp(num,den)Gs1=zpk(z,p,k)%W1S=tf([0.2 3],[1 0.001])%W2S=tf([0.002],[1])%W3S=tf([0.002 0.00001],[1]) %W1S=tf([0.2 4],[1 0.002])%W2S=tf([0.003],[1])%W3S=tf([0.001 0.00001],[1]) %加权函数W1S=tf([0.2 10],[1 0.001])W2S=tf([0.001],[1])W3S=tf([0.001 0.00001],[1]) W3NUM=W3S.num{1,1};W3DEN=W3S.den{1,1};i=0;j=0;a=length(W3NUM);b=length(W3DEN);for c=1:1:aif W3NUM(c-i)~=0break;elseW3NUM(c-i)=[];i=i+1;endendfor d=1:1:bif W3DEN(d-j)~=0break;elseW3DEN(d-j)=[];j=j+1;endendGg=ss(Ag,Bg,Cg,Dg);[Aw1,Bw1,Cw1,Dw1]=tf2ss(W1S.num{1,1},W1S.den{1,1});[Aw2,Bw2,Cw2,Dw2]=tf2ss(W2S.num{1,1},W2S.den{1,1});if Polyorder(W3NUM)>Polyorder(W3DEN)[Q,R]=deconv(W3NUM,W3DEN)W3poly=Q[Aw3,Bw3,Cw3,Dw3]=tf2ss(R,W3DEN)else[Aw3,Bw3,Cw3,Dw3]=tf2ss(W3NUM,W3DEN)W3poly=[]endW1=[Aw1,Bw1;Cw1,Dw1]W2=[Aw2,Bw2;Cw2,Dw2]W3=[Aw3,Bw3;Cw3,Dw3][A,B1,B2,C1,C2,D11,D12,D21,D22]=augss(Ag,Bg,Cg,Dg,Aw1,Bw1,Cw1,Dw1,Aw2,Bw2,Cw2,Dw2,A w3,Bw3,Cw3,Dw3,W3poly)G=[A,B1,B2;C1,D11,D12;C2,D21,D22];[Acp,Bcp,Ccp,Dcp,Ac1,Bc1,Cc1,Dc1]=hinf(A,B1,B2,C1,C2,D11,D12,D21,D22)[knum,kden]=ss2tf(Acp,Bcp,Ccp,Dcp) %控制器的状态空间形式转换成传递函数形式K=tf(knum,kden) %控制器的传递函数[z,p,k]=tf2zp(knum,kden) %控制器的传递函数转换成零极点增益形式Kzpk=zpk(z,p,k)%{[gnum,gden]=ss2tf(Ag,Bg,Cg,Dg)Gs=tf(gnum,gden) %标称系统的传递函数P=nd2sys(gnum,gden,1) %标称系统的系统矩阵形式K1=nd2sys(knum,kden,1) %鲁棒控制器的系统矩阵形式[type,out,in,n]=minfo(P)I=eye(out)S=minv(madd(I,mmult(P,K1))) %灵敏度函数(1+GK)^-1 ，其中P和K1都是系统矩阵的形式T=msub(I,S)W1S1=nd2sys(W1S.num{1,1},W1S.den{1,1},1)INVW1S1=minv(W1S1) %加权函数W1的逆W1^-1INVW3S1=nd2sys(W3S.den{1,1},W3S.num{1,1},1) %加权函数W3的逆W3^-1w=logspace(-5,5,500)Sw=vsvd(frsp(S,w)) %灵敏度函数的频域响应的奇异值Tw=vsvd(frsp(T,w)) %补灵敏度函数的频域响应的奇异值INVW1S1w=vsvd(frsp(INVW1S1,w))INVW3S1w=vsvd(frsp(INVW3S1,w))figure(1)vplot('liv,lm',Sw,'r-',INVW1S1w,'b--')title('Singular values of sensivity function S and W1^{-1}')set(gca,'color','w')xlabel('Frequency(rad/sec)')ylabel('Amplitude')grid onfigure(2)vplot('liv,lm',Tw,'r-',INVW3S1w,'b--')title('Singular values of complementary sensivity function T and W3^{-1}')set(gca,'color','w')xlabel('Frequency(rad/sec)')ylabel('Amplitude')grid on%}w=logspace(-5,5,500);figure(3)bode(W1S,'-b',W3S,'-r',w)legend('W1','W3')title('Bode Diagram Of Weighed Function W1 and W3')grid on[Acg,Bcg,Ccg,Dcg]=series(Acp,Bcp,Ccp,Dcp,Ag,Bg,Cg,Dg)[As,Bs,Cs,Ds]=feedbk(Acg,Bcg,Ccg,Dcg,1) %灵敏度函数的状态空间(1+GK)^-1 [At,Bt,Ct,Dt]=feedbk(Acg,Bcg,Ccg,Dcg,2) %补灵敏度函数的状态空间GK(1+GK)^-1svs=sigma(As,Bs,Cs,Ds,w);svs=20*log10(svs);[Aw1i,Bw1i,Cw1i,Dw1i]=unpck(minv(pck(Aw1,Bw1,Cw1,Dw1)));svw1i=sigma(Aw1i,Bw1i,Cw1i,Dw1i,w);svw1i= 20*log10(svw1i);figure(4)semilogx(w,svw1i,'b--',w,svs,'r-')title('Singular values of sensivity function S and W1^{-1}')xlabel('Frequency(rad/sec)')ylabel('Amplitude(db)')legend('W1^{-1}','S');grid on[Aw3i,Bw3i,Cw3i,Dw3i]=tf2ss(W3S.den{1,1},W3S.num{1,1})svw3i=sigma(Aw3i,Bw3i,Cw3i,Dw3i,w);svw3i=20*log10(svw3i);svt=sigma(At,Bt,Ct,Dt,w);svt=20*log10(svt);figure(5)semilogx(w,svw3i,'b--',w,svt,'r-')title('Singular values of complementary sensivity function T and W3^{-1}')xlabel('Frequency(rad/sec)')ylabel('Amplitude(db)')legend('W3^{-1}','T');grid onsvtt=sigma(Ac1,Bc1,Cc1,Dc1,1,w);svtt=20*log10(svtt);figure(6)semilogx(w,svtt,'b-')xlabel('Frequency(rad/sec)')ylabel('Amplitude(db)')title('Frequency Characteristic of Closed-Loop System Twz') grid onsvcp=sigma(Acp,Bcp,Ccp,Dcp,1,w);svcp=20*log10(svcp);figure(7)semilogx(w,svcp,'b-')xlabel('Frequency(rad/sec)')ylabel('Amplitude(db)')title('Frequency Characteristic of H∞Controller')grid ont=[0:0.01:15];setvalue=8y=setvalue*step(At,Bt,Ct,Dt,1,t);figure(8)plot(t,y,'-b')axis([0,10,0,10])xlabel('Time(s)')ylabel('L/Min)')grid onfigure(9)bode(Ag,Bg,Cg,Dg,1,w)title('Bode Diagram Of The Plant')grid on[cnum,cden]=ss2tf(At,Bt,Ct,Dt,1)figure(10)step(Ag,Bg,Cg,Dg,1)grid on。