CN105759609A - 一种基于显式非线性模型预测控制的舰载机自主着舰方法 - Google Patents

Abstract

一种基于显式非线性模型预测控制的舰载机自主着舰方法，包括如下几个步骤：步骤一：由飞机所处的飞行情况简化纵向运动非线性方程；步骤二：确定控制量和被控量；步骤三：确定每个输出量的相对阶数；步骤四：K是由K_i组成的对角阵，由矩阵的第一行求得K_i的值。步骤五：对指令信号y_Di进行泰勒展开，得到指令信号的各阶导函数；步骤六：进行积分即可得到控制量。步骤七：得到了一种基于显式非线性模型预测控制的舰载机自主着舰纵向控制律。本发明可实现舰载机自主对于指令信号的快速跟踪，并且避免在线优化过程，节约了大量时间，节省了硬件资源。显式非线性模型预测控制也可以广泛应用于机器人、航空、航天、工业生产等涉及非线性控制的领域。

Description

一种基于显式非线性模型预测控制的舰载机自主着舰方法

【技术领域】

本发明是一种基于显式非线性模型预测控制的舰载机自主着舰方法，属于航空航天控制领域。

【背景技术】

舰载机着舰引导技术是舰载机在较差的海况、甲板运动、大气紊流、舰尾流、夜间或者低能见度下能安全、精准着陆的保障，是航母舰面系统关键技术之一。由于驾驶员对于飞机的熟练度、飞行时间不同引起的差异，着舰时天气、云层、海况、能见度的不同，人工着舰会受到很多因素的影响，为了解决这些问题，产生了自动着舰系统(AutomaticCarrierLandingSystem,ACLS)，一种基于显式非线性模型预测控制的舰载机自主着舰方法能较好的解决着舰时ACLS中指令输入问题，能在干扰下实现较高的着舰精度和在着舰过程中满足较高的指标要求。

模型预测控制可以实现对被控对象的精确控制，但有些模型预测控制依赖在线的优化过程，比较耗时。对于有些线性模型预测控制问题来说，可以求取其最优控制的显式解，但对于非线性问题，求取显式解比较困难。显示非线性模型预测控制(ExplicitNonlinearModelPredictiveControl，ENMPC)应用泰勒展开原理，可以求取非线性模型预测控制的显式解，避免在线优化过程。本发明将ENMPC应用在舰载机自主着舰的控制问题中，在控制过程中，求出指令信号的导函数信息，即可实现舰载机自主对姿态和速度指令的快速准确跟踪。

ENMPC方法应用在舰载机自主着舰过程的控制中。推导ENMPC控制律时需要大量的计算，而舰载机六自由度运动方程比较复杂，给控制律的推导带来很大麻烦。由于舰载机着舰时的主要运动为纵向的运动，横侧向运动幅度相对较小，因此相对纵向状态量而言，横侧向状态量的变化不太大。针对这一特点，本发明将舰载机纵向控制器和横侧向控制器分开来设计，在设计纵向控制器时，将横侧向状态量近似为零，在设计横侧向控制器时，将纵向状态量当做已知量，这样就简化了控制器的设计，降低了设计难度。由于横侧向运动滚转和偏航是耦合的，运动方程比纵向运动方程复杂得多，为了叙述简单，本发明只叙述纵向控制律的设计过程，显示非线性模型预测控制方法可以通过同样的方法应用于舰载机的横侧向控制律设计。

显式非线性模型预测控制方法对系统定义了相对阶数和控制阶数。

考虑下面的一般非线性系统，系统输入为N维，系统输出为P维。

$\left\{\begin{array}{c}\stackrel{\·}{x}=f(x,u)\\ y=h\left(x\right)\end{array}\right.---\left(1\right)$

其中，u为系统输入，x为系统状态变量，y为系统输出

在一个非线性系统中，假设输入和输出量的各阶导数都存在，对于每个系统输出都可以定义一个相对阶数ρ，和一个控制阶数r。其中相对阶数是指，对一个非线性系统的某一输出量求ρ阶导数之后，控制量就会首次出现在该输出的ρ阶导数之中，那么ρ就成为该输出量的相对阶数。当控制量出现之后，如果继续对其求r阶导数，在得到的结果中就会出现控制量r阶导数，这时定义ρ+r为该输出量的控制阶数。选择控制阶数时，还需要与相对阶数结合起来考虑，对于某一特定的相对阶数，其对应的控制阶数必须大于一定的值才能使得控制系统是稳定的，这一点在理论上可以被证明，国外学者给出了保证系统稳定的相对阶数与控制阶数的对应表格,如表1所示。

表1相对阶数、控制阶数与稳定性的关系

针对公式(1)中的非线性系统，设y_i为系统的第i个输出，u_i为系统的第i个输入,ρ_i是第i个输出的相对阶数，y_i的各阶导数如下式所示：

$\left\{\begin{array}{c}\frac{d^j{y}_i}{{\mathrm{dt}}^j}={\mathrm{\&alpha;}}_j^i\left(x\right),0\&le;j<{\mathrm{\&rho;}}_i\\ \frac{d^j{y}_i}{{\mathrm{dt}}^j}={\mathrm{\&alpha;}}_j^i(x,u_{1~N}),j={\mathrm{\&rho;}}_i\\ \frac{d^{{\mathrm{\&rho;}}_i+1}{y}_i}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+1}}={\mathrm{\&alpha;}}_j^i(x,u_{1~N})+{\mathrm{\&beta;}}_{i1}\left(x\right){\stackrel{\&CenterDot;}{u}}_1+\mathrm{...}+{\mathrm{\&beta;}}_{iN}\left(x\right){\stackrel{\&CenterDot;}{u}}_N,j={\mathrm{\&rho;}}_i+1\\ \frac{d^{{\mathrm{\&rho;}}_i+1}{y}_i}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+1}}={\mathrm{\&alpha;}}_j^i(x,u_{1~N},\frac{{\mathrm{du}}_{1~N}}{dt},\mathrm{...},\frac{d^{j-{\mathrm{\&rho;}}_i-1}{u}_{1~N}}{{\mathrm{dt}}^{j-{\mathrm{\&rho;}}_i-1}})+{\mathrm{\&beta;}}_{i1}\left(x\right)\frac{d^{j-{\mathrm{\&rho;}}_i}{u}_1}{{\mathrm{dt}}^{j-{\mathrm{\&rho;}}_i}}+\mathrm{...}+{\mathrm{\&beta;}}_{iN}\left(x\right)\frac{d^{j-{\mathrm{\&rho;}}_i}{u}_N}{{\mathrm{dt}}^{j-{\mathrm{\&rho;}}_i}},j>{\mathrm{\&rho;}}_i+1\end{array}\right.---\left(2\right)$

