CN103760905B - 基于模糊前馈单旋翼无人直升机姿态非线性鲁棒控制方法 - Google Patents

Abstract

本发明属于微小型旋翼式无人飞行器自主飞行控制研究领域，为提供一种无人机控制方法，使无人直升机实现快速、准确的镇定控制，并且该控制器对模型先验知识依赖较低，对系统的不确定性具有良好的鲁棒性，为此，本发明采用的技术方案是，基于模糊前馈单旋翼无人直升机姿态非线性鲁棒控制方法，包括如下步骤：一、单旋翼无人直升机动力学模型单旋翼无人直升机姿态动力学模型形式如下：<maths num="0001"></maths>二、基于模糊前馈的非线性半连续鲁棒控制器设计。本发明主要应用于微小型旋翼式无人飞行器的设计制造。

Description

基于模糊前馈单旋翼无人直升机姿态非线性鲁棒控制方法

技术领域

本发明属于微小型旋翼式无人飞行器自主飞行控制研究领域，主要针对一种单旋翼无人飞行器的控制算法设计，包括无人直升机姿态动力学模型介绍、基于模糊前馈非线性半连续鲁棒控制器设计以及姿态飞行控制实验，具体讲，涉及基于模糊前馈的单旋翼无人直升机姿态非线性鲁棒控制方法。

背景技术

单旋翼无人直升机是旋翼式飞行器的一种，一般搭载自主飞行控制器、传感器以及无线通信系统等，能够实现自主飞行以及自主执行多任务。单旋翼无人飞行器能够实现垂直起降，在空中能够实现悬停、前飞、后飞、侧飞等飞行动作，具有很强的机动性和灵活性。因此，单旋翼无人直升机近年来被大量应用于军事和民用领域，受到了大批研究人员的关注。但单旋翼无人直升机有极强的静不稳定性，且在动力学上具备欠驱动、强耦合与非线性等特点，使得单旋翼无人直升机飞行控制器设计难度大大提升。

目前国内已经有很多高校和科研机构都在进行单旋翼无人直升机方面的研究。如南京航空航天大学研究小组，针对单旋翼无人机非线性非仿射系统，设计了一种模糊自适应控制器以及输出反馈控制器，实现了系统的输出反馈控制。通过无人机飞行控制数值仿真验证了算法有效性，并没有进行相应实验验证(期刊：南京航空航天大学学报；著者：常勇，卢广山，姜长生；出版年月：2013年；文章题目：无人机非线性非仿射飞控系统的自适应模糊H∞输出反馈控制及其应用；页码：第45卷第1期99-103)。又如北京航空航天大学的研究小组，根据已辨识模型，提出了一种改进的无人直升机H∞鲁棒控制器设计方法，使系统具有更宽的鲁棒稳定裕度。使用该鲁棒控制器进行了数值仿真验证，并没有进行相应的实验验证(期刊：航空学报；著者：刘鹏，王强，蒙志君；出版年月：2012年；文章题目：基于飞行品质评估的无人直升机鲁棒控制器设计；页码：第33卷第9期1587-1597)。

另一方面，国外研究人员在单旋翼无人直升机的控制方面也取得了一定的成果。如澳大利亚国防学院的研究小组，针对无人直升机真实飞行的阵风干扰，基于无人直升机线性模型，利用反步法设计了滚转、俯仰通道的鲁棒控制器，并通过与LQR控制器对比验证了控制器有效性。论文仅完成了数值仿真验证，并未进行相应实际飞行实验(会议：the31stChineseControlConference；著者：RoyTK，PotaHR，GarrattM；出版年月：2012年；文章题目：Robustcontrolforlongitudinalandlateraldynamicsofsmallscalehelicopter；页码：2607-2612)。韩国建国大学的研究小组，针对小型单旋翼无人直升机姿态，利用加权函数以及H∞控制方法，设计了无人机姿态鲁棒控制器，使用无人机非线性模型进行了数值仿真，验证了控制器抗风性能，并未进行相关飞行实验验证(期刊：JournalofAerospaceEngineering；著者：JeongDY，KangT，DharmayandaHR；出版年月：2011年；文章题目：H-InfinityAttitudeControlSystemDesignforaSmall-ScaleAutonomousHelicopterwithNonlinearDynamicsandUncertainties；页码：第25卷第4期501-518)。

上述研究机构都针对单旋翼无人直升机姿态控制提出了较好的解决方案，对于控制器抗干扰能力及鲁棒性进行了改进。但是由于部分控制器设计对模型依赖较强并且算法复杂度较高，大多停留在数值仿真实验中，在实际飞行实验中的控制精度与鲁棒性仍然未知。

发明内容

本发明旨在解决克服现有技术的不足，为提供一种无人机控制方法，使无人直升机实现快速、准确的镇定控制，并且该控制器对模型先验知识依赖较低，对系统的不确定性具有良好的鲁棒性，为此，本发明采用的技术方案是，基于模糊前馈单旋翼无人直升机姿态非线性鲁棒控制方法，包括如下步骤：

一、单旋翼无人直升机动力学模型

单旋翼无人直升机姿态动力学模型形式如下：

$M\left(x\right)\stackrel{\&CenterDot;\&CenterDot;}{x}+C(x,\stackrel{\&CenterDot;}{x})\stackrel{\&CenterDot;}{x}+G\left(x\right)=\mathrm{\&tau;}---\left(1\right)$

其中，x＝[φθψ]^T为列向量，φ、θ、ψ分别表示滚转角、俯仰角、偏航角，表示对x的一阶导，表示对x的二阶导，M(x)∈R^3×3表示惯性矩阵，表示科里奥利力矩阵，G(x)∈R³表示保守力向量，τ∈R³表示转矩输入向量，其中R表示全体实数；

式(1)具有如下三个性质：

性质1.惯性矩阵M(x)是对称，正定的，并满足下面的不等式

$m_1\left|\right|\mathrm{\&xi;}|{|}^2\&le;{\mathrm{\&xi;}}^{T}M\left(x\right)\mathrm{\&xi;}\&le;m_2|\left|\mathrm{\&xi;}\right|{|}^2,\&ForAll;\mathrm{\&xi;}\&Element;R^3---\left(2\right)$

其中m₁，m₂是正常数；

性质2.科里奥利力和保守力满足如下不等式

其中，是有界正常数；

性质3.科里奥利力矩阵满足下面的关系

$C(x,\mathrm{\&xi;})\mathrm{\&upsi;}=C(x,\mathrm{\&upsi;})\mathrm{\&xi;},\&ForAll;\mathrm{\&xi;},\mathrm{\&upsi;}\&Element;R^3---\left(4\right)$

式(1)中，转矩输入τ可通过纵向挥舞角a、横向挥舞角b和尾桨推力T_T表示，即：

τ＝S^-T(A(T_M)υ_c+B(T_M)),(5)

其中υ_c＝[abT_T]^T，T_M∈R表示主桨推力，A(T_M)∈R^3×3是可逆矩阵，B(T_M)∈R³是主桨推力的分力向量,S表示从体坐标系到惯性坐标系的平移矩阵，其表达式如下：

$S=\left[\begin{array}{ccc}1& \frac{s_{\mathrm{\&phi;}}{s}_{\mathrm{\&theta;}}}{c_{\mathrm{\&theta;}}}& \frac{c_{\mathrm{\&phi;}}{s}_{\mathrm{\&theta;}}}{c_{\mathrm{\&theta;}}}\\ 0& c_{\mathrm{\&phi;}}& -s_{\mathrm{\&phi;}}\\ 0& \frac{s_{\mathrm{\&phi;}}}{c_{\mathrm{\&theta;}}}& \frac{c_{\mathrm{\&phi;}}}{c_{\mathrm{\&theta;}}}\end{array}\right]---\left(6\right)$

其中，s_φ表示sinφ，s_θ表示sinθ，c_φ表示cosφ，c_θ表示cosθ；

单旋翼无人直升机旋翼动力学模型形式如下：

