自适应控制——仿真实验一 用李雅普诺夫稳定性理论设计自适应规律

一、问题描述

设控制对象的状态方程为
$${\dot{\mathbf{x}}}_p={\mathbf{A}}_p(t)x_p+{\mathbf{b}}_p(t)u\text{\hspace{1em}(1)}$$
式中
$${\mathbf{A}}_p=\left[\begin{array}{cc}0& 1\\ -6& -7\end{array}\right],{\mathbf{b}}_p=\left[\begin{array}{l}2\\ 4\end{array}\right]\text{\hspace{1em}(2)}$$
参考模型的状态方程为
$${\dot{\mathbf{x}}}_m={\mathbf{A}}_m{x}_m+{\mathbf{b}}_{m}r\text{\hspace{1em}(3)}$$
式中
$${\mathbf{A}}_m=\left[\begin{array}{cc}0& 1\\ -10& -5\end{array}\right],{\mathbf{b}}_m=\left[\begin{array}{l}1\\ 2\end{array}\right]\text{\hspace{1em}(4)}$$
用李雅普诺夫稳定性理论设计自适应规律。

二、问题建模

由于控制对象的参数（状态矩阵 ${\mathbf{A}}_p$ 和控制矩阵 ${\mathbf{b}}_p$ ）一般是未知的，且无法直接调整。所以为改变控制对象的动态特性，需采用前馈控制加反馈控制。

控制信号 $u$ 由前馈信号 $Kr$ 和反馈信号 $F{x}_p$ 组成，即
$$u=Kr+F{\mathbf{x}}_p\text{\hspace{1em}(5)}$$
式中，$r$ 为 $m$ 维输入向量，${\mathbf{x}}_p$ 为 $n$ 维状态向量，$K$ 为 $m\times m$ 前馈增益矩阵，$F$ 为 $m\times n$ 反馈增益矩阵；具体在本次仿真实验中，输入向量维度 $m=1$，状态向量维度 $n=2$。

将(5)式代入控制对象的状态方程，可得
$${\dot{\mathbf{x}}}_p=\left[{\mathbf{A}}_p(t)+{\mathbf{b}}_p(t)F\right]{\mathbf{x}}_p+{\mathbf{b}}_p(t)Kr\text{\hspace{1em}(6)}$$
设系统的广义状态误差向量为
$$\mathbf{e}={\mathbf{x}}_m-{\mathbf{x}}_p\text{\hspace{1em}(7)}$$
由参考模型的状态方程，结合(6)式及(7)式，可得：
$$\dot{\mathbf{e}}={\mathbf{A}}_m\mathbf{e}+\left({\mathbf{A}}_m-{\mathbf{A}}_p-{\mathbf{b}}_{p}F\right){\mathbf{x}}_p+\left({\mathbf{b}}_m-{\mathbf{b}}_{p}K\right)r\text{\hspace{1em}(8)}$$
在理想情况，即 $e\to 0$ 的情况下，(8)式等号右端后两项应等于0。设前馈增益矩阵 $K$ 和反馈增益矩阵 $F$ 的理想值分别为 $\bar{K}$ 和 $\bar{F}$。

则最终可将(8)式写成
$$\dot{\mathbf{e}}={\mathbf{A}}_m\mathbf{e}+{\mathbf{b}}_m{\bar{K}}^{-1}\Phi {\mathbf{x}}_p+{\mathbf{b}}_m{\bar{K}}^{-1}\Psi r\text{\hspace{1em}(9)}$$
式中，$\Phi =\bar{F}-F$ 为 $m\times n$ 矩阵，$\Psi =\bar{K}-K$ 为 $m\times m$ 矩阵。