其中，j为自然数，x为系统状态变量，u_1～N为系统第1～N个输入量，为系统第1～N个输入量的一阶导数，0≤j＜ρ_i时，只是系统状态的函数；j＝ρ_i时，是系统状态和输入的函数；j＞ρ_i时，是系统状态、输入以及输入的1～(j-ρ_i-1)阶导数的函数；β_i1～β_iN都是系统状态的函数。

假设系统输出有足够阶的导数，那么就可以通过泰勒展开对输出量进行近似。所有的输出量的控制阶数都设置为r，输出量y_i利用(ρ_i+r)阶泰勒展开近似，输入量u_i利用r阶泰勒展开近似：

$\left\{\begin{array}{c}{y}_i(t+\mathrm{\&tau;})\&ap;y_i\left(t\right)+\mathrm{\&tau;}\frac{{\mathrm{dy}}_i}{dt}+...+\frac{{\mathrm{\&tau;}}^{{\mathrm{\&rho;}}_i+r}}{({\mathrm{\&rho;}}_i+r)!}\frac{d^{{\mathrm{\&rho;}}_i+r}{y}_i}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+r}}\\ u_i(t+\mathrm{\&tau;})\&ap;u_i\left(t\right)+\mathrm{\&tau;}\frac{{\mathrm{du}}_i}{dt}+...+\frac{{\mathrm{\&tau;}}^r}{r!}\frac{d^r{u}_i}{{\mathrm{dt}}^r}\end{array}\right.---\left(3\right)$

其中，τ为任意自然数，但满足(3)式近似相等关系，为τ的(ρ_i+r)次方，(ρ_i+r)！为(ρ_i+r)的阶乘。

对公式(3)中u_i(t+τ)等式两边先对τ求导再积分可以得到：

$u_i(t+\mathrm{\&tau;})\&ap;\&Integral;(\frac{{\mathrm{du}}_i}{dt}+\mathrm{\&tau;}\frac{d^2{u}_i}{{\mathrm{dt}}^2}+...+\frac{{\mathrm{\&tau;}}^{r-1}}{r!}\frac{d^r{u}_i}{{\mathrm{dt}}^r})d\mathrm{\&tau;}---\left(4\right)$

从上式可以看出，τ时间段内的连续输入量u_i(t+τ)可以用当前输入量的1到r阶导数表示，这就把模型预测控制中的连续优化问题转化为了离散的优化问题。由于在模型需要控制中，需要求解的只是当前τ＝0时刻的控制量，因此τ＝0时公式(4)可以简化为：

$u_i\left(t\right)\&ap;\&Integral;\left(\frac{{\mathrm{du}}_i}{dt}\right)d\mathrm{\&tau;}---\left(5\right)$

由此可以看出，在优化问题中只要求解出u_i的一阶导数即可，控制量u_i可以通过其导数的积分得到。

为书写方便起见，将y_i(t+τ)改写成下面的形式：

y_i(t+τ)≈T_iY_i(6)

其中T_i为行向量，Y_i为列向量，T_i和Y_i如下所示：

$\left\{\begin{array}{c}{T}_i=\left[\begin{array}{ccccc}1& \frac{\mathrm{\&tau;}}{1!}& \frac{{\mathrm{\&tau;}}^2}{2!}& ...& \frac{{\mathrm{\&tau;}}^{{\mathrm{\&rho;}}_i+r}}{({\mathrm{\&rho;}}_i+r)!}\end{array}\right]\\ Y_i={\left[\begin{array}{ccccc}{y}_i& \frac{{\mathrm{dy}}_i}{dt}& \frac{{\mathrm{dy}}_i^2}{{\mathrm{dt}}^2}& ...& \frac{{\mathrm{dy}}_i^{{\mathrm{\&rho;}}_i+r}}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+r}}\end{array}\right]}^T\end{array}\right.---\left(7\right)$

针对多输出的情况，可以将Y_i写在一个列向量里，公式(6)和公式(7)就可以被扩充成：

$y(t+\mathrm{\&tau;})\&ap;\stackrel{\&OverBar;}{T}\stackrel{\&OverBar;}{Y}---\left(8\right)$

为了计算误差函数，需要把指令信号y_D(t+τ)写成跟输出信号y(t+τ)一致的形式：

$y_D(t+\mathrm{\&tau;})\&ap;\stackrel{\&OverBar;}{T}{\stackrel{\&OverBar;}{Y}}_D---\left(10\right)$

将非线性系统输入、输出以及指令信号用泰勒展开近似后，可以定义指令跟踪误差e(t+τ)如下式所示：

$e(t+\mathrm{\&tau;})=y_D(t+\mathrm{\&tau;})-y(t+\mathrm{\&tau;})=\stackrel{\&OverBar;}{T}({\stackrel{\&OverBar;}{Y}}_D-\stackrel{\&OverBar;}{Y})---\left(11\right)$

在模型预测控制中，代价函数J定义为T₁～T₂时间区间上跟踪误差e(t+τ)平方的积分，其中Q是权重矩阵，选择不同的权值可以使得某个输出量在代价函数中占有的比重发生变化。

$J=\frac{1}{2}{({\stackrel{\&OverBar;}{Y}}_D-\stackrel{\&OverBar;}{Y})}^T\left({\&Integral;}_{T_1}^{T_2}\&lsqb;{\stackrel{\&OverBar;}{T}}^{T}Q\stackrel{\&OverBar;}{T}\&rsqb;d\mathrm{\&tau;}\right)({\stackrel{\&OverBar;}{Y}}_D-\stackrel{\&OverBar;}{Y})---\left(12\right)$

在代价函数J的积分计算中，和都是与τ无关的变量，因此可以提到积分号外面，进而代价函数中的积分项只由权重矩阵、控制阶数、相对阶数以及积分时间区间(T₁,T₂)决定，而与系统的状态无关。令(ρ_i+r)＝S_i如下式所示：

前面已经阐述过，将输入量u_i(t+τ)用r阶泰勒展开近似后，u_i(t+τ)就可以用当前τ＝0时刻输入量的1致r阶导数来描述，因此优化问题需要优化的输入量可以写成如下形式：

