H Infinity contrl |
您所在的位置：网站首页 › mat2dec函数 › H Infinity contrl |
H Infinity contrl
|
H Infinity contrl ---LMI方法（Yalmip 求解）和Riccati 方法
参考书籍：Yu, Hai-Hua_ Duan, Guangren - LMIs in control systems _ analysis, design and applications (2013, CRC Press) - libgen.lc(一书的附录有详细介绍LMI的编写方式-有很多例子) 参考网站：http://www.doc88.com/p-8718494012962.html YALMIP Download - YALMIP 有问题可以上google论坛--在matlab 问答模块可以搜到网址--可以提问相关问题， 关于LMI 本身求解参考： 线性矩阵不等式(LMI)的_MATLAB求解- http://www.doc88.com/p-7874379836838.html Matlab中LMI(线性矩阵不等式)工具箱使用教程 https://www.doc88.com/p-304944065894.html 线性矩阵不等式的LMI工具箱求解-LMI三个函数的应用实例； https://www.doc88.com/p-3925042294217.html https://www.doc88.com/p-660150804881.html https://blog.csdn.net/qq_28093585/article/details/69358180 ceshi_Hinf_lmi  
clc;clear; close all;
A = [2 1 -2;1 -1 -3; 4 0 -1];
B = [1 0;0 3;3 1];
E1 = [1 0.2 -0.5]';
C = [2 1 -0.5];
D = [1 1];
E2 = [0.05]
P = 10e+7;
X0 = [8.7047 -6.5321 4.2500;
-6.5321 8.8601 -4.2821;
4.2500 -4.2821 5.0227];
W0 = [-5.2786 2.9364 -7.7991;
-3.4736 -0.8733 6.0925];
W =P*W0; X = X0*P;
K1 =W*inv(X)
Ac =A-B*K1;
Bc = B;
Cc = C;
Dc = D;
sys_cl = ss(Ac,Bc,Cc,Dc);
step(sys_cl)
% figure(1)
% t = 0:0.01:15;
% r =[0.5*ones(size(t));
% 0.5*ones(size(t))];
% [y,t,x]=lsim(sys_cl,r,t);
% plot(t,y)
%
%
% K2 = [1.2850 0.5978 -2.0283;
% -1.4101 -0.7807 1.4861]
% Ac =A-B*K2;
% Bc = B;
% Cc = C;
% Dc = D;
% sys_cl = ss(Ac,Bc,Cc,Dc);
%
% figure(2)
% t = 0:0.01:15;
% r =[0.5*ones(size(t));
% 0.5*ones(size(t))];
% [y,t,x]=lsim(sys_cl,r,t);
% plot(t,y)
View Code
%连续H2 optimal control
clc;clear;close all;
%% 书籍：LMIs in control systems _ analysis, design and applications中428页中例题的结果--
% 连续H2 state feedback control
% Step 1. Describe this LMIs in MATLAB 网址：http://www.doc88.com/p-8718494012962.html
%{
% Step 1. Describe this LMIs in MATLAB
% specify parameter matrices
A = [-3 -2 1;1 2 1;1 -1 -1]; B1 = [1; 0; 1];
B2 = [2 0;0 2;0 1]; C = [1 0 1;0 1 0]; D = eye(2);
% define the LMI system
setlmis ([]) % initialing
% declare the LMI variable
X = lmivar(1, [3 1]); W = lmivar(2,[2,3]);
Z = lmivar(1, [2 1]); rou = lmivar(1, [1 1]);
% #1 LMI, left-hand side
lmiterm ( [1 1 1 X],A, 1, 's');
lmiterm ( [1 1 1 W],B2, 1, 's');
lmiterm ( [1 1 1 0],B1*B1');
% #2 LMI, left-hand side
lmiterm ( [2 1 1 Z],1,-1); % the (1,1) block
lmiterm ( [2 1 2 W],D,1); % the (1,2) block
lmiterm ( [2 1 2 X],C,1); % the (1,2) block
lmiterm ( [2 2 2 X],1,-1); % the (2,2) block
% #3 LMI, left-hand side
lmiterm ( [3 1 1 Z],[1 0],[1 0]');
lmiterm ( [3 1 1 Z],[0 1],[0 1]');
% #3 LMI, right-hand side
lmiterm ( [3 1 1 rou],1,-1);
lmisys=getlmis; % finishing description
% Step 2. Solve this LMI problem and use dec2mat to get the
% corresponding matrix X
% c = mat2dec(lmisys,0,0,ones(2)); % obtain vector c
n=decnbr(lmisys);%获取LMI中决策变量的个数
c=zeros(n,1);
% for jj=1:n
% [Zj]=defcx(lmisys,jj,Z);
% c(jj)=trace(Zj)
% end
for jj=1:n
[pj]=defcx(lmisys,jj,rou);
c(jj)=(pj);
end
options=[1e-5 0 0 0 1];
%[copt,xopt] = mincx(lmisys,c);
[copt,xopt] = mincx(lmisys,c,options); % solve the optimization problem
X = dec2mat(lmisys,xopt,1); % extract matrix X
W = dec2mat(lmisys,xopt,2); % extract matrix W
Z = dec2mat(lmisys,xopt,3); % extract matrix Z
K = W*inv(X); % get the feedback gain
%}
%注：求出的结果为：[-1.00000000000000,5.44822811876415e-18,-1.00000000000000;
% 1.75680133525483e-16,-1.00000000000000,-7.35314080057690e-18]
% 书中的结果是[-1 0 -1;0 -1 0]
%%
%{
%% 书籍：LMIs in control systems _ analysis, design and applications中259页中例题的结果
% Step 1. Describe this LMIs in MATLAB 网址：http://www.doc88.com/p-8718494012962.html
% specify parameter matrices
A = [-3 -2 1;1 2 1;1 -1 -1]; B1 = [3; 0; 1];%干扰矩阵
B2 = [2 0;0 2;0 1]; %控制矩阵
C = [1 0 1;0 1 1]; D = [1 1;0 1];
% define the LMI system
setlmis ([]) % initialing
% declare the LMI variable
X = lmivar(1, [3 1]); W = lmivar(2,[2,3]);
Z = lmivar(1, [2 1]); rou = lmivar(1, [1 1]);
% #1 LMI, left-hand side
lmiterm ( [1 1 1 X],A, 1, 's');
lmiterm ( [1 1 1 W],B2, 1, 's');
lmiterm ( [1 1 1 0],B1*B1');
% #2 LMI, left-hand side
lmiterm ( [2 1 1 Z],1,-1); % the (1,1) block
lmiterm ( [2 1 2 W],D,1); % the (1,2) block
lmiterm ( [2 1 2 X],C,1); % the (1,2) block
lmiterm ( [2 2 2 X],1,-1); % the (2,2) block
% #3 LMI, left-hand side
lmiterm ( [3 1 1 Z],[1 0],[1 0]');
lmiterm ( [3 1 1 Z],[0 1],[0 1]');
% #3 LMI, right-hand side
lmiterm ( [3 1 1 rou],1,-1);
lmisys=getlmis; % finishing description
% Step 2. Solve this LMI problem and use dec2mat to get the
% corresponding matrix X
% c = mat2dec(lmisys,0,0,ones(2)); % obtain vector c
n=decnbr(lmisys);%获取LMI中决策变量的个数
c=zeros(n,1);
% for jj=1:n
% [Zj]=defcx(lmisys,jj,Z);
% c(jj)=trace(Zj)
% end
for jj=1:n
pj=defcx(lmisys,jj,rou);
c(jj)=(pj);
end
%options=[1e-5 0 0 0 1];
[copt,xopt] = mincx(lmisys,c);
%[copt,xopt] = mincx(lmisys,c,options) % solve the optimization problem
X = dec2mat(lmisys,xopt,1) % extract matrix X
W = dec2mat(lmisys,xopt,2) % extract matrix W
Z = dec2mat(lmisys,xopt,3) % extract matrix Z
K = W*inv(X) % get the feedback gain
%}
%注：求出的结果为：K =[-0.911192567348860,1.75193380502826,-0.249064838067320;
% -0.106341894142934,-1.90039834212258,-0.701758897247713];
% pou =3.349135e-04;
%书中的结果：[-0.9112 1.7519 -0.2491; -0.1063 -1.9004 -0.7018];
%% 采用yample优化工具箱的求出结果；
%示例 https://github.com/eoskowro/LMI
% Arizona State University
% MAE 598 LMIs in Control Systems
%{
clear;close all; clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LMI for Opimal FSF H_2 Control
% Bring in Data
A = [-1 1 0 1 0 1;...
-1 -2 -1 0 0 1;...
1 0 -2 -1 1 1;...
-1 1 -1 -2 0 0;...
-1 -1 1 1 -2 -1;...
0 -1 0 0 -1 -3];
B = [0 -1 -1;...
0 0 0;...
-1 1 1;...
-1 0 0;...
0 0 1;...
-1 1 1];
C = [0 1 0 -1 -1 -1;...
0 0 0 -1 0 0;...
1 0 0 0 -1 0];
D = zeros(3);
% Determine sizes
ns = size(A,1); % number of states
nai = size(B,2); % number of actuator inputs
nmo = size(C,1); % number of measured outputs
nd = nai+nmo; % number of disturbances
nro = nmo+nai; % number of regulated outputs
eta = 0.0001;
% 9 Matrix Representation
B1 = [B zeros(ns,nmo)];
B2 = B;
C1 = [C;...
zeros(nai,ns)];
C2 = C;
D11 = [D zeros(nai);...
zeros(nmo) zeros(nmo)];
D12 = [D;...
eye(nmo)];
D21 = [D eye(nmo)];
D22 = D;
% Settings
opt = sdpsettings('verbose',0,'solver','sedumi'); %Create/alter solver options structure.
% Define Variables
gam= sdpvar(1);
X = sdpvar(6);
W = sdpvar(6);
Z = sdpvar(3,6);
% Define matricies
mat1 = [A B2]*[X;Z] + [X Z']*[A';B2']+ B1*B1';
mat2 = [X (C1*X+D12*Z)';C1*X+D12*Z W];
mat3 = trace(W);
% Define Constraints
F = [];
F = [F, X>=eta*eye(6)];
F = [F, mat1=eta*eye(12)];
F = [F, mat3=0; Z>=0 ;trace(Z)= 0; P>=0];
optimize(F, gam);
Kd = value(Fd)*inv(value(P))
Hinf_norm = (value(gam))
L1 =eig(A-Bd2*Kd)
%% simulation
% sysd = ss(A,Bd2,Cd2,Dd22);
% feedin = [1 2];
% feedout = [1 2 3 4];
% sys_cl = feedback(sysd,Kd,feedin,feedout)
%
% Ac = A-Bd2*Kd;
% Bc = Bd2;
% Cc = Cd2;
% Dc = Dd22;
% sys_cl = ss(Ac,Bc,Cc,Dc)
% T = 0:0.1:10;
% step(sys_cl,T)
% clear;
% clc;
A = [0.1 0.2 0.3;
0.4 0.5 0.6;
0.7 0.8 0.9];
Bd1 = ones(3,1);
Bd2 = ones(3,1);
Cd1 = ones(1,3);
Dd12 = 1;
Dd11 = 0;
P = sdpvar(3);
Z = sdpvar(1);
Fd = sdpvar(1,3);
gam = sdpvar(1);
mat1 = [P A*P+Bd2*Fd Bd1 zeros(3,1);
(A*P+Bd2*Fd)' P zeros(3,1) P*Cd1'-Fd'*Dd12';
Bd1' zeros(1,3) gam*eye(1) Dd11';
zeros(1,3) (P*Cd1'-Fd'*Dd12')' Dd11 gam*eye(1)];
F = [mat1 >= 0; P>=0];
optimize(F, gam);
Kd = value(Fd)*inv(value(P))
Hinf_norm = (value(gam))
L2 =eig(A-Bd2*Kd)
View Code
H2 状态反馈控制--模型参数例子来于书本：robust control design _269页--三个例子分开运行
clc;clear;close all; %文中有三个例子，需要分开运行才可以
%% 模型参数例子来于书本：robust control design _269页
%{
A=[-5 2 -4 ;0 -3 0; 0 7 -1];
B1=[7 -3 1]'; B2=[6 8 -5]';
C1=[-2 9 4]; C2=[6 3 -1];
D11=[0]; D12=[1]; D21=[2]; D22=[0];
G=ss(A,B2,C2,D22);
P =ltisys(A,[B1 B2],[C1;C2],[D11 D12;D21 D22]);
P1 =pck(A,[B1 B2],[C1;C2],[D11 D12;D21 D22]);
[Aa,Bb,Cc,Dd]=branch(P);
sys=ss(Aa,Bb,Cc,Dd);
po=pole(sys)
nmeas = 1;
ncont = 1;
[K,CL,gamma] = h2syn(sys,nmeas,ncont);
[Ak,Bk,Ck,Dk]=branch(K);
GK=ss(Ak,Bk,Ck,Dk);
K_tf=zpk();
[Knum,Kden]=ss2tf(Ak,Bk,Ck,Dk);
K_tf=tf(Knum,Kden)
% [Aa,Bb,Cc,Dd]=branch(P);
% sys=ss(Aa,Bb,Cc,Dd);
% P_tf=zpk(ss(Aa,Bb,Cc,Dd))
step(feedback(G*GK,-1),30)
% clsys =slft(P,K);
% splot(clsys,'st')
%}
%{
%% matlab帮助中的h2syn函数的例子
A = [5 6 -6
6 0 5
-6 5 4];
B = [0 4 0 0
1 1 -2 -2
4 0 0 -3];
C = [-6 0 8
0 5 0
-2 1 -4
4 -6 -5
0 -15 7];
D = [0 0 0 0
0 0 0 1
0 0 0 0
0 0 3 6
8 0 -7 0];
P = ss(A,B,C,D);
P1 = pck(A,B,C,D);
pole(P)
nmeas = 2;
ncont = 1;
[K,CL,gamma] = h2syn(P,nmeas,ncont);
pole(CL)
step(CL)
clsys =slft(P,K);
splot(clsys,'st')
[A,B1,B2,C1,C2,D11,D12,D21,D22]=branch(P)
%}
%% 原来旧函数的应用--https://wenku.baidu.com/view/d79bb3cb0722192e44
A =[0 1 0 0;
-5000 -100/3 500 100/3;
0 -1 0 1;
0 100/3 -4 -60];
B =[0 25/3 0 -1]';
C=[0 0 1 0];
D=0;
G =ss(A,B,C,D);
W1=[0 100; 1 1];
W2=1e-5;
W3=[1 0; 0 1000];
T_ss=augtf(G,W1,W2,W3);
[Aa,Bb,Cc,Dd]=branch(T_ss);
[Aa,Bb,Cc,Dd]=ssdata(T_ss);
[a,b1,b2,c1,c2,d11,d12,d21,d22]=branch(T_ss);%当使用augtf函数获取增广系统矩阵表达式后，
%可用branch函数获得各个增光系统矩阵的参数
P =rt2smat(T_ss) %变换成系统矩阵P
K_h2=h2lqg(T_ss) %获得H2最优控制器
[AK0,BK0,CK0,DK0]=branch(K_h2);
[AK0,BK0,CK0,DK0]=ssdata(K_h2);%两个函数均可实现这个功能
K_H2_tf=zpk(ss(AK0,BK0,CK0,DK0))
K_Hinf=hinf(T_ss) %获得H00优控制器
[AK1,BK1,CK1,DK1]=branch(K_Hinf);
K_Hinf_tf=zpk(ss(AK1,BK1,CK1,DK1))
[gamma,K_opt]=hinfopt(T_ss); %获得H00最优控制器
[AK2,BK2,CK2,DK2]=branch(K_opt);
K_opt_tf=zpk(ss(AK2,BK2,CK2,DK2))
% 绘制闭环节约响应下系统的开环bode图
bode(G*K_H2_tf,'-',G*K_Hinf_tf,'--',G*K_opt_tf,':')
% 闭环阶跃响应曲线
figure(1)
step(feedback(G,K_H2_tf,-1),'r-')
figure(2)
step(feedback(G*K_H2_tf,1),'r-')
% step(feedback(G*K_H2_tf,1),'r-',feedback(G*K_Hinf_tf,1),'g--',feedback(G*K_opt_tf,1),':')
%系统矩阵变换函数 -参见：https://wenku.baidu.com/view/d79bb3cb0722192e44
function P =rt2smat(A,B1,B2,C1,C2,D11,D12,D21,D22)
if nargin ==1,
[A,B1,B2,C1,C2,D11,D12,D21,D22]=branch(A);
end
n =length(A);
P=[A,B1,B2;C1,D11,D12;C2,D21,D22];
P(size(P,1)+1,size(P,2)+1)=-inf;
m=size(P,2);
P(1,m)=n;
end
%}
View Code
H_infinity_looping_shape_control---参考文章：BACHELOR’S FINAL PROJECT（有代码有几种控制方法lqgLQRHOOPID）
clear;clc;
close all;
%% State-space model
A=[-0.38 0.60 -0.36 -9.80;
-0.98 -10.65 16.74 -0.21;
0.18 -5.39 -16.55 0;
0 0 1 0];
B=[-0.36;
-3.62;
-141.57;
0];
C4=[0 0 0 1]; %longitudinal speed u
D=0;
sys=ss(A,B,C4,D);
%% Define system with uncertainties
a11=ureal('a11',1,'Percentage',10);
a12=ureal('a12',1,'Percentage',10);
a13=ureal('a13',1,'Percentage',10);
a14=ureal('a14',1,'Percentage',10);
a21=ureal('a21',1,'Percentage',10);
a22=ureal('a22',1,'Percentage',5);
a23=ureal('a23',1,'Percentage',5);
a24=ureal('a24',1,'Percentage',10);
a31=ureal('a31',1,'Percentage',10);
a32=ureal('a32',1,'Percentage',5);
a33=ureal('a33',1,'Percentage',5);
a34=ureal('a34',1,'Percentage',10);
a41=ureal('a41',1,'Percentage',10);
a42=ureal('a42',1,'Percentage',10);
a43=ureal('a43',1,'Percentage',10);
a44=ureal('a44',1,'Percentage',10);
uA=[-0.38*a11 0.60*a12 -0.36*a13 -9.80*a14;
-0.98*a21 -10.65*a22 16.74*a23 -0.21*a24;
0.18*a31 -5.39*a32 -16.55*a33 0;
0 0 1*a43 0];
b1=ureal('b1',1,'Percentage',10);
b2=ureal('b2',1,'Percentage',5);
b3=ureal('b3',1,'Percentage',5);
b4=ureal('b4',1,'Percentage',10);
uB=[-0.36*b1;
-3.62*b2;
-141.57*b3;
0];
usys=ss(uA,uB,C4,D);
%% H-infinity loop shaping
s=tf('s'); %Define s as the Laplace variable
%Singular value diagram to obtain the crossover frequency (Wc)
figure(1)
sigma(sys) %Wc=6.11 rad/s (Value seen on figure 1)
Wc=6.11;
G=Wc/(s); %desired loopshape
% Compute the optimal loop shaping controller K
[K,CL,GAM] = loopsyn(sys,G);
[K.A,K.B,K.C,K.D]=ssdata(K); %得出控制器状态空间表达式参数
%% Closed-loop system
L=K*sys; %Adding disturbance to the system (这句注释的不对)
L_close=feedback(L,1); %Closed-loop
figure(2)
step(L_close);
str=sprintf('Step response of the pitch angle \\theta (H-infinity loopshaping)');
title(str);
ylabel('\theta (rad)');
%System information
S_final = stepinfo(L_close)
[y,t]=step(100*L_close);
sserror_final=100-y(end);
%Display gain and phase margins
figure(3)
margin(L_close);
%% Controller simplification
%Hankel singular values of K
figure(4)
hsv = hankelsv(K);
semilogy(hsv,'*--'), grid
title('Hankel singular values of K (\theta)'), xlabel('Order')
%Reduce controller from 7 to 5 states
Kr = reduce(K,5);
order(Kr) %check simplification
%Closed-loop responses comparison (K and Kr)
Lr=Kr*sys;
Lr_close=feedback(Lr,1);
figure(5)
step(L_close,'b',Lr_close,'r-.');
str=sprintf('Original and simplified controller (K & Kr) comparison for \\theta');
title(str);
ylabel('\theta (rad)');
legend('K','Kr','Location','Southeast')
%% Plot step response with uncertainties
uLr=Kr*usys;
uLr_close=feedback(uLr,1); %Closed loop system with uncertainties
figure(6);
step(uLr_close);
str=sprintf('Step response of the pitch angle \\theta (H-infinity loopshaping, uncertain plant)');
title(str);
ylabel('\theta (rad)')
View Code
H2_optimal_control_fixed_wing-----基于Ｈ２最优控制的小型无人机飞行姿态控制器设计
clc;clear;close all;
%{
%% State-space system
A=[-0.38 0.60 -0.36 -9.80;
-0.98 -10.65 16.74 -0.21;
0.18 -5.39 -16.55 0;
0 0 1 0];
B2=[-0.36; -3.62; -141.57; 0];
B1 =[0 0.6 0 0.2]'; %干扰矩阵
C1 =[0.3 0 0 0; %这个C1矩阵怎样选择
0 0.4 0 0;
0 0 0.1 0;
0 0 0 0.02]
C2= [0 0 0 1];
D11=zeros(4,1);
D12=[0 0 0 1]';
D21=[1];
D22 =[0];
G =ss(A,B2,C2,D22);
P =ltisys(A,[B1 B2],[C1;C2],[D11 D12;D21 D22])
[A1,B1,C1,D1]=branch(P)
nmeas =1; %测量输出个数
ncont =1; %控制输入个数
gmin =0.1;
gmax=18;
tol =0.0001
% [K,CL,gamma] = hinfsyn(P,nmeas,ncont) %这里为什么不能使用这个函数会出错？
[K,CL,gamma] = hinfsyn(P,nmeas,ncont,gmin,gmax,tol);
%% 画出阶跃响应信号
figure(1)
clsys =slft(P,K); %得到闭环系统
spol(clsys); %该函数返回闭环系统的极点
splot(clsys,'st',1);%画出阶跃响应曲线
% %% 画出控制信号
% figure(2)
% clsys =slft(K,P); %得到闭环系统
% spol(clsys); %该函数返回闭环系统的极点
% splot(clsys,'st');%画出阶跃响应曲线
[Ak,Bk,Ck,Dk]=branch(K)
sys=ss(Ak,Bk,Ck,Dk)
K_tf=zpk(ss(Ak,Bk,Ck,Dk))
[Knum2,Kden2]=ss2tf(Ak,Bk,Ck,Dk)
K_tf1=tf(Knum2,Kden2)
figure(3)
step(feedback(G*K_tf,-1),40,'-')
%%
G =ss(A,B2,C2,D22);
s=tf('s')
W1=(100)/(10*s+1);
W2=1e-6;
W3=(10*s)/(s+0.000008);
T_ss=augtf(G,W1,W2,W3);
[Aa,Bb,Cc,Dd]=branch(T_ss);
[Aa,Bb,Cc,Dd]=ssdata(T_ss);
K_Hinf=hinf(T_ss) %获得H00优控制器
[AK1,BK1,CK1,DK1]=branch(K_Hinf);
[a,b1,b2,c1,c2,d11,d12,d21,d22]=branch(T_ss)
figure(4)
step(feedback(G,K_Hinf,-1),40,'-')
%}
%% 例子二： 系统参数数据-来自论文：基于Ｈ２最优控制的小型无人机飞行姿态控制器设计
A =[-0.136 11.52 -9.8 0;
-0.0551 -0.2059 0 0.98;
0 0 0 1;
0.026 -0.66 0 -0.5534];
B =[0.34 0.001;
-0.00167 -0.0063;
0 0;
-0.000304 -0.00985];
C =[ 0 0 0 1];
% C=eye(4);
D=[0 0];
% D =[0 0 0 0;0 0 0 0]';
G =ss(A,B,C,D);
num=[0.008 5]; den =[1 0.0005]; w1 =tf(num,den);
p1=w1;
W1s=w1*eye(4)
W2s=1e-5*eye(2);
W2s=tf(W2s)
p2=W2s;
num1=[8 0.1]; den1 =[0.00001 40]; w2 =tf(num1,den1);
p3=w2;
% num2=[1 1.667e-5]; den2 =[0.00001 40]; w3 =tf(num2,den2);
% W3s=[0 0 0 0;
% 0 w2 0 0;
% 0 0 0 0;
% 0 0 0 w3];
% T_ss=augtf(G,W1s,W2s,W3s);
T_ss=augtf(G,p1,p2,p3);
[Aa,Bb,Cc,Dd]=branch(T_ss);
[Aa,Bb,Cc,Dd]=ssdata(T_ss)
[a,b1,b2,c1,c2,d11,d12,d21,d22]=branch(T_ss);%当使用augtf函数获取增广系统矩阵表达式后，
%可用branch函数获得各个增光系统矩阵的参数
P =rt2smat(T_ss) %变换成系统矩阵P
K_h2=h2lqg(T_ss) %用H2lqg函数获得H2最优控制器
nmeas = 1;
ncont = 2;
[K,CL,gamma] = h2syn(T_ss,nmeas,ncont) %用H2syn函数获得H2最优控制器
% [] = h2syn(T_ss,nmeas,ncont) %用H2syn函数获得H2最优控制器
[AK0,BK0,CK0,DK0]=branch(K_h2);
[AK0,BK0,CK0,DK0]=ssdata(K_h2);%两个函数均可实现这个功能
K_H2_tf=zpk(ss(AK0,BK0,CK0,DK0))
figure(1)
bode(G*K_H2_tf,'-');
figure(2)
% feedin=[2];
% feedout=[ 4];
% step(feedback(G*K_H2_tf,1,feedin,feedout),'r-')
figure(1)
step(feedback(G*K_H2_tf,1),'r-',6)
figure(2)
step(feedback(G*K,1),'r-',6)
%% 问题 采用H2syn函数时 求出gamma 是无穷代表啥意思
%% 系统进行离散化
Ts=0.1; %采样时间
Gd=c2d(G,Ts,'z') %采用零阶保持器的方法对系统离散化
T_ssd=c2d(T_ss,Ts,'z') %广义系统的离散化
P1 =rt2smat(T_ssd)
[Aad,Bbd,Ccd,Ddd]=branch(T_ss);
[ad,bd1,bd2,cd1,cd2,dd11,dd12,dd21,dd22]=branch(T_ssd);
% dd11=zeros(4,1);
% Pd =ltisys(ad,[bd1,bd2],[cd1,cd2],[dd11,dd12,dd21,dd22])
Pd =mksys(ad,bd1,bd2,cd1,cd2,dd11,dd12,dd21,dd22,'tss')
Kd_h2=dh2lqg(Pd) %获得离散的H2最优控制器，这里使用的是命令dh2lqg
[AK0d,BK0d,CK0d,DK0d]=branch(Kd_h2);
Kd_h2_tf=zpk(ss(AK0d,BK0d,CK0d,DK0d,Ts)) %将控制器转换成零极点白表达式
figure(3)
bode(Gd*Kd_h2_tf,'-'); %画出系统伯德图
figure(4)
step(feedback(Gd*Kd_h2_tf,1),'r-',6)
%}
%系统矩阵变换函数 -参见：https://wenku.baidu.com/view/d79bb3cb0722192e44
function P =rt2smat(A,B1,B2,C1,C2,D11,D12,D21,D22)
if nargin ==1,
[A,B1,B2,C1,C2,D11,D12,D21,D22]=branch(A);
end
n =length(A);
P=[A,B1,B2;C1,D11,D12;C2,D21,D22];
P(size(P,1)+1,size(P,2)+1)=-inf;
m=size(P,2);
P(1,m)=n;
end
View Code
H infinty 状态反馈-基于Riccati方程方法
%小型固定翼无人机俯仰控制——系统参数参见文章：BACHELOR’S FINAL PROJECT(好done)
% H infinty 状态反馈-基于Riccati方程方法
clear
close all;
clc
%% State-space system
A=[-0.38 0.60 -0.36 -9.80;
-0.98 -10.65 16.74 -0.21;
0.18 -5.39 -16.55 0;
0 0 1 0];
B2=[-0.36; -3.62; -141.57; 0];
B1 =[0 0.86 0 4.2]'; %干扰矩阵
C1 =[0.1 0 0 0; %这个C1矩阵怎样选择
0 0.04 0 0;
0 0 0.02 0;
0 0 0 0.01;
0 0 0 0]
C2= [0 0 0 1];% 测量矩阵
% D11 =[0];
D11 =[0;0;0;0;0];
D12 =[0;0;0;0;1];
D21=0;
D22=0;
P =ltisys(A,[B1 B2],[C1;C2],[D11 D12;D21 D22])
%% 判断系统是否满足基于Riccati方程的假设条件
sys =ss(A,B2,C2,D22);
P0 = pole(sys)
Cc =rank(ctrb(A,B2))
Ob =rank(obsv(A,C2))
S =D12'*[C1 D12]
%% 利用Riccati方程求解H infinity 在增益K, A'*X+X*A+X*(B1*B1'-B2*B2')*X+C'*C=0
% A'*X+X*A-X*[B1 B2]*[-I 0;0 I]*[B1';B2']*X+C1'*C1=0
B =[B1 B2];
R =[-1 0;0 1];
C =C1;
X =care(A,B,C'*C,R)
K =-B2'*X
root =eig(A-K) %判断闭环系统是否稳定；看其特征值是否都有负实部
[Aa,Bb,Cc,Dd]=branch(P)
sys=ss(Aa,Bb,Cc,Dd)
P_tf=zpk(ss(Aa,Bb,Cc,Dd))
%
feedin=[1];
feedout=[1 2 3 4];
step(feedback(sys,K,feedin,feedout,-1));
% clsys1 =slft(P,K); %得到闭环系统
% spol(clsys1,'st') %验证闭环系统极点
View Code
|
【本文地址】
今日新闻 |
推荐新闻 |