【控制】自适应控制，模型参考自适应控制，公式推导，有程序有结果图

在模型参考自适应控制框架中，可以人为构造一个所谓的“参考模型”，用以表征希望的闭环系统控制性能。而模型参考自适应控制则是希望求得一种动态调整的反馈控制律，使得系统的闭环控制性能与参考模型的性能可以保持一致。

控制流程如下：

1. 外部指令 $r$ 同时给到参考模型 reference model 和被控对象 plant 的控制器 controller，此时控制器的参数是一个初始的参数，控制器根据这个初始的参数给出一个控制量 $u$ 输入到被控对象，被控对象作出对应的响应 $x$，同时参考模型也根据指令作响应 $x_r$；
2. 由于被控对象存在不确定性，所以其输出与参考模型存在一定的误差 $e=x-x_r$，这个误差会输入到自适应律中，而自适应律正是根据这个误差来调整控制器参数的；
3. 自适应律更新控制参数之后会传递给控制器，在下一个采样时刻，控制器用新的参数计算和新的状态反馈、指令计算控制量，并输出给被控对象；
4. 整个参数更新过程一直循环，知道被控对象的输出与参考模型一致，即两者误差为零的时候停止更新。

1 问题描述

1.1 被控对象

假设存在一个单输入单输出的标量系统，其模型表达式为式 (1)
$$\dot{x}=ax+bu\text{\hspace{1em}(1)}$$

其中 $x$ 表示系统状态，$u$ 表示系统控制输入，$a,b$ 是常数。

1.2 参考模型

关于参考模型，我们有
$${\dot{x}}_r=a_r{x}_r+b_{r}r\text{\hspace{1em}(2)}$$

其中 $a_r<0$。因为只有 $a_r$ 是负定的，系统存在平衡状态。即当外部指令 $r=0$ 时，参考模型能够自动收敛到平衡状态，这样的参考系统才有构建和研究的必要。

1.3 控制目标

要求设计一个合理的控制信号 $u$，使得系统状态 $x(t)$ 跟踪上参考模型的状态 $x_r(t)$，假设 $x_r(t)$ 是解析并有界的，且其微分 ${\dot{x}}_r(t)$ 也是连续且有界的。这个假设在实际中是可以满足的，因为被跟踪信号 $x_r(t)$ 往往是人为设计的。

2 常规解决过程

首先看下当被控系统的参数 $a$ 和 $b$ 均已知时的常规解决办法。

2.1 误差动力学系统

令误差 $e=x-x_r$，则误差的动力学方程为

e ˙ = x ˙ − x ˙ r = a x + b u − ( a r x r + b r r ) = a x + b u − a r x r + a r x − a r x − b r r = a x + b u + a r e − a r x − b r r = a r e + ( a − a r ) x + b u − b r r (3)

\tag{3} e˙​=x˙−x˙r​=ax+bu−(ar​xr​+br​r)=ax+bu−ar​xr​+ar​x−ar​x−br​r=ax+bu+ar​e−ar​x−br​r=ar​e+(a−ar​)x+bu−br​r​(3)

2.2 控制器设计

根据一个事实，即针对一个满足微分方程 $\dot{e}=a\cdot e$ 的变量 $e$，其中 $a$ 是常数。当 $a<0$ 时，对变量 $e$ 的任意初态，均有 $e\to 0$。

那么回到我们的误差动力学方程式 (3)。若想让误差 $e\to 0$，那么就需要式 (3) 右边的所有项小于0。同时在前文已经令 $a_r<0$，那么等式 (3) 右边除 $a_{r}e$ 的其他项咱们也拼凑成小于0，则满足整体小于0。

拼凑形式可以为 $\dot{e}=a_{r}e-Ke$ 的形式，其中要求 $K>0$，整体还是小于 0 的。

因此，可以设计控制器为
( a − a r ) x + b u − b r r = − K e b u = ( a r − a ) x + b r r − K e u = a r − a b x + b r b r − K b e (4)

\tag{4} (a−ar​)x+bu−br​rbuu​=−Ke=(ar​−a)x+br​r−Ke=bar​−a​x+bbr​​r−bK​e​(4)

将控制器式 (4) 代入误差动力学方程 (3) 有
e ˙ = a r e + ( a − a r ) x + b u − b r r = a r e + ( a − a r ) x + b ( a r − a b x + b r b r − K b e ) − b r r = a r e − K e = ( a r − K ) e (5)

\tag{5} e˙​=ar​e+(a−ar​)x+bu−br​r=ar​e+(a−ar​)x+b(bar​−a​x+bbr​​r−bK​e)−br​r=ar​e−Ke=(ar​−K)e​(5)

