一个非奇异快速终端滑模控制(NTSM)实例及仿真

一、被控对象

考虑这么一个被控对象
J θ ¨ ( t ) = u ( t ) + d ( t ) J \ddot\theta(t) = u(t) + d(t) Jθ¨(t)=u(t)+d(t)
其中, J J J 为转动惯量, θ \theta θ 为角度, u u u 为控制量, d d d 为扰动,且 d ( t ) < D d(t) < D d(t)<D 扰动有界

二、设计滑模面

2.1、传统快速Terminal滑模面

s = x ˙ + β x q / p = 0 s = \dot x + \beta x ^ {q/p} = 0 s=x˙+βxq/p=0

其中 x ∈ R 1 x \in R ^ 1 xR1 属于状态变量, β > 0 , p , q ( p > q ) \beta > 0,p,q(p>q) β>0,p,q(p>q) 为正奇数。由上式可得
d x d t = − β x q / p \frac{dx}{dt} = - \beta x ^ {q/p} dtdx=βxq/p

d t = − 1 β x q / p d x dt = - \frac{1}{\beta} x ^ {q/p}dx dt=β1xq/pdx

∫ 0 t d t = − ∫ x 0 0 1 β x q / p d x \int_{0}^{t} dt = -\int_{x_0}^{0}\frac{1}{\beta} x ^ {q/p}dx 0tdt=x00β1xq/pdx

从而得到从任意初始状态 x ( 0 ) ≠ 0 x(0) \ne 0 x(0)=0 沿滑动模态到达平衡状态 x = 0 x=0 x=0 的时间为
t s = p β ( p − q ) ∣ x ( 0 ) ∣ ( p − q ) / p t_s = \frac{p}{\beta(p-q)}|x(0)|^{(p-q)/p} ts=β(pq)px(0)(pq)/p

2.2、非奇异Terminal滑模面

为了解决普通Terminal滑模控制中的奇异性问题,提出了非奇异的Terminal滑模控制,
s = x 1 + 1 β x 2 p / q s = x_1 + \frac{1}{\beta} x_2 ^ {p/q} s=x1+β1x2p/q

一些思考

