利用matlab实现Hinfinity鲁棒控制.docx
《利用matlab实现Hinfinity鲁棒控制.docx》由会员分享，可在线阅读，更多相关《利用matlab实现Hinfinity鲁棒控制.docx（14页珍藏版）》请在冰豆网上搜索。
利用matlab实现Hinfinity鲁棒控制
利用Matlab实现H∞控制
Prof.Dr.-Ing.F.Allgwer
InstituteforSystemsTheoryandAutomaticControl
http:
//www.ist.uni-stuttgart.de/education/courses/robust
1引言
H∞控制器设计原理容易理解，难点在于编程。
这里简单介绍Matlab里面几个相关函数的用法，希望能帮助你设计第一个H∞控制器。
Matlab提供了很多H∞设计函数，与H∞设计相关的几个重要的工具箱有：
ControlSystemToolbox，mu-AnalysisandSynthesisToolbox（mu-tools），RobustControlToolbox（RCT）和LMIControlToolbox。
Matlab7.0之后的版本中，LMI和mu-tools都包含在RCTv3.0.1中，Matlab7.0之前的版本中这些工具箱是独立的。
本文中用到的函数都写在了一个m文件中（见附录），也可以从网站下载。
利用混合S/KS问题说明H∞相关函数的用法。
首先回顾这个问题。
2回路成形传递函数
混合S/KS问题可用图1来说明。
从w到z的闭环传递函数
可以表示为
广义过程模型P（s）（见图2）可以表示为
假设上面几个状态空间变量具有如下的形式：
于是可以得到P（s）的一个可能的状态空间实现形式：
WS和WKS为调整参数。
一种选择方法为：
其中A<1为允许的最大稳态误差，
为期望带宽，M为灵敏度峰值（一般情况下A=0.01，M=2）。
从控制器设计的方面来说，
的倒数为回路成形期望灵敏度的上限，
影响控制器的输出u。
有些情况下希望对补灵敏度函数
进行回路成形设计（在图1中增加一个输出
）。
一种选择方法为：
此函数与
函数成
轴对称，见图3所示。
图中参数设置为A=0.01（=-40dB），M=2（=6dB），
。
3子系统的实现
在Matlab中有几种方式得到G，WS和WKS。
例如ControlSystemToolbox提供的ss，tf和zpk等函数。
Mu-tools也提供了诸如pck，nd2sys，zp2sys等函数，也可以用mksys和tree等方法。
需要注意的是，Mu-tools提供了一种与ControlSystemToolbox不一样的表达方式：
系统矩阵（systemmatrix）。
ControlSystemToolbox里面可以写成
，对mu-tools则不适用。
Mu-tools：
4广义系统P的实现
广义系统P也有多种产生方式。
下面我列举了五种：
（1）直接写出传递函数矩阵，见2节。
我倾向于利用mu-tools。
如果后面需要转化为状态空间模型，需要利用minreal之类的函数得到最小化实现。
重要的函数有：
sbs（side-by-side），abv（above），mmult（multiply），minv（inverse）。
（2）写出状态空间系数矩阵A,B,C,D，然后利用3节中的指令：
P=pck（A,B,C,D）。
（3）利用sysic函数（systeminterconnect）。
先在一个m文件中利用mu-tools设置了子系统，然后利用该函数将子系统互连起来。
（4）利用sconnect函数，这是LMI-tools提供的函数。
子系统、输入、输出以参数的形式传递，sconnect返回互连的系统。
（5）利用RCTv3.0.1提供的iconnect函数，该函数功能与sysic类似。
这些方法中，我倾向于sysic和iconnect函数，因为它们使用灵活，能够搭建方法1和方法2难于搭建的复杂系统。
一般情况下，得到均衡实现形式以避免数值计算问题是一个很好的想法，均衡实现形式可以利用ControlSystemToolbox提供的balreal实现。
5控制器设计
设计问题是寻找一控制器K，使之稳定系统G，并使下列的H∞范数最小：
有很多种得到H∞控制器的方法，例如hinfsyn，hinfric，hinflmi，这些函数将P作为输入，并以系统矩阵（mu-tools）方式表达。
RCTv3.0.1提供了mixsyn函数，将
作为输入（
为补灵敏度函数权重），并不需要事先得到广义系统模型P。
这些方法的主要区别在于是否用到Riccati方程和γ迭代或者线性矩阵不等式来求解最优化问题。
LMI方法不需要求解Riccati时设定的假设条件。
另外还有些指令：
ncfsyn和loopsyn（对开环传递函数L=GK进行H∞回路成形设计），hinfmix和msfsyn（多目标）。
具体请查阅手册。
6结果分析
当控制器设计好后，需要对结果进行分析，这时可以利用ControlSystemToolbox提供的函数，例如利用lsim，step（阶跃响应），bode（伯德图），sigma（奇异值），freqresp（频域响应）等函数对传递函数S，KS，T，K，GK进行分析。
Mu-tools提供的函数有：
trsp（时域响应），frsp（频域响应），vsvd（奇异值），vplot。
7结论
从前面可见，有很多种进行H∞控制器设计的方法。
为了避免混淆，我建议尽可能使用mu-tools和RCT。
如果你知道ControlSystemToolbox存在某个函数，那么mu-tools也很可能存在具有相同功能的函数，只不过名称稍有些不同。
如果你清楚自己的目的，但不知道函数的名字，那么建议你浏览帮助手册中的函数索引。
希望此处简短的介绍能够帮助你进行H∞控制器的设计。
好运！
附录
%下面的代码展示了如何在Matlab中进行H-infinity控制器的设计。
此处举的例子与混合
%S/KS问题有些不同。
此处用到的模型和权重函数见Skogestad和Postlethwaite，1996，
%ed.1,p.60.权重函数并不是“最优”的。
%
%大部分函数来自mu-tools，一些来自lmi-tools。
mu-tools和lmi-tools均包含在RCT
%v3.0.1中。
%-JorgenJohnsen14.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.051]））,200）;
%方法2：
controlsystemtoolbox
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/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）;
[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,1）;[q,p]=size（D）;%原文此处为[p,q]=size（D）;
%全系统的ss实现
Ap=[A1zeros（n1,n2）-B1*C;
zeros（n2,n1）A2zeros（n2,n）;
zeros（n,n1）zeros（n,n2）A];
Bp=[B1-B1*D;
zeros（n2,p）B2;
zeros（n,p）B];
Cp=[C1zeros（q1,n2）-D1*C;
zeros（q2,n1）C2zeros（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='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='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函数仅适用于RobustControlToolboxv3.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）;
%利用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,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（'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））;
%画奇异值
%注意：
如果愿意的话，可以将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（'Stepresponse'）;ylabel（'y'）
subplot（2,1,2）
vplot（'iv,d',u）
ylabel（'y'）;xlabel（'time'）
升级会员 