1 给出倒立摆系统简化模型[1]

Matlab+simulink H无穷控制器的设计(作业向)

控制器封装库(五)鲁棒控制器
Chenglin Li 的视频
 · 4821 播放

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——

参考

  1. ^https://zhuanlan.zhihu.com/p/113053634