$\left\{\begin{array}{c}{U}_i={\left[\begin{array}{cccc}\frac{{\mathrm{du}}_i}{dt}& \frac{d^2{u}_i}{{\mathrm{dt}}^2}& ...& \frac{d^r{u}_i}{{\mathrm{dt}}^r}\end{array}\right]}^T\\ \stackrel{\&OverBar;}{U}=\&lsqb;{\begin{array}{cccc}{U}_1& U_2& ...& U_N\end{array}}^T\&rsqb;\end{array}\right.---\left(15\right)$

模型预测控制中，需要解决的优化问题可以表述为：选取一个向量使得在T₁～T₂时间区间上的代价函数J达到最小值。最优解需要满足的条件为：代价函数J对的偏导数为零。

$\frac{\&part;J}{\&part;\stackrel{\&OverBar;}{U}}=-{({\stackrel{\&OverBar;}{Y}}_D-\stackrel{\&OverBar;}{Y})}^T\&Pi;\frac{\&part;\stackrel{\&OverBar;}{Y}}{\&part;\stackrel{\&OverBar;}{U}}=0---\left(16\right)$

通过推导证明可以得到多输入多输出非线性系统ENMPC控制量的显示解的形式为：

${\stackrel{\&CenterDot;}{U}}_c=B^{-1}\&lsqb;K(Y_D^U-Y^U)+Y_D^{L1}-A\&rsqb;---\left(17\right)$

其中，K是由K_i组成的对角阵，K_i为矩阵的第一行，Y^U是系统实际的输出，为指令信号y_Di的(ρ_i+1)阶导数。

其中，

【发明内容】

1、发明目的：

本发明提出了一种基于显式非线性模型预测控制的舰载机自主着舰方法，其目的是提供一种在舰载机着舰时具有较高可靠性的纵向控制方法，也可以应用于横侧向的控制。

该方法采用基于显式非线性模型预测控制的方法。显式非线性模型预测控制应用泰勒展开原理，可以通过求取非线性模型预测控制的显式解，来解决非线性问题，并且避免在线优化过程，节约了大量时间，节省了硬件资源。

2、技术方案：

基于显式非线性模型预测控制的舰载机着舰纵向控制律设计如下：

步骤一：由飞机所处的飞行情况简化纵向运动非线性方程。

由于着舰过程中航迹角μ较小，所以其正弦值可以用弧度值近似；着舰过程中迎角α也比较小，因此这里将推力的方向近似为与速度轴平行；由于升降舵产生的升力很小，因此在升力中忽略升降舵的作用。纵向控制律设计中考虑的系统方程如下：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{h}=V\mathrm{\&mu;}\\ \stackrel{\&CenterDot;}{V}=\&lsqb;T-(C_{D0}+C_{D\mathrm{\&alpha;}}\mathrm{\&alpha;})\stackrel{\&OverBar;}{q}S-mgsin\mathrm{\&mu;}\&rsqb;/m\\ \stackrel{\&CenterDot;}{\mathrm{\&mu;}}=\&lsqb;(C_{L0}+C_{L\mathrm{\&alpha;}}\mathrm{\&alpha;})\stackrel{\&OverBar;}{q}S-mg\cos \mathrm{\&mu;}\&rsqb;/\left(mV\right)\\ \stackrel{\&CenterDot;}{\mathrm{\&alpha;}}=q-\stackrel{\&CenterDot;}{\mathrm{\&mu;}}\\ \stackrel{\&CenterDot;}{q}=\&lsqb;C_{m0}+C_{m\mathrm{\&alpha;}}\mathrm{\&alpha;}+C_{mq}\&CenterDot;c/\left(2V\right)\&CenterDot;q+C_{{\mathrm{m\&delta;}}_e}{\mathrm{\&delta;}}_e\&rsqb;/I_y\end{array}\right.---\left(18\right)$

其中，m为飞机质量，g为重力加速度，T为飞机发动机推力，为飞机所在高度的一阶微分，V为飞机速度，为飞机速度的一阶微分，α为飞机攻角，为飞机攻角的一阶导数，q为飞机俯仰角的角速度，为飞机俯仰角的角加速度，δ_e为飞机的升降舵偏角，为动压，S为飞机机翼的参考面积，c为飞机机翼的平均几何弦长，为飞机航迹倾斜角的一阶导数，C_D0和C_Dα为飞机阻力系数，C_L0和C_Lα为飞机升力系数，C_m0，C_mα，C_mq和为飞机的俯仰力矩系数，I_y为飞机以y轴为转轴的转动惯量。

步骤二：确定控制量和被控量。

纵向控制中的控制量为油门和升降舵，被控量为舰载机速度和高度，横侧向控制中控制量是副翼和方向舵，被控量是滚转角和偏航角。

步骤三：确定每个输出量的相对阶数，并由表1中的数据确定所需的控制阶数，求取被控量y_i的(ρ_i+1)阶对t的导数(也可求ρ_i阶导数，最后通过式(17)直接求得输入量)，并使用纵向运动方程(18)进行变量替换，使y_i只由输入量及其导数和确定的参数表示。

步骤四：K是由K_i组成的对角阵，由矩阵的第一行求得K_i的值。

对于每一个输出量，由相对阶数，控制阶数及预测时间T₁,T₂可以求得：

${\mathrm{\&Pi;}}_i^{12}=\left[\begin{array}{ccc}\frac{T_2^{{\mathrm{\&rho;}}_i+2}-T_1^{{\mathrm{\&rho;}}_i+2}}{({\mathrm{\&rho;}}_i+1)!({\mathrm{\&rho;}}_i+2)}& ...& \frac{T_2^{S_i+1}-T_1^{S_i+1}}{(S_i+1)!}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{T_2^{2{\mathrm{\&rho;}}_i+2}-T_1^{2{\mathrm{\&rho;}}_i+2}}{{\mathrm{\&rho;}}_i!({\mathrm{\&rho;}}_i+1)!(2{\mathrm{\&rho;}}_i+2)}& ...& \frac{T_2^{{\mathrm{\&rho;}}_i+S_i+1}-T_1^{{\mathrm{\&rho;}}_i+S_i+1}}{{\mathrm{\&rho;}}_i!S_i!({\mathrm{\&rho;}}_i+S_i+1)}\end{array}\right]---\left(19\right)$

${\mathrm{\&Pi;}}_i^{22}=\left[\begin{array}{ccc}\frac{T_2^{2{\mathrm{\&rho;}}_i+3}-T_1^{2{\mathrm{\&rho;}}_i+3}}{({\mathrm{\&rho;}}_i+1)!({\mathrm{\&rho;}}_i+1)!(2{\mathrm{\&rho;}}_i+3)}& ...& \frac{T_2^{{\mathrm{\&rho;}}_i+S_i+2}-T_1^{{\mathrm{\&rho;}}_i+S_i+2}}{({\mathrm{\&rho;}}_i+1)!S_i!({\mathrm{\&rho;}}_i+S_i+2)}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{T_2^{{\mathrm{\&rho;}}_i+S_i+2}-T_1^{{\mathrm{\&rho;}}_i+S_i+2}}{({\mathrm{\&rho;}}_i+1)!S_i!({\mathrm{\&rho;}}_i+S_i+2)}& ...& \frac{T_2^{2{\mathrm{\&rho;}}_i+1}-T_1^{2{\mathrm{\&rho;}}_i+1}}{{\mathrm{\&rho;}}_i!{\mathrm{\&rho;}}_i!(2{\mathrm{\&rho;}}_i+1)}\end{array}\right]---\left(20\right)$

(若在步骤三中对被控量y_i求ρ_i阶导数，则会发生变化，需根据之前对Π_i的分块求得)因此可以求得K_i。

步骤五：对指令信号y_Di进行泰勒展开，得到指令信号的各阶导函数。

当系统存在多输出的情况时，公式(10)中的为：

${\stackrel{\&OverBar;}{Y}}_D={\left[\begin{array}{cccc}{Y}_{D1}& Y_{D2}& ...& Y_{Dp}\end{array}\right]}^T$

$Y_{Di}={\left[\begin{array}{ccccc}{y}_{Di}& \frac{{\mathrm{dy}}_{Di}}{dt}& \frac{{\mathrm{dy}}_{Di}^2}{{\mathrm{dt}}^2}& ...& \frac{{\mathrm{dy}}_{Di}^{{\mathrm{\&rho;}}_i+r}}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+r}}\end{array}\right]}^T,$ 其中i＝1,2,...,p。