$\begin{array}{c}\stackrel{\&CenterDot;}{a}=-\frac{{\mathrm{\&tau;}}_{mr}+K_{sb}{\mathrm{\&tau;}}_{sb}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}q-\frac{1}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}a+\frac{{\mathrm{\&tau;}}_{mr}{A}_b}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}b+\frac{A_{lon}+K_{sb}{C}_{lon}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}{\mathrm{\&delta;}}_{lon}\\ \stackrel{\&CenterDot;}{b}=-\frac{{\mathrm{\&tau;}}_{mr}+K_{sb}{\mathrm{\&tau;}}_{sb}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}p+\frac{{\mathrm{\&tau;}}_{mr}{B}_a}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}a-\frac{1}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}b+\frac{B_{lat}+K_{sb}{D}_{lat}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}{\mathrm{\&delta;}}_{lat}\\ {\stackrel{\&CenterDot;}{T}}_T=B_{ped}+K_{ped}{\mathrm{\&delta;}}_{ped},\end{array}---\left(7\right)$

其中，a表示纵向挥舞角，表示a的一阶导数，b表示横向挥舞角，表示b的一阶导数，T_T表示尾桨推力，表示T_T的一阶导数，p表示俯仰角速度，q表示滚转角速度，δ_lon表示控制输入纵向周期变矩，δ_lat表示控制输入横向周期变矩，δ_ped表示控制输入尾桨矩，τ_mr表示主旋翼挥舞时间常数，τ_sb表示副翼挥舞时间常数，A_b表示主旋翼纵向伺服输入比例系数，B_a表示主旋翼横向伺服输入比例系数，C_lon表示副翼纵向伺服输入比例系数，D_lat表示副翼横向伺服输入比例系数，K_sb表示主旋翼与副翼伺服输入比值，B_ped表示尾桨输入常数，K_ped表示尾桨伺服输入比例系数,A_lon表示主旋翼纵向指令输入比例系数，B_lat表示主旋翼横向指令输入比例系数；

在悬停状态下，对上述模型进行简化，得到如下形式：

$\begin{array}{c}a=A_{b}b-A_{lon0}{\mathrm{\&delta;}}_{lon}\\ b=-B_{a}a+B_{lat0}{\mathrm{\&delta;}}_{lat}\\ T_T=K_{ped0}{\mathrm{\&delta;}}_{ped},\end{array}---\left(8\right)$

其中，A_lon0、B_lat0和K_ped0为简化后常数，将式(8)带入式(5)，可得简化的单旋翼无人直升机旋翼动力学模型形式如下：

τ＝S^-T(A(T_M)Cδ+B(T_M)),(9)

其中δ＝[δ_lonδ_latδ_ped]^T表示控制输入，常数阵C∈R^3×3定义如下：

$C=\left[\begin{array}{ccc}-\frac{A_{lon0}}{A_b{B}_a+1}& \frac{A_b{B}_{lat0}}{A_b{B}_a+1}& 0\\ \frac{B_{lat0}}{A_b{B}_a+1}& \frac{B_a{A}_{lon0}}{A_b{B}_a+1}& 0\\ 0& 0& K_{ped0}\end{array}\right]---\left(10\right)$

二、基于模糊前馈的非线性半连续鲁棒控制器设计

为实现控制目标，首先定义姿态跟踪误差：

e₁＝x_d-x,(11)

其中x_d∈R³表示期望姿态向量，x∈R³表示实际姿态向量；

定义滤波误差e₂和r，表达式如下：

$e_2={\stackrel{\&CenterDot;}{e}}_1+e_1---\left(12\right)$

$r={\stackrel{\&CenterDot;}{e}}_2+{\mathrm{\&alpha;e}}_2---\left(13\right)$

其中，α∈R^3×3是正定对角常数阵，e₂＝[e₂₁e₂₂e₂₃]^T，e₂₁，e₂₂和e₂₃表示e₂向量在滚转、俯仰、偏航三个方向的分量，r＝[r₁r₂r₃]^T，r₁、r₂和r₃表示r向量在滚转、俯仰、偏航三个方向的分量，表示e₁的一阶导数，表示e₂的一阶导数；

对式(13)求导，等式的两边同时左乘M(x)，然后带入式(1)，可得

$M\left(x\right)\stackrel{\&CenterDot;}{r}=-\frac{1}{2}\stackrel{\&CenterDot;}{M}\left(x\right)r-e_2+N-S^{-T}AC\stackrel{\&CenterDot;}{\mathrm{\&delta;}},---\left(14\right)$

其中，A同公式(5)中A(T_M)，表示控制输入δ的一阶导数，辅助函数N定义如下：

$\begin{array}{c}N=M\left(x\right)({\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{x}}_d+{\stackrel{\&CenterDot;\&CenterDot;}{e}}_1+\mathrm{\&alpha;}{\stackrel{\&CenterDot;}{e}}_2)+\stackrel{\&CenterDot;}{M}\left(x\right)(\stackrel{\&CenterDot;\&CenterDot;}{x}+\frac{1}{2}r)+e_2+\stackrel{\&CenterDot;}{C}\stackrel{\&CenterDot;}{x}+C\stackrel{\&CenterDot;\&CenterDot;}{x}+\stackrel{\&CenterDot;}{G}-\\ {\stackrel{\&CenterDot;}{S}}^{-T}{S}^T\left(M\right(x)\stackrel{\&CenterDot;\&CenterDot;}{x}+C\stackrel{\&CenterDot;}{x}+G-S^{-T}B)\end{array}---\left(15\right)$

其中，B同公式(5)中B(T_M)，G同公式(1)中G(x)，引入有||N_d||，L_∞表示无穷范数，并且，N_d可以通过模糊推理方法来逼近，模糊输出：

$\hat{f}={\widehat{\mathrm{\&omega;}}}^T\mathrm{\&sigma;}(q_d,{\stackrel{\&CenterDot;}{q}}_d,{\stackrel{\&CenterDot;\&CenterDot;}{q}}_d,{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{q}}_d),---\left(16\right)$

其中为模糊推理权值，σ∈R^27×1，q_d＝[φ_d,θ_d,ψ_d]^T，φ_d，θ_d，ψ_d分别表示滚转角、俯仰角、偏航角参考；模糊系统输入共有12个变量，分别来自3个通道，每个通道有4个状态，由于输入变量过多，通过如下变换，即可使用3个变量用于模糊推理；

$\begin{array}{c}{\stackrel{\&OverBar;}{q}}_1={\mathrm{\&phi;}}_d+{\mathrm{\&theta;}}_d+{\mathrm{\&psi;}}_d+{\stackrel{\&CenterDot;}{\mathrm{\&phi;}}}_d+{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_d+{\stackrel{\&CenterDot;}{\mathrm{\&psi;}}}_d\\ {\stackrel{\&OverBar;}{q}}_2={\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&psi;}}}_d\\ {\stackrel{\&OverBar;}{q}}_3={\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&psi;}}}_d\end{array}---\left(17\right)$

在模糊推理中，隶属度函数选择高斯函数，表达式如下：

${\mathrm{\&mu;}}_{\stackrel{\&OverBar;}{q}}=e^{\frac{-0.5\&times;{(\stackrel{\&OverBar;}{q}-c)}^2}{{\mathrm{\&omega;}}^2}},---\left(18\right)$

其中，c和ω为高斯函数中实数，高斯函数将状态输入分为3个区间：正、零、负；模糊推理选择乘积运算，由于输入包含3个变量，每个变量分为3个区间，故可产生27条规则；解模糊方法选择中心平均法，表达式如下：

${\mathrm{\&sigma;}}_i=\frac{{\mathrm{\&mu;}}_{F_1^i\&times;F_2^i\&times;F_3^i}}{\underset{i=1}{\overset{p}{\&Sigma;}}{\mathrm{\&mu;}}_{F_1^i\&times;F_2^i\&times;F_3^i}},---\left(19\right)$

其中p＝27，表示第i条规则时模糊推理乘积运算结果；

设计自适应更新率：