在描述参考模型时，已经分析过参数 $a_r$ 应该小于0。这里我们再要求 $K>0$，那么就能保证 $a_r-K<0$。显而易见的是此时的误差动力学微分方程，无论误差的初态在哪里，最终都能收敛到 0。

2.3 实验验证

clear
clc

%% Initial states
% controlled system
a = 1.6;
b = 0.5;
x(1,1) = rand;
u(1,1) = rand;

% reference model
a_r = -4;
b_r =  4;
x(1,1) = rand;
x_r(1,1) = 0;
r(1,1) = 0;

% adaptive parameters
K = 1;

%% Time parameters
tBegin = 0;
tFinal = 50;
dT = 0.1;
times = (tFinal-tBegin)/dT;
t(1,1) = tBegin;

% Iteration
for k=1:times
    % record time
    t(:,k+1) = t(:,k) + dT;
    
    % update x_r
    if mod(t(:,k),4) >= 2
        % non-control
        r(:,k+1) = 2;  
    else
        r(:,k+1) = 0;
    end
    dot_x_r = a_r * x_r(:,k) + b_r * r(:,k+1);
    x_r(:,k+1) = x_r(:,k) + dT * dot_x_r;
    
    % calculate error
    e = x(:,k) - x_r(:,k);

    % calculate input
    u(:,k+1) = (a_r-a)/b * x(:,k) + b_r/b * r(:,k) - K/b * e;
    
    % update controlled system
    dot_x = a * x(:,k) + b * u(:,k+1);
    x(:,k+1) = x(:,k) + dT * dot_x;
end

%% 
subplot(2,2,1)
plot(t,r, t,x_r, 'linewidth',1.5)
legend('reference signal $r$', 'reference model state $x_r$', 'interpreter','latex')
grid on

subplot(2,2,2)
plot(t,r, t,x_r, t,x, 'linewidth',1.5)
legend('reference signal $r$', 'reference model state $x_r$', 'controlled state $x$', 'interpreter','latex')
grid on

subplot(2,2,3)
plot(t,x_r-x, 'linewidth',1.5)
legend('error $e$', 'interpreter','latex')
grid on

subplot(2,2,4)
plot(t,u, 'linewidth',1.5)
legend('control $u$', 'interpreter','latex')
grid on
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

3 自适应解决过程，针对参数 $a$ 未知

当然本节的重点不在于常规求解，而是需要自适应求解。也就是当被控系统的参数 $a$ 或 $b$ 未知时，如何使用自适应控制对其进行估计。

3.1 误差动力学系统

令误差 $e=x-x_r$，则误差的动力学方程为

e ˙ = x ˙ − x ˙ r = a x + b u − ( a r x r + b r r ) = a x + b u − a r x r + a r x − a r x − b r r = a x + b u + a r e − a r x − b r r = a r e + ( a − a r ) x + b u − b r r (6)

\tag{6} e˙​=x˙−x˙r​=ax+bu−(ar​xr​+br​r)=ax+bu−ar​xr​+ar​x−ar​x−br​r=ax+bu+ar​e−ar​x−br​r=ar​e+(a−ar​)x+bu−br​r​(6)

3.2 控制器设计

同样的，可以设计控制器为
$$u=\frac{a_r-\hat{a}}{b}x+\frac{b_r}{b}r-\frac{K}{b}e\text{\hspace{1em}(7)}$$

其中 $\hat{a}$ 表示 $a$ 的估计值。因为我们对被控系统的实际参数并不知道，因此控制器也不可能获取到真实的 $a$，所以需要使用估计值来设计控制器。但是控制器的结构不会发生变换，还可常规的结构相同。

将控制器式 (7) 代入误差动力学方程 (6) 有
e ˙ = a r e + ( a − a r ) x + b u − b r r = a r e + ( a − a r ) x + b ( a r − a ^ b x + b r b r − K b e ) − b r r = a r e + ( a − a r ) x + ( a r − a ^ ) x + b r r − K e − b r r = a r e + ( a − a ^ ) x − K e (8)

\tag{8} e˙​=ar​e+(a−ar​)x+bu−br​r=ar​e+(a−ar​)x+b(bar​−a^​x+bbr​​r−bK​e)−br​r=ar​e+(a−ar​)x+(ar​−a^)x+br​r−Ke−br​r=ar​e+(a−a^)x−Ke​(8)

---

3.3 更新估计参数