我们可以看到上面的滑模面,有的是 x ˙ \dot x x˙ x x x ,有的是 x 1 x_1 x1 x 2 x_2 x2 ,为什么呢?这里我有一些我的思考,其实他们都是一样的,他们的控制对象是下面这个不确定非线性动态系统,这个不确定非线性动态系统的状态 x 1 x_1 x1 就是 x x x x 2 x_2 x2 就是 x ˙ \dot x x˙
{ x ˙ 1 = x 2 x ˙ 2 = f ( x ) + g ( x ) u + d ( x , t ) \begin{cases} \dot x_1 = x_2\\ \dot x_2 = f(x) + g(x) u + d(x,t) \end{cases} { x˙1=x2x˙2=f(x)+g(x)u+d(x,t)
其中跟踪误差及其导数为, θ d \theta_d θd 为理想的角度信号
{ e ( t ) = θ ( t ) − θ d ( t ) e ˙ ( t ) = θ ˙ ( t ) − θ ˙ d ( t ) e ¨ ( t ) = θ ¨ ( t ) − θ ¨ d ( t ) \begin{cases} e(t) = \theta (t) - \theta_d(t) \\ \dot e(t) = \dot\theta (t) - \dot\theta_d(t) \\ \ddot e(t) = \ddot\theta (t) - \ddot\theta_d(t) \\ \end{cases} e(t)=θ(t)θd(t)e˙(t)=θ˙(t)θ˙d(t)e¨(t)=θ¨(t)θ¨d(t)

所以放在这个真实的系统中,状态 x 1 x_1 x1 就是 e ( t ) e(t) e(t) x 2 x_2 x2 就是 e ˙ ( t ) \dot e(t) e˙(t)

三、获得控制量

传统快速Terminal滑模的控制量为:
u = − g − 1 ( x ) ( f ( x ) + β q p x 1 q / p − 1 x 2 + ( D + ε ) s g n ( s ) ) u = - g^{-1}(x)(f(x) + \beta \frac{q}{p}x_1^{q/p-1}x_2 + (D + \varepsilon)sgn(s)) u=g1(x)(f(x)+βpqx1q/p1x2+(D+ε)sgn(s))
非奇异Terminal滑模的控制量为:
u = − g − 1 ( x ) ( f ( x ) + β q p x 2 2 − p / q x 2 + ( D + ε ) s g n ( s ) ) u = - g^{-1}(x)(f(x) + \beta \frac{q}{p}x_2^{2-p/q}x_2 + (D + \varepsilon)sgn(s)) u=g1(x)(f(x)+βpqx22p/qx2+(D+ε)sgn(s))

四、证明稳定性

4.1、传统快速Terminal滑模面

设LYapunov函数 V = 1 2 S 2 V=\frac{1}{2}S^2 V=21S2,则有
V ˙ = s s ˙ = s ( x ˙ 2 + β q p x 1 q p − 1 x ˙ 1 ) = s ( f ( x ) + g ( x ) u + d ( x , t ) + β q p x 1 q p − 1 x ˙ 1 ) = s ( d ( x , t ) − ( D + ε ) s g n ( s ) ) = s d ( x , t ) − ( D + ε ) ∣ s ∣ ≤ − ε ∣ s ∣ \begin{align} \dot V &= s \dot s \\ &= s(\dot x_2 + \beta\frac{q}{p}x_1^{\frac{q}{p}-1}\dot x_1) \\ &= s(f(x) + g(x)u + d(x,t) + \beta\frac{q}{p}x_1^{\frac{q}{p}-1}\dot x_1) \\ &= s(d(x,t)-(D+\varepsilon) sgn(s)) \\ &=sd(x,t) - (D+\varepsilon)|s| \\ &\le -\varepsilon|s| \end{align} V˙=ss˙=s(x˙2+βpqx1pq1x˙1)=s(f(x)+g(x)u+d(x,t)+βpqx1pq1x˙1)=s(d(x,t)(D+ε)sgn(s))=sd(x,t)(D+ε)sεs
可见,当干扰$d(t) $ 有上界时,即 d ( t ) < D d(t) < D d(t)<D 扰动有界,若 ε \varepsilon ε 取值合理,则系统稳定。

4.2、非奇异Terminal滑模面

设LYapunov函数 V = 1 2 S 2 V=\frac{1}{2}S^2 V=21S2,则有
V ˙ = s s ˙ = s ( x ˙ 1 + 1 β p q x 2 p q − 1 x ˙ 2 ) = s 1 β p q x 2 p q − 1 x ˙ 2 ( d ( x , t ) − ( D + ε ) s g n ( s ) ) = 1 β p q x 2 p q − 1 x ˙ 2 ( s d ( x , t ) − ( D + ε ) ∣ s ∣ ) \begin{align} \dot V &= s \dot s \\ &= s(\dot x_1 + \frac{1}{\beta}\frac{p}{q}x_2^{\frac{p}{q}-1}\dot x_2) \\ &= s \frac{1}{\beta}\frac{p}{q}x_2^{\frac{p}{q}-1}\dot x_2 (d(x,t)-(D+\varepsilon) sgn(s)) \\ &=\frac{1}{\beta}\frac{p}{q}x_2^{\frac{p}{q}-1}\dot x_2(sd(x,t) - (D+\varepsilon)|s|) \\ \end{align} V˙=ss˙=s(x˙1+β1qpx2qp1x˙2)=sβ1qpx2qp1x˙2(d(x,t)(D+ε)sgn(s))=β1qpx2qp1x˙2(sd(x,t)(D+ε)s)
由于 1 < p q < 2 1<\frac{p}{q}<2 1<qp<2 所以 0 < p q − 1 < 1 0 < \frac{p}{q} - 1 < 1 0<qp1<1 β > 0 \beta > 0 β>0 ,所以 x 2 p q − 1 > 0 ( x 2 ≠ 0 ) x_2 ^ {\frac{p}{q} - 1} > 0(x_2 \ne 0) x2qp1>0(x2=0) ,所以 V ˙ = − η ′ ∣ s ∣ ≤ 0 \dot V = -\eta ' |s| \le 0 V˙=ηs0 ,根据当 V ˙ ≡ 0 \dot V \equiv 0 V˙0 时, s ≡ 0 s \equiv 0 s0,LaSalle不变形原理, t → ∞ t \to \infty t 时, s → 0 s \to 0 s0 ,根据快速滑动模态特性, t → ∞ t \to \infty t x → 0 x \to 0 x0

可见,当 x 2 ≠ 0 x_2 \ne 0 x2=0 时,控制器满足LYapunov稳定条件。

五、实验

这里我们仅验证非奇异Terminal滑模面,其实其他的是一样的

需要注意一件事,在编程时, x 5 3 x ^ {\frac{5}{3}} x35 这种分数形式,不能直接写成

x ** (5 / 3)

这样写的话会出现虚数,需要写成下面的格式

np.abs(x) ** (5 / 3) * np.sign(x)

这里我的示例是python,但是其他的编程语言也一样,包括Matlab,C++等。

代码如下

import numpy as np
import matplotlib.pyplot as plt


class ControlModel():
    def __init__(self):
        self.T = 10
        self.times = 10 * 1000
        self.time_T = self.T / self.times

        self.pos = 0

        self.theta = np.zeros(self.times, dtype='float64')
        self.dot_theta = np.zeros(self.times, dtype='float64')
        self.u = np.zeros(self.times, dtype='float64')
        self.d = np.zeros(self.times, dtype='float64')
        self.J = 10

        self.Max_u = 30

    def reset(self):
        self.pos = 0
        self.dot_theta[0] = np.deg2rad(0)
        self.theta[0] = np.deg2rad(0)
        self.d[0] = 0

    def step(self, u):
        if self.pos >= self.times:
            return True

        if np.abs(u) > self.Max_u:
            u = np.sign(u) * self.Max_u

        self.u[self.pos] = u
        self.d[self.pos] = 2
        self.stepOnce()
        self.pos += 1

        return False

    def stepOnce(self):
        data = np.array([self.u[self.pos],
                         self.d[self.pos],
                         self.dot_theta[self.pos],
                         self.theta[self.pos],
                         self.J])

        k1 = self.time_T * self.iterateOnce(data)
        k2 = self.time_T * self.iterateOnce(data + 0.5 * k1)
        k3 = self.time_T * self.iterateOnce(data + 0.5 * k2)
        k4 = self.time_T * self.iterateOnce(data + k3)

        data = data + (k1 + 2 * k2 + 2 * k3 + k4) / 6

        self.dot_theta[self.pos + 1] = data[2]
        self.theta[self.pos + 1] = data[3]

    def iterateOnce(self, data):
        u = data[0]
        d = data[1]
        dot_theta = data[2]
        theta = data[3]
        J = data[4]

        _dot_theta = (u + d) / J
        _theta = dot_theta

        return np.array([0, 0, _dot_theta, _theta, 0])

    def get_theta(self):
        return self.theta[self.pos - 1]

    def get_dot_theta(self):
        return self.dot_theta[self.pos - 1]


class Control:
    def __init__(self):
        self.e = 0
        self.dot_e = 0
        self.ddot_e = 0
        self.T = 10
        self.times = 10 * 1000
        self.time_T = self.T / self.times

    def insert(self, e):
        self.ddot_e = (e - self.e - self.dot_e)
        self.dot_e = (e - self.e)
        self.e = e

    def get_e(self):
        return self.e

    def get_dot_e(self):
        return self.dot_e

    def get_ddot_e(self):
        return self.ddot_e


M = ControlModel()
M.reset()
C = Control()
T = Control()

theta_list = []
for i in range(M.times - 1):
    theta_d = np.sin(i / 1000)
    theta_list.append(theta_d)
    e = (M.get_theta() - theta_d) / M.time_T
    C.insert(e)
    T.insert(theta_d)

    c = 10
    varesplion = 3
    beta = 3
    q = 3
    p = 5
    T1 = np.abs(C.get_dot_e()) ** (p / q) * np.sign(C.get_dot_e())
    T2 = np.abs(C.get_dot_e()) ** (2 - p / q) * np.sign(C.get_dot_e())
    s = C.get_e() + 1 / beta * T1
    u = - M.J * (beta * q / p * T2 + varesplion * np.sign(s))

    M.step(u)

plt.figure(2)
plt.plot(M.theta)
plt.plot(theta_list)
plt.show()

plt.figure(3)
plt.plot(M.u)
plt.show()

跟踪效果如图所示

请添加图片描述

控制量

请添加图片描述

我们发现跟踪效果不错,但是控制量是存在抖振的,而这种抖振是需要我们在控制过程中去避免的。

猜你喜欢

转载自blog.csdn.net/weixin_43903639/article/details/130693828