%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')