令真实参数 $a$ 与其估计值 $\hat{a}$ 的误差为
a ~ = a ^ − a (9)

\tag{9} a~​=a^−a​(9)

我们将估计值 $\hat{a}$ 的和其与真实值的误差 $\tilde{a}$ 代入到式 (8) 中，这样做的目的是消除未知的真实参数 $a$，得到一个误差仅和估计值有关的式子。

e ˙ = a r e + ( a − a ^ ) x − K e = a r e − K e − a ~ x (10)

\tag{10} e˙​=ar​e+(a−a^)x−Ke=ar​e−Ke−a~x​(10)

---

定义李雅普诺夫函数为
V ( e , a ~ ) = 1 2 e 2 + 1 2 η a a ~ 2 (11)

\tag{11} V(e,a~)​=21​e2+2ηa​1​a~2​(11)

李雅普诺夫函数式 (11) 关于时间求导有
V ˙ ( e , a ~ ) = e e ˙ + 1 η a a ~ a ~ ˙ = e e ˙ + 1 η a a ~ ( a ^ ˙ − a ˙ ) = e e ˙ + 1 η a a ~ a ^ ˙ = e ( a r e − K e − a ~ x ) + 1 η a a ~ a ^ ˙ = a r e 2 − K e 2 − a ~ x e + 1 η a a ~ a ^ ˙ = ( a r − K ) e 2 + a ~ ( 1 η a a ^ ˙ − x e ) (12)

\tag{12} V˙(e,a~)​=ee˙+ηa​1​a~a~˙=ee˙+ηa​1​a~(a^˙−a˙)=ee˙+ηa​1​a~a^˙=e(ar​e−Ke−a~x)+ηa​1​a~a^˙=ar​e2−Ke2−a~xe+ηa​1​a~a^˙=(ar​−K)e2+a~(ηa​1​a^˙−xe)​(12)

在描述参考模型时，已经分析过参数 $a_r$ 应该小于0。这里我们再要求 $K>0$，那么就能保证 $a_r-K<0$。

关于等式 (12) 右侧的第二项，我们就可以设计自适应律为
( 1 η a a ^ ˙ − x e ) = 0 1 η a a ^ ˙ = x e a ^ ˙ = η a ⋅ x ⋅ e (13)

\tag{13} (ηa​1​a^˙−xe)ηa​1​a^˙a^˙​=0=xe=ηa​⋅x⋅e​(13)

显而易见的是此时的误差动力学微分方程，无论初态在哪里，最终都能收敛到0。

3.4 实验验证

clear
clc

%% Initial states
% controlled system
a = 1.6;
b = 0.5;
x(1,1) = rand;
u(1,1) = rand;

% reference model
a_r = -4;
b_r =  4;
x(1,1) = rand;
x_r(1,1) = 0;
r(1,1) = 0;

% adaptive parameters
K = 1;
eta_a = 0.4;
hat_a(1,1) = rand;

%% Time parameters
tBegin = 0;
tFinal = 50;
dT = 0.1;
times = (tFinal-tBegin)/dT;
t(1,1) = tBegin;

% Iteration
for k=1:times
    % record time
    t(:,k+1) = t(:,k) + dT;
    
    % update x_r
    if mod(t(:,k),4) >= 2
        % non-control
        r(:,k+1) = 2;  
    else
        r(:,k+1) = 0;
    end
    dot_x_r = a_r * x_r(:,k) + b_r * r(:,k+1);
    x_r(:,k+1) = x_r(:,k) + dT * dot_x_r;
    
    % calculate error
    e = x(:,k) - x_r(:,k);
    
    % calculate adaptive parameters
    dot_hat_a = -eta_a * x(:,k) * e;
    hat_a(:,k+1) = hat_a(:,k) + dT * dot_hat_a;

    % calculate input
    u(:,k+1) = (a_r-hat_a(:,k+1))/b * x(:,k) + b_r/b * r(:,k) - K/b * e;
    
    % update controlled system
    dot_x = hat_a(:,k+1) * x(:,k) + b * u(:,k+1);
    x(:,k+1) = x(:,k) + dT * dot_x;
end

%% 
subplot(2,2,1)
plot(t,r, t,x_r, 'linewidth',1.5)
legend('reference signal $r$', 'reference model state $x_r$', 'interpreter','latex')
grid on

subplot(2,2,2)
plot(t,r, t,x_r, t,x, 'linewidth',1.5)
legend('reference signal $r$', 'reference model state $x_r$', 'controlled state $x$', 'interpreter','latex')
grid on

subplot(2,2,3)
plot(t,u, 'linewidth',1.5)
legend('control $u$', 'interpreter','latex')
grid on

