一、问题描述

设控制对象的状态方程为
$${\dot{\mathbf{x}}}_p={\mathbf{A}}_p(t)x_p+{\mathbf{b}}_p(t)u\text{\hspace{1em}(1)}$$
式中
A p = [ 0 1 − 6 − 7 ] , b p = [ 2 4 ] (2) \boldsymbol{A}_{p}=\left[

\right], \quad \boldsymbol{b}_{p}=\left[\right] \tag{2} Ap​=[0−6​1−7​],bp​=[24​](2)
参考模型的状态方程为
$${\dot{\mathbf{x}}}_m={\mathbf{A}}_m{x}_m+{\mathbf{b}}_{m}r\text{\hspace{1em}(3)}$$
式中
A m = [ 0 1 − 10 − 5 ] , b m = [ 1 2 ] (4) \boldsymbol{A}_{m}=\left[

01−10−5" role="presentation" style="position: relative;">0−101−5$$\begin{array}{cc}0& 1\\ -10& -5\end{array}$$

\right], \quad \boldsymbol{b}_{m}=\left[

12" role="presentation" style="position: relative;">12$$\begin{array}{l}1\\ 2\end{array}$$

\right] \tag{4} Am​=[0−10​1−5​],bm​=[12​](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)式，再根据矩阵迹的性质，于是有
V ˙ = 1 2 e T ( P A m + A m T P ) e + tr ⁡ ( Φ ˙ T Γ 1 − 1 Φ + x p e T P b m K ˉ − 1 Φ ) + tr ⁡ ( Ψ ˙ T Γ 2 − 1 Ψ + r e T P b m K ˉ − 1 Ψ ) (12)

\tag{12} V˙=​21​eT(PAm​+AmT​P)e+tr(Φ˙TΓ1−1​Φ+xp​eTPbm​Kˉ−1Φ)+tr(Ψ˙TΓ2−1​Ψ+reTPbm​Kˉ−1Ψ)​(12)
为满足李雅普诺夫第二法，需保证(12)式是负定的，对应的情况为(12)式第一项是负定的，后两项都为零。

因为 ${\mathbf{A}}_m$ 为稳定矩阵，则可选定正定对称阵 $Q$，使 $\mathbf{P}{\mathbf{A}}_m+{\mathbf{A}}_m^{\mathbf{T}}\mathbf{P}=-\mathbf{Q}$ 成立。同时根据上述对应情况，$\Phi$ 和 $\Psi$ 的选择如下：
Φ ˙ = − Γ 1 ( b m K ˉ − 1 ) T P e x p T Ψ ˙ = − Γ 2 ( b m K ˉ − 1 ) T P e r T (13)

\tag{13} Φ˙Ψ˙​=−Γ1​(bm​Kˉ−1)TPexpT​=−Γ2​(bm​Kˉ−1)TPerT​(13)
当 ${\mathbf{A}}_p$ 和 ${\mathbf{b}}_p$ 为常值或缓慢变化时，可得自适应调节规律：
F ( t ) = ∫ 0 t Γ 1 ( b m K ˉ − 1 ) T P e x p T d τ + F ( 0 ) K ( t ) = ∫ 0 t Γ 2 ( b m K ˉ − 1 ) T P e r d τ + K ( 0 ) (14)

F(t)=∫0tΓ1(bmK¯−1)TPexpTdτ+F(0)K(t)=∫0tΓ2(bmK¯−1)TPerdτ+K(0)" role="presentation" style="position: relative;">F(t)K(t)=∫t0Γ1(bmK¯−1)TPexTpdτ+F(0)=∫t0Γ2(bmK¯−1)TPerdτ+K(0)$$\begin{array}{rl}F(t)& ={\int }_0^t{\mathrm{\Gamma }}_1{\left({\mathit{b}}_m{\overline{K}}^{-1}\right)}^T\mathit{P}\mathit{e}{\mathit{x}}_p^{T}d\tau +F(0)\\ K(t)& ={\int }_0^t{\mathrm{\Gamma }}_2{\left({\mathit{b}}_m{\overline{K}}^{-1}\right)}^T\mathit{P}\mathit{e}rd\tau +K(0)\end{array}$$

