1 给出倒立摆系统简化模型[1]
Matlab+simulink H无穷控制器的设计(作业向)
控制器封装库(五)鲁棒控制器
2 LMI矩阵不等式
- 可行性问题
- 特征值问题
3 绘制Simulink仿真图
4 在扰动情况下的求解结果
最优H∞控制决策变量
5 参数求解程序
(1)主程序
function [A,B1,B2,K]=Parameters()
clear,clc
M=1; %小车质量
m=0.1; %倒立摆质量
L=0.5; %摆的长度
g=9.8 ;%重力加速度
J=m*L^2/3; %转动惯量
k=J*(M+m)+m*M*L^2 ;
k1=(M+m)*m*g*L/k ;
k2=-m^2*g*L^2/k ;
k3=-m*L/k ;
k4=( J+m*L^2 ) /k ; %中间变量
%------PlantA------------
A=[0,0,1,0 ;0,0,0,1 ;k1,0,0,0 ;k2,0,0,0];
% B1=[0; 0; 0.09; 0.11 ];
B1=[0;0;0.2;0.2] ;
B2=[0; 0; k3; k4];
C1=[1,0,0,0; 0,1,0,0; 0,0,1,0; 0,0,0,1; 0,0,0,0] ;
% C=[1,0,0,0;0,0,0,0] ;
N=length(A);
D11=[0;0;0];
D12=[0; 0; 0; 0; 1]; %状态方程各个矩阵
% D12=0;
%--------直接求出K-------------
[X,W,gamma]=LMID(A,B1,B2,C1,D11,D12,N)
K=W*inv(X)
%--------------------------
% gamma=1.5;
%-----根据所得的次优gamma重新求解K--------
K=LMIC(A,B1,B2,C1,D11,D12,gamma,N)
% A=[0,1,0,0;10.78,0,0,0;0,0,0,1;-0.98,0,0,0];
% B1=[0;0;0.09;0.11];
% B2=[0;-1;0;0.5];
% N=length(A);
% D12=[0];
% C=[1,0,0,0];
% D11=[];
% gamma=4;
%-----LMI求解------------
% gamma=2;
% K=LMIC(A,B1,B2,C,D11,D12,gamma,N) %LMI求解
%-------------------
% K=Riccati(A,B1,B2,C1);
% K=[51.1721, 3.6803, 9.2785, 7.9564];
%------PlantB------------
% A=[0 1 0 0;
% 0 -0.0883 0.6293 0;
% 0 0 0 1;
% 0 -0.2357 27.8285 0] ;
% B1=[0 2.3566 0 104.2027]';
% B2=[0 0.8832 0 2.3566]';
% C=[0.064 0 0 0;
% 0 1e-3 0 0;
% 0 0 0.11 0;
% 0 0 0 0.01;
% 0 0 0 0];
% D12=[0;0;0;0;1];
%-----Raccati方程求解---可行--------
% B=[B1,B2];
% R=[-1, 0; 0, 1];
% X=care(A, B, C'*C, R) ;
% K=-B2'*X %控制器
% K=[ 112.1078 99.9032 -364.5474 -72.8227];
% X=[0.1253,0.0042,-0.6266,-0.0266;
% 0.0042,0.565,0.0156,-0.3558;
% -0.6266,0.0156,3.6427,-0.1747;
% -0.0266,-0.3558,-0.1747,0.5286];
% W=[0.4583,0.0057,-0.1069,-0.8101];
% K=W*inv(X)
end
(2)函数调用(特征值问题,求最小决策变量)
function [X,W,gamma]=LMID(A,B1,B2,C1,D11,D12,N)
%{
程序功能:
1、利用LMI求解倒立摆
2、状态反馈控制,求控制器
3、\gamam-次优控制
参考文献:[1]吕 申,武俊峰.基于 LMI 优化的鲁棒控制器设计[J].工业仪表与自动化装置,2017,3:123-128.
%}
% clear,clc
% A=[0,1,0,0; 0,-0.0618,-0.7167,0; 0,0,0,1; 0,0.2684,31.6926,0];
% B1=[0; -2.6838; 0; 118.6765];
% B2=[0 ;0.8906;0;-2.6838];
% C1=[1e-5,0,0,0; 0,-1e-4,0,0; 0,0,-0.01,0; 0,0,0,-0.01];
% D11=[0;0;0];
% D12=[0;0;0];% 参数初始化
% gamma=1;
%构建矩阵变量
setlmis([]);
X=lmivar(1,[N,1]); %4阶的对称满块
W=lmivar(2,[1,N]); %的矩阵
g2=lmivar(1,[1,1]);
%描述LMI
lmiterm([1,1,1,X],A,1,'s'); %AX+(AX)'
lmiterm([1,1,1,W],B2,1,'s'); %B2*W+(B2*W)'
lmiterm([1,1,2,0],B1); %B1
lmiterm([1,1,3,-W],1,D12');%W*D2'
lmiterm([1,1,3,X],1,C1');%X*C1'
lmiterm([1,2,2,0],-1); %-I
lmiterm([1,2,3,0],D11' );%D11'
lmiterm([1,3,3,g2], -1,1);%-gamma^2*I
lmiterm([-2,1,1,X],1,1); %X正定
lmisys=getlmis ;%获取lmi信息
%-------------------------------
% lmiterm([1,1,1,X],A,1,'s'); %AX+(AX)'
% lmiterm([1,1,1,W],B2,1,'s'); %B2*W+(B2*W)'
% lmiterm([1,1,1,0],gamma^(-2)*(B1*B1') );
% lmiterm([1,2,1,X],C1,1 ); %C1*X
% lmiterm([1,2,1,W],D12,1 ); %D12*W
% lmiterm([1,2,2,0], -1 ); %-I
% lmisys=getlmis ;%获取lmi信息
%求解可行解X,W
% [tmin, xfeas]=feasp(lmisys);
% P=ltisys( A,[B1,B2],[C1;C2],[D11,D12;D21,D22] );
% gopt=ltiflmi(P,[1,1]);
% c=mat2dec(lmisys, eye(3) ); %目标值:最优X的trace
n=decnbr(lmisys);
c=zeros(n,1);
for j=1:n
bj=defcx(lmisys, j, g2);
c(j)=bj;
end
% options=[1e-5, 0, 0, 0, 0]; %计算精度
[copt,xopt]=mincx(lmisys,c );
% if(tmin<0)
% sys=[];
% disp('Feasible');
% else
% return
%
% end
X=dec2mat(lmisys, xopt, X);
W=dec2mat(lmisys, xopt, W);
gamma=sqrt(copt);
% W=dec2mat(lmisys, xfeas, W);
%
% %获取反馈控制器
% K=W*inv(X);
% K=W/X;
% sys=P;
end
(3)可行性问题(求解给定决策变量下的可行性矩阵)
function sys=LMIC(A,B1,B2,C1,D11,D12,gamma,N)
%{
程序功能:
1、利用LMI求解倒立摆
2、状态反馈控制,求控制器
3、
参考文献:[1]吕 申,武俊峰.基于 LMI 优化的鲁棒控制器设计[J].工业仪表与自动化装置,2017,3:123-128.
%}
% clear,clc
% A=[0,1,0,0; 0,-0.0618,-0.7167,0; 0,0,0,1; 0,0.2684,31.6926,0];
% B1=[0; -2.6838; 0; 118.6765];
% B2=[0 ;0.8906;0;-2.6838];
% C1=[1e-5,0,0,0; 0,-1e-4,0,0; 0,0,-0.01,0; 0,0,0,-0.01];
% D11=[0;0;0];
% D12=[0;0;0];% 参数初始化
% gamma=1;
%构建矩阵变量
setlmis([]);
[X,n,Sx]=lmivar(1,[N,1]); %4阶的对称满块
[W,n,Sw]=lmivar(2,[1,N]); %4*1的矩阵
%描述LMI
%{
lmiterm([1,1,1,X],A,1,'s'); %AX+(AX)'
lmiterm([1,1,1,W],B2,1,'s'); %B2*W+(B2*W)'
lmiterm([1,1,2,0],B1); %B1
lmiterm([1,1,3,-W],1,D12');%W*D2'
lmiterm([1,1,3,X],1,C1');%X*C1'
lmiterm([1,2,2,0],-1); %-I
lmiterm([1,2,3,0],D11' );%D11'
lmiterm([1,3,3,0], -gamma^2);%-gamma^2*I
%}
%-------------------------------
lmiterm([1,1,1,X],A,1,'s'); %AX+(AX)'
lmiterm([1,1,1,W],B2,1,'s'); %B2*W+(B2*W)'
lmiterm([1,1,2,0], B1); %C1*X
lmiterm([1,1,3,X],1,C1' ); %D12*W
lmiterm([1,1,3,-W],1,D12' ); %W'*D12'
lmiterm([1,2,2,0],-gamma^2);%-gamma
lmiterm([1,2,3,0],D11'); %D12'
lmiterm([1,3,3,0], -1 ); %-I
%------------------------------------------
lmiterm([-2,1,1,X],1,1); %X正定
%------------------------------------------
lmisys=getlmis ;%获取lmi信息
%求解可行解X,W
[tmin, xfeas]=feasp(lmisys);
if(tmin<0)
disp('Feasible');
else
sys=[];
return
end
X=dec2mat(lmisys, xfeas, X);
W=dec2mat(lmisys, xfeas, W);
%获取反馈控制器
% K=W*inv(X);
K=W/X;
sys=K;
end
——2020.03.14——
评论(0)
您还未登录,请登录后发表或查看评论