$\stackrel{\&CenterDot;}{\widehat{\mathrm{\&omega;}}}\left(t\right)=-k\widehat{\mathrm{\&omega;}}\left(t\right)+\mathrm{\&tau;}\mathrm{\&sigma;}{\&lsqb;sat(e_{2\left(t\right)}+p_{2\left(t\right)})\&rsqb;}^T$

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{\widehat{\mathrm{\&omega;}}}=-k\widehat{\mathrm{\&omega;}}+\mathrm{\&tau;}\mathrm{\&sigma;}{\&lsqb;sat(e_{2\left(t\right)}+p_{2\left(t\right)})\&rsqb;}^T\\ p_{2\left(t\right)}=\frac{1}{\mathrm{\&epsiv;}}{e}_{2\left(t\right)}-\frac{1}{\mathrm{\&epsiv;}}{y}_{\left(t\right)}\\ {\stackrel{\&CenterDot;}{y}}_{\left(t\right)}=\frac{1}{\mathrm{\&epsiv;}}\&lsqb;-y_{\left(t\right)}+e_{2\left(t\right)}\&rsqb;\end{array}\right.---\left(20\right)$

其中，e_2(t)同式(12)中e₂，p_2(t)和y_(t)∈R³是辅助滤波信号，k和ε是正常数，τ∈R^27×27是常数对角阵，σ∈R²⁷，sat(·)∈R³是由三个标准的饱和函数组成的向量,可以看出，和是有界的，并且和是有界的，饱和函数表达式如下：

$sat\left(x\right)=\left\{\begin{array}{c}M,x>M\\ x,N\&le;x\&le;M\\ N,x<N\end{array}\right.,---\left(21\right)$

其中，M为饱和函数输出上限，N为饱和函数输出下限；

式(14)的右边加上再减去N_d，得到开环误差系统：

$M\left(x\right)\stackrel{\&CenterDot;}{r}=-\frac{1}{2}\stackrel{\&CenterDot;}{M}\left(x\right)r-e_2+\tilde{N}+N_d-S^{-T}AC\stackrel{\&CenterDot;}{\mathrm{\&delta;}}---\left(22\right)$

其中由中值定理可知其中，z＝[e₁e₂r]^T；

根据开环误差系统，设计控制输入

$\stackrel{\&CenterDot;}{\mathrm{\&delta;}}\left(t\right)=C^{-1}{A}^{-1}{S}^T(u+\hat{f})---\left(23\right)$

其中，

u＝(K_s+I_3×3)r+βSgn(e₂)(24)

其中，e₂＝[e₂₁e₂₂e₂₃]^T，e₂₁，e₂₂和e₂₃表示e₂向量在滚转、俯仰、偏航三个方向的分量；并且，K_s和β∈R^3×3是对角正定矩阵，I_3×3∈R^3×3是单位矩阵,符号函数为列向量，即Sgn(e₂)＝[Sgn(e₂₁),Sgn(e₂₂),Sgn(e₂₃)]^T；

对式(23)积分，可得控制输入δ(t)表达式如下：

$\mathrm{\&delta;}\left(t\right)=C^{-1}{A}^{-1}\underset{0}{\overset{t}{\&Integral;}}{S}^T\left(\right(K_s+I_{3\&times;3})r+\mathrm{\&beta;}Sgn(e_2)+\hat{f})d\mathrm{\&tau;},---\left(25\right)$

带入r可得：

$\mathrm{\&delta;}\left(t\right)=C^{-1}{A}^{-1}\underset{0}{\overset{t}{\&Integral;}}{S}^T\left(\right(K_s+I_{3\&times;3}\left)\right({\stackrel{\&CenterDot;}{e}}_2+{\mathrm{\&alpha;e}}_2)+\mathrm{\&beta;}Sgn(e_2)+\hat{f})d\mathrm{\&tau;},---\left(26\right)$

其中，作为控制器前馈项，用于改善跟踪性能，减小控制增益β；

将式(16)、式(23)带入开环误差系统(22)，得到闭环误差系统：

$M\left(x\right)\stackrel{\&CenterDot;}{r}=-\frac{1}{2}\stackrel{\&CenterDot;}{M}\left(x\right)r-e_2+\tilde{N}+N_d-(K_s+I_{3\&times;3})r-\mathrm{\&beta;}Sgn\left(e_2\right)---\left(27\right).$

本发明针对单旋翼无人直升机的姿态控制问题，提出了一套可行的实施方案。其所具有的优点和有益效果如下：

1、根据单旋翼无人直升机姿态动力学模型特性，设计了一种新的非线性半连续鲁棒控制算法，从而消除了使用符号函数引起的抖震现象，改善了控制性能；

2、在此基础上，设计了模糊前馈控制器，对系统的不确定性以及外界扰动进行补偿，从而也减小了主控制器的控制增益，改善了控制性能；

3、考虑到以往设计的控制算法往往局限于数值仿真，本实验组自主开发了三自由度飞行平台，进行了相应飞行实验，在阵风扰动条件下通过与LQR控制方法的对比实验，验证了本发明较高的控制精度和较强鲁棒性。

附图说明

图1：单旋翼无人直升机硬件在环飞行实验平台

图2：LQR无风扰实验结果，姿态角：滚转角、俯仰角、偏航角。

图3：LQR无风扰实验结果，控制量：滚转通道、俯仰通道、偏航通道。

图4：本发明提出控制器无风扰实验结果，姿态角：滚转角、俯仰角、偏航角。

图5：本发明提出控制器无风扰实验结果，控制量：滚转通道、俯仰通道、偏航通道。

图6：LQR风扰实验结果，姿态角：滚转角、俯仰角、偏航角。

图7：LQR风扰实验结果，控制量：滚转通道、俯仰通道、偏航通道。

图8：本发明提出控制器风扰实验结果，姿态角：滚转角、俯仰角、偏航角。

图9：本发明提出控制器风扰实验结果，控制量：滚转通道、俯仰通道、偏航通道。

图10：本发明工作原理框图。

具体实施方式

本发明属于微小型旋翼式无人飞行器自主飞行控制研究领域，主要针对一种单旋翼无人飞行器的控制算法设计，包括无人直升机姿态动力学模型介绍、基于模糊前馈非线性半连续鲁棒控制器设计以及姿态飞行控制实验。

本发明针对单旋翼无人直升机的姿态控制问题，设计了基于模糊前馈补偿的非线性半连续鲁棒控制算法，并进行了基于Lyapunov方法的稳定性分析，证明了设计的控制器能够实现无人直升机姿态的半全局渐进跟踪控制。姿态控制飞行实验结果表明，本发明可以使无人直升机实现快速、准确的镇定控制，并且该控制器对模型先验知识依赖较低，对系统的不确定性具有良好的鲁棒性。

本发明提出了一种新颖的基于模糊前馈补偿的非线性半连续鲁棒控制方法。该方法对于单旋翼无人直升机系统模型的不确定性以及环境的干扰具有很强的适应性，可显著提高单旋翼无人直升机的姿态控制精度，增强系统鲁棒性。基于本实验组自主开发的单旋翼无人直升机三自由度实验平台进行了算法实验验证，利用实验平台主控制器完成本发明控制算法运算，通过数据链路将控制命令发送至底层控制器，最后由底层控制器产生舵机驱动信号，通过控制舵机位置改变单旋翼无人直升机平衡盘角度，从而改变主旋翼分力方向，进而完成单旋翼无人直升机的姿态控制。

一、单旋翼无人直升机动力学模型

单旋翼无人直升机姿态动力学模型形式如下：

$M\left(x\right)\stackrel{\&CenterDot;\&CenterDot;}{x}+C(x,\stackrel{\&CenterDot;}{x})\stackrel{\&CenterDot;}{x}+G\left(x\right)={\mathrm{\&tau;}}^I---\left(1\right)$

其中，x＝[φθψ]^T为列向量，分别表示滚转角、俯仰角、偏航角，表示对x的一阶导，表示对x的二阶导，M(x)∈R^3×3表示惯性矩阵，表示科里奥利力矩阵，G(x)∈R³表示保守力向量，τ^I(t)∈R³表示转矩输入向量，其中R表示全体实数。

式(1)具有如下三个性质：

性质4.惯性矩阵M(x)是对称，正定的，并满足下面的不等式

$m_1\left|\right|\mathrm{\&xi;}|{|}^2\&le;{\mathrm{\&xi;}}^{T}M\left(x\right)\mathrm{\&xi;}\&le;m_2|\left|\mathrm{\&xi;}\right|{|}^2,\&ForAll;\mathrm{\&xi;}\&Element;R^3---\left(2\right)$

其中m₁，m₂是正常数。

性质5.科里奥利力和保守力满足如下不等式

$\left|\right|C(x,\stackrel{\&CenterDot;}{x})\left|\right|\&le;{\mathrm{\&zeta;}}_c\left|\right|\stackrel{\&CenterDot;}{x}\left|\right|,\left|\right|G\left(x\right)\left|\right|\&le;{\mathrm{\&zeta;}}_g---\left(3\right)$

其中，是有界正常数。

性质6.科里奥利力矩阵满足下面的关系

$C(x,\mathrm{\&xi;})\mathrm{\&upsi;}=C(x,\mathrm{\&upsi;})\mathrm{\&xi;},\&ForAll;\mathrm{\&xi;},\mathrm{\&upsi;}\&Element;R^3---\left(4\right)$

式(1)中，转矩输入τ^I(t)可通过纵向挥舞角a(t)、横向挥舞角b(t)和尾桨推力T_T表示，即：

τ^I＝S^-T(A(T_M)υ_c+B(T_M)),(5)

其中υ_c＝[abT_T]^T，T_M∈R表示主桨推力，A(T_M)∈R^3×3是可逆矩阵，B(T_M)∈R³是主桨推力的分力向量,S表示从体坐标系到惯性坐标系的平移矩阵，其表达式如下：

$S=\left[\begin{array}{ccc}1& \frac{s_{\mathrm{\&phi;}}{s}_{\mathrm{\&theta;}}}{c_{\mathrm{\&theta;}}}& \frac{c_{\mathrm{\&phi;}}{s}_{\mathrm{\&theta;}}}{c_{\mathrm{\&theta;}}}\\ 0& c_{\mathrm{\&phi;}}& -s_{\mathrm{\&phi;}}\\ 0& \frac{s_{\mathrm{\&phi;}}}{c_{\mathrm{\&theta;}}}& \frac{c_{\mathrm{\&phi;}}}{c_{\mathrm{\&theta;}}}\end{array}\right].---\left(6\right)$

其中，s_φ表示sinφ，s_θ表示sinθ，c_φ表示cosφ，c_θ表示cosθ。

单旋翼无人直升机旋翼动力学模型形式如下：

$\begin{array}{c}\stackrel{\&CenterDot;}{a}=-\frac{{\mathrm{\&tau;}}_{mr}+K_{sb}{\mathrm{\&tau;}}_{sb}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}q-\frac{1}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}a+\frac{{\mathrm{\&tau;}}_{mr}{A}_b}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}b+\frac{A_{lon}+K_{sb}{C}_{lon}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}{\mathrm{\&delta;}}_{lon}\\ \stackrel{\&CenterDot;}{b}=-\frac{{\mathrm{\&tau;}}_{mr}+K_{sb}{\mathrm{\&tau;}}_{sb}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}p+\frac{{\mathrm{\&tau;}}_{mr}{B}_a}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}a-\frac{1}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}b+\frac{B_{lon}+K_{sb}{D}_{lat}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}{\mathrm{\&delta;}}_{lat}\\ {\stackrel{\&CenterDot;}{T}}_T=B_{ped}+K_{ped}{\mathrm{\&delta;}}_{ped},\end{array}---\left(7\right)$