选取李雅普诺夫函数为：
$$V=\frac{1}{2}\left[{\mathbf{e}}^T\mathbf{P}\mathbf{e}+\mathrm{tr}\left({\Phi }^T{\Gamma }_1^{-1}\Phi +{\Psi }^T{\Gamma }_2^{-1}\Psi \right)\right]\text{\hspace{1em}(10)}$$
式中，$\mathbf{P}$ 为 $n\times n$ 维正定对称阵，${\Gamma }_1$ 和 ${\Gamma }_2$ 均为 $m\times m$ 维正定对称阵；符号 $\mathrm{tr}$ 表示矩阵的迹。

求(10)式对时间的导数，得
$$\dot{V}=\frac{1}{2}\left[\dot{\mathbf{e}}\mathbf{P}\mathbf{e}+{\mathbf{e}}^T\mathbf{P}\dot{\mathbf{e}}+\mathrm{tr}\left({\dot{\Phi }}^T{\Gamma }_1^{-1}\Phi +{\Phi }^T{\Gamma }_1^{-1}\dot{\Phi }+{\dot{\Psi }}^T{\Gamma }_2^{-1}\Psi +{\Psi }^T{\Gamma }_2^{-1}\dot{\Psi }\right)\right]\text{\hspace{1em}(11)}$$
将(9)式代入(11)式，再根据矩阵迹的性质，于是有
$$\begin{array}{rl}\dot{V}=& \frac{1}{2}{\mathbf{e}}^T\left(\mathbf{P}{\mathbf{A}}_m+{\mathbf{A}}_m^{\mathbf{T}}\mathbf{P}\right)\mathbf{e}+\mathrm{tr}\left({\dot{\Phi }}^T{\Gamma }_1^{-1}\Phi +{\mathbf{x}}_p{\mathbf{e}}^T\mathbf{P}{\mathbf{b}}_m{\bar{K}}^{-1}\Phi \right)\\ & +\mathrm{tr}\left({\dot{\Psi }}^T{\Gamma }_2^{-1}\Psi +r{\mathbf{e}}^T\mathbf{P}{\mathbf{b}}_m{\bar{K}}^{-1}\Psi \right)\end{array}\text{\hspace{1em}(12)}$$
为满足李雅普诺夫第二法，需保证(12)式是负定的，对应的情况为(12)式第一项是负定的，后两项都为零。

因为 ${\mathbf{A}}_m$ 为稳定矩阵，则可选定正定对称阵 $Q$，使 $\mathbf{P}{\mathbf{A}}_m+{\mathbf{A}}_m^{\mathbf{T}}\mathbf{P}=-\mathbf{Q}$ 成立。同时根据上述对应情况，$\Phi$ 和 $\Psi$ 的选择如下：
$$\begin{array}{rl}\dot{\Phi }& =-{\Gamma }_1{\left({\mathbf{b}}_m{\bar{K}}^{-1}\right)}^T\mathbf{P}\mathbf{e}{\mathbf{x}}_p^T\\ \dot{\Psi }& =-{\Gamma }_2{\left({\mathbf{b}}_m{\bar{K}}^{-1}\right)}^T\mathbf{P}\mathbf{e}{r}^T\end{array}\text{\hspace{1em}(13)}$$
当 ${\mathbf{A}}_p$ 和 ${\mathbf{b}}_p$ 为常值或缓慢变化时，可得自适应调节规律：
$$\begin{array}{rl}F(t)& ={\int }_0^t{\Gamma }_1{\left({\mathbf{b}}_m{\bar{K}}^{-1}\right)}^T\mathbf{Pe}{\mathbf{x}}_p^{T}d\tau +F(0)\\ K(t)& ={\int }_0^t{\Gamma }_2{\left({\mathbf{b}}_m{\bar{K}}^{-1}\right)}^T\mathbf{Pe}rd\tau +K(0)\end{array}\text{\hspace{1em}(14)}$$
需额外说明的一点是，按上述步骤推导得到的自适应调节规律要求 ${\mathbf{x}}_p$ 与 $r$ 线性独立。两者独立的条件是 $r(t)$ 为具有一定频率的方波信号或为 $q$ 个不同频率的正弦信号组成的分段连续信号，其中 $q>n/2$ 或 $q>(n-1)/2$。