步骤六：将技术方案中步骤三、四、五所计算得到的被控量y_i的表达式，K的值和指令信号的各阶导函数带入公式(17)，即可解出控制量的一阶导数，进行积分即可得到控制量(若在步骤三中对被控量y_i求ρ_i阶导数，则直接可得到控制量)。

步骤七：至此，得到了一种基于显式非线性模型预测控制的舰载机自主着舰纵向控制律，将该控制律加入飞机的非线性模型进行仿真验证，通过给出着舰时舰载机的下滑高度指令信号，观察在着舰过程中舰载机对于高度指令信号的跟踪情况，进而验证该控制律的效果及稳定性。

3、优点及效果：

本发明提出了一种基于显式非线性模型预测控制的舰载机自主着舰方法，其中显式非线性模型预测控制应用泰勒展开原理，在控制过程中，求出指令信号的导函数信息，可以通过求取非线性模型预测控制的显式解，来解决非线性问题，即可实现舰载机自主对于指令信号的快速跟踪，并且避免在线优化过程，节约了大量时间，节省了硬件资源。显式非线性模型预测控制也可以广泛应用于机器人、航空、航天、工业生产等涉及非线性控制的领域。

【附图说明】

图1为纵向ENMPC算法控制框图。

图2a和图2b为考虑发动机模型时的速度响应曲线和油门推力。

图3a和图3b为斜坡信号下高度响应和升降舵舵偏。

图4为基于显式非线性模型预测控制的舰载机自主着舰方法流程图。

图中序号、符号及代号说明如下：

h——舰载机高度；V——舰载机速度；μ——舰载机航迹倾斜角；

α——舰载机攻角；q——舰载机俯仰角速度；δ_e——舰载机升降舵偏角；

T——舰载机发动机推力；

【具体实施方式】

下面通过一种舰载机F/A-18A为例，利用显式非线性模型预测控制和F/A-18A线性方程来设计其纵向着舰控制律。之后搭建simulink模型，利用F/A-18A模型进行仿真，MATLABR2014a版本。

基于显式非线性模型预测控制的舰载机自主着舰纵向控制律设计如下：

由于油门对速度的控制较为简单，在此不做叙述，只叙述用升降舵对于高度的控制，在设计对高度的控制过程中假设速度由油门控制且已经为恒定，在接下来推导利用升降舵δ_e控制高度h的控制算法时，假设V为常值，即

步骤一：由飞机所处的飞行情况简化纵向运动非线性方程。

如前所述，利用着舰时特点将方程简化为(18)。

步骤二：确定控制量和被控量。

由于油门响应延迟较大，因此使用油门控制速度，使用升降舵控制高度。设第一个控制量为y₁＝V，由推力T控制，系统的第二个输出y₂＝h，由δ_e控制。

步骤三：对于系统的第二个输出y₂＝h，控制量δ_e出现在y₂的四阶导数中，因此定义输出量y₂的相对阶数ρ₂＝4，并由表1中的数据确定所需的控制阶数为r₂＝7，根据公式(18)中的系统状态方程可以得到下式：

$\left\{\begin{array}{c}{y}_2=h\\ {\stackrel{\&CenterDot;}{y}}_2=\stackrel{\&CenterDot;}{h}=V\mathrm{\&mu;}\\ y_2^{\left(2\right)}=V\stackrel{\&CenterDot;}{\mathrm{\&mu;}}=\&lsqb;(C_{L0}+C_{L\mathrm{\&alpha;}}\mathrm{\&alpha;})\stackrel{\&OverBar;}{q}S-mgcos\mathrm{\&mu;}\&rsqb;/m\\ y_2^{\left(3\right)}=(C_{L\mathrm{\&alpha;}}\stackrel{\&OverBar;}{q}S\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}+mgsin\mathrm{\&mu;}\&CenterDot;\stackrel{\&CenterDot;}{\mathrm{\&mu;}})/m\\ y_2^{\left(4\right)}=\&lsqb;C_{L\mathrm{\&alpha;}}\stackrel{\&OverBar;}{q}S(\stackrel{\&CenterDot;}{q}-\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}})+mg\cos \mathrm{\&mu;}\&CenterDot;{\stackrel{\&CenterDot;}{\mathrm{\&mu;}}}^2+mg\sin \mathrm{\&mu;}\&CenterDot;\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}\&rsqb;/m\\ y_2^{\left(5\right)}=\&lsqb;C_{L\mathrm{\&alpha;}}\stackrel{\&OverBar;}{q}S(\stackrel{\&CenterDot;\&CenterDot;}{q}-\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}})-mg\sin \mathrm{\&mu;}\&CenterDot;{\stackrel{\&CenterDot;}{\mathrm{\&mu;}}}^3+2mg\cos \mathrm{\&mu;}\&CenterDot;\stackrel{\&CenterDot;}{\mathrm{\&mu;}}\&CenterDot;\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}+mgcos\mathrm{\&mu;}\&CenterDot;\stackrel{\&CenterDot;}{\mathrm{\&mu;}}\&CenterDot;\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}+mgsin\mathrm{\&mu;}\&CenterDot;\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}\&rsqb;/m\end{array}\right.---\left(21\right)$

