PSO优化滑模控制器参数
PSO优化滑模控制器参数
1.粒子群算法
1.1粒子群算法原理
粒子群优化算法(Particle Swarm Optimization,简称PSO), 是1995年Eberhart博士和Kennedy博士一起提出的,它是源于对鸟群捕食行为的研究。粒子群优化算法的基本核心是利用群体中的个体对信息的共享从而使得整个群体的运动在问题求解空间中产生从无序到有序的演化过程,从而获得问题的最优解。
简单来说,PSO就是通过初始化一群随机粒子,然后进行迭代来逐渐靠近最优解,最终收敛至最优解。算法的流程图和伪代码见下图。
算法中最为关键的步骤是每轮迭代时更新粒子的速度和位置。
其中速度更新的计算主要依赖于上一时刻的速度值、当前位置与粒子历史最优位置的差值和当前位置与粒子群历史最优位置的差值,其具体计算公式如下:
v
i
(
t
+
1
)
=
ω
v
i
(
t
)
+
c
1
r
1
(
p
b
e
s
t
i
−
x
i
(
t
)
)
+
c
2
r
2
(
g
b
e
s
t
−
x
i
(
t
)
)
\large v_{i}(t+1) = \omega v_{i}(t) + c_{1} r_{1}(pbest_{i} - x_{i}(t)) + c_{2} r_{2}(gbest - x_{i}(t))
vi(t+1)=ωvi(t)+c1r1(pbesti−xi(t))+c2r2(gbest−xi(t))
其中
v
i
(
t
)
v_{i}(t)
vi(t)代表第
i
i
i个粒子在
t
t
t时刻的速度值;
x
i
(
t
)
x_{i}(t)
xi(t)代表第
i
i
i个粒子在
t
t
t时刻的位置;
p
b
e
s
t
i
pbest_{i}
pbesti代表第
i
i
i个粒子的历史最优位置;
g
b
e
s
t
gbest
gbest代表第整个粒子群的历史最优位置。系数
ω
\omega
ω是速度惯性权重,用于调节对解空间的搜索范围;
c
1
c_{1}
c1是自我学习因子,
c
2
c_{2}
c2是群体学习因子,共同用于调节最大步长;
r
1
r_{1}
r1和
r
2
r_{2}
r2是两个随机数,范围介于0到1之间,主要用于增加搜索的随机性。
而位置的更新则较为简单,只需要上一时刻的位置值和速度值,计算公式如下:
x
i
(
t
+
1
)
=
x
i
(
t
)
+
v
i
\large x_{i}(t+1) = x_{i}(t) + v_{i}
xi(t+1)=xi(t)+vi
1.2粒子群算法求函数最值
基于上述原理,将粒子群算法用于求解函数的最值。以自变量的值作为粒子位置值,以函数值作为种群适应度,在matlab中编写代码粒子群算法,实现寻找函数极值点与极值的功能,
为测试代码,选取测试函数如下:
f
(
x
)
=
x
s
i
n
x
c
o
s
2
x
−
2
x
s
i
n
3
x
f(x) = xsinxcos2x-2xsin3x
f(x)=xsinxcos2x−2xsin3x
对上述函数求解
[
0
,
20
]
[0,20]
[0,20]范围内的最小值,其结果如下。其中下图1为函数曲线,红点为PSO寻优寻找到的最小值标记点;下图2为PSO寻优过程中的收敛曲线。
最终得到函数的最值点为 x 0 = 19.4113 x_{0}=19.4113 x0=19.4113,函数的最值 f ( x 0 ) = − 30.0963 f(x_{0})=-30.0963 f(x0)=−30.0963。测试PSO寻优最值成功。
2.滑模变结构控制
2.1被控对象
选取简单的电机作为被控对象,其数学模型如下:
θ
¨
(
t
)
=
−
f
(
θ
,
t
)
+
b
u
(
t
)
\ddot \theta (t) = -f(\theta,t) + b u(t) \\
θ¨(t)=−f(θ,t)+bu(t)
其中
θ
(
t
)
\theta(t)
θ(t)为电机转动的角度;
u
(
t
)
u(t)
u(t)为输入给电机的控制量。
f
f
f为模型中的非线性项,其表达式为
f
(
θ
,
t
)
=
25
θ
˙
f(\theta,t) = 25 \dot \theta
f(θ,t)=25θ˙。
2.2滑模控制器设计
对于被控电机,设其期望跟踪角度为 θ d ( t ) = s i n t \theta_d (t) = sint θd(t)=sint。现设计滑模控制器实现闭环控制,使得电机实际角度 θ ( t ) \theta (t) θ(t)的变化可以跟踪上期望角度$\theta_{d} $的变化。
定义系统的跟踪误差为
e
(
t
)
=
θ
d
−
θ
e(t) = \theta_{d} - \theta
e(t)=θd−θ,系统跟踪误差的微分为
e
˙
=
θ
˙
d
−
θ
˙
\dot e = \dot \theta_{d} - \dot \theta
e˙=θ˙d−θ˙。取线性滑模面
s
=
c
e
+
e
˙
s=ce+\dot e
s=ce+e˙,并且采用等速+指数趋近率,最后设计得到滑模控制器表达式如下:
u
(
t
)
=
1
b
(
ϵ
s
g
n
s
+
k
s
+
c
e
˙
+
θ
¨
d
+
f
(
θ
,
t
)
)
\large u(t) = \frac{1}{b}(\epsilon sgns + ks + c \dot e + \ddot \theta_{d} + f(\theta,t))
u(t)=b1(ϵsgns+ks+ce˙+θ¨d+f(θ,t))
控制器设计过程与稳定性证明见附录。
2.3实验仿真
在simulink中搭建相应的仿真模型如下:
部分重要参数选取如下:
名称 | 符号 | 数值 |
---|---|---|
滑模面系数 | c c c | 15 |
指数趋近率 | k k k | 10 |
等速趋近率 | ϵ \epsilon ϵ | 5 |
输入系数 | b b b | 133 |
仿真结果见下图,可见系统输出在
t
=
0.568
s
t = 0.568s
t=0.568s时跟踪上期望输入。
3.粒子群算法优化滑模控制参数
3.1 优化目标
借鉴最优控制的思想:以一控制系统作为研究对象,考虑其中某些参数可变或者不确定的情况,寻找一组最优数,使得控制系统满足某些特殊条件或约束。
本文以上述控制系统作为对象,为滑模面系数 c c c和等速趋近率 ϵ \epsilon ϵ寻找一组最优参数,使得系统的能量消耗和跟踪误差最少。
基于上述要求,以系统跟踪误差
e
e
e和系统控制量
u
t
ut
ut作为指标,提出了优化控制的目标函数如下:
J
=
∫
0
∞
a
∣
e
(
t
)
∣
+
b
∣
u
t
(
t
)
∣
d
t
\large J = \int ^{\infty} _{0} a|e(t)| +b| ut(t)| \space dt
J=∫0∞a∣e(t)∣+b∣ut(t)∣ dt
其中系数
a
a
a和
b
b
b分别为跟踪误差与控制量的权重,代表了对能量消耗要求和误差跟踪要求的重视程度。
3.2 实验仿真
本文实验平台为matlab2020b版本。实验思路为利用simulink仿真计算系统的能量消耗和跟踪误差值,再导入workspace,通过matlab计算目标函数值并进行粒子群寻优,最终获得最优参数解。
由于数值计算具有离散型与有限性,这与上述目标函数相悖,因此在实际实验与计算时,将目标函数的无限积分更改为有限累加,更改后的目标函数如下:
J
=
∑
i
=
s
t
a
r
t
i
=
e
n
d
(
a
∣
e
(
i
)
∣
∗
T
s
t
e
p
+
b
∣
u
t
(
i
)
∣
∗
T
s
t
e
p
)
J = \sum _{i=start} ^{i=end} (a\space |e(i)| * T_{step} \space + b \space |ut(i)| * T_{step} \space)
J=i=start∑i=end(a ∣e(i)∣∗Tstep +b ∣ut(i)∣∗Tstep )
其中
T
s
t
e
p
T_{step}
Tstep为仿真步长,本次实验统一使用固定步长仿真,取
T
s
t
e
p
=
0.001
s
T_{step}=0.001s
Tstep=0.001s。
3.3 实验结果
以下为粒子群位置分布图,其中下图1为粒子群初始时的位置,其满足随机分布;经过迭代后,得到粒子群末态分布位置图,粒子基本分布与最优参数解 c = 2.2007 , ϵ = 8.8920 c=2.2007,\epsilon=8.8920 c=2.2007,ϵ=8.8920的附近,此时目标函数取最小值 J = 1.1980 J=1.1980 J=1.1980。
再统计历代粒子最优解可得进化过程曲线,如下图。并且发现当迭代至
N
=
10
N=10
N=10代时,即可取得最优值。
附录
[A]滑模控制器设计与稳定性证明
已知系统跟踪误差为
e
(
t
)
=
θ
d
−
θ
e(t) = \theta_{d} - \theta
e(t)=θd−θ,系统跟踪误差的微分为
e
˙
=
θ
˙
d
−
θ
˙
\dot e = \dot \theta_{d} - \dot \theta
e˙=θ˙d−θ˙,对
s
=
c
e
+
e
˙
s=ce+\dot e
s=ce+e˙取微分得
s
˙
=
c
e
˙
+
e
¨
=
c
e
˙
+
θ
¨
d
−
θ
¨
=
c
e
˙
+
θ
¨
d
+
f
(
θ
,
t
)
−
b
u
(
t
)
\dot s = c \dot e + \ddot e \\ =c \dot e + \ddot \theta_{d} - \ddot \theta \\ =c \dot e + \ddot \theta_{d} + f(\theta,t) - bu(t) \\
s˙=ce˙+e¨=ce˙+θ¨d−θ¨=ce˙+θ¨d+f(θ,t)−bu(t)
令其等于指数+等速趋近率
s
˙
=
−
k
s
g
n
(
s
)
−
ϵ
\dot s = -ksgn(s) - \epsilon
s˙=−ksgn(s)−ϵ推到得控制器表达式。
s
˙
=
c
e
˙
+
θ
¨
d
+
f
(
θ
,
t
)
−
b
u
(
t
)
=
−
k
s
g
n
(
s
)
−
ϵ
⇓
u
(
t
)
=
1
b
(
ϵ
s
g
n
s
+
k
s
+
c
e
˙
+
θ
¨
d
+
f
(
θ
,
t
)
)
\dot s = c \dot e + \ddot \theta_{d} + f(\theta,t) - bu(t) = -ksgn(s) - \epsilon \\ \Downarrow \\ u(t) = \frac{1}{b}(\epsilon sgns + k s + c \dot e + \ddot \theta_{d} + f(\theta,t))
s˙=ce˙+θ¨d+f(θ,t)−bu(t)=−ksgn(s)−ϵ⇓u(t)=b1(ϵsgns+ks+ce˙+θ¨d+f(θ,t))
利用李雅普诺夫稳定理论证明稳定性,构造备选李雅普诺夫函数
V
=
1
2
s
2
V = \frac{1}{2}s^2
V=21s2,对其求微分:
V
˙
=
s
s
˙
=
s
∗
(
−
ϵ
s
g
n
s
−
k
s
)
=
−
ϵ
s
−
k
s
2
=
−
ϵ
s
−
2
k
V
\dot V = s \dot s \\ = s * (- \epsilon sgns - k s) \\ = - \epsilon s - ks^2 \\ = - \epsilon s - 2k V \\
V˙=ss˙=s∗(−ϵsgns−ks)=−ϵs−ks2=−ϵs−2kV
再利用放缩可得:
V
˙
≦
−
2
k
V
⇒
V
(
t
)
≤
e
−
2
k
(
t
−
t
0
)
V
(
t
0
)
\dot V \leqq -2k V \space \Rightarrow \space V(t) \leq e^{-2k(t-t_{0})} V(t_{0})
V˙≦−2kV ⇒ V(t)≤e−2k(t−t0)V(t0)
证明
V
(
t
)
V(t)
V(t)以指数形式收敛至0,证得系统稳定。
[B]仿真代码
源代码已上传github,链接如下:粒子群算法优化滑模控制