一种基于卡尔曼滤波器的预测滑模控制方法

1.本发明属于系统控制技术领域，具体涉及一种基于卡尔曼滤波器的预测滑模控制方法。


背景技术：

2.在实际的控制系统中，如飞行控制系统、倒立摆系统、机械系统等，由于受工作环境，工作任务等一系列问题的影响，在建立其数学模型时通常都面临着外部干扰随机、传感器测量误差等一系列复杂性问题。以高超声速飞行器姿态控制系统为例，当其以超高的速度在大气层中飞行时，存在着气动参数的剧烈变化、随机阵风干扰、以及传感器对姿态角的测量误差等一系列不确定因素，这些都是在控制器设计时不得不考虑的问题。因此在实际系统的控制当中，减弱或消除外部随机干扰和传感器测量误差对系统的影响，是非常有必要的。
3.在现有的研究中，针对非线性系统的控制，提出了一系列的控制方法。例如，模型预测控制方法、基于干扰观测器的控制方法、滑模控制方法等，均取得的了一定的成果。但是在现实的控制系统中，外部干扰的随机性以及传感器的测量误差，在一定程度上限制了这些方法的使用。


技术实现要素：

4.解决的技术问题：针对上述技术问题，本发明提供了一种基于卡尔曼滤波器的预测滑模控制方法，能够使被控系统在随机外部干扰和传感器测量误差的情况下保持稳定，并提高被控系统的控制精度。
5.技术方案：一种基于卡尔曼滤波器的预测滑模控制方法，包括步骤如下：
6.步骤1：在不考虑噪声干扰的情况下，以标称系统状态空间模型为预测模型，采用预测滑模控制方法，设计预测滑模控制器；
7.步骤2：将外部随机干扰等效为系统的过程噪声，加入系统的输入端；将传感器的测量误差等效为系统的观测噪声，加入系统的输出端；并设计卡尔曼滤波器，得到滤波后的新状态；
8.步骤3：用滤波后的新状态代替之前的无干扰状态，进而求出基于卡尔曼滤波器的预测滑模控制器。
9.优选的，所述步骤2中过程噪声和观测噪声均为白噪声。
10.优选的，所述步骤1具体包括步骤如下：
11.1.1以k时刻为当前时刻，计算k+p时刻系统预测输出；
12.1.2设计滑模面，并预测k+p时刻滑模面及趋近律；
13.1.3设计无干扰情况下的预测滑模控制器。
14.优选的，所述步骤1.1具体包括步骤如下：
15.1.1.1选取预测模型；
16.1.1.2选取基函数；
17.1.1.3设计参考轨迹；
18.1.1.4计算k+p时刻预测模型输出；
19.1.1.5设计误差补偿，计算k+p时刻系统预测输出。
20.优选的，所述步骤1.1.2中基函数为阶跃函数。
21.优选的，所述步骤2具体包括步骤如下：
22.2.1选取实际受扰系统的状态方程；
23.2.2设计卡尔曼滤波器。
24.优选的，所述步骤2.2具体包括步骤如下：
25.2.2.1选取无干扰及初始输入为0的系统状态空间表达式；
26.2.2.2计算输出信号估计误差；
27.2.2.3计算未校正状态估计误差均方值；
28.2.2.4计算校正后状态估计误差均方值；
29.2.2.5计算最佳卡尔曼增益和最小校正后状态估计误差均方值；
30.2.2.6总结卡尔曼递推公式。
31.有益效果：本发明能够保持被控系统稳定，并在白噪声干扰的作用下仍然能够精确的跟踪参考轨迹，具有较强的鲁棒性。
32.本发明从随机干扰产生的机理入手，在预测滑模控制的基础之上，来研究系统在白噪声干扰下的控制问题；并将外部随机干扰等效为系统的过程噪声，将传感器的测量误差等效为系统的观测噪声，提出一种基于卡尔曼滤波器的预测滑模控制方法，以解决预测滑模控制方法在白噪声干扰时，控制效果不佳的问题，从而实现白噪声干扰下对系统的有效控制。
附图说明
33.图1：本发明控制系统原理图；
34.图2：本发明设计流程示意图；
35.图3：卡尔曼滤波器递推示意图；
36.图4：一级倒立摆系统示意图；
37.图5：(a)白噪声干扰下滤波前小车位置仿真示意图，(b)白噪声干扰下滤波后小车位置仿真示意图；
38.图6：(a)白噪声干扰下滤波前杆的摆角仿真示意图，(b)白噪声干扰下滤波后杆的摆角仿真示意图；
39.图7：(a)白噪声干扰下滤波前小车速度仿真示意图，(b)白噪声干扰下滤波后小车速度仿真示意图；
40.图8：(a)白噪声干扰下滤波前杆的角速度仿真示意图，(b)白噪声干扰下滤波后杆的角速度仿真示意图；
41.图9：(a)白噪声干扰下滤波前系统控制输入仿真示意图，(b)白噪声干扰下滤波后系统控制输入仿真示意图。
具体实施方式
42.下面结合附图和具体实施例对本发明作进一步描述。
43.实施例1
44.一种基于卡尔曼滤波器的预测滑模控制方法，具体的控制系统原理图如图1所示，将受扰输入和输出经过卡尔曼滤波器后，估计新的状态变量，再作用于预测滑模控制器，从而得到新的基于卡尔曼滤波器的预测滑模控制器。
45.具体的控制器设计方法如图2所示，包括以下步骤：
46.步骤1：在不考虑噪声干扰的情况下，以系统的状态方程为预测模型，以k时刻为当前时刻，采用预测滑模控制方法，设计预测滑模控制器。具体如下：
47.1.1计算k+p时刻系统预测输出
48.(1)选取预测模型
49.选取标称系统离散状态空间方程为预测模型如下：
50.x(k)＝ax(k-1)+bu(k-1)
51.y(k)＝cx(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
52.(2)选取基函数
53.本发明选取阶跃函数作为输入的基函数，因此由阶跃函数的性质可得如下关系：
54.u(k+i)＝u(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
55.(3)设计参考轨迹
56.参考轨迹可根据预测模型输出和实际输出设计，可使系统输出平稳趋近于设定值，则在预测时域i∈{1,2,...,p}内参考轨迹为：
57.yr(k+i)＝cr(k)-λi(cr(k)-y
p
(k))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
58.其中p为预测步长，yr(k+i)为k+i时刻的参考输出，cr(k)为设定值，(k)为设定值，ts为采样时间，tr为期望相应时间，y
p
(k)为k时刻系统实际输出。
59.(4)计算k+p时刻预测模型输出
60.由预测模型式(1)和式(2)可知：
61.x(k+1)＝ax(k)+bu(k)
62.x(k+2)＝ax(k+1)+bu(k+1)＝a2x(k)+(ab+b)u(k)
63....
64.x(k+p)＝a
p
x(k)+(a
p-1
b+a
p-2
b+...+ab+b)u(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
65.因此k+p时刻预测模型输出为：
66.y(k+p)＝cx(k+p)
67.ꢀꢀꢀꢀꢀꢀ
＝ca
p
x(k)+c(a
p-1
b+a
p-2
b+...+ab+b)u(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
68.令：
69.g＝c(a
p-1
b+a
p-2
b+...+ab+b)
70.则：
71.y(k+p)＝ca
p
x(k)+gu(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
72.(5)设计误差补偿，计算k+p时刻系统预测输出
73.选取k时刻实际输出与预测模型输出之间的误差作为整个预测时域i∈{1,2,...,
p}内的误差补偿，具体补偿方式如下：
74.e(k+i)＝y
p
(k)-y(k)＝y
p
(k)-cx(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
75.用e(k+i)对预测p步时的预测模型输出进行修正，得到k+p时刻的系统预测输出为：
76.ym(k+p)＝y(k+p)+e(k+i)
77.ꢀꢀꢀꢀꢀꢀꢀ
＝(ca
p-c)x(k)+gu(k)+y
p
(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
78.1.2设计滑模面，并预测k+p时刻滑模面及趋近律
79.(1)设计滑模面
80.在实际系统中，初始状态可处于任意位置，因此为使系统平稳地趋近稳定状态，设计一个滑模面，使系统的运动不断地趋近于滑模面稳定。
81.设计滑模面为：
82.s(k)＝ce(ym(k)-yr(k))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
83.计算k+p时刻滑模面为：
84.s(k+p)＝ce(ym(k+p)-yr(k+p))
85.ꢀꢀꢀꢀꢀꢀ
＝ce((ca
p-c)x(k)+(i-λ
p
)y
p
(k)-(i-λ
p
)cr(k)+gu(k))
86.式中i为单位阵，令：
87.k0＝i-λ
p
88.k1＝c-ca
p
89.则：
90.s(k+p)＝ce(-k1x(k)+k0y
p
(k)-k0cr(k)+gu(k))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
91.(2)设计趋近律
92.选取指数趋近律如下所示：
[0093][0094]
其中ε，q均是大于零的常数。
[0095]
则在预测时域i∈{1,2,...,p}内的指数趋近律为：
[0096][0097]
sr(k+i)＝-εit
s sgn(s(k))+(1-qitss(k))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0098]
则在预测p步的滑模趋近律为：
[0099]
sr(k+p)＝-εpt
s sgn(s(k))+(1-qptss(k))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0100]
1.3计算预测滑模控制器
[0101]
由式(10)和式(13)可知，预测滑模控制器为：
[0102]
u(k)＝(ceg)-1
((1-qptss(k))-εptssgn(s(k))-ce(-k1x(k)+k0y
p
(k)-k0cr(k)))
ꢀꢀ
(14)
[0103]
步骤2：分析随机干扰产生原因，将外部随机干扰等效为系统的过程噪声，将传感器的测量误差等效为系统的观测噪声，且均为白噪声，分别在系统输入端加入过程噪声，在系统的输出端加入测量噪声，并设计卡尔曼滤波器，得到滤波后的新状态。具体如下：
[0104]
2.1选取实际受扰系统的状态方程
[0105]
选取实际受扰系统如下所示：
[0106]
x(k)＝ax(k-1)+bu(k-1)+w(k-1)
[0107]
y(k)＝cx(k)+v(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0108]
其中过程噪声w(k)和测量噪声v(k)均为正态白噪声。
[0109]
2.2设计卡尔曼滤波器
[0110]
卡尔曼滤波器是采用递推算法实现的，其基本思想是，先不考虑输入信号、过程噪声以及测量噪声对系统的影响，得到状态变量和输出信号的估计值，再用输出信号的估计误差加权后来校正状态变量的估计值。
[0111]
(1)选取无干扰及初始输入为0的系统状态空间表达式：
[0112]
当不考虑输入信号、过程噪声和测量噪声对系统的影响时，系统状态方程如下所示：
[0113][0114][0115]
其中为未经校正的状态变量估计值，为前一时刻校正后的状态变量估计值，为输出信号估计值。
[0116]
(2)计算输出信号估计误差ey(k)：
[0117][0118]
则误差补偿后系统的状态变量估计值为：
[0119][0120]
其中h为卡尔曼增益，具体为加权矩阵，卡尔曼滤波器的关键就在于调节h的大小，使得校正后状态的估计误差均方值最小。
[0121]
(3)计算未校正状态估计误差均方值
[0122]
设p
′
(k)为当前时刻未校正状态估计误差均方值，且为对称矩阵，则：
[0123][0124]
其中p(k-1)为前一时刻校正后状态的估计误差均方值，q(k)为过程噪声方差。
[0125]
(4)计算校正后状态估计误差均方值
[0126]
设e
x
(k)为当前时刻校正后状态估计误差且u(k-1)＝0，则：
[0127][0128]
由式(18)知：
[0129][0130]
则：
[0131][0132]
令：
[0133][0134]
b＝(i-hc)w(k-1)
[0135]
c＝hv(k)
[0136]
则：
[0137]ex
(k)＝a+b-c
[0138][0139]
其中假设状态变量的增益矩阵a不随时间发生变化，可知x(k-1)与w(k-1)、v(k)互不相关，与w(k-1)、v(k)互不相关。可得：
[0140]
p(k)＝e[aa
t
]+e[bb
t
]+e[cc
t
]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(24)
[0141]
化简可得：
[0142][0143]
e[bb
t
]＝[(i-hc)w(k-1)w
t
(k-1)(i-hc)
t
]
[0144]
ꢀꢀꢀꢀꢀꢀ
＝(i-hc)q(k-1)(i-hc)
t
[0145]
e[cc
t
]＝e[hv(k)v
t
(k)h
t
]＝hr(k)h
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(25)
[0146]
将式(19)、式(25)带入式(24)得：
[0147]
p(k)＝(i-hc)p
′
(k)(i-hc)
t
+hr(k)h
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(26)
[0148]
(5)计算最佳卡尔曼增益，和最小校正后状态估计误差均方值
[0149]
由式(26)知：
[0150]
p(k)＝p
′
(k)-hcp
′
(k)-p
′
(k)c
tht
+hcp
′
(k)c
tht
+hr(k)h
t
[0151]
ꢀꢀꢀꢀ
＝p
′
(k)-hcp
′
(k)-p
′
(k)c
tht
+h(cp
′
(k)c
t
+r(k))h
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(27)
[0152]
其中cp
′
(k)c
t
+r(k)为正定矩阵，记：
[0153]
cp
′
(k)c
t
+r(k)＝ss
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(28)
[0154]
令：
[0155]ft
＝(p
′
(k)c
t
)
t
＝cp
′
(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(29)
[0156]
将式(29)带入式(27)中得：
[0157]
p(k)＝p
′
(k)-hf
t-fh
t
+h(ss
t
)h
t
[0158]
ꢀꢀꢀꢀ
＝(hs-f(s
t
)-1
)(hs-f(s
t
)-1
)
t
+p
′
(k)-f(ss
t
)-1ft
[0159]
ꢀꢀꢀꢀ
＝(hs-f(s
t
)-1
)(hs-f(s
t
)-1
)
t
+p
′
(k)-p
′
(k)c
t
(cp
′
(k)c
t
+r(k))-1
cp
′
(k)
ꢀꢀ
(30)
[0160]
式中后两项与h无关，要使p(k)最小，则：
[0161]
hs-f(s
t
)-1
＝0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(31)
[0162]hopt
＝u(s
t
)-1
s-＝u(ss
t
)-1
[0163]
ꢀꢀꢀ
＝p
′
(k)c
t
(cp
′
(k)c
t
+r(k))-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(32)
[0164]
则p(k)最小为：
[0165]
p(k)＝p
′
(k)-p
′
(k)c
t
(cp
′
(k)c
t
+r(k))-1
cp
′
(k)
[0166]
ꢀꢀꢀꢀ
＝(i-h
opt
c)p
′
(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(33)
[0167]
(1)总结卡尔曼递推公式
[0168]
由以上分析可知，卡尔曼递推公式如下所示：
[0169]
p
′
(k)＝ap(k-1)a
t
+q(k-1)
[0170]
h＝p
′
(k)c
t
(cp
′
(k)c
t
+r(k))-1
[0171][0172]
p(k)＝(i-hc)p
′
(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(34)
[0173]
卡尔曼滤波器递推流程图如图3所示。
[0174]
步骤3：用滤波后的新状态代替之前的无干扰状态，进而求出基于卡尔曼滤波器的预测滑模控制器。具体如下：
[0175]
3.1选取系统新的状态空间表达式
[0176]
由上一步可知，经过卡尔曼滤波后，系统新的状态空间表达式为：
[0177][0178][0179]
3.2设计新的控制器
[0180]
结合前两步，将经过卡尔曼滤波后的新的状态变量运用到预测滑模控制器当中，得到新的基于卡尔曼滤波器的预测滑模控制器，具体如下所示：
[0181][0182]
本发明从随机干扰产生的机理入手，在预测滑模控制的基础之上，来研究系统在白噪声干扰下的控制问题；并将外部随机干扰等效为系统的过程噪声，将传感器的测量误差等效为系统的观测噪声，提出一种基于卡尔曼滤波器的预测滑模控制方法，以解决预测滑模控制方法在白噪声干扰时，控制效果不佳的问题，从而实现白噪声干扰下对系统的有效控制。
[0183]
下面以一级倒立摆模型(如图4所示)为被控对象，验证本发明的有效性和稳定性。
[0184]
1选取预测模型
[0185]
选取标称系统为预测模型，设计预测滑模控制器，具体如下：
[0186]
x(k)＝ax(k-1)+bu(k-1)
[0187]
y(k)＝cx(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(37)
[0188][0189][0190]a32
＝-3(c
1-mga)/[a(4m+m)]
[0191]a42
＝-3(m+m)(c
1-mga)/[a2m(4m+m)]
[0192]
b3＝4/(4m+m)
[0193]
b4＝3/(4m+m)
[0194]
其中c1为弹性硬度，g＝9.8m/s2，m为小车质量，m为摆的质量，杆长为2a，x为小车位置，α为杆的摆角。
[0195]
2选取实际受扰系统
[0196]
模拟实际受扰系统，将外部随机干扰等效为系统的过程噪声，将传感器测量误差等效为观测噪声，在系统的输出端加上幅值为1的随机白噪声干扰信号，在输入端加上幅值为0.015的随机白噪声干扰信号，以此来设计卡尔曼滤波器，具体如下：
[0197]
x(k)＝ax(k-1)+bu(k-1)+w(k-1)
[0198]
y(k)＝cx(k)+v(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(38)
[0199]
其中w(k)为幅值为1的正态白噪声，v(k)为幅值为0.015的正态白噪声。
[0200]
3仿真结果分析
[0201]
为验证本发明的有效性，给出了该控制方法在matlab仿真平台上的仿真结果。参数设计如下：g＝9.8m/s2(重力加速度)，m＝5.0kg(小车质量)，m＝1.0kg(杆的质量)，a＝0.5m(杆的半长)，c1＝1(弹性硬度)，w(k)为幅值为1的正态白噪声，v(k)为幅值为0.015的正态白噪声。
[0202]
采样时ts＝0.02s，仿真时间30s，初始条件为：x(0)＝0.5，α(0)＝0.3，x
′
(0)＝0，α
′
(0)＝0。期望状态为：x＝0，α＝0，x
′
＝0，α
′
＝0。其中摆动角度单位为弧度。
[0203]
通过分析仿真结果，具体如图5～9所示。可知，虽然噪声幅值很小，但是噪声带来的干扰对系统的稳定性有着很大的影响。在经过卡尔曼滤波后，实际小车位移在[-2,2]cm之间，摆角在[-0.01,0.01]rad之间，系统控制量在[-0.5,0.5]m/s2之间，同时小车移动速度和杆的摆动角速度的幅度也缩小10倍之多，较之前未经滤波的结果均有大幅度改善，摆动范围很小，可认为是稳定的，从而验证了本发明的有效性。