等价输入干扰EID（四）：“Yalmip与LMI实战”

等价输入干扰EID（四）：“Yalmip与LMI实战”

reference：**“Improve Disturbance-Rejection Performance for an Equivalent-Input-Disturbance-Based Control System by Incorporating a ProportionalIntegral Observer” (pdf) (Wu 等, 2020, p. 1254)

A.镇稳目标

使EID系统镇稳的LMI推导详细请见上篇。

本篇论文中，使整个系统稳定的LMI条件为：

$$\begin{gathered} [\begin{array}{cccc}{\Psi }_{11}& W_{2}C+B{M}^{\mathrm{T}}{X}_2& {\Psi }_{13}& 0\\ \star & {\Psi }_{22}& {\Psi }_{23}& {\Psi }_{24}\\ \star & \star & {\Psi }_{33}& {\Psi }_{34}\\ \star & \star & \star & {\Psi }_{44}\end{array}]<0 \\ (1) \end{gathered}$$

$$\begin{gathered} \left[\begin{array}{cccc}{X}_1& 0& 0& 0\\ \star & X_2& X_{2}M& 0\\ \star & \star & X_3& 0\\ \star & \star & \star & X_4\end{array}\right]>0 \\ (2) \end{gathered}$$

待定LMI变量为$X_1,X_{11},X_{22},X_3,X_4,W_1,W_2,W_3$，其中前四个为对称半正定矩阵，后三个为合适矩阵即可

(1)式和(2)式内各参数为：

$$C=U[\begin{array}{cc}S& 0\end{array}]{\hat{T}}^{\mathrm{T}}$$

$$X_2=\hat{T}\left[\begin{array}{cc}{X}_{11}& 0\\ 0& X_{22}\end{array}\right]{\hat{T}}^{\mathrm{T}}$$

$$\begin{array}{rl}& {\Psi }_{11}=(A+{\alpha }_{1}I)X_1+X_1(A+{\alpha }_{1}I{)}^{\mathrm{T}}+B{W}_1+W_1^{\mathrm{T}}{B}^{\mathrm{T}}\\ & {\Psi }_{13}=W_{2}CM+B{X}_3\\ & {\Psi }_{22}=(A+{\alpha }_{2}I)X_2+X_2(A+{\alpha }_{2}I{)}^{\mathrm{T}}-W_{2}C-C^{\mathrm{T}}{W}_2^{\mathrm{T}}\\ & -B{M}^{\mathrm{T}}{X}_2-X_{2}M{B}^{\mathrm{T}}\\ & {\Psi }_{23}=-B{X}_3+C^{\mathrm{T}}{W}_3^{\mathrm{T}}+(A+{\alpha }_{2}I)X_{2}M-W_{2}CM\\ & {\Psi }_{24}=-B{C}_F{X}_4+(B_F{B}^{+}{W}_{2}C{)}^{\mathrm{T}}+X_{2}M{B}_F\\ & {\Psi }_{33}=W_{3}CM+(W_{3}CM{)}^{\mathrm{T}}\\ & {\Psi }_{34}=X_3{B}_F^{\mathrm{T}}+(B_F{B}^{+}{W}_{2}CM{)}^{\mathrm{T}}\\ & {\Psi }_{44}=(A_F+B_F{C}_F)X_4+X_4(A_F+B_F{C}_F{)}^{\mathrm{T}}\end{array}$$

求解后可得参数：