subplot(2,2,4)
plot(t,hat_a, t,a*ones(times+1,1), 'linewidth',1.5)
legend('estimator $\hat{a}$', 'unknown param $a$', 'interpreter','latex')
grid on
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
76
77
78
79

4 自适应解决过程，针对参数 $b$ 未知

当然本节的重点不在于常规求解，而是需要自适应求解。也就是当被控系统的参数 $a$ 或 $b$ 未知时，如何使用自适应控制对其进行估计。

4.1 误差动力学系统

令误差 $e=x-x_r$，则误差的动力学方程为

e ˙ = a r e + ( a − a r ) x + b u − b r r (14)

\tag{14} e˙​=ar​e+(a−ar​)x+bu−br​r​(14)

4.2 控制器设计

同样的，可以设计控制器为
$$u=\frac{a_r-a}{\hat{b}}x+\frac{b_r}{\hat{b}}r-\frac{K}{\hat{b}}e\text{\hspace{1em}(15)}$$

其中 $\hat{b}$ 表示 $b$ 的估计值。因为我们对被控系统的实际参数并不知道，因此控制器也不可能获取到真实的 $b$，所以需要使用估计值来设计控制器。但是控制器的结构不会发生变换，还可常规的结构相同。

将控制器式 (15) 代入误差动力学方程 (14) 有
e ˙ = a r e + ( a − a r ) x + b u − b r r = a r e + ( a − a r ) x + b ( a r − a b ^ x + b r b ^ r − K b ^ e ) − b r r = a r e + a x − a r x + b a r b ^ x − b a b ^ x + b b r b ^ r − b K b ^ e − b r r (16)

\tag{16} e˙​=ar​e+(a−ar​)x+bu−br​r=ar​e+(a−ar​)x+b(b^ar​−a​x+b^br​​r−b^K​e)−br​r=ar​e+ax−ar​x+b^bar​​x−b^ba​x+b^bbr​​r−b^bK​e−br​r​(16)

---

4.3 更新估计参数

令真实参数 $b$ 与其估计值 $\hat{b}$ 的误差为
b ~ = b ^ − b (17)

\tag{17} b~​=b^−b​(17)

我们将估计值 $\hat{b}$ 的和其与真实值的误差 $\tilde{b}$ 代入到式 (16) 中，这样做的目的是消除未知的真实参数 $b$，得到一个误差仅和估计值有关的式子。

e ˙ = a r e + a x − a r x + b a r b ^ x − b a b ^ x + b b r b ^ r − b K b ^ e − b r r = a r e + a x − a r x − b r r + ( b ^ − b ~ ) a r b ^ x − ( b ^ − b ~ ) a b ^ x + ( b ^ − b ~ ) b r b ^ r − ( b ^ − b ~ ) K b ^ e = a r e + a x − a r x − b r r + a r x − b ~ a r b ^ x − a x + b ~ a b ^ x + b r r − b ~ b r b ^ r − K e + b ~ K b ^ e = a r e − K e − b ~ a r b ^ x + b ~ a b ^ x − b ~ b r b ^ r + b ~ K b ^ e = a r e − K e + − b ~ a r x + b ~ a x − b ~ b r r + b ~ K e b ^ (18)

\tag{18} e˙​=ar​e+ax−ar​x+b^bar​​x−b^ba​x+b^bbr​​r−b^bK​e−br​r=ar​e+ax−ar​x−br​r+b^(b^−b~)ar​​x−b^(b^−b~)a​x+b^(b^−b~)br​​r−b^(b^−b~)K​e=ar​e+ax−ar​x−br​r+ar​x−b^b~ar​​x−ax+b^b~a​x+br​r−b^b~br​​r−Ke+b^b~K​e=ar​e−Ke−b^b~ar​​x+b^b~a​x−b^b~br​​r+b^b~K​e=ar​e−Ke+b^−b~ar​x+b~ax−b~br​r+b~Ke​​(18)

---

定义李雅普诺夫函数为
V ( e , b ~ ) = 1 2 e 2 + 1 2 η b b ~ 2 (19)

\tag{19} V(e,b~)​=21​e2+2ηb​1​b~2​(19)

