预设性能控制简单算例仿真——双弹簧阻尼（附matlab代码）

突然发现研一的时候写的预设性能控制的博客一直有人看，估计是PPC又火起来了，虽然现在不做这个方向了，但想了想把之前做的仿真也加上，这样理论和仿真都有了更完整一点，可以帮助大家更快的入门PPC。

本仿真复现裕魏才盛学长博士论文第二章《一类非线性系统低复杂度静态预设性能控制方法》

有兴趣的同学也可以去学习一下这篇博士论文

---

1. 问题描述

(1) 动力学模型

大型非线性系统由N 个互联的子系统组成，具体动力学模型如下所示：
$$\{\begin{array}{l}{\dot{x}}_{i,j}=x_{i,j+1},j=1,...,n_i-1\\ {\dot{x}}_{i,n_i}=f_i\left(x_i,u_i\right)+h_i\left(x_1,...,x_N\right)+d_i\left(t\right)\\ y_i=x_{i,1}\end{array}$$
${x_i=[\begin{array}{c}{x}_{i,1},...,x_{i,n_i}\end{array}]}^T\in {\mathbb{R}}^n,y_i\in \mathbb{R}\left(i=1,...,N\right)$表示第$i$个子系统的状态向量和输出。

(2) 控制目标

a. 在不依赖于神经网络和模糊系统逼近的情况下，设计一个鲁棒的分散式控制器，使得每个子系统能够追踪上期望指令。同时追踪误差系统（2-10）中的所有闭环变量都是一致最终有界的；

b. 追踪误差系统的瞬态与稳态性能能够先验地设计并在全时间域上得以保证。

2. 低复杂度静态预设性能控制器设计

根据预设性能控制的基本原理，低复杂度预设性能控制器设计包含三个步骤。

首先，设计预设性能函数；其次，基于设计的预设性能函数，设计预设性能控制器；最后，在设计的控制器下，给出受控系统的稳定性证明。

(1) 预设函数设计

a. 滑模面设计

$$s_i=c_{i,1}{e}_{i,1}+c_{i,2}{e}_{i,2}+|\cdots +c_{i,n_i-1}{e}_{i,n_i-1}+e_{i,n_i}$$

其中，${c_i=[\begin{array}{c}{c}_{i,1},...,c_{i,n_i-1}\end{array}]}^T\in {\mathbb{R}}^{n_i}$ 是待设计的参数向量并且使得多项式$C_i\left(\lambda \right)=c_{i,1}+c_{i,2}\lambda +\cdots +c_{i,n,-1}{\lambda }^{n_i-1}+{\lambda }^{n_i}$ 是 Hurwitz 稳定的（$\lambda$ 是Laplace 算子）。同时，参
数$c_i$的选择应该保证当多项式$C_i(\lambda )=0$时有$n_i$个不同的解。值得注意的是，$s_i$的大小可以看作当前追踪误差到时滑模面的距离。基于上述滑模面可得：
$$\frac{e_{i,1}\left(\lambda \right)}{s_i\left(\lambda \right)}=\frac{1}{c_{i,1}+c_{i,2}\lambda +\cdots +c_{i,n_i-1}{\lambda }^{n_i-1}+{\lambda }^{n_i}}=\frac{1}{{\Pi }_{j^{\ast }=1}^{n_i}\left(\lambda +p_{i,j^{\ast }}\right)}$$
其中$:p_{i,j},>0$是多项式$C_i(\lambda )=0$时的第$j^{\ast }$个解$(i=1,...,N,j^{\ast }=1,...,n_t)$。

b. 预设函数设计

定义性能包络：
$$-{\delta }_{i,1}{\mu }_i(t)<s_i(t)<{\delta }_{i,2}{\mu }_i(t)$$
其中：${\delta }_{i,1},{\delta }_{i,2}$是正的设计参数。${\mu }_i(t)>0$ 是严格连续递减性能函数。不失一般性，性能函数选择为
$${\mu }_i(t)=\left({\mu }_{i0}-{\mu }_{i\infty }\right){\mathrm{e}}^{-{\kappa }_{i}t}+{\mu }_{i\sigma }$$
其中 ，${\mu }_{i0}>{\mu }_{i\infty }>0,{\kappa }_i>0$ 是待设计参数。

为了减轻控制系统设计的复杂度，定义如下约束无关连续转换函数：
$$s_i(t)\triangleq \rho (z_i){\mu }_i$$
其中：严格单调函数 $\rho (\bullet )$ 满足 $\rho (0)\ne 0,{\lim }_{z_i\to +\infty }\rho (z_i)={\delta }_{i,2},{\lim }_{z_i\to -\infty }\rho (z_i)=-{\delta }_{i,1}$。不失一般性，函数$\rho (\bullet )$选为：

$$\rho (z_i)=\frac{{\delta }_{i,2}\exp (z_i)-{\delta }_{i,1}\exp (-z_i)}{\exp (z_i)+\exp (-z_i)}$$
其中：$z_i\in \mathbb{R}$ 是新定义的误差转化变量。基于式 上式, 可得转换误差变量等于

$$z_i={\rho }^{-1}\left(s_i\right)=\frac{1}{2}\ln \left[\left({\delta }_{i,1}+{\Lambda }_i\right)/\left({\delta }_{i,2}-{\Lambda }_i\right)\right]\text{(}{\Lambda }_i=s_i/{\mu }_i)$$
将$z_i$对时间求导数可得：
$$\{\begin{array}{ll}{\dot{z}}_i& =\frac{1}{2{\mu }_i}\frac{{\delta }_{i,1}+{\delta }_{i,2}}{\left({\delta }_{i,1}+{\Lambda }_i\right)\left({\delta }_{i,2}-{\Lambda }_i\right)}({\dot{s}}_i-{\Lambda }_i{\dot{\mu }}_i)\\ & ={\xi }_i[f_{u_i,\theta }{u}_i+f_i(x_i,u_i^{\ast })-{\int }_{u_i,\theta }{u}_i^{\ast }+h_i(x_1,...,x_N)\\ & +d_i\left(t\right)-y_{di}^{(n_i)}+{\sum }_{j=1}^{n_i-1}{c}_{i,j}{e}_{i,j+1}-{\Lambda }_i{\dot{\mu }}_i]\\ {\xi }_i& =\frac{1}{2{\mu }_i}\frac{{\delta }_{i,1}+{\delta }_{i,2}}{\left({\delta }_{i,1}+{\Lambda }_i\right)({\delta }_{i,2}-{\Lambda }_i)}\end{array}$$
简写为${\dot{z}}_i={\xi }_i\left(f_{u_{i,\theta }}{u}_i+g_i\right)$。

(2) 预设性能控制器设计

基于$\rho (\bullet )$设计的无约束转换函数，设计的低复杂度预设性能控制器如下所示：
$$u_i=-k_i{\xi }_i{z}_i=-k_i{\xi }_i\ln \left(\frac{{\delta }_{i,1}+{\Lambda }_i}{{\delta }_{i,2}-{\Lambda }_i}\right)$$
其中：$k_i$是正的控制增益，其他参数含义见上文。

给出的控制器形式类似传统PID 控制中的静态比例（P 型）控制器，形式简单，复杂度低，容易在线实现。

(3) 稳定性分析

略

3.算例仿真验证

(1) 仿真算例

(2) 仿真结果

4. 源代码

clc;
clear all;
close all;

M1 = 0.25; M2 = 0.2;          %物块质量
k1 = 0.1; k2 = 0.1;           
K1 = 1; K2 = 1;
deta11 = 0.5; deta12 = 1;
deta21 = 1; deta22 = 0.7;
mu10 = 2; mu11 = 0.1;
mu20 = 3; mu21 = 0.1;
circle = 3000;

a21 = 10; a22 = 10;
c11 = 1; c21 = 1;
X = zeros(4,circle);
X(:,1) = [0.5; 0; -1; 0];

h = 0.01;
e = zeros(4,circle);

x11 = 0.5;
x12 = 0;
x21 = -1;
x22 = 0;

for i = 1:circle
    
    yd1(i) = 0.5*(sin(1.5*h*i)+sin(0.5*h*i));
    yd2(i) = sin(h*i);
    dyd1(i) = 0.5*(1.5*cos(1.5*h*i)+0.5*cos(0.5*h*i));
    dyd2(i) = cos(h*i);

    delta1 = 0.2*sin(3*h*i)*exp(-0.2*h*i);
    delta2 = 0.2*cos(3*h*i)*exp(-0.1*h*i);

    e(:,i) = [x11; x12; x21; x22]-[yd1(i);dyd1(i);yd2(i);dyd2(i)];


    s1(i) = c11*e(1,i)+e(2,i);       %滑模面
    s2(i) = c21*e(3,i)+e(4,i);

    mu1(i) = (mu10-mu11)*exp(-K1*h*i)+mu11;
    mu2(i) = (mu20-mu21)*exp(-K2*h*i)+mu21;

    gama1(i) = s1(i)/mu1(i); gama2 = s2(i)/mu2(i);

    ksi1(i) = (1/(2*mu1(i)))*((deta11+deta12)/((deta11+gama1(i))*(deta12-gama1(i))));
    ksi2(i) = (1/(2*mu2(i)))*((deta21+deta22)/((deta21+gama2)*(deta22-gama2)));

    u1(i) = -k1*ksi1(i)*log((deta11+gama1(i))/(deta12-gama1(i)));%位控制输出
    u2(i) = -k2*ksi2(i)*log((deta21+gama2)/(deta22-gama2));

    sat_u1(i) = sign(u1(i))*min(abs(u1(i)),a21);%饱和输出
    sat_u2(i) = sign(u2(i))*min(abs(u2(i)),a22);%
    
    fs1 = X(1,i);
    fs2 = 2*(X(3,i)-X(1,i))+0.12*(X(3,i)-X(1,i))^3;
    fd1 = 2*X(2,i)+0.2*(X(2,i))^2;
    fd2 = 2.2*(X(4,i)-X(2,i))+0.15*(X(4,i)-X(2,i))^2;
    fc1 = 0.02*sign(X(2,i));
    fc2 = 0.02*sign(X(4,i)-X(2,i));
    
    
    dx11(i) = x12;
    dx12(i) = sat_u1(i)/M1 + (1/M1)*(-fs1-fd1+fs2+fd2-fc1+fc2+delta1);
    dx21(i) = x22;
    dx22(i) = sat_u2(i)/M2 + (1/M2)*(-fs2-fd2-fc2+delta2);
    
    dX = [dx11; dx12; dx21; dx22];
    
    X(:,i) = h*dX(:,i) + [x11; x12; x21; x22];
    
    x11 = X(1,i); 
    x12 = X(2,i);
    x21 = X(3,i);
    x22 = X(4,i);
    

end
x=linspace(0,circle,circle);
figure(1);
subplot(2,1,1);
plot(x*h,X(1,:),'r',x*h,yd1,'-.b','linewidth',1.2);grid;
ylabel('x_{11}(m)');xlabel('time(s)');
legend('x_{11}','y_{d1}');
subplot(2,1,2);
plot(x*h,X(3,:),'r',x*h,yd2,'-.b','linewidth',1.2);grid;
ylabel('x_{21}(m)');xlabel('time(s)');
legend('x_{21}','y_{d2}');
sgtitle('质量块位移响应图','FontSize',12);

figure(2)
subplot(2,1,1);
plot(x*h,sat_u1,'r','linewidth',1);grid;
ylabel('u_{1}(N)');xlabel('time(s)');
subplot(2,1,2);
plot(x*h,sat_u2,'b','linewidth',1);grid;
ylabel('u_{2}(N)');xlabel('time(s)');
sgtitle('质量块位控制输入图','FontSize',12);

figure(3)
plot(x*h,s1,'r',x*h,deta12*mu1,':k',x*h,-deta11*mu1,':k','linewidth',1);grid;
ylabel('s');xlabel('time(s)');
legend('s1','PPC');
sgtitle('滑模面s_1随时间响应图','FontSize',12);
figure(4)
plot(x*h,s2,'b',x*h,deta22*mu2,':k',x*h,-deta21*mu2,':k','linewidth',1);grid;
ylabel('s');xlabel('time(s)');
legend('s2','PPC');
sgtitle('滑模面s_2随时间响应图','FontSize',12);

figure(5);
subplot(2,1,1);
plot(x*h,X(2,:),'r',x*h,dyd1,'-.b','linewidth',1.2);grid;
ylabel('x_{12}(m/s)');xlabel('time(s)');
legend('x_{12}','dy_{d1}');
subplot(2,1,2);
plot(x*h,X(4,:),'r',x*h,dyd2,'-.b','linewidth',1.2);grid;
ylabel('x_{22}(m/s)');xlabel('time(s)');
legend('x_{22}','dy_{d2}');
sgtitle('质量块速度响应图','FontSize',12);

本文来源于魏才盛学长博士论文《航天器姿态预设性能控制方法研究》，证明过程省略，可能也省略了一些其他的重要部分；
如有疏漏以论文为主，有兴趣强烈推荐阅读原文。
航天器姿态预设性能控制方法研究 - 中国知网 (cnki.net)