$$\{\begin{array}{rl}{K}_P& =W_1{X}_1^{-1}\\ L_1& =W_{2}US{X}_{11}^{-1}{S}^{-1}{U}^{\mathrm{T}}\\ L_2& =W_{3}US{X}_{11}^{-1}{S}^{-1}{U}^{\mathrm{T}}\end{array}$$

B.基于 Yalmip 的 Matlab 程序编写

1.初始化准备

%清空环境
clear;  clc;

%系统矩阵
A = [-2 3; 4 -5];
B = [1;1]; 
C = [2 1];
Bd = [2;1];
M = [-1;0];

%系统参数
AF = -101;  
BF = 100;
CF = 1; 

%计算伪逆BZ
BZ = pinv(B);

%svd分解
[U,S,V] = svd(C);
V1 = V(:,1);
V2 = V(:,2);
S = S(1); 

%权重，不晓得干嘛的
w = 0.5;

%其他参数
a=1;
b=1;
alfa_1 = 1;
alfa_2 = 1;

2.定义sdpvar变量

% 半正定矩阵变量X
X1 = sdpvar(2);
X11 = sdpvar(1);
X22 = sdpvar(1);
X3 = sdpvar(1);
X4 = sdpvar(1);

%X2表达式
%V1*X11*V1'+V2*X22*V2'
X2=V1*X11*V1'+V2*X22*V2';

% 合适矩阵变量W
W1 = sdpvar(1,2);
W2 = sdpvar(2,1);
W3 = sdpvar(1);

3.构建LMI不等式

%不等式1
pusai_11 = (A*X1+alfa_1*X1)+(A*X1+alfa_2*X1)'+B*W1+(B*W1)';
pusai_12 = W2*C+B*M'*X2;
pusai_13 = W2*C*M+B*X3;
pusai_14 = 0*ones(size(pusai_13));%适当维度的全零矩阵
pusai_22 = (A*X2+alfa_2*X2)+(A*X2+alfa_2*X2)'-W2*C-(W2*C)'-B*M'*X2-(B*M'*X2)';
pusai_23 = -B*X3+C'*W3'+A*X2*M+alfa_2*X2*M-W2*C*M;
pusai_24 = -B*CF*X4+(BF*BZ*W2*C)'+X2*M*BF;
pusai_33 = W3*C*M+(W3*C*M)';
pusai_34 = X3*BF'+(BF*BZ*W2*C*M)';
pusai_44 = (AF+BF*CF)*X4+((AF+BF*CF)*X4)';

eq1 = [pusai_11  pusai_12  pusai_13  pusai_14;
       pusai_12' pusai_22  pusai_23  pusai_24;
       pusai_13' pusai_23' pusai_33  pusai_34;
       pusai_14' pusai_24' pusai_34' pusai_44;
      ];

%不等式2
Z_22 = zeros(size(X1));          %2*2全零矩阵
Z_21 = [0;0];
eq2 = [X1    Z_22    Z_21  Z_21;
       Z_22' X2      X2*M  Z_21;
       Z_21' (X2*M)' X3    0;
       Z_21' Z_21'   0     X4;
      ]

注意上式的全零阵的维度，也就是观察整个矩阵保持方阵的维度正确

4.LMI约束编程

% LMI约束
F = [X1 >= 0];
F = F + [X11 >= 0];
F = F + [X22 >= 0];
F = F + [X3 >= 0];
F = F + [X4 >= 0];
F = F + [eq1 <= 0];
F = F + [eq2 >= 0];

5.求解LMI

% 求解LMI 
sol = optimize(F);
if sol.problem == 0
    disp("求解成功");
elseif sol.problem == 1
        disp( "求解失败，无可行解");
    else
        disp( "求解失败，其他错误");
end

6.提取解

% 提取解
X1 = value(X1);
X11 = value(X11);
X22 = value(X22);
X3 = value(X3);
X4 = value(X4);
W1 = value(W1);
W2 = value(W2);
W3 = value(W3);

这里我都提取了，其实只需要些许解就ok了，不过无伤大雅

7.计算反馈控制器

% 计算反馈控制器
L1 = W2*U*S*inv(X11)*inv(S)*U'
L2 = W3*U*S*inv(X11)*inv(S)*U'
L3 = BZ*L1
Kp = W1*inv(X1)
B1=BZ;

C.仿真实战

1.仿真模型（s-function搭建）

2.干扰信号图

3.论文中LMI解得参数抑制后的输出效果

4.本文yalmip工具解得参数抑制输出效果

结论：用另外一种LMI工具一样可以解出可行参数使系统达到镇定效果

NOTICE！！！

本文使用的仿真文件的参数都在设置的Parameters.m文件，写入之后运行到工作区才能进行仿真。

上传的仿真文件是一个星期之前的版本，要自行改动些许，仿真自己琢磨去，重点是LMI求解方式！！！