其中，a表示纵向挥舞角，表示a的一阶导数，b表示横向挥舞角，表示b的一阶导数，T_T表示尾桨推力，表示T_T的一阶导数，p表示俯仰角速度，表示p的一阶导数，q表示滚转角速度，表示q的一阶导数，δ_lon表示控制输入纵向周期变矩，δ_lat表示控制输入横向周期变矩，δ_ped表示控制输入尾桨矩，τ_mr表示主旋翼挥舞时间常数，τ_sb表示副翼挥舞时间常数，A_b表示主旋翼纵向伺服输入比例系数，B_a表示主旋翼横向伺服输入比例系数，C_lon表示副翼纵向伺服输入比例系数，D_lat表示副翼横向伺服输入比例系数，K_sb表示主旋翼与副翼伺服输入比值，B_ped表示尾桨输入常数，K_ped表示尾桨伺服输入比例系数，A_lon表示主旋翼纵向指令输入比例系数，B_lat表示主旋翼横向指令输入比例系数。

在悬停状态下，可对上述模型进行简化，得到如下形式：

$\begin{array}{c}a=A_{b}b-A_{lon0}{\mathrm{\&delta;}}_{lon}\\ b=-B_{a}a+B_{lat0}{\mathrm{\&delta;}}_{lat}\\ T_T=K_{ped0}{\mathrm{\&delta;}}_{ped},\end{array}---\left(8\right)$

其中，A_lon0、B_lat0和K_ped0为简化后常数，将式(8)带入式(5)，可得简化的单旋翼无人直升

机旋翼动力学模型形式如下：

τ^I＝S^-T(A(T_M)Cδ+B(T_M)),(9)

其中δ＝[δ_lonδ_latδ_ped]^T表示控制输入，常数阵C∈R^3×3定义如下：

$C=\left[\begin{array}{ccc}-\frac{A_{lon}}{A_b{B}_a+1}& \frac{A_b{B}_{lat}}{A_b{B}_a+1}& 0\\ \frac{B_{lat}}{A_b{B}_a+1}& \frac{B_a{A}_{lon}}{A_b{B}_a+1}& 0\\ 0& 0& K_{ped0}\end{array}\right].---\left(10\right)$

二、基于模糊前馈的非线性半连续鲁棒控制器设计

为实现控制目标，首先定义姿态跟踪误差：

e₁＝x_d-x,(11)

其中x_d∈R³表示期望姿态向量，x∈R³表示实际姿态向量。

定义滤波误差：

$e_2={\stackrel{\&CenterDot;}{e}}_1+e_1---\left(12\right)$

$r={\stackrel{\&CenterDot;}{e}}_2+{\mathrm{\&alpha;e}}_2---\left(13\right)$

其中，α∈R^3×3是正定对角常数阵，e₂＝[e₂₁e₂₂e₂₃]^T，e₂₁，e₂₂和e₂₃表示e₂向量在滚转、俯仰、偏航三个方向的分量，r＝[r₁r₂r₃]^T，r₁、r₂和r₃表示r向量在滚转、俯仰、偏航三个方向的分量，表示e₁的一阶导数，表示e₂的一阶导数；

对式(13)求导，等式的两边同时左乘M(x)，然后带入式(1)，可得

$M\left(x\right)\stackrel{\&CenterDot;}{r}=-\frac{1}{2}\stackrel{\&CenterDot;}{M}\left(x\right)r-e_2+N-S^{-T}AC\stackrel{\&CenterDot;}{\mathrm{\&delta;}},---\left(14\right)$

其中，A同公式(5)中A(T_M)，表示控制输入δ的一阶导数。辅助函数N定义如下：

$\begin{array}{c}N=M\left(x\right)({\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{x}}_d+{\stackrel{\&CenterDot;\&CenterDot;}{e}}_1+\mathrm{\&alpha;}{\stackrel{\&CenterDot;}{e}}_2)+\stackrel{\&CenterDot;}{M}\left(x\right)(\stackrel{\&CenterDot;\&CenterDot;}{x}+\frac{1}{2}r)+e_2+\stackrel{\&CenterDot;}{C}\stackrel{\&CenterDot;}{x}+C\stackrel{\&CenterDot;\&CenterDot;}{x}+\stackrel{\&CenterDot;}{G}-\\ {\stackrel{\&CenterDot;}{S}}^{-T}{S}^T\left(M\right(x)\stackrel{\&CenterDot;\&CenterDot;}{x}+C\stackrel{\&CenterDot;}{x}+G-S^{-T}B)\end{array}---\left(15\right)$

其中，B同公式(5)中B(T_M)，G同公式(1)中G(x)。引入有||N_d(t)||，L_∞表示无穷范数，并且，N_d可以通过模糊推理方法来逼近，模糊输出：

${\hat{f}}_{\left(t\right)}={\widehat{\mathrm{\&omega;}}}^T\mathrm{\&sigma;}(q_d,{\stackrel{\&CenterDot;}{q}}_d,{\stackrel{\&CenterDot;\&CenterDot;}{q}}_d,{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{q}}_d),---\left(16\right)$

其中σ∈R^27×1，q_d＝[φ_d,θ_d,ψ_d]^T，φ_d，θ_d，ψ_d分别表示滚转角、俯仰角、偏航角参考。模糊系统输入共有12个变量，分别来自3个通道，每个通道有4个状态。由于输入变量过多，通过如下变换，即可使用3个变量用于模糊推理。