其中，控制量δ_e存在于中，控制量的导数存在于中。

$\stackrel{\&CenterDot;\&CenterDot;}{q}=\stackrel{\&OverBar;}{q}S\stackrel{\&OverBar;}{c}\&CenterDot;\&lsqb;C_{m\mathrm{\&alpha;}}\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}+C_{mq}(c/2V)\stackrel{\&CenterDot;}{q}+C_{{\mathrm{m\&delta;}}_e}{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_e\&rsqb;/I_y---\left(22\right)$

航迹角μ的高阶导数可以按照下面的公式计算。

$\left\{\begin{array}{c}\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}=(C_{L\mathrm{\&alpha;}}\stackrel{\&OverBar;}{q}S\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}+mg\sin \mathrm{\&mu;}\&CenterDot;\stackrel{\&CenterDot;}{\mathrm{\&mu;}})/m\\ \stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}=\&lsqb;C_{L\mathrm{\&alpha;}}\stackrel{\&OverBar;}{q}S(\stackrel{\&CenterDot;}{q}-\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}})+mg\cos \mathrm{\&mu;}\&CenterDot;{\stackrel{\&CenterDot;}{\mathrm{\&mu;}}}^2+mg\sin \mathrm{\&mu;}\&CenterDot;\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}\&rsqb;/m\end{array}\right.---\left(23\right)$

步骤四：由矩阵的第一行求得K_i的值。

确定相对阶数，控制阶数及预测时间T₁,T₂可以由公式(19)(20)可以求得K_i,i＝1,2,3,4,5，参数及K_i如下表所示：

步骤五：对指令信号y_Di进行泰勒展开，得到指令信号的各阶导函数。

着舰过程中的高度指令y_D2是一条随时间匀速下降的直线，其一阶导函数 ${\stackrel{\&CenterDot;}{y}}_{D2}=-69.96sin(3.5/57.3)=-4.2731m/s$ 为常值，其二阶及以上导函数均为零。

步骤六：将步骤三、四、五所计算得到的公式或参数带入公式(17)，可得：

$y_2^{\left(5\right)}=\&lsqb;\begin{array}{ccccc}{k}_1& k_2& k_3& k_4& k_5\end{array}\&rsqb;\left[\begin{array}{c}{y}_{D2}-h\\ {\stackrel{\&CenterDot;}{y}}_{D2}-\stackrel{\&CenterDot;}{h}\\ y_{D2}^{\left(2\right)}-h^{\left(2\right)}\\ y_{D2}^{\left(3\right)}-h^{\left(3\right)}\\ y_{D2}^{\left(4\right)}-h^{\left(4\right)}\end{array}\right]+y_{D2}^{\left(5\right)}---\left(24\right)$

求解公式(24)中的方程可以求得升降舵的导数如下式所示：

${\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_e=\frac{\frac{\&lsqb;m\left(K_2\right(Y_{D2}^U-Y_2^U)+g({\stackrel{\&CenterDot;}{\mathrm{\&mu;}}}^{3}sin\mathrm{\&mu;}-3\stackrel{\&CenterDot;}{\mathrm{\&mu;}}\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}cos\mathrm{\&mu;}-\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}sin\mathrm{\&mu;}\left)\right)+C_{L\mathrm{\&alpha;}}\stackrel{\&OverBar;}{q}S\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{\mathrm{\&mu;}}\&rsqb;I_y}{C_{L\mathrm{\&alpha;}}{\stackrel{\&OverBar;}{q}}^2{S}^2\stackrel{\&OverBar;}{c}}-C_{m\mathrm{\&alpha;}}\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}-\frac{C_{mq}\stackrel{\&OverBar;}{c}\stackrel{\&CenterDot;}{q}}{2V}}{C_{{\mathrm{m\&delta;}}_e}}---\left(25\right)$

对(25)的结果加入积分环节，即可得到δ_e。

对推力T的设计过程和上述步骤相同，而且更加简单。

步骤七：将得到的控制律加入未化简的模型进行仿真验证。

Simulink模型框图如图1所示。

由于设计控制算法时使用了简化方法，因此之前存在建模误差，希望证明在存在简化计算的情况下，ENMPC控制方法依然能够满足控制需求。首先对速度控制进行验证，在加入发动机模型的情况下进行仿真。由于设计控制律时没有考虑发动机动态模型，而且发动机存在较大的延迟，因此加入该模型会对速度响应造成一定的影响，但从图2a和图2b的仿真结果中可以看出，加入发动机模型的速度响应依然可以满足控制需求。

对于无人机高度控制进行仿真时，由于在推导高度控制的控制律时，假设速度为常值，因此在应用ENMPC控制器对无人机高度进行控制时，需要应用上面介绍的速度控制器将速度保持为常值。在本仿真中，无人机空速被设定为69.96m/s，为了验证高度控制算法的有效性，本仿真对控制器输入高度斜坡信号，控制器采用了之前表中的控制阶数、预测时间以及控制增益，仿真结果如图3a和图3b所示。

由图3a和图3b可知，利用ENMPC方法设计的无人机着舰高度控制器对阶跃指令和斜坡指令都有很好的控制效果，而且升降舵偏转的幅值和速率也在就可以接受的范围内。上述仿真中考虑的升降舵的动态特性、饱和以及速率限制，同时也考虑了发动机模型的动态特性和饱和限制，仿真结果满足设计要求。

Claims (2)

1.一种基于显式非线性模型预测控制的舰载机自主着舰方法，其中，显式非线性模型预测控制方法对非线性系统定义了相对阶数和控制阶数；

