目录
无人直升机的建模与控制
长期以来,直升机执行着诸如人员运输、货物运输、信息传输与监测等任务。在将这些任务转移到无人机的过程中,从数学建模的角度来看,无人直升机是一个过渡性技术。
无人直升机的分类
单旋翼直升机由四部分组成 —— 机身、尾翼系统、主旋翼和稳定器。同轴旋翼直升机也由四部分组成 —— 机身、下旋翼、上旋翼和稳定器。每架直升机都有一个称为斜盘的特殊动作机构。斜盘安装在直升机的主桅杆上,用于驱动旋翼桨叶俯仰角。
各种角度可参见 航向角的解释
单旋翼直升机的建模
单旋翼直升机由大量刚体组成。因此,首先需要导出每个刚体的模型;然后,将所有刚体的模型组合在一起,并推导出单旋翼直升机的完整数学模型。
伺服电机动力学 (Dynamics of Servomotors)
在单旋翼直升机中,有五个伺服电机用于驱动 —— 副翼伺服电机、升降舵伺服电机、集合俯仰伺服电机、方向舵伺服电机和油门伺服电机。伺服电机的传递公式如下:
G
s
(
s
)
=
ω
n
s
2
s
2
+
2
ζ
s
ω
n
s
s
+
ω
n
s
2
.
(1)
G_{s}(s) = \frac{\omega_{ns}^{2}}{s^{2} + 2 \zeta_{s} \omega_{ns} s + \omega_{ns}^{2}}. \tag{1}
Gs(s)=s2+2ζsωnss+ωns2ωns2.(1)
其中,
- ω n s \omega_{ns} ωns 是自然角频率;
- ζ s \zeta_{s} ζs 是阻尼系数。
俯仰与滚动 (Pitching and Rolling Motions)
当副翼和升降舵伺服电机的旋转角度改变时,主叶片的循环俯仰角也会改变。因此,转子盘旋转并产生陀螺力矩。俯仰和滚动力矩由主转子叶片的这些陀螺力矩产生。
我们假设陀螺运动对俯仰和滚动方向的影响是相同的;因此,由于机身惯性矩的不同,模型增益存在差异。在这种情况下,从伺服电机旋转角度到机身旋转角度的传递函数定义如下:
G
(
s
)
=
K
(
T
s
+
1
)
s
,
(2)
G(s) = \frac{K}{(T s + 1) s}, \tag{2}
G(s)=(Ts+1)sK,(2)
其中,
T
T
T 是系统的时间常数,
K
K
K 是模型增益。
将 (1) 与 (2) 结合,得到从脉冲输入到机身旋转角度的传递函数:
G
θ
(
s
)
=
ω
n
s
2
K
θ
(
s
2
+
2
ζ
s
ω
n
s
s
+
ω
n
s
2
)
(
T
θ
s
+
1
)
s
,
G
ϕ
(
s
)
=
ω
n
s
2
K
ϕ
(
s
2
+
2
ζ
s
ω
n
s
s
+
ω
n
s
2
)
(
T
ϕ
s
+
1
)
s
.
G_{\theta}(s) = \frac{\omega_{ns}^{2} K_{\theta}}{\left(s^{2} + 2 \zeta_{s} \omega_{ns} s + \omega_{ns}^{2}\right) (T_{\theta} s + 1) s}, \\ G_{\phi}(s) = \frac{\omega_{ns}^{2} K_{\phi}}{\left(s^{2} + 2 \zeta_{s} \omega_{ns} s + \omega_{ns}^{2}\right) (T_{\phi} s + 1) s}.
Gθ(s)=(s2+2ζsωnss+ωns2)(Tθs+1)sωns2Kθ,Gϕ(s)=(s2+2ζsωnss+ωns2)(Tϕs+1)sωns2Kϕ.
其中,
θ
\theta
θ 和
ϕ
\phi
ϕ 分别表示俯仰方向和滚动方向。进一步地,在实际设计控制器的过程中,考虑时滞因素,因而传递函数如下:
G
θ
(
s
)
=
e
−
L
s
ω
n
s
2
K
θ
(
s
2
+
2
ζ
s
ω
n
s
s
+
ω
n
s
2
)
(
T
θ
s
+
1
)
s
,
G
ϕ
(
s
)
=
e
−
L
s
ω
n
s
2
K
ϕ
(
s
2
+
2
ζ
s
ω
n
s
s
+
ω
n
s
2
)
(
T
ϕ
s
+
1
)
s
.
G_{\theta}(s) = e^{- L s} \frac{\omega_{ns}^{2} K_{\theta}}{\left(s^{2} + 2 \zeta_{s} \omega_{ns} s + \omega_{ns}^{2}\right) (T_{\theta} s + 1) s}, \\ G_{\phi}(s) = e^{- L s} \frac{\omega_{ns}^{2} K_{\phi}}{\left(s^{2} + 2 \zeta_{s} \omega_{ns} s + \omega_{ns}^{2}\right) (T_{\phi} s + 1) s}.
Gθ(s)=e−Ls(s2+2ζsωnss+ωns2)(Tθs+1)sωns2Kθ,Gϕ(s)=e−Ls(s2+2ζsωnss+ωns2)(Tϕs+1)sωns2Kϕ.
偏航运动 (Dynamics of Yawing Motion)
当方向舵伺服电机的旋转角度改变时,尾旋翼叶片俯仰角也会改变。因此,尾桨的推力也会发生变化。偏航力矩由尾桨推力产生,同时直升机沿偏航方向旋转。应该注意的是,小型直升机的偏航运动系统不是一个简单的开环系统。如果偏航系统是一个开环系统,直升机将通过主旋翼的反力矩旋转。然而,在这种情况下,控制直升机的航向是极其困难的。
出于上述原因,所有小型直升机都有一个称为速率陀螺的局部反馈系统。速率陀螺仪反馈由陀螺仪传感器测量的偏航速率,以抵消主转子的反力矩。下文中使用了一个主动速度控制系统 (AVCS) 陀螺,其中实现了比例积分 (PI) 控制系统。
假设带有 AVCS 陀螺仪、执行器和偏航动力学的闭环系统为二阶系统。此外,我们考虑上述的积分元素和时滞。那么,从方向舵脉冲输入到偏航方向旋转角
ψ
\psi
ψ 的传递函数如下:
G
ψ
(
s
)
=
e
−
L
s
ω
n
s
2
K
ψ
(
s
2
+
2
ζ
s
ω
n
s
s
+
ω
n
s
2
)
s
.
G_{\psi}(s) = e^{- L s} \frac{\omega_{ns}^{2} K_{\psi}}{\left(s^{2} + 2 \zeta_{s} \omega_{ns} s + \omega_{ns}^{2}\right) s}.
Gψ(s)=e−Ls(s2+2ζsωnss+ωns2)sωns2Kψ.
垂荡方向 (Dynamics of Heave Direction)
由于主旋翼推力的变化,直升机上下移动。根据 blade element theory,主转子产生的推力计算如下:
T
=
b
4
ρ
a
c
Ω
2
R
3
(
θ
t
+
ϕ
t
)
,
T = \frac{b}{4} \rho a c \Omega^{2} R^{3} \left(\theta_{t} + \phi_{t}\right),
T=4bρacΩ2R3(θt+ϕt),
此处,
- T T T 是主转子的推力;
- b b b 是转子叶片的数量;
- ρ \rho ρ 是空气密度;
- a a a 是二维升力曲线斜率;
- c c c 是叶片弦长;
- Ω \Omega Ω 是主转子的转速;
- R R R 是主转子的半径;
- θ t \theta_{t} θt 是主转子的总变桨角;
- ϕ t \phi_{t} ϕt 是流入角,定义为 ϕ t = v d + V z 2 R Ω , \phi_{t} = \frac{v_{d} + V_{z}}{2 R \Omega}, ϕt=2RΩvd+Vz, 其中 , v d v_{d} vd 是直升机的诱导速度, V z V_{z} Vz 是直升机的上升速度。
尽管此处参数众多,在实际的直升机飞行中可变的参数只有三个 Ω \Omega Ω, θ t \theta_{t} θt 和 ϕ t \phi_{t} ϕt 。
- 特别是在悬停飞行中, V z V_{z} Vz 为零, v d v_{d} vd 为常数;因此, ϕ t \phi_{t} ϕt 仅随 Ω \Omega Ω 变化。
- Ω \Omega Ω, θ t \theta_{t} θt 的变化分别与油门伺服电机和集合伺服电机的旋转角度成比例。此时,使用一种称为调速器的发动机转速控制器。调速器保持恒定的转子转速。
- 此外,假设 v d ≪ Ω v_{d} \ll \Omega vd≪Ω,因此, ϕ t ≈ 0 \phi_{t} \approx 0 ϕt≈0.
在这些假设下,可以说主转子推力
T
T
T 仅随总变桨角
θ
t
\theta_{t}
θt 而变化,即
T
=
K
t
θ
t
,
T = K_{t} \theta_{t},
T=Ktθt,
其中
K
t
K_{t}
Kt 是依赖于任何常参数的增益。
假设直升机自身重量为
M
M
M,则垂荡方向的运动方程为
M
V
˙
z
=
−
K
v
V
z
+
K
t
θ
t
,
M \dot{V}_{z} = - K_{v} V_{z} + K_{t} \theta_{t},
MV˙z=−KvVz+Ktθt,
新增加的项(即右侧第一项)是空气阻力,
K
v
K_{v}
Kv 是空气阻力的平均系数。进一步,可得到聚合脉冲到直升机高度的传递函数为
G
u
d
(
s
)
=
K
t
K
θ
t
(
M
s
+
K
v
)
s
,
G_{u d} (s) = \frac{K_{t} K_{\theta_{t}}}{\left(M s + K_{v}\right) s},
Gud(s)=(Ms+Kv)sKtKθt,
此处,
K
θ
t
K_{\theta_{t}}
Kθt 是聚合脉冲输入与悬停飞行中中主旋翼桨叶集体俯仰角的比率。
水平速度与位置 (Horizontal Velocity and Position)
主转子盘因俯仰方向的旋转运动而倾斜。据此,主转子推力分解为水平方向(向前)和垂直方向。此外,主旋翼推力的水平分量使直升机向前移动。现假设直升机保持静止或以恒定速度上下移动。主旋翼推力的垂直分量等于直升机的重量。记
- M M M 是直升机的质量;
- g g g 是重力引起的加速度;
- θ \theta θ 是直升机的俯仰旋转角度;
- V x V_{x} Vx 是直升机的前进速度。
假设
θ
\theta
θ 足够小,并获得前进方向的运动方程,
M
V
˙
x
=
−
M
g
tan
θ
≈
−
M
g
θ
.
(3)
M \dot{V}_{x} = - M g \tan \theta \approx - M g \theta. \tag{3}
MV˙x=−Mgtanθ≈−Mgθ.(3)
以类似的方式,记
ϕ
\phi
ϕ 作为直升机的旋转转角,
V
y
V_{y}
Vy 作为直升机的右旋速度;假设
ϕ
\phi
ϕ 足够小,右向运动方程为
M
V
˙
y
=
M
g
tan
ϕ
≈
M
g
ϕ
.
(4)
M \dot{V}_{y} = M g \tan \phi \approx M g \phi. \tag{4}
MV˙y=Mgtanϕ≈Mgϕ.(4)
对 (3) 和 (4) 做 Lapace 变换,得到
G
v
x
=
V
x
θ
=
−
g
s
,
G
v
y
=
V
y
ϕ
=
g
s
.
G_{v x} = \frac{V_{x}}{\theta} = - \frac{g}{s}, \\ G_{v y} = \frac{V_{y}}{\phi} = \frac{g}{s}.
Gvx=θVx=−sg,Gvy=ϕVy=sg.
还可以再次模型上进行修正,得到更符合需求的传递方程。
同轴旋翼直升机的建模
假设每个轴的耦合项可以忽略(与单旋翼直升机的情况相同)。此外,垂荡和航向运动的动力学与单旋翼直升机的动力学相似;因此,仅介绍伺服电机的动力学、横摇和俯仰运动的动力学以及水平运动的动力学。
伺服电机的动力学
此处的同轴旋翼直升机,无需驱动总俯仰角和尾桨俯仰角。因此,只有伺服电机用于驱动俯仰和滚动运动。伺服电机本身较小,其运动范围不大;因此,高阶动力学可以忽略。假设伺服电机的动力学可以是一个近似的一阶方程
G
m
(
s
)
=
K
m
T
m
s
+
1
,
G_{m} (s) = \frac{K_{m}}{T_{m} s + 1},
Gm(s)=Tms+1Km,
此处,
T
m
T_{m}
Tm 表示伺服电机的时间常数,
K
m
K_{m}
Km 表示适当的标量增益。
滚动和俯仰运动的动力学
对于同轴旋翼直升机,
- 副翼和升降舵伺服电机的旋转角度改变时,斜盘移动;下部转子的循环变桨角度发生变化。然后,下转子的旋转平面倾斜,并获得扑动方向。此处,伺服电机旋转角度和下转子倾斜角度之间的系统可近似为一阶传递函数,
G φ l ( s ) = K φ l T φ l s + 1 . G_{\varphi l} (s) = \frac{K_{\varphi l}}{T_{\varphi l} s + 1}. Gφl(s)=Tφls+1Kφl.
其中 φ l \varphi l φl 表示这是下转子倾斜角的传递函数。 - 力矩随后由下旋翼的倾斜产生,机身的旋转角度发生变化。然而,安装在车身顶部的稳定杆利用陀螺效应保持旋转平面。这意味着主桅杆和稳定杆之间的旋转角发生变化,该系统可近似为以下传递函数
G φ s ( s ) = K φ s ( 1 − 1 T φ s s + 1 ) . G_{\varphi s} (s) = K_{\varphi s} \left(1 - \frac{1}{T_{\varphi s} s + 1}\right). Gφs(s)=Kφs(1−Tφss+11).
其中 φ s \varphi s φs 表示这是稳定杆倾斜角度的传递函数。 - 接下来,当稳定杆的倾斜角改变时,连接到稳定杆的上转子的循环俯仰角改变,并且上转子的旋转平面倾斜。这里,稳定杆倾斜角度和上转子之间的传递函数可近似如下
G φ u ( s ) = K φ u T φ u s + 1 . G_{\varphi u} (s) = \frac{K_{\varphi u}}{T_{\varphi u} s + 1}. Gφu(s)=Tφus+1Kφu.
其中 φ l \varphi l φl 表示这是上转子倾斜角的传递函数。
现通过上述三个传递函数和伺服电机的传递函数的组合来描述表示滚动和俯仰运动动力学的传递函数。外部力矩对机身姿态的变化有重大影响,它是由上下转子倾角的变化产生的。现在,考虑
- T u T_{u} Tu 和 T l T_{l} Tl 分别为上、下转子产生的推力;
- L u L_{u} Lu 和 L l L_{l} Ll 分别为机身重心与上下转子之间的距离;
- φ u \varphi_{u} φu 和 φ l \varphi_{l} φl 分别是上下转子的倾斜角。
假设
φ
u
\varphi_{u}
φu 和
φ
l
\varphi_{l}
φl 足够小,则上下转子围绕机身重心分别产生的力矩
M
u
M_{u}
Mu 和
M
l
M_{l}
Ml 如下所示
M
u
=
T
u
L
u
sin
φ
u
≈
T
u
L
u
φ
u
,
M
l
=
T
l
L
l
sin
φ
l
≈
T
l
L
l
φ
l
.
M_{u} = T_{u} L_{u} \sin \varphi_{u} \approx T_{u} L_{u} \varphi_{u}, \\ M_{l} = T_{l} L_{l} \sin \varphi_{l} \approx T_{l} L_{l} \varphi_{l}.
Mu=TuLusinφu≈TuLuφu,Ml=TlLlsinφl≈TlLlφl.
假设
θ
\theta
θ 是俯仰或滚转姿态角,
J
J
J 表示各轴的惯性矩;关于机身姿态的运动方程的一般形式
J
θ
¨
=
M
u
+
M
l
.
J \ddot{\theta} = M_{u} + M_{l}.
Jθ¨=Mu+Ml.
水平运动的动力学
转子产生的推力的水平分量主要影响水平运动的动力学。现考虑由上、下转子产生的推力水平分量分别为
F
u
h
F_{u h}
Fuh 和
F
l
h
F_{l h}
Flh ;这些力与上下转子倾角之间的关系可表示为
F
u
h
=
T
u
sin
(
θ
+
φ
u
)
≈
T
u
(
θ
+
φ
u
)
,
F
l
h
=
T
l
sin
(
θ
+
φ
l
)
≈
T
l
(
θ
+
φ
l
)
.
F_{u h} = T_{u} \sin \left(\theta + \varphi_{u}\right) \approx T_{u} \left(\theta + \varphi_{u}\right), \\ F_{l h} = T_{l} \sin \left(\theta + \varphi_{l}\right) \approx T_{l} \left(\theta + \varphi_{l}\right).
Fuh=Tusin(θ+φu)≈Tu(θ+φu),Flh=Tlsin(θ+φl)≈Tl(θ+φl).
假设每个转子的姿态角和倾斜角足够小。
x
x
x 被认为是
X
X
X 轴和
Y
Y
Y 轴方向上的位移,
M
M
M 是物体的质量,运动方程的一般形式表示为
M
x
¨
=
−
D
x
x
˙
+
F
u
h
+
F
l
h
.
M \ddot{x} = - D_{x} \dot{x} + F_{u h} + F_{l h}.
Mx¨=−Dxx˙+Fuh+Flh.
小型无人直升机控制系统设计
基于之前的数学模型,进行控制系统的设计。具体而言,介绍了一种基于最优控制理论的方法,作为悬停等静止飞行控制典型设计方法的实例。
最优控制
在进行控制系统设计之前,有必要将前一节推导的传递函数转换为状态空间方程,因为最优控制设计需要状态空间方程。将传递函数转换为状态空间方程称为实现 (Realization)。一般来说,众所周知,传递函数存在多个实现。在下面的讨论中,希望系统是可控的和可观察的;可控和可观的实现称为最小实现。
基于之前的传递函数的实现为
x
˙
=
A
x
(
t
)
+
B
u
(
t
)
,
x
∈
R
n
,
u
∈
R
,
y
=
C
x
(
t
)
,
y
∈
R
,
\dot{\mathbf x} = {\mathbf A}{\mathbf x}(t) + {\mathbf B}{u}(t), \quad {\mathbf x} \in \mathbb{R}^{n}, u \in \mathbb{R}, \\ y = {\mathbf C}{\mathbf x}(t), \quad y \in \mathbb{R},
x˙=Ax(t)+Bu(t),x∈Rn,u∈R,y=Cx(t),y∈R,
其中,
(
A
,
B
)
({\mathbf A}, {\mathbf B})
(A,B) 是可控的,
(
A
,
C
)
({\mathbf A}, {\mathbf C})
(A,C) 是可观的。
考虑一个恒定的参考值
r
r
r,则将参考与误差之间的积分值定义为
e
(
t
)
=
∫
0
t
r
−
y
(
τ
)
d
τ
.
e(t) = \int_{0}^{t} r - y(\tau) \; d \tau.
e(t)=∫0tr−y(τ)dτ.
定义
x
a
=
[
x
T
,
e
]
T
{\mathbf x}_{a} = [\mathbf{x}^{\rm T}, e]^{\rm T}
xa=[xT,e]T,则增广的状态空间方程为
x
˙
a
=
A
a
x
a
+
B
a
u
+
E
r
.
\dot{\mathbf x}_{a} = {\mathbf A}_{a} {\mathbf x}_{a} + {\mathbf B}_{a} {u} + {\mathbf E} r.
x˙a=Aaxa+Bau+Er.
现基于增广的状态空间方程设计最优控制。最优控制是控制系统设计方法的一个目标。在该方法中,通过使用能够最小化某些准则的状态反馈来计算最优控制输入。
现考虑如下所示的二次准则
J
=
∫
0
∞
x
a
T
(
t
)
Q
x
a
(
t
)
+
R
u
2
(
t
)
d
t
,
(5)
J = \int_{0}^{\infty} {\mathbf x}_{a}^{\rm T}(t) {\mathbf Q} {\mathbf x}_{a}(t) + R u^{2}(t) \; dt, \tag{5}
J=∫0∞xaT(t)Qxa(t)+Ru2(t)dt,(5)
控制目标的加权矩阵
Q
{\mathbf Q}
Q 是一个非负定对称度量矩阵,控制输入的权重
R
R
R 是正标量。最小化 (5) 的最优控制输入为
u
(
t
)
=
−
F
x
a
(
t
)
,
F
=
R
−
1
B
a
T
P
.
u(t) = - {\mathbf F} {\mathbf x}_{a}(t), \quad {\mathbf F} = R^{- 1} {\mathbf B}_{a}^{\rm T}{\mathbf P}.
u(t)=−Fxa(t),F=R−1BaTP.
此处的
P
{\mathbf P}
P 是 Riccati 方程的解,且为正定对称矩阵。