李雅普诺夫函数式 (19) 关于时间求导有
V ˙ ( e , b ~ ) = e e ˙ + 1 η b b ~ b ~ ˙ = e e ˙ + 1 η b b ~ ( b ^ ˙ − b ˙ ) = e e ˙ + 1 η b b ~ b ^ ˙ = e ( a r e − K e + − b ~ a r x + b ~ a x − b ~ b r r + b ~ K e b ^ ) + 1 η b b ~ b ^ ˙ = a r e 2 − K e 2 + − b ~ a r x e + b ~ a x e − b ~ b r r e + b ~ K e 2 b ^ + 1 η b b ~ b ^ ˙ = a r e 2 − K e 2 + b ~ ( − a r x e + a x e − b r r e + K e 2 b ^ + 1 η b b ^ ˙ ) (20)

\tag{20} V˙(e,b~)​=ee˙+ηb​1​b~b~˙=ee˙+ηb​1​b~(b^˙−b˙)=ee˙+ηb​1​b~b^˙=e(ar​e−Ke+b^−b~ar​x+b~ax−b~br​r+b~Ke​)+ηb​1​b~b^˙=ar​e2−Ke2+b^−b~ar​xe+b~axe−b~br​re+b~Ke2​+ηb​1​b~b^˙=ar​e2−Ke2+b~(b^−ar​xe+axe−br​re+Ke2​+ηb​1​b^˙)​(20)

在描述参考模型时，已经分析过参数 $a_r$ 应该小于0。这里我们再要求 $K>0$，那么就能保证 $a_r-K<0$。

关于等式 (12) 右侧的第二项，我们就可以设计自适应律为
− a r x e + a x e − b r r e + K e 2 b ^ + 1 η b b ^ ˙ = 0 1 η b b ^ ˙ = a r x e − a x e + b r r e − K e 2 b ^ b ^ ˙ = η b a r x e − a x e + b r r e − K e 2 b ^ (13)

\tag{13} b^−ar​xe+axe−br​re+Ke2​+ηb​1​b^˙ηb​1​b^˙b^˙​=0=b^ar​xe−axe+br​re−Ke2​=ηb​b^ar​xe−axe+br​re−Ke2​​(13)

显而易见的是此时的误差动力学微分方程，无论初态在哪里，最终都能收敛到0。

4.4 实验验证

clear
clc

%% Initial states
% controlled system
a = 1.6;
b = 0.5;
x(1,1) = rand;
u(1,1) = rand;

% reference model
a_r = -4;
b_r =  4;
x(1,1) = rand;
x_r(1,1) = 0;
r(1,1) = 0;

% adaptive parameters
K = 1;
eta_b = 1;
hat_b(1,1) = rand;

%% Time parameters
tBegin = 0;
tFinal = 50;
dT = 0.1;
times = (tFinal-tBegin)/dT;
t(1,1) = tBegin;

% Iteration
for k=1:times
    % record time
    t(:,k+1) = t(:,k) + dT;
    
    % update x_r
    if mod(t(:,k),4) >= 2
        % non-control
        r(:,k+1) = 2;  
    else
        r(:,k+1) = 0;
    end
    dot_x_r = a_r * x_r(:,k) + b_r * r(:,k+1);
    x_r(:,k+1) = x_r(:,k) + dT * dot_x_r;
    
    % calculate error
    e = x(:,k) - x_r(:,k);
    
    % calculate adaptive parameters
    dot_hat_b = eta_b * (a_r*x(:,k)*e - a*x(:,k)*e + b_r*r(:,k)*e - K*e^2) / hat_b(:,k);
    hat_b(:,k+1) = hat_b(:,k) + dT * dot_hat_b;

    % calculate input
    u(:,k+1) = (a_r-a)/hat_b(:,k+1) * x(:,k) + b_r/hat_b(:,k+1) * r(:,k) - K/hat_b(:,k+1) * e;
    
    % update controlled system
    dot_x = a * x(:,k) + hat_b(:,k+1) * u(:,k+1);
    x(:,k+1) = x(:,k) + dT * dot_x;
end

%% 
subplot(2,2,1)
plot(t,r, t,x_r, 'linewidth',1.5)
legend('reference signal $r$', 'reference model state $x_r$', 'interpreter','latex')
grid on

subplot(2,2,2)
plot(t,r, t,x_r, t,x, 'linewidth',1.5)
legend('reference signal $r$', 'reference model state $x_r$', 'controlled state $x$', 'interpreter','latex')
grid on

subplot(2,2,3)
plot(t,u, 'linewidth',1.5)
legend('control $u$', 'interpreter','latex')
grid on

subplot(2,2,4)
plot(t,hat_b, t,b*ones(times+1,1), 'linewidth',1.5)
legend('estimator $\hat{b}$', 'unknown param $b$', 'interpreter','latex')
grid on
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
76
77
78
79

5 总结

Ref.