非线性系统输入为N维，非线性系统输出为P维；

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=f(x,u)\\ y=h\left(x\right)\end{array}\right.---\left(1\right)$

其中，u为非线性系统输入，x为非线性系统状态变量，y为非线性系统输出；

在一个非线性系统中，假设输入和输出量的各阶导数都存在，对于每个非线性系统输出都定义一个相对阶数ρ，和一个控制阶数r；其中相对阶数是指，对一个非线性系统的某一输出量求ρ阶导数之后，控制量就会首次出现在该输出的ρ阶导数之中，那么ρ就成为该输出量的相对阶数；当控制量出现之后，如果继续对其求r阶导数，在得到的结果中就会出现控制量r阶导数，定义ρ+r为该输出量的控制阶数；选择控制阶数时，还需要与相对阶数结合起来考虑，对于某一特定的相对阶数，其对应的控制阶数必须大于一定的值才能使得控制非线性系统是稳定的；保证非线性系统稳定的相对阶数与控制阶数的对应表格,如表1所示；表1为相对阶数、控制阶数与稳定性的关系；

针对公式(1)中的非线性系统，设y_i为非线性系统的第i个输出，u_i为非线性系统的第i个输入,ρ_i是第i个输出的相对阶数，y_i的各阶导数如下式所示：

$\left\{\begin{array}{c}\frac{d^j{y}_i}{{\mathrm{dt}}^j}={\mathrm{\&alpha;}}_j^i\left(x\right),0\&le;j<{\mathrm{\&rho;}}_i\\ \frac{d^j{y}_i}{{\mathrm{dt}}^j}={\mathrm{\&alpha;}}_j^i(x,u_{1~N}),j={\mathrm{\&rho;}}_i\\ \frac{d^{{\mathrm{\&rho;}}_i+1}{y}_i}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+1}}={\mathrm{\&alpha;}}_j^i(x,u_{1~N})+{\mathrm{\&beta;}}_{i1}\left(x\right){\stackrel{\&CenterDot;}{u}}_1+\mathrm{...}+{\mathrm{\&beta;}}_{iN}\left(x\right){\stackrel{\&CenterDot;}{u}}_N,j={\mathrm{\&rho;}}_i+1\\ \frac{d^{{\mathrm{\&rho;}}_i+1}{y}_i}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+1}}={\mathrm{\&alpha;}}_j^i(x,u_{1~N},\frac{{\mathrm{du}}_{1~N}}{dt},\mathrm{...},\frac{d^{j-{\mathrm{\&rho;}}_i-1}{u}_{1~N}}{{\mathrm{dt}}^{j-{\mathrm{\&rho;}}_i-1}})+{\mathrm{\&beta;}}_{i1}\left(x\right)\frac{d^{j-{\mathrm{\&rho;}}_i}{u}_1}{{\mathrm{dt}}^{j-{\mathrm{\&rho;}}_i}}+\mathrm{...}+{\mathrm{\&beta;}}_{iN}\left(x\right)\frac{d^{j-{\mathrm{\&rho;}}_i}{u}_N}{{\mathrm{dt}}^{j-{\mathrm{\&rho;}}_i}},j>{\mathrm{\&rho;}}_i+1\end{array}\right.---\left(2\right)$

其中，j为自然数，x为非线性系统状态变量，u_1～N为非线性系统第1～N个输入量，为非线性系统第1～N个输入量的一阶导数，0≤j＜ρ_i时，只是非线性系统状态的函数；j＝ρ_i时，是非线性系统状态和输入的函数；j＞ρ_i时，是非线性系统状态、输入以及输入的1～(j-ρ_i-1)阶导数的函数；β_i1～β_iN都是非线性系统状态的函数；

假设非线性系统输出有足够阶的导数，那么就通过泰勒展开对输出量进行近似；所有的输出量的控制阶数都设置为r，输出量y_i利用(ρ_i+r)阶泰勒展开近似，输入量u_i利用r阶泰勒展开近似：

$\left\{\begin{array}{c}{y}_i(t+\mathrm{\&tau;})\&ap;y_i\left(t\right)+\mathrm{\&tau;}\frac{{\mathrm{dy}}_i}{dt}+...+\frac{{\mathrm{\&tau;}}^{{\mathrm{\&rho;}}_i+r}}{({\mathrm{\&rho;}}_i+r)!}\frac{d^{{\mathrm{\&rho;}}_i+r}{y}_i}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+r}}\\ u_i(t+\mathrm{\&tau;})\&ap;u_i\left(t\right)+\mathrm{\&tau;}\frac{{\mathrm{du}}_i}{dt}+...+\frac{{\mathrm{\&tau;}}^r}{r!}\frac{d^r{u}_i}{{\mathrm{dt}}^r}\end{array}\right.---\left(3\right)$

其中，τ为任意自然数，但满足(3)式近似相等关系，为τ的(ρ_i+r)次方，(ρ_i+r)！为(ρ_i+r)的阶乘；

对公式(3)中u_i(t+τ)等式两边先对τ求导再积分得到：

$u_i(t+\mathrm{\&tau;})\&ap;\&Integral;(\frac{{\mathrm{du}}_i}{dt}+\mathrm{\&tau;}\frac{d^2{u}_i}{{\mathrm{dt}}^2}+...+\frac{{\mathrm{\&tau;}}^{r-1}}{r!}\frac{d^r{u}_i}{{\mathrm{dt}}^r})d\mathrm{\&tau;}---\left(4\right)$

从上式看出，τ时间段内的连续输入量u_i(t+τ)用当前输入量的1到r阶导数表示，这就把模型预测控制中的连续优化问题转化为了离散的优化问题；由于在模型需要控制中，需要求解的只是当前τ＝0时刻的控制量，因此τ＝0时公式(4)简化为：

$u_i\left(t\right)\&ap;\&Integral;\left(\frac{{\mathrm{du}}_i}{dt}\right)d\mathrm{\&tau;}---\left(5\right)$

在优化问题中只要求解出u_i的一阶导数即可，控制量u_i通过其导数的积分得到；

为了简便，将y_i(t+τ)改写成下面的形式：

y_i(t+τ)≈T_iY_i(6)

其中T_i为行向量，Y_i为列向量，T_i和Y_i如下所示：