三、问题求解

由上述推导可知，为采取李雅普诺夫稳定性理论设计该MRACS，需引入前馈增益矩阵 $K$ 和反馈增益矩阵 $F$，设计的目标是确定 $K$ 和 $F$ 的系数。

在引入两个增益矩阵进行自适应控制后，可调系统的状态方程变为：
$${\dot{\mathbf{x}}}_p=\left[{\mathbf{A}}_p(t)+{\mathbf{b}}_p(t)F\right]{\mathbf{x}}_p+{\mathbf{b}}_p(t)Kr\text{\hspace{1em}(15)}$$
由之前的推导可知，(14)式中的 ${\mathbf{b}}_m{\bar{K}}^{-1}$ 与 ${\mathbf{b}}_p$ 的关系如下：
$${\mathbf{b}}_m{\bar{K}}^{-1}={\mathbf{b}}_p=\left[\begin{array}{l}2\\ 4\end{array}\right]\text{\hspace{1em}(16)}$$
选取(14)式中的部分自适应参数如下：
$$\mathbf{P}=\left[\begin{array}{ll}3& 1\\ 1& 1\end{array}\right],{\Gamma }_1={\Gamma }_2=1\text{\hspace{1em}(17)}$$
所以可得最终的自适应规律：
$$\begin{array}{rl}F(t)& ={\int }_0^t\left[\begin{array}{ll}2& 4\end{array}\right]\left[\begin{array}{ll}3& 1\\ 1& 1\end{array}\right]\mathbf{e}{\mathbf{x}}_p^{T}d\tau +F(0)\\ K(t)& ={\int }_0^t\left[\begin{array}{ll}2& 4\end{array}\right]\left[\begin{array}{ll}3& 1\\ 1& 1\end{array}\right]\mathbf{erd}\tau +K(0)\end{array}\text{\hspace{1em}(18)}$$
下将上述连续自适应规律进行离散化，用于实际的数值仿真实验。设数值积分步长为 $h$，则各时刻的参考模型状态向量及控制对象状态向量如下：
$$\begin{array}{rl}{\mathbf{x}}_m(k+1)& ={\mathbf{x}}_m(k)+h\left[{\mathbf{A}}_m(k){\mathbf{x}}_m(k)+{\mathbf{B}}_m(k)r(k)\right]\\ {\mathbf{x}}_p(k+1)& ={\mathbf{x}}_p(k)+h\left[{\mathbf{A}}_p(k){\mathbf{x}}_p(k)+{\mathbf{B}}_p(k)u(k)\right]\end{array}\text{\hspace{1em}(19)}$$
由于上述推导得到的自适应控制规律要求 ${\mathbf{x}}_p$ 与 $r$ 线性独立，即要求 $r(t)$ 为具有一定频率的方波信号或为 $q$ 个不同频率的正弦信号组成的分段连续信号，其中$q>n/2$ 或 $q>(n-1)/2$。在本次实验中，$n=2$，对应就要求 $q>1$，所以本次实验中选取由3个不同频率的正弦信号组成的分段连续信号，具体的输入信号的形式如下：
$$r(k)=\sin (0.01\pi k)+4\sin (0.2\pi k)+\sin (\pi k)\text{\hspace{1em}(20)}$$
我们设计自适应规律时引入的控制信号 $u$ 的离散化形式如下：
$$u(k)=K(k)r(k)+F(k){\mathbf{x}}_p(k)\text{\hspace{1em}(21)}$$
最终，还需将自适应规律离散化：
$$\begin{array}{rl}F(k)& =h\cdot \sum _{j=0}^k{\mathbf{b}}_p^T\mathbf{P}\mathbf{e}(k){\left({\mathbf{x}}_p(k)\right)}^T+F(0)\\ K(k)& =h\cdot \sum _{j=0}^k{\mathbf{b}}_p^T\mathbf{P}\mathbf{e}(k)r(k)+K(0)\end{array}\text{\hspace{1em}(22)}$$
在推导出全部的自适应规律并对相应规律进行离散化后，通过MATLAB进行了相关的仿真实验。