$\begin{array}{c}{\stackrel{\&OverBar;}{q}}_1={\mathrm{\&phi;}}_d+{\mathrm{\&theta;}}_d+{\mathrm{\&psi;}}_d+{\stackrel{\&CenterDot;}{\mathrm{\&phi;}}}_d+{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_d+{\stackrel{\&CenterDot;}{\mathrm{\&psi;}}}_d\\ {\stackrel{\&OverBar;}{q}}_2={\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&psi;}}}_d\\ {\stackrel{\&OverBar;}{q}}_3={\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&psi;}}}_d.\end{array}---\left(17\right)$

在模糊推理中，隶属度函数选择高斯函数，表达式如下：

${\mathrm{\&mu;}}_{\stackrel{\&OverBar;}{q}}=e^{\frac{-0.5\&times;{(\stackrel{\&OverBar;}{q}-c)}^2}{{\mathrm{\&omega;}}^2}},---\left(18\right)$

其中，c和ω为高斯函数中实数，高斯函数将状态输入分为3个区间：正、零、负；模糊推理选择乘积运算，由于输入包含3个变量，每个变量分为3个区间，故可产生27条规则；解模糊方法选择中心平均法，表达式如下：

$\mathrm{\&alpha;}=\frac{{\mathrm{\&mu;}}_{F_1^i\&times;F_2^i\&times;F_3^i}}{\underset{i=1}{\overset{p}{\&Sigma;}}{\mathrm{\&mu;}}_{F_1^i\&times;F_2^i\&times;F_3^i}},---\left(19\right)$

其中p＝27，表示第i条规则时模糊推理乘积运算结果。

设计自适应更新率：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{\widehat{\mathrm{\&omega;}}}=-k\widehat{\mathrm{\&omega;}}+\mathrm{\&tau;}\mathrm{\&sigma;}{\&lsqb;sat(e_{2\left(t\right)}+p_{2\left(t\right)})\&rsqb;}^T\\ p_{2\left(t\right)}=\frac{1}{\mathrm{\&epsiv;}}{e}_{2\left(t\right)}-\frac{1}{\mathrm{\&epsiv;}}{y}_{\left(t\right)}\\ {\stackrel{\&CenterDot;}{y}}_{\left(t\right)}=\frac{1}{\mathrm{\&epsiv;}}\&lsqb;-y_{\left(t\right)}+e_{2\left(t\right)}\&rsqb;\end{array}\right.---\left(20\right)$

其中，e_2(t)同式(12)中e₂，p_2(t)和y_(t)∈R³是辅助滤波信号，k和ε是正常数，τ∈R^27×27是常数对角阵，σ∈R²⁷，sat(·)∈R³是由三个标准的饱和函数组成的向量,可以看出，和是有界的，并且和是有界的。饱和函数表达式如下：

$sat\left(x\right)=\left\{\begin{array}{c}M,x>M\\ x,N\&le;x\&le;M\\ N,x<N\end{array}\right.,---\left(21\right)$

其中，M为饱和函数输出上限，N为饱和函数输出下限。

式(14)的右边加上再减去N_d(t)，得到开环误差系统：

$M\left(x\right)\stackrel{\&CenterDot;}{r}=-\frac{1}{2}\stackrel{\&CenterDot;}{M}\left(x\right)r-e_2+\tilde{N}+N_d-S^{-T}AC\stackrel{\&CenterDot;}{\mathrm{\&delta;}}---\left(22\right)$

其中由中值定理可知其中z＝[e₁e₂r]^T。

根据开环误差系统，设计控制输入

$\stackrel{\&CenterDot;}{\mathrm{\&delta;}}\left(t\right)=C^{-1}{A}^{-1}{S}^T(u+\hat{f})---\left(23\right)$

其中，

u＝(K_s+I_3×3)r+βSgn(e₂)(24)

其中，e₂＝[e₂₁e₂₂e₂₃]^T，e₂₁，e₂₂和e₂₃表示e₂向量在滚转、俯仰、偏航三个方向的分量。

并且，K_s和β∈R^3×3是对角正定矩阵，I_3×3∈R^3×3是单位矩阵,符号函数为列向量，即Sgn(e₂)＝[Sgn(e₂₁),Sgn(e₂₂),Sgn(e₂₃)]^T。

对式(23)积分，可得控制输入δ(t)表达式如下：

$\mathrm{\&delta;}\left(t\right)=C^{-1}{A}^{-1}\underset{0}{\overset{t}{\&Integral;}}{S}^T\left(\right(K_s+I_{3\&times;3})r+\mathrm{\&beta;}Sgn(e_2)+\hat{f})d\mathrm{\&tau;},---\left(25\right)$

带入r可得：

$\mathrm{\&delta;}\left(t\right)=C^{-1}{A}^{-1}\underset{0}{\overset{t}{\&Integral;}}{S}^T\left(\right(K_s+I_{3\&times;3}\left)\right({\stackrel{\&CenterDot;}{e}}_2+{\mathrm{\&alpha;e}}_2)+\mathrm{\&beta;}Sgn(e_2)+\hat{f})d\mathrm{\&tau;},---\left(26\right)$

其中，作为控制器前馈项，用于改善跟踪性能，减小控制增益β。

将式(16)、式(23)带入开环误差系统(22)，得到闭环误差系统：

$M\left(x\right)\stackrel{\&CenterDot;}{r}=-\frac{1}{2}\stackrel{\&CenterDot;}{M}\left(x\right)r-e_2+\tilde{N}+N_d-(K_s+I_{3\&times;3})r-\mathrm{\&beta;}Sgn\left(e_2\right)---\left(27\right)$

三、本发明的理论支持

在进行稳定性分析之前，首先介绍两个引理。

引理1：定义辅助函数L＝r^T(N_d-βSgn(e₁))，控制增益α，β如式(13)、(24)中定义，并满足如下条件：

$\mathrm{\&alpha;}\mathrm{\&beta;}\&GreaterEqual;\left|\right|\mathrm{\&alpha;}\left|\right|\left|\right|N_d\left(t\right)\left|\right|+\left|\right|{\stackrel{\&CenterDot;}{N}}_d\left(t\right)\left|\right|,---\left(28\right)$

则有如下不等式成立：

${\&Integral;}_0^{t}L\left(\mathrm{\&tau;}\right)\mathrm{d\&tau;}\&le;{\mathrm{\&zeta;}}_b,$

其中正常数定义如下：

${\mathrm{\&zeta;}}_b=\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}\left({\mathrm{\&beta;}}_i\right|e_{2i}\left(0\right)|-|e_{2i}\left(0\right)\left|N_{\mathrm{di}}\left(0\right)\right).---\left(30\right)$

引理2：对于系统f:R^m×R_≥0→R^m，解是存在的。定义区域D＝{y∈R^m|||y||<ε}，ε是正常数，令V:D×R_≥0→R^m连续可微，并满足如下条件：

W₁(y)≤V(y,t)≤W₂(y)(31)

$\stackrel{\&CenterDot;}{V}(y,t)\&le;-W\left(y\right),---\left(32\right)$

其中W₁(y)，W₂(y)是连续正定函数，W(y)是一致连续半正定函数。只需上两式成立，且y(0)∈S，则下式成立：

W(y(t))→0,t→∞(33)

其中域S定义如下：

$S=\{y\&Element;D|W_2\left(y\right)\&le;\mathrm{\&delta;}\},\mathrm{\&delta;}<\underset{\left|\right|\mathrm{\&xi;}\left|\right|=\mathrm{\&epsiv;}}{min}{W}_1\left(y\right)---\left(34\right)$

其中δ是正常数。

本发明的稳定性分析主要结果可由下列定理给出。

定理1：对于单旋翼无人机系统，式(23)设计控制器可以使飞行器姿态获得半全局渐进稳定的跟踪效果，即：

当t→∞时，有e₁(t)→0,y(0)∈S.(35)

证明:

首先定义辅助函数P(t)∈R，表达式如下：

$P\left(t\right)={\mathrm{\&zeta;}}_b-{\&Integral;}_0^{t}L\left(\mathrm{\&tau;}\right)\mathrm{d\&tau;},---\left(36\right)$