$\left\{\begin{array}{c}{T}_i=\&lsqb;\begin{array}{ccccc}1& \frac{\mathrm{\&tau;}}{1!}& \frac{{\mathrm{\&tau;}}^2}{2!}& \mathrm{...}& \frac{{\mathrm{\&tau;}}^{{\mathrm{\&rho;}}_i+r}}{({\mathrm{\&rho;}}_i+r)!}\end{array}\&rsqb;\\ Y_i={\left[\begin{array}{ccccc}{y}_i& \frac{{\mathrm{dy}}_i}{dt}& \frac{{\mathrm{dy}}_i^2}{{\mathrm{dt}}^2}& \mathrm{...}& \frac{{\mathrm{dy}}_i^{{\mathrm{\&rho;}}_i+r}}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+r}}\end{array}\right]}^T\end{array}\right.---\left(7\right)$

针对多输出的情况，将Y_i写在一个列向量里，公式(6)和公式(7)就能被扩充成：

$y(t+\mathrm{\&tau;})\&ap;\stackrel{\&OverBar;}{T}\stackrel{\&OverBar;}{Y}---\left(8\right)$

为了计算误差函数，需要把指令信号y_D(t+τ)写成跟输出信号y(t+τ)一致的形式：

$y_D(t+\mathrm{\&tau;})\&ap;\stackrel{\&OverBar;}{T}{\stackrel{\&OverBar;}{Y}}_D---\left(10\right)$

将非线性系统输入、输出以及指令信号用泰勒展开近似后，定义指令跟踪误差e(t+τ)如下式所示：

$e(t+\mathrm{\&tau;})=y_D(t+\mathrm{\&tau;})-y(t+\mathrm{\&tau;})=\stackrel{\&OverBar;}{T}({\stackrel{\&OverBar;}{Y}}_D-\stackrel{\&OverBar;}{Y})---\left(11\right)$

在模型预测控制中，代价函数J定义为T₁～T₂时间区间上跟踪误差e(t+τ)平方的积分，其中Q是权重矩阵，选择不同的权值使得某个输出量在代价函数中占有的比重发生变化；

$J=\frac{1}{2}{({\stackrel{\&OverBar;}{Y}}_D-\stackrel{\&OverBar;}{Y})}^T\left({\&Integral;}_{T_1}^{T_2}\&lsqb;{\stackrel{\&OverBar;}{T}}^{T}Q\stackrel{\&OverBar;}{T}\&rsqb;d\mathrm{\&tau;}\right)({\stackrel{\&OverBar;}{Y}}_D-\stackrel{\&OverBar;}{Y})---\left(12\right)$

在代价函数J的积分计算中，和都是与τ无关的变量，因此提到积分号外面，进而代价函数中的积分项只由权重矩阵、控制阶数、相对阶数以及积分时间区间(T₁,T₂)决定，而与非线性系统的状态无关；令(ρ_i+r)＝S_i如下式所示：

将输入量u_i(t+τ)用r阶泰勒展开近似后，u_i(t+τ)用当前τ＝0时刻输入量的1致r阶导数来描述，因此优化问题需要优化的输入量写成如下形式：