\tag{14} F(t)K(t)​=∫0t​Γ1​(bm​Kˉ−1)TPexpT​dτ+F(0)=∫0t​Γ2​(bm​Kˉ−1)TPerdτ+K(0)​(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$ 的关系如下：
b m K ˉ − 1 = b p = [ 2 4 ] (16) \boldsymbol{b}_{m} \bar{K}^{-1}=\boldsymbol{b}_{p}=\left[

\right] \tag{16} bm​Kˉ−1=bp​=[24​](16)
选取(14)式中的部分自适应参数如下：
P = [ 3 1 1 1 ] , Γ 1 = Γ 2 = 1 (17) \boldsymbol{P}=\left[

3111" role="presentation" style="position: relative;">3111$$\begin{array}{ll}3& 1\\ 1& 1\end{array}$$

\right], \quad \Gamma_{1}=\Gamma_{2}=1 \tag{17} P=[31​11​],Γ1​=Γ2​=1(17)
所以可得最终的自适应规律：
F ( t ) = ∫ 0 t [ 2 4 ] [ 3 1 1 1 ] e x p T d τ + F ( 0 ) K ( t ) = ∫ 0 t [ 2 4 ] [ 3 1 1 1 ] e r d τ + K ( 0 ) (18)

F(t)=∫0t[24][3111]expTdτ+F(0)K(t)=∫0t[24][3111]erdτ+K(0)" role="presentation" style="position: relative;">F(t)K(t)=∫t0[24][3111]exTpdτ+F(0)=∫t0[24][3111]erdτ+K(0)$$\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]\mathit{e}{\mathit{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]\mathit{e}\mathit{r}\mathit{d}\tau +K(0)\end{array}$$

\tag{18} F(t)K(t)​=∫0t​[2​4​][31​11​]expT​dτ+F(0)=∫0t​[2​4​][31​11​]erdτ+K(0)​(18)
下将上述连续自适应规律进行离散化，用于实际的数值仿真实验。设数值积分步长为 $h$，则各时刻的参考模型状态向量及控制对象状态向量如下：
x m ( k + 1 ) = x m ( k ) + h [ A m ( k ) x m ( k ) + B m ( k ) r ( k ) ] x p ( k + 1 ) = x p ( k ) + h [ A p ( k ) x p ( k ) + B p ( k ) u ( k ) ] (19)

xm(k+1)=xm(k)+h[Am(k)xm(k)+Bm(k)r(k)]xp(k+1)=xp(k)+h[Ap(k)xp(k)+Bp(k)u(k)]" role="presentation" style="position: relative;">xm(k+1)xp(k+1)=xm(k)+h[Am(k)xm(k)+Bm(k)r(k)]=xp(k)+h[Ap(k)xp(k)+Bp(k)u(k)]$$\begin{array}{rl}{\mathit{x}}_m(k+1)& ={\mathit{x}}_m(k)+h\left[{\mathit{A}}_m(k){\mathit{x}}_m(k)+{\mathit{B}}_m(k)r(k)\right]\\ {\mathit{x}}_p(k+1)& ={\mathit{x}}_p(k)+h\left[{\mathit{A}}_p(k){\mathit{x}}_p(k)+{\mathit{B}}_p(k)u(k)\right]\end{array}$$

\tag{19} xm​(k+1)xp​(k+1)​=xm​(k)+h[Am​(k)xm​(k)+Bm​(k)r(k)]=xp​(k)+h[Ap​(k)xp​(k)+Bp​(k)u(k)]​(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)}$$
最终，还需将自适应规律离散化：
F ( k ) = h ⋅ ∑ j = 0 k b p T P e ( k ) ( x p ( k ) ) T + F ( 0 ) K ( k ) = h ⋅ ∑ j = 0 k b p T P e ( k ) r ( k ) + K ( 0 ) (22)

F(k)=h⋅∑j=0kbpTPe(k)(xp(k))T+F(0)K(k)=h⋅∑j=0kbpTPe(k)r(k)+K(0)" role="presentation" style="position: relative;">F(k)K(k)=h⋅∑j=0kbTpPe(k)(xp(k))T+F(0)=h⋅∑j=0kbTpPe(k)r(k)+K(0)$$\begin{array}{rl}F(k)& =h\cdot \sum _{j=0}^k{\mathit{b}}_p^T\mathit{P}\mathit{e}(k){\left({\mathit{x}}_p(k)\right)}^T+F(0)\\ K(k)& =h\cdot \sum _{j=0}^k{\mathit{b}}_p^T\mathit{P}\mathit{e}(k)r(k)+K(0)\end{array}$$

\tag{22} F(k)K(k)​=h⋅j=0∑k​bpT​Pe(k)(xp​(k))T+F(0)=h⋅j=0∑k​bpT​Pe(k)r(k)+K(0)​(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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

参考书目

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