其中，L(t)如引理1定义。由引理1可以得到P(t)≥0。

定义V:R³×R_≥0×R_≥0→R_≥0如下：

$V(y,t)=\frac{1}{2}(e_1^T{e}_1+e_2^T{e}_2)+\frac{1}{2}{r}^{T}M\left(x\right)r+P---\left(37\right)$

其中y(t)定义如下：

$y\left(t\right)={\left[\begin{array}{cc}{z}^T& \sqrt{P}\end{array}\right]}^T---\left(38\right)$

其中，z定义如式(22)。由模型的性质1，对式(37)进行放缩，可得：

λ₁||y||²≤V≤λ₂||y||²,(39)

其中

将式(37)对时间求导数，带入式(12)，(13)，(25)可得：

$\stackrel{\&CenterDot;}{V}=-e_1^T{e}_1-e_2^T{\mathrm{\&alpha;}}^T{e}_2+e_1^T{e}_2-r^{T}r+r^T\tilde{N}-r^T{K}_{s}r,---\left(40\right)$

利用对上式进行放缩可得：

$\stackrel{\&CenterDot;}{V}\&le;-{\mathrm{\&lambda;}}_3\left|\right|z|{|}^2+|\left|r\right|\left|\mathrm{\&rho;}\left(\right|\left|z\right|\left|\right)\right|\left|z\right||-K_s|\left|r\right|{|}^2,---\left(41\right)$

进一步放缩最终得到：

$\stackrel{\&CenterDot;}{V}\&le;-({\mathrm{\&lambda;}}_3-\frac{{\mathrm{\&rho;}}^2\left(\right|\left|z\right|\left|\right)}{4K_s})\left|\right|z|{|}^2,---\left(42\right)$

其中λ₃＝min{1/2,α-1/2}，α>1/2，由式(42)，可得只要满足下面条件：

$K_s>\frac{1}{4{\mathrm{\&lambda;}}_3}{\mathrm{\&rho;}}^2\left(\right|\left|z\right|\left|\right),\left|\right|z\left|\right|<{\mathrm{\&rho;}}^{-1}\left(2\sqrt{{\mathrm{\&lambda;}}_3{K}_s}\right).---\left(43\right)$

下面，对式(39)至式(43)应用引理2。可以得到式(37)的上下界函数和式(37)导数的上界函数如下:

W₁(y)＝λ₁||y||²,W₂(y)＝λ₂||y||²,W(y)＝γ||z||²(44)

利用式(43)定义区域D，表达式如下：

$D=\{y\&Element;R^3\&times;R_{\&GreaterEqual;0}|\left|\right|y\left|\right|<{\mathrm{\&rho;}}^{-1}\left(2\sqrt{{\mathrm{\&lambda;}}_3{K}_s}\right)\},---\left(45\right)$

由式(37)至式(43)，可知V(y(t),t)∈L_∞，因此e₁(t),e₂(t),r(t)∈L_∞。由式(13)可知由式(25)和式(43)可知则有W(y(t))一致连续。定义区域S，表达式如下：

$S=\{y\&Element;D|W_2\left(y\right)<{\mathrm{\&lambda;}}_1{\left({\mathrm{\&rho;}}^{-1}\right(2\sqrt{{\mathrm{\&lambda;}}_3{K}_s}\left)\right)}^2\},---\left(46\right)$

则由引理2可得，当t→∞时，有

||z||²→0,y(0)∈S,(47)

进一步可以得到，当t→∞时，有

e₁(t)→0,y(0)∈S.(48)

当控制增益K_s足够大，吸引域式(46)可以包含更多的初始状态，这样就可以得到一个半全局渐进稳定的结论。由式(44)和式(46)可以算得如下吸引域：

$W_2\left(y\right(0\left)\right)<{\mathrm{\&lambda;}}_1{\left({\mathrm{\&rho;}}^{-1}\right(2\sqrt{{\mathrm{\&lambda;}}_3{K}_s}\left)\right)}^2\&DoubleRightArrow;\left|\right|y\left(0\right)\left|\right|<\sqrt{\frac{{\mathrm{\&lambda;}}_1}{{\mathrm{\&lambda;}}_2\left(\right|\left|y\right(0)\left|\right|)}}{\mathrm{\&rho;}}^{-1}\left(2\sqrt{{\mathrm{\&lambda;}}_3{K}_s}\right)---\left(49\right)$

即

下面结合附图和具体实施方式进一步说明本发明。

一、系统硬件连接及配置

图1展示了单旋翼无人直升机硬件在环飞行实验平台的硬件组成和连接关系。机身本体选用TREX450小型无线电操控直升机，该小型航模直升机机身长640mm，主桨长度为710mm，飞机总重约638g，有效负载约为500g。机载传感器选用Xsens公司生产的MTI姿态航向参考系统，该传感器最高更新频率为120Hz，提供三轴角速度及三轴姿态角，其中俯仰角和滚转角精度为±0.5°，偏航角精度为±1°。该平台由宿主机生成仿真代码，下载到PC/104目标机，由目标机进行复杂控制算法的计算，其采样频率最高可达100kHZ，保证了控制系统的实时性。主控制器PC/104主要分为三个模块:数据采集模块，该模块负责惯性导航单元的数据采集与处理；飞行控制模块，该模块负责控制器算法的运行；数据通讯模块，该模块负责主控制器与底层控制器之间的数据传输，如控制量等。

此外，本研究组自主设计开发了基于DSP(型号为TMS320F28335)处理器的底层控制器，该底层控制器配有主控模块、数据采集模块、通讯模块及手自动切换模块。其中，主控模块负责控制算法的运算，数据采集模块负责传感器MTI的数据采集，通讯模块负责DSP与上位机信息交互，手自动切换模块负责接收机PPM信号捕捉和舵机PWM信号输出。

二、无风扰对比实验

在设计硬件在环飞行实验时，模糊前馈控制器参数选择：对应隶属度函数参数为c₁＝-5，c₂＝0，c₃＝5，w₁＝10，对应隶属度函数参数为c₁＝-2，c₂＝0，c₃＝2，w₁＝4，对应隶属度函数参数为c₁＝-5，c₂＝0，c₃＝5，w₁＝10，k＝0.001，ε＝0.01，sat(·)上下界±100，τ＝diag[600,600,500,400,600,800,800,600,500,800,800,800,500,500,500,600,500,800,500,600,500,600,500,800,800,500,600]。主控制器参数选择：α＝diag[0.0002,0.0002，0.00001],β＝diag[0.06,0.06,0.00002]以及K_s＝diag[1.5,1.8,0.005-1]。通过编写控制算法，由宿主机将生成代码下载到PC/104控制器中，进行小型无人直升机实物飞行实验。

首先，利用PC/104进行本发明控制算法δ(t)的运算，该向量三个分量分别对应单旋翼无人直升机滚转、俯仰、偏航通道的控制输入；然后，通过PC/104数据通讯模块将控制输入δ(t)发送至DSP底层控制器；最后，由DSP底层控制器完成单旋翼无人直升机舵机驱动信号的产生，从而完成单旋翼无人直升机的姿态控制。

图2和图3分别给出了传统LQR控制器无风干扰条件下姿态镇定实验的姿态角及控制量，在37s切换为自动，实验进行约100s，整个过程中，三通道控制精度均为2度，控制量为归一化后的数据，均在正常范围以内。

图4和图5分别给出了本发明设计控制器无风干扰条件下姿态镇定实验的姿态角以及控制量，在52s时切换为自动控制，约3s后无人机达到稳定状态。在125s时再次进行手动自动切换，相当于跟踪15度阶跃信号，大约在8s后无人机达到稳定状态。在180s和190s时进行了偏航通道50度阶跃信号跟踪，大约3s后达到稳定状态。整个过程，滚转、俯仰通道稳态精度达1度，偏航通道稳态精度达1.5度，控制量为归一化后的数据，均在正常范围以内。

通过对比实验验证了本发明设计控制器对于单旋翼无人直升机姿态具有较高的控制精度，验证了控制算法有效性。

三、抗风扰对比实验

在阵风条件下进行单旋翼无人直升机姿态镇定实验，参数选择如前所述，加入阵风风速大小为4m/s-6.5m/s。