$\left\{\begin{array}{c}{U}_i={\left[\begin{array}{cccc}\frac{{\mathrm{du}}_i}{dt}& \frac{d^2{u}_i}{{\mathrm{dt}}^2}& \mathrm{...}& \frac{d^r{u}_i}{{\mathrm{dt}}^r}\end{array}\right]}^T\\ \stackrel{\&OverBar;}{U}={\left[\begin{array}{cccc}{U}_1& U_2& \mathrm{...}& U_N\end{array}\right]}^T\end{array}\right.---\left(15\right)$

模型预测控制中，需要解决的优化问题表述为：选取一个向量使得在T₁～T₂时间区间上的代价函数J达到最小值；最优解需要满足的条件为：代价函数J对的偏导数为零；

$\frac{\&part;J}{\&part;\stackrel{\&OverBar;}{U}}=-{({\stackrel{\&OverBar;}{Y}}_D-\stackrel{\&OverBar;}{Y})}^T\mathrm{\&Pi;}\frac{\&part;\stackrel{\&OverBar;}{Y}}{\&part;\stackrel{\&OverBar;}{U}}=0---\left(16\right)$

通过推导证明得到多输入多输出非线性系统ENMPC控制量的显示解的形式为：

${\stackrel{\&CenterDot;}{U}}_c=B^{-1}\&lsqb;K(Y_D^U-Y^U)+Y_D^{L1}-A\&rsqb;---\left(17\right)$

其中，K是由K_i组成的对角阵，K_i为矩阵的第一行，Y^U是非线性系统实际的输出，为指令信号y_Di的(ρ_i+1)阶导数；

其中，

$Y_D^U=\left[\begin{array}{c}{Y}_{D1}^U\\ .\\ .\\ .\\ Y_{\mathrm{DP}}^U\end{array}\right],Y^U=\left[\begin{array}{c}{Y}_1^U\\ .\\ .\\ .\\ Y_P^U\end{array}\right],Y_D^{L1}=\left[\begin{array}{c}{y}_{D1}^{({\mathrm{\&rho;}}_1+1)}\\ .\\ .\\ .\\ y_{D_P}^{({\mathrm{\&rho;}}_P+1)}\end{array}\right],A=\left[\begin{array}{c}{\mathrm{\&alpha;}}_{{\mathrm{\&rho;}}_1+1}^1\\ .\\ .\\ .\\ {\mathrm{\&alpha;}}_{{\mathrm{\&rho;}}_P+1}^P\end{array}\right],{\stackrel{.}{U}}_c=\left[\begin{array}{c}{\stackrel{.}{u}}_1\\ .\\ .\\ .\\ {\stackrel{.}{u}}_N\end{array}\right]$

其特征在于，具体步骤如下：

步骤一：由飞机所处的飞行情况简化纵向运动非线性方程；

由于着舰过程中航迹角μ很小，所以其正弦值用弧度值近似；着舰过程中迎角α也很小，因此将推力的方向近似为与速度轴平行；由于升降舵产生的升力很小，因此在升力中忽略升降舵的作用；纵向控制律非线性系统方程如下：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{h}=V\mathrm{\&mu;}\\ \stackrel{\&CenterDot;}{V}=\&lsqb;T-(C_{D0}+C_{D\mathrm{\&alpha;}}\mathrm{\&alpha;})\stackrel{\&OverBar;}{q}S-mgsin\mathrm{\&mu;}\&rsqb;/m\\ \stackrel{\&CenterDot;}{\mathrm{\&mu;}}=\&lsqb;(C_{L0}+C_{L\mathrm{\&alpha;}}\mathrm{\&alpha;})\stackrel{\&OverBar;}{q}S-mg\cos \mathrm{\&mu;}\&rsqb;/\left(mV\right)\\ \stackrel{\&CenterDot;}{\mathrm{\&alpha;}}=q-\stackrel{\&CenterDot;}{\mathrm{\&mu;}}\\ \stackrel{\&CenterDot;}{q}=\&lsqb;C_{m0}+C_{m\mathrm{\&alpha;}}\mathrm{\&alpha;}+C_{mq}\&CenterDot;c/\left(2V\right)\&CenterDot;q+C_{{\mathrm{m\&delta;}}_e}{\mathrm{\&delta;}}_e\&rsqb;/I_y\end{array}\right.---\left(18\right)$

其中，m为飞机质量，g为重力加速度，T为飞机发动机推力，为飞机所在高度的一阶微分，V为飞机速度，为飞机速度的一阶微分，α为飞机攻角，为飞机攻角的一阶导数，q为飞机俯仰角的角速度，为飞机俯仰角的角加速度，δ_e为飞机的升降舵偏角，为动压，S为飞机机翼的参考面积，c为飞机机翼的平均几何弦长，为飞机航迹倾斜角的一阶导数，C_D0和C_Dα为飞机阻力系数，C_L0和C_Lα为飞机升力系数，C_m0，C_mα，C_mq和为飞机的俯仰力矩系数，I_y为飞机以y轴为转轴的转动惯量；

步骤二：确定控制量和被控量；

纵向控制中的控制量为油门和升降舵，被控量为舰载机速度和高度，横侧向控制中控制量是副翼和方向舵，被控量是滚转角和偏航角；

步骤三：确定每个输出量的相对阶数，并由表1中的数据确定所需的控制阶数，求取被控量y_i的(ρ_i+1)阶对t的导数，并使用纵向运动方程(18)进行变量替换，使y_i只由输入量及其导数和确定的参数表示；

步骤四：K是由K_i组成的对角阵，由矩阵的第一行求得K_i的值；

对于每一个输出量，由相对阶数，控制阶数及预测时间T₁,T₂求得：

${\mathrm{\&Pi;}}_i^{12}=\left[\begin{array}{ccc}\frac{T_2^{{\mathrm{\&rho;}}_i+2}-T_1^{{\mathrm{\&rho;}}_i+2}}{({\mathrm{\&rho;}}_i+1)!({\mathrm{\&rho;}}_i+2)!}& \mathrm{...}& \frac{T_2^{S_i+1}-T_1^{S_i+1}}{(S_i+1)!}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{T_2^{2{\mathrm{\&rho;}}_i+2}-T_1^{2{\mathrm{\&rho;}}_i+2}}{{\mathrm{\&rho;}}_i!({\mathrm{\&rho;}}_i+1)!(2{\mathrm{\&rho;}}_i+2)}& \mathrm{...}& \frac{T_2^{{\mathrm{\&rho;}}_i+S_i+1}-T_1^{{\mathrm{\&rho;}}_i+S_i+1}}{{\mathrm{\&rho;}}_i!S_i!({\mathrm{\&rho;}}_i+S_i+1)}\end{array}\right]---\left(19\right)$

${\mathrm{\&Pi;}}_i^{22}=\left[\begin{array}{ccc}\frac{T_2^{2{\mathrm{\&rho;}}_i+3}-T_1^{2{\mathrm{\&rho;}}_i+3}}{({\mathrm{\&rho;}}_i+1)!({\mathrm{\&rho;}}_i+1)!(2{\mathrm{\&rho;}}_i+3)}& \mathrm{...}& \frac{T_2^{{\mathrm{\&rho;}}_i+S_i+2}-T_1^{{\mathrm{\&rho;}}_i+S_i+2}}{({\mathrm{\&rho;}}_i+1)!S_i!({\mathrm{\&rho;}}_i+S_i+2)}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{T_2^{{\mathrm{\&rho;}}_i+S_i+2}-T_1^{{\mathrm{\&rho;}}_i+S_i+2}}{({\mathrm{\&rho;}}_i+1)!S_i!({\mathrm{\&rho;}}_i+S_i+2)}& \mathrm{...}& \frac{T_2^{2{\mathrm{\&rho;}}_i+1}-T_1^{2{\mathrm{\&rho;}}_i+1}}{{\mathrm{\&rho;}}_i!{\mathrm{\&rho;}}_i!(2{\mathrm{\&rho;}}_i+1)}\end{array}\right]---\left(20\right)$

步骤五：对指令信号y_Di进行泰勒展开，得到指令信号的各阶导函数；

当非线性系统存在多输出的情况时，公式(10)中的为：

${\stackrel{\&OverBar;}{Y}}_D={\left[\begin{array}{cccc}{Y}_{D1}& Y_{D2}& \mathrm{...}& Y_{Dp}\end{array}\right]}^T$

$Y_{Di}={\left[\begin{array}{ccccc}{y}_{Di}& \frac{{\mathrm{dy}}_{Di}}{dt}& \frac{{\mathrm{dy}}_{Di}^2}{{\mathrm{dt}}^2}& \mathrm{...}& \frac{{\mathrm{dy}}_{Di}^{{\mathrm{\&rho;}}_i+r}}{{\mathrm{dt}}^{{\mathrm{\&rho;}}_i+r}}\end{array}\right]}^T,$ 其中i＝1,2,...,p；

步骤六：将技术方案中步骤三、步骤四、步骤五所计算得到的被控量y_i的表达式，K的值和指令信号的各阶导函数带入公式(17)，即解出控制量的一阶导数，进行积分即得到控制量；

步骤七：由此，得到了一种基于显式非线性模型预测控制的舰载机自主着舰纵向控制律，将该控制律加入飞机的非线性模型进行仿真验证，通过给出着舰时舰载机的下滑高度指令信号，观察在着舰过程中舰载机对于高度指令信号的跟踪情况，进而验证该控制律的效果及稳定性。

2.根据权利要求1所述的一种基于显式非线性模型预测控制的舰载机自主着舰方法，其特征在于：在步骤三中，确定每个输出量的相对阶数，由表1中的数据确定所需的控制阶数，求ρ_i阶导数，最后通过式(17)直接求得输入量，使y_i只由输入量及其导数和确定的参数表示。