可以得到2个维度的状态向量的参考模型值与可调系统值的情况如下：

图1. 状态向量的参考模型值与可调系统值

可以看到，可调系统并没有很好的跟踪参考模型，这是由于在该例中不存在最优匹配。

附录：实现MATLAB代码

% 课本习题3.4-用李雅普诺夫稳定性理论设计自适应规律
clear, clc;
close all;

h=0.01;L=100/h;     % 数值积分步长和仿真步数
% 可调系统的系数矩阵
Ap = [0 1;-6 -7];
Bp = [2; 4];
% 参考模型的系数矩阵
Am = [0 1;-10 -5];
Bm = [1; 2];
% n为行向量维数、m为列向量维数，Bp是n*m的矩阵
n = size(Bp, 1);
m = size(Bp, 2);

P = [3 1;1 1];              % 经计算得到的用于自适应规律的正定对称矩阵

% 设定所有参数的初始值
yr0 = zeros(m, 1);
xp0 = zeros(n, 1);
xm0 = zeros(n, 1);
u0 = zeros(m, 1);
e0 = zeros(n, 1);
F0 = zeros(m, n);           % 反馈增益矩阵初始值
K0 = zeros(m, m);           % 前馈增益矩阵初始值

% 初始分配参数空间
time = zeros(1, L);         % 用于记录仿真的时刻，对应绘图的横轴
yr = zeros(m, L);           % 输入信号(L个m维向量)
xp = zeros(n, L);           % 可调系统的状态向量(L个n维向量)
xm = zeros(n, L);           % 参考模型的状态向量(L个n维向量)
u = zeros(m, L);            % 控制信号(L个m维向量)
e = zeros(n, L);            % 系统的广义状态误差向量(L个n维向量)

for k = 1:L
    time(k) = k*h;
    % 输入信号
    yr(k) = 1*sin(0.01*pi*time(k))+4*sin(0.2*pi*time(k))+sin(1*pi*time(k));
    xp(:,k) = xp0+h*(Ap*xp0+Bp*u0);     % 计算xp
    xm(:,k) = xm0+h*(Am*xm0+Bm*yr0);    % 计算xm
    e(:,k) = xm(:,k)-xp(:,k);           % e=xm-xp
    
    % 代入F和K的自适应控制规律
    F = F0+h*(Bp'*P*e0*xp0');
    K = K0+h*(Bp'*P*e0*yr0);

    % 控制信号u=K*r+F*xp（K是前馈增益矩阵，F是反馈增益矩阵）
    u(:,k) = K*yr(k)+F*xp(:,k);
    
    % 将本轮求解得到的参数赋值给参数初始值，方便下一轮迭代使用
    yr0 = yr(:,k);
    u0 = u(:,k);
    e0 = e(:,k);
    xp0 = xp(:,k);
    xm0 = xm(:,k);
    F0 = F;
    K0 = K;
end

subplot(2,1,1);
plot(time, xm(1,:), 'Color', 'b', 'LineWidth', 0.9);
hold on
plot(time, xp(1,:), 'Color', 'r', 'LineStyle', '--', 'LineWidth', 1.1);
xlabel('t');
ylabel('x_m_1(t)、x_p_1(t)');
legend('x_m_1(t)','x_p_1(t)');
hold off
subplot(2,1,2);
plot(time, xm(2,:), 'Color', 'b', 'LineWidth', 0.9)
hold on
plot(time, xp(2,:), 'Color', 'r', 'LineStyle', '--', 'LineWidth', 1.1)
xlabel('t');
ylabel('x_m_2(t)、x_p_2(t)');
legend('x_m_2(t)', 'x_p_2(t)');
hold off

参考书目

李言俊, 张科. 自适应控制理论及应用[M]. 西北工业大学出版社, 2005.