首先，利用PC/104进行姿态传感器MTI的数据采集并完成本发明中控制算法δ(t)的运算，控制量δ(t)中三个分量分别对应单旋翼无人直升机滚转、俯仰、偏航通道的控制输入；然后，通过PC/104数据通讯模块将控制量δ(t)发送至DSP底层控制器；最后，由DSP底层控制器完成单旋翼无人直升机舵机驱动信号的产生，从而驱动舵机完成单旋翼无人直升机的姿态控制。

图6和图7分别给出了传统LQR控制器阵风干扰条件下姿态镇定实验的姿态角及控制量，在28s切换为自动，三通道稳态精度均为2度。在无人机稳定后，70s时加入阵风扰动，随即发散，失去控制。

图8和图9分别给出了本发明设计控制器阵风干扰条件下姿态镇定实验的姿态角以及控制量，在10s时切换自动控制，约4s后无人机达到稳定状态，在80s时加入阵风干扰，实验进行约180s，整个过程，滚转、俯仰通道稳态精度达2度，偏航通道稳态精度达2度，控制量为归一化后的数据，均在正常范围以内。

通过对比实验验证了本发明设计控制器对于单旋翼无人直升机姿态具有较高的抗干扰能力，验证了控制算法的鲁棒性。

Claims (1)

1.一种基于模糊前馈单旋翼无人直升机姿态非线性鲁棒控制方法，其特征是，包括如下步骤：

一、单旋翼无人直升机动力学模型

单旋翼无人直升机姿态动力学模型形式如下：

$M\left(x\right)\stackrel{\&CenterDot;\&CenterDot;}{x}+C(x,\stackrel{\&CenterDot;}{x})\stackrel{\&CenterDot;}{x}+G\left(x\right)=\mathrm{\&tau;}---\left(1\right)$

其中，x＝[φθψ]^T为列向量，φ、θ、ψ分别表示滚转角、俯仰角、偏航角，表示对x的一阶导，表示对x的二阶导，M(x)∈R^3×3表示惯性矩阵，表示科里奥利力矩阵，G(x)∈R³表示保守力向量，τ∈R³表示转矩输入向量，其中R表示全体实数；

式(1)具有如下三个性质：

性质1.惯性矩阵M(x)是对称，正定的，并满足下面的不等式

$m_1\left|\right|\mathrm{\&xi;}|{|}^2\&le;{\mathrm{\&xi;}}^{T}M\left(x\right)\mathrm{\&xi;}\&le;m_2|\left|\mathrm{\&xi;}\right|{|}^2,\&ForAll;\mathrm{\&xi;}\&Element;R^3---\left(2\right)$

其中m₁，m₂是正常数；

性质2.科里奥利力和保守力满足如下不等式

其中是有界正常数；

性质3.科里奥利力矩阵满足下面的关系

$C(x,\mathrm{\&xi;})\mathrm{\&upsi;}=C(x,\mathrm{\&upsi;})\mathrm{\&xi;},\&ForAll;\mathrm{\&xi;},\mathrm{\&upsi;}\&Element;R^3---\left(4\right)$

式(1)中，转矩输入τ可通过纵向挥舞角a、横向挥舞角b和尾桨推力T_T表示，即：

τ＝S^-T(A(T_M)υ_c+B(T_M)),(5)

其中υ_c＝[abT_T]^T，T_M∈R表示主桨推力，A(T_M)∈R^3×3是可逆矩阵，B(T_M)∈R³是主桨推力的分力向量,S表示从体坐标系到惯性坐标系的平移矩阵，其表达式如下：

$S=\left[\begin{array}{ccc}1& \frac{s_{\mathrm{\&phi;}}{s}_{\mathrm{\&theta;}}}{c_{\mathrm{\&theta;}}}& \frac{c_{\mathrm{\&phi;}}{s}_{\mathrm{\&theta;}}}{c_{\mathrm{\&theta;}}}\\ 0& c_{\mathrm{\&phi;}}& -s_{\mathrm{\&phi;}}\\ 0& \frac{s_{\mathrm{\&phi;}}}{c_{\mathrm{\&theta;}}}& \frac{c_{\mathrm{\&phi;}}}{c_{\mathrm{\&theta;}}}\end{array}\right]---\left(6\right)$

其中，s_φ表示sinφ，s_θ表示sinθ，c_φ表示cosφ，c_θ表示cosθ；

单旋翼无人直升机旋翼动力学模型形式如下：

$\begin{array}{c}\stackrel{\&CenterDot;}{a}=-\frac{{\mathrm{\&tau;}}_{mr}+K_{sb}{\mathrm{\&tau;}}_{sb}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}q-\frac{1}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}a+\frac{{\mathrm{\&tau;}}_{mr}{A}_b}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}b+\frac{A_{lon}+K_{sb}{C}_{lon}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}{\mathrm{\&delta;}}_{lon}\\ \stackrel{\&CenterDot;}{b}=-\frac{{\mathrm{\&tau;}}_{mr}+K_{sb}{\mathrm{\&tau;}}_{sb}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}p+\frac{{\mathrm{\&tau;}}_{mr}{B}_a}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}a-\frac{1}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}b+\frac{A_{lat}+K_{sb}{C}_{lat}}{{\mathrm{\&tau;}}_{mr}+{\mathrm{\&tau;}}_{sb}}{\mathrm{\&delta;}}_{lat}\\ {\stackrel{\&CenterDot;}{T}}_T=B_{ped}+K_{ped}{\mathrm{\&delta;}}_{ped},\end{array}---\left(7\right)$

其中，a表示纵向挥舞角，表示a的一阶导数，b表示横向挥舞角，表示b的一阶导数，T_T表示尾桨推力，表示T_T的一阶导数，p表示俯仰角速度，q表示滚转角速度，δ_lon表示控制输入纵向周期变矩，δ_lat表示控制输入横向周期变矩，δ_ped表示控制输入尾桨矩，τ_mr表示主旋翼挥舞时间常数，τ_sb表示副翼挥舞时间常数，A_b表示主旋翼纵向伺服输入比例系数，B_a表示主旋翼横向伺服输入比例系数，C_lon表示副翼纵向伺服输入比例系数，D_lat表示副翼横向伺服输入比例系数，K_sb表示主旋翼与副翼伺服输入比值，B_ped表示尾桨输入常数，K_ped表示尾桨伺服输入比例系数,A_lon表示主旋翼纵向指令输入比例系数，B_lat表示主旋翼横向指令输入比例系数；

在悬停状态下，对上述模型进行简化，得到如下形式：

$\begin{array}{c}a=A_{b}b-A_{lon0}{\mathrm{\&delta;}}_{lon}\\ b=-B_{a}a+B_{lat0}{\mathrm{\&delta;}}_{lat}\\ T_T=K_{ped0}{\mathrm{\&delta;}}_{ped},\end{array}---\left(8\right)$

其中，A_lon0、B_lat0和K_ped0为简化后常数，将式(8)带入式(5)，可得简化的单旋翼无人直升机旋翼动力学模型形式如下：

τ＝S^-T(A(T_M)Cδ+B(T_M)),(9)

其中δ＝[δ_lonδ_latδ_ped]^T表示控制输入，常数阵C∈R^3×3定义如下：

$C=\left[\begin{array}{ccc}-\frac{A_{lon0}}{A_b{B}_a+1}& \frac{A_b{B}_{lat0}}{A_b{B}_a+1}& 0\\ \frac{B_{lat0}}{A_b{B}_a+1}& \frac{B_a{A}_{lon0}}{A_b{B}_a+1}& 0\\ 0& 0& K_{ped0}\end{array}\right]---\left(10\right)$

二、基于模糊前馈的非线性半连续鲁棒控制器设计

为实现控制目标，首先定义姿态跟踪误差：

e₁＝x_d-x,(11)

其中x_d∈R³表示期望姿态向量，x∈R³表示实际姿态向量；

定义滤波误差e₂和r，表达式如下：

$e_2={\stackrel{\&CenterDot;}{e}}_1+e_1---\left(12\right)$

$r={\stackrel{\&CenterDot;}{e}}_2+{\mathrm{\&alpha;e}}_2---\left(13\right)$

其中，α∈R^3×3是正定对角常数阵，e₂＝[e₂₁e₂₂e₂₃]^T，e₂₁，e₂₂和e₂₃表示e₂向量在滚转、俯仰、偏航三个方向的分量，r＝[r₁r₂r₃]^T，r₁、r₂和r₃表示r向量在滚转、俯仰、偏航三个方向的分量，表示e₁的一阶导数，表示e₂的一阶导数；

对式(13)求导，等式的两边同时左乘M(x)，然后带入式(1)，可得

$M\left(x\right)\stackrel{\&CenterDot;}{r}=-\frac{1}{2}\stackrel{\&CenterDot;}{M}\left(x\right)r-e_2+N-S^{-T}AC\stackrel{\&CenterDot;}{\mathrm{\&delta;}},---\left(14\right)$

其中，A同公式(5)中A(T_M)，表示控制输入δ的一阶导数，辅助函数N定义如下：

$\begin{array}{c}N=M\left(x\right)({\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{x}}_d+{\stackrel{\&CenterDot;\&CenterDot;}{e}}_1+\mathrm{\&alpha;}{\stackrel{\&CenterDot;}{e}}_2)+\stackrel{\&CenterDot;}{M}\left(x\right)(\stackrel{\&CenterDot;\&CenterDot;}{x}+\frac{1}{2}r)+e_2+\stackrel{\&CenterDot;}{C}\stackrel{\&CenterDot;}{x}+C\stackrel{\&CenterDot;\&CenterDot;}{x}+\stackrel{\&CenterDot;}{G}-\\ {\stackrel{\&CenterDot;}{S}}^{-T}{S}^T\left(M\right(x)\stackrel{\&CenterDot;\&CenterDot;}{x}+C\stackrel{\&CenterDot;}{x}+G-S^{-T}B)\end{array}---\left(15\right)$

其中，B同公式(5)中B(T_M)，G同公式(1)中G(x)，引入有||N_d||，L_∞表示无穷范数，并且，N_d可以通过模糊推理方法来逼近，模糊输出：

$\hat{f}={\widehat{\mathrm{\&omega;}}}^T\mathrm{\&sigma;}\left(q_d,{\stackrel{\&CenterDot;}{q}}_d,{\stackrel{\&CenterDot;\&CenterDot;}{q}}_d,{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{q}}_d\right),---\left(16\right)$

其中为模糊推理权值，σ∈R^27×1，q_d＝[φ_d,θ_d,ψ_d]^T，φ_d，θ_d，ψ_d分别表示滚转角、俯仰角、偏航角参考；模糊系统输入共有12个变量，分别来自3个通道，每个通道有4个状态，由于输入变量过多，通过如下变换，即可使用3个变量用于模糊推理；

$\begin{array}{c}{\stackrel{\&OverBar;}{q}}_1={\mathrm{\&phi;}}_d+{\mathrm{\&theta;}}_d+{\mathrm{\&psi;}}_d+{\stackrel{\&CenterDot;}{\mathrm{\&phi;}}}_d+{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_d+{\stackrel{\&CenterDot;}{\mathrm{\&psi;}}}_d\\ {\stackrel{\&OverBar;}{q}}_2={\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&psi;}}}_d\\ {\stackrel{\&OverBar;}{q}}_3={\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&phi;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_d+{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&psi;}}}_d\end{array}---\left(17\right)$

在模糊推理中，隶属度函数选择高斯函数，表达式如下：

${\mathrm{\&mu;}}_{\stackrel{\&OverBar;}{q}}=e^{\frac{-0.5\&times;{(\stackrel{\&OverBar;}{q}-c)}^2}{{\mathrm{\&omega;}}^2}},---\left(18\right)$

其中，c和ω为高斯函数中实数，高斯函数将状态输入分为3个区间：正、零、负；模糊推理选择乘积运算，由于输入包含3个变量，每个变量分为3个区间，故可产生27条规则；解模糊方法选择中心平均法，表达式如下：

${\mathrm{\&sigma;}}_i=\frac{{\mathrm{\&mu;}}_{{F_1}^i\&times;{F_2}^i\&times;{F_3}^i}}{\underset{i=1}{\overset{p}{\&Sigma;}}{\mathrm{\&mu;}}_{{F_1}^i\&times;{F_2}^i\&times;{F_3}^i}},---\left(19\right)$

其中p＝27，表示第i条规则时模糊推理乘积运算结果；

设计自适应更新率：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{\widehat{\mathrm{\&omega;}}}\left(t\right)=-k\widehat{\mathrm{\&omega;}}\left(t\right)+\mathrm{\&tau;}\mathrm{\&sigma;}{\&lsqb;sat(e_{2\left(t\right)}+p_{2\left(t\right)})\&rsqb;}^T\\ p_{2\left(t\right)}=\frac{1}{\mathrm{\&epsiv;}}{e}_{2\left(t\right)}-\frac{1}{\mathrm{\&epsiv;}}{y}_{\left(t\right)}\\ {\stackrel{\&CenterDot;}{y}}_{\left(t\right)}=\frac{1}{\mathrm{\&epsiv;}}\&lsqb;-y_{\left(t\right)}+e_{2\left(t\right)}\&rsqb;\end{array}\right.---\left(20\right)$

其中，e_2(t)同式(12)中e₂，p_2(t)和y_(t)∈R³是辅助滤波信号，k和ε是正常数，τ∈R^27×27是常数对角阵，σ∈R²⁷，sat(·)∈R³是由三个标准的饱和函数组成的向量,可以看出，和是有界的，并且和是有界的，饱和函数表达式如下：

$sat\left(x\right)=\left\{\begin{array}{c}M,x>M\\ x,N\&le;x\&le;M\\ N,x<N\end{array}\right.,---\left(21\right)$

其中，M为饱和函数输出上限，N为饱和函数输出下限；

式(14)的右边加上再减去N_d，得到开环误差系统：

$M\left(x\right)\stackrel{\&CenterDot;}{r}=-\frac{1}{2}\stackrel{\&CenterDot;}{M}\left(x\right)r-e_2+\tilde{N}+N_d-S^{-T}AC\stackrel{\&CenterDot;}{\mathrm{\&delta;}}---\left(22\right)$

其中由中值定理可知其中，z＝[e₁e₂r]^T；

根据开环误差系统，设计控制输入

$\stackrel{\&CenterDot;}{\mathrm{\&delta;}}\left(t\right)=C^{-1}{A}^{-1}{S}^T(u+\hat{f})---\left(23\right)$

其中，

u＝(K_s+I_3×3)r+βSgn(e₂)(24)

其中，e₂＝[e₂₁e₂₂e₂₃]^T，e₂₁，e₂₂和e₂₃表示e₂向量在滚转、俯仰、偏航三个方向的分量；并且，K_s和β∈R^3×3是对角正定矩阵，I_3×3∈R^3×3是单位矩阵,符号函数为列向量，即Sgn(e₂)＝[Sgn(e₂₁),Sgn(e₂₂),Sgn(e₂₃)]^T；

对式(23)积分，可得控制输入δ(t)表达式如下：

$\mathrm{\&delta;}\left(t\right)=C^{-1}{A}^{-1}\underset{0}{\overset{t}{\&Integral;}}{S}^T\left(\right(K_s+I_{3\&times;3})r+\mathrm{\&beta;}Sgn(e_2)+\hat{f})d\mathrm{\&tau;},---\left(25\right)$

带入r可得：

$\mathrm{\&delta;}\left(t\right)=C^{-1}{A}^{-1}\underset{0}{\overset{t}{\&Integral;}}{S}^T\left(\right(K_s+I_{3\&times;3}\left)\right({\stackrel{\&CenterDot;}{e}}_2+{\mathrm{\&alpha;e}}_2)+\mathrm{\&beta;}Sgn(e_2)+\hat{f})d\mathrm{\&tau;},---\left(26\right)$

其中，作为控制器前馈项，用于改善跟踪性能，减小控制增益β；

将式(16)、式(23)带入开环误差系统(22)，得到闭环误差系统：

$M\left(x\right)\stackrel{\&CenterDot;}{r}=-\frac{1}{2}\stackrel{\&CenterDot;}{M}\left(x\right)r-e_2+\tilde{N}+N_d-(K_s+I_{3\&times;3})r-\mathrm{\&beta;}Sgn\left(e_2\right)---\left(27\right).$