专利名称：一种空天飞机高空爬升纵向控制方法
技术领域：
本发明属航空航天控制领域，涉及一种空天飞机在高空爬升时的纵向控制方法。
背景技术：
空天飞机的概念于20世纪80年代初提出。随着空间活动的增加，特别是载人航天的发展，一次性使用火箭、飞船和航天飞机的高额发射费用日益成为大规模开展空间活动的“瓶颈”，迫切需要一种既能像普通飞机一样起降又能往返于天地之间的经济、安全的飞行器，这就是空天飞机。它既能完成民用航空航天任务，又能执行多种军事航空航天任务，将是21世纪控制空间、争夺制天权的关键武器装备之一，是一种具有广阔发展前景的载人航天兵器，将为未来战争带来重大变革。
空天飞机由载机携带到一定高度进行投放，投放后，火箭发动机点火，空天飞机迅速拉起，进行迅速爬升，当燃油消耗完后，发动机停机，进入无动力爬升阶段，达到预定高度和轨道。在大气层内爬升飞行时，由于巨大的发动机推力，空天飞机在爬升过程中爬升迅速，飞行速度和高度急剧变化，从亚音速到高超音速，使空天飞机的气动特性产生急剧变化，空天飞机的升力特性、焦点及操纵面效率都会产生急剧变化。为了达到预定的高度，空天飞机需携带大量的燃料，这些燃料在爬升过程中快速消耗掉。因此会引起爬升过程中空天飞机的重量，重心和惯矩的急剧而大幅度的变化。这些变化会对空天飞机的操纵稳定特性和动态响应特性产生较大的影响。
对于高空爬升段空天飞机的控制律设计问题，通过对各项数据进行详细分析，可以看出爬升段设计中主要面临如下困难(1)高度和速度协调控制问题。由于巨大的发动机推力，空天飞机速度上升很快。爬升段控制目标为飞行高度超过预先制定的高度范围，飞行速度则应小于机体结构限制的最大马赫数要求。因此如何实现高度和速度协调控制是设计的一个难点。(2)长周期不可控问题。在高空高速的范围内，由于空气密度急剧减小，导致速度变量一直是快速增大的，且爬升角难以控制。即在此范围内，外环长周期不能控的问题。(3)如何针对短周期变量，设计合适的控制目标，如攻角和俯仰角速率跟踪曲线，以实现爬升段的控制，也是面临的问题之一。
从上面的分析可以看出，空天飞机高空爬升段与传统飞机控制模式是不同的，但公开发表的文献尚没有相关控制方法。因此，有必要寻找一种新的方法来完成空天飞机高空爬升控制，解决上述三个难点，即设计攻角和俯仰角速率跟踪曲线和适于工程应用的控制方法，以实现速度和高度的协调控制，并且满足俯仰角和过载限制，克服飞机控制中的调参问题，从而降低设计的复杂程度。
发明内容
本发明的技术解决问题是克服现有技术的不足，提供了一种可协调控制性好、满足俯仰角和过载限制、可进行六自由度控制的空天飞机高空爬升纵向控制方法。
本发明的技术解决方案是一种空天飞机高空爬升纵向控制方法，步骤如下 (1)按照空天飞机动力学方程和高空气动特点，确定空天飞机高空爬升速度上升曲线； (2)根据步骤(1)得到的速度上升曲线，得到爬升角γ的跟踪曲线； (3)根据步骤(2)得到的爬升角γ的跟踪曲线确定攻角α的跟踪曲线； (4)确定俯仰角速率q的跟踪曲线； (5)根据步骤(3)和(4)的结果，确定攻角α和俯仰角速率q的跟踪控制律。
所述步骤(1)中得到的空天飞机高空爬升速度上升曲线方程为 其中，V为空天飞机的爬升速度，P为空天分机受到的推力，m为空天飞机的质量，μ为地球引力常数，r为空天飞机距地心的距离，γ为空天飞机的爬升角。
所述步骤(2)中得到的爬升角γ的跟踪曲线方程为 其中，γ为空天飞机的爬升角，P为空天分机受到的推力，m为空天飞机的质量，μ为地球引力常数，r为空天飞机距地心的距离，V为空天飞机的爬升速度，a为声速。
所述步骤(3)中得到的攻角α的跟踪曲线方程为 其中，α为攻角，γ为空天飞机的爬升角。
所述步骤(4)中俯仰角速率q的跟踪曲线方程为 q＝0。
所述步骤(5)中攻角α和俯仰角速率q的跟踪控制律采用基于特征模型的全系数自适应控制方法得到，步骤为 首先分别建立攻角α和俯仰角速率q的特征模型如下 α(k+2)＝f11(k)α(k+1)+f12(k)α(k)+g1(k)u(k) q(k+2)＝f21(k)q(k+1)+f22(k)q(k)+g2(k)u(k) 式中，fi1(k)∈[1.4331，1.9974]，i＝1，2；fi2(k)∈[-0.9999，-0.5134]，i＝1，2； 然后确定gi(k)，i＝1，2的范围，使用最小二乘方法确定式中的6个参数，和 最后，利用得到的和设计全系数自适应控制律uj＝ugj+uij+udj，j＝1，2， 式中， uij(k)＝uij(k-1)-kijej(k)，j＝1，2 udj(k)＝-kdj(ej(k)-ej(k-1))，j＝1，2 l1＝0.382，l2＝0.618 e1(k)＝α(k)-αr(k)，α(k)为状态变量，αr(k)表示攻角跟踪曲线； e2(k)＝q(k)-qr(k)，q(k)为状态变量，qr(k)表示俯仰角速度跟踪曲线； kijl，j，l＝1，2，cdj，ldj，j＝1，2为需调试参数。
本发明与现有技术相比的优点在于 (1)本发明的控制方法包括跟踪曲线设计和跟踪控制方法设计两部分，充分考虑了被控对象本身的动力学特点和高空气动特性，不仅实现了协调控制的目标，而且满足俯仰角和过载限制，为进一步六自由度高空爬升控制设计提供了思路； (2)本发明方法通过仔细分析高空被控对象动力学方程的各项数据，明确了飞行器长周期变量高空不可控的问题，揭示了高空飞行器飞行的本质。根据被控对象自身特点得到了速度上升曲线，为空天飞机协调控制奠定了基础； (3)本发明方法根据空天飞机高空被控对象动力学和协调控制特点，给出了爬升角计算公式，为跟踪曲线设计奠定了理论依据； (4)本发明方法基于空天飞机动力学方程和飞机控制的本质特点，给出了攻角计算公式，为攻角跟踪曲线设计提供了理论依据； (5)本发明方法根据飞机控制本质要求，明确给出了俯仰角速率跟踪曲线，为在高空空天飞机内环控制奠定了基础； (6)本发明设计采用基于特征模型的全系数自适应控制方法，克服了现有自适应理论在应用中的不足，具有较强的鲁棒性和对初始参数的适应性，适于工程应用，克服了现有飞机控制中的调参问题，降低了设计的复杂性。
图1为本发明方法的流程框图； 图2为本发明实施例中空天飞机重心的变化规律图； 图3为本发明实施例中空天飞机质量的变化规律图； 图4为本发明实施例中空天飞机惯矩的变化规律图； 图5为本发明实施例中依本发明方法所设计的速度曲线与直线上升曲线的比较示意图； 图6为本发明实施例中依本发明方法所得到的爬升角γ跟踪曲线； 图7为本发明实施例中依本发明方法所得到的攻角α跟踪曲线； 图8为采用图7的攻角α跟踪曲线得到的跟踪结果； 图9为本发明实施例中俯仰角速率q的跟踪结果； 图10为本发明实施例中的高度仿真结果； 图11为本发明实施例中的马赫数仿真结果； 图12为本发明实施例中的俯仰角仿真结果； 图13为本发明实施例中的爬升角仿真结果； 图14为本发明实施例中的控制输入仿真结果； 图15为本发明实施例中的过载仿真结果； 图16为本发明实施例中的攻角跟踪结果； 图17为本发明实施例中的俯仰角速度跟踪结果。
具体实施例方式 如图1所示，为本发明方法的流程框图，步骤为(1)按照空天飞机动力学方程和高空气动特点，确定空天飞机高空爬升速度上升曲线；(2)根据速度上升曲线，计算爬升角γ的跟踪曲线；(3)根据爬升角γ的跟踪曲线确定攻角α的跟踪曲线；(4)设计俯仰角速率q的跟踪曲线；(5)根据步骤(3)和(4)的结果，设计攻角α和俯仰角速率q的跟踪控制律。
本发明中，速度相关动力学方程采用其中V，P，r，m，γ，μ，α，D分别表示空天飞机的速度、推力、距地心距离、质量、爬升角、地球引力常数、攻角和外部阻力。由于高空空气密度非常小，导致
的第三项很小，可忽略不计，而cosα约为1，因此
可简化为 综合上述推导，可以得到高空速度上升曲线计算公式为其中γ按照工作点角度计算。
考虑到声速a的变化率远小于速度的变化率，由其中M表示马赫数，可得由于空天飞机爬升控制目标要求马赫数和高度具有近似相同的上升斜率，因此有由以上两式可得又根据动力学方程可以得到考虑到上述的空天飞机高空爬升速度上升曲线方程可以得到 综合上述推导，可以得到γ角计算公式为 按照飞机控制的特点，设计俯仰角速率跟踪曲线q＝0。进而由动力学方程以及上述γ的跟踪曲线可以得到攻角α的跟踪曲线。
本发明中，攻角α和俯仰角速率q的跟踪控制律采用基于特征模型的全系数自适应控制方法得到，首先分别建立攻角和俯仰角速度的特征模型如下 α(k+2)＝f11(k)α(k+1)+f12(k)α(k)+g1(k)u(k) q(k+2)＝f21(k)q(k+1)+f22(k)q(k)+g2(k)u(k) 其中fi1(k)∈[1.4331，1.9974]，i＝1，2；fi2(k)∈[-0.9999，-0.5134]，i＝1，2；适当选取gi(k)，i＝1，2的范围，使用最小二乘方法辨识上述6个参数，辨识结果分别记为和利用上述参数设计全系数自适应控制律 uj＝ugj+uij+udj，j＝1，2 其中， uij(k)＝uij(k-1)-kijej(k)，j＝1，2 udj(k)＝-kdj(ej(k)-ej(k-1))，j＝1，2 l1＝0.382，l2＝0.618 e1(k)＝α(k)-αr(k)，αr(k)表示攻角跟踪曲线； e2(k)＝q(k)-qr(k)，qr(k)表示俯仰角速度跟踪曲线； kijl，j，l＝1，2，cdj，ldj，j＝1，2，为需调试参数。(详细推导过程可参见吴宏鑫，全系数自适应控制理论及其方法，北京国防工业出版社，1990) 实施例 下面以空间飞行器X-34的高空纵向爬升控制为例详细说明本发明方法。
空天飞行器X-34由载机携带到高度7公里(0.7M)进行投放，飞行器投放后，火箭发动机点火，飞行器迅速拉起，进行迅速爬升，到72公里(7.2M)，燃油消耗完后，发动机停机，进入无动力爬升阶段，达到预定高度和轨道，最后进行无动力下降。这里考虑高空爬升段25公里(2.5M)到72公里(7.2M)的控制律设计问题。对于上述动力学方程中的各项气动参数CLB，CDB，CMB，CLDE，CDDE，CMDE，CMQ，采用相关参考文献中的气动数据(此处给出两篇文献1.Aerodynamic Characteristics and Development of theAerodynamic Database of the X-34 Reusable Launch Vehicle.2.Aerodynamic Characteristics，Database Development and Flight Simulation of the X-34 Vehicle.)。在爬升过程中，由于巨大的发动机推力，飞行器重心、质量和惯矩随时间剧烈变化，其规律如图2～4所示。空天飞机控制舵面为升降舵和襟翼。
跨大气层爬升段的任务目标是在经历动力上升和无动力上升后，因热防护材料的因素规定(这里体现协调控制的必要性，若飞行器在较低的高度以较大速度飞行，则气动热较高，若在较高的高度，则气动效应差)，飞行高度应能超过预先制定的高度范围，飞行速度则应小于机体结构限制的最大马赫数要求。除此之外，要求在200秒时空天飞机到达72公里(7.2M)，且满足过载要求。
在圆形旋转地球情形下，飞行器纵向动力学方程为 r＝h+RE CL＝CLB+CLDE CD＝CDB+CDDE CZ＝-CLcosα-CDsinα 上述各个变量的物理意义为 V航迹速度(飞行器质心相对大地的速度) h飞行器高度 γ爬升角 m飞行器质量 Iy绕俯仰轴转动惯量 P推力 q俯仰角速率 α攻角 μ地球引力常数(3.986005×1014m3/s2) ωE地球旋转角速度(7.292116×10-5rad/s) RE地球半径，6371386m M马赫数 L，D，My分别表示气动升力、气动阻力和俯仰力矩 ρ空气密度 S飞行器参考面积(33.2m2) c飞行器参考弦长(4.43m) CL，CD，CM分别表示升力、阻力和俯仰力矩系数 CZ本体坐标系Z方向总的气动力系数 CLB、CDB、CMB分别表示升力系数、阻力系数和俯仰力矩系数基本量 CLDE、CDDE、CMDE分别表示升降舵产生的气动升力、阻力以及俯仰力矩系数 CMQ俯仰阻尼系数 XCG重心位置 CGREF参考重心(10.668m) qs动压 通过分析动力学方程(1)中的各项数据μ/r2、P/m和qs/m可以看出，在高空时速度是不可控的，为了按照设计时间到达既定速度和高度，必须按照飞行器自身特点从后向前进行设计，所设计的速度上升曲线如图5所示。图5中圆圈组成的线为所设计的速度上升曲线。
由动力学方程(3)可知，如果速度已按照要求达到控制目标，那么通过设计合适的γ角，即可实现高度的设计。因此实现高度和速度协调控制的关键是γ角的控制问题。下面通过分析动力学方程及其各项数据，沿着上述速度上升曲线得到了γ角的跟踪曲线，计算公式如下 在40公里以上爬升段 在25到40公里的爬升段 利用上述公式(4)和(5)，所设计的γ跟踪曲线如图6所示。
由于在高空时，大气密度急剧下降。通过分析动力学方程(3)的各项数据可知γ是不可控的，在166-200s时，
约为-0.12度/秒。但在图6中166秒以后要求γ呈上升趋势，这在实际中是无法实现的，必须按照被控对象的特点对上述曲线进行修正。仿真结果表明，按照修正的γ曲线进行控制设计，能够达到最终控制目标。
通过上述分析可知，在高空高速时，飞行器长周期变量是不可控的，需要将其控制问题转化为短周期变量攻角α和俯仰角速率q的控制问题。这里我们按照下述规则设计α和q的跟踪曲线。通过对工作点数据进行分析，限制攻角α在5～10度，按照公式 和实际γ跟踪曲线进行α的跟踪曲线设计。
按照飞机控制的特点，设计俯仰角速率跟踪曲线q＝0。进而由动力学方程以及上述γ的跟踪曲线可以得到攻角α的跟踪曲线，如图7所示。从上面的分析可知，跨大气层飞行器的重量和气动参数等在爬升过程中发生巨大的变化。这里采用基于特征模型的全系数自适应控制律，减少了控制器设计的复杂程度。
选择采样周期为50ms，首先分别建立攻角和俯仰角速度的特征模型如下 α(k+2)＝f11(k)α(k+1)+f12(k)α(k)+g1(k)u(k) q(k+2)＝f21(k)q(k+1)+f22(k)q(k)+g2(k)u(k) 其中fi1(k)∈[1.4331，1.9974]，i＝1，2；fi2(k)∈[-0.9999，-0.5134]，i＝1，2；适当选取gi(k)∈
，i＝1，2，使用最小二乘方法辨识上述6个参数，辨识结果分别记为和并把其限制在上述范围内。利用上述参数设计全系数自适应控制律 uj＝ugj+uij+udj，j＝1，2 (7) 其中， uij(k)＝uij(k-1)-kijej(k)，j＝1，2 udj(k)＝-kdj(ej(k)-ej(k-1))，j＝1，2 l1＝0.382，l2＝0.618 e1(k)＝α(k)-αr(k)，αr(k)表示攻角跟踪曲线； e2(k)＝q(k)-qr(k)，qr(k)表示俯仰角速度跟踪曲线； kijl，j，l＝1，2，cdj，ldj，j＝1，2，为需调试参数。
这里需设计的控制舵面为升降舵偏差和襟翼。通过对俯仰气动参数进行考察，可知襟翼为-15°时，可以提供最大俯仰力矩，而小于-15°时出现了反效现象。因此这里首先固定襟翼为-15°，以提供俯仰力矩。仅设计自适应控制输入升降舵偏差为u1+u2，u1，u2见(7)式。
在初始参数发动机工作初始时间为ENGWT＝110秒，h＝25千米，M＝2.5马赫，α＝5°，γ＝28°，DE＝-12°时，仿真结果如图8-15所示，图8和图9为攻角α和俯仰角速率q的跟踪结果；图10-15分别为高度仿真结果、马赫数、俯仰角、爬升角、控制输入和过载仿真结果。从仿真结果可见，上述设计不仅实现了高度和速度的协调控制，而且满足俯仰角和过载等的限制要求，并且减少了设计的复杂程度。
下面通过仿真考察上述控制律的对初始参数的鲁棒性。上述仿真结果初始参数分别为发动机工作初始时间ENGWT，高度初值h，马赫数初值M，攻角初值α，爬升角初值γ，升降舵偏差初值DE，下面在上述初值变化的下列几种情况下进行了仿真 1、ENGWT＝110秒，M＝2.5马赫，α＝5°，γ＝28°，DE＝-12°，高度初值h在23～25.3千米范围内； 2、ENGWT＝110秒，h＝25千米，α＝5°，γ＝28°，DE＝-12°，马赫数初值M在2～2.9马赫范围内； 3、ENGWT＝110秒，h＝25千米，M＝2.5，α＝5°，DE＝-12°，爬升角初值γ在20°～29°范围内； 4、ENGWT＝110秒，h＝25千米，M＝2.5，α＝5°，γ＝28°，升降舵偏差DE初值在-9°～-14°范围内； 5、ENGWT＝110秒，h＝25千米，M＝2.5，γ＝28°，DE＝-10°，攻角初值α在4.8°～5.2°范围内； 6、h＝25千米，M＝2.5，α＝5°，γ＝28°，DE＝-10°，发动机工作初始时间ENGWT在109.1秒～110.4秒范围内。
仿真结果表明，在上述初值条件下，对于攻角和俯仰角速度都有很好的跟踪效果。由于在时间为90秒时，推力从240千牛变为0，发生了很大变化。延长仿真时间，仿真结果如图16和17所示，分别为攻角α和俯仰角速率q的跟踪结果。可见控制律在90秒后工作正常。
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。
权利要求
1、一种空天飞机高空爬升纵向控制方法，其特征在于步骤如下
(1)按照空天飞机动力学方程和高空气动特点，确定空天飞机高空爬升速度上升曲线；
(2)根据步骤(1)得到的速度上升曲线，得到爬升角γ的跟踪曲线；
(3)根据步骤(2)得到的爬升角γ的跟踪曲线确定攻角α的跟踪曲线；
(4)确定俯仰角速率q的跟踪曲线；
(5)根据步骤(3)和(4)的结果，确定攻角α和俯仰角速率q的跟踪控制律。
2、根据权利要求1所述的一种空天飞机高空爬升纵向控制方法，其特征在于所述步骤(1)中得到的空天飞机高空爬升速度上升曲线方程为
其中，V为空天飞机的爬升速度，P为空天分机受到的推力，m为空天飞机的质量，μ为地球引力常数，r为空天飞机距地心的距离，γ为空天飞机的爬升角。
3、根据权利要求2所述的一种空天飞机高空爬升纵向控制方法，其特征在于所述步骤(2)中得到的爬升角γ的跟踪曲线方程为
其中，γ为空天飞机的爬升角，P为空天分机受到的推力，m为空天飞机的质量，μ为地球引力常数，r为空天飞机距地心的距离，V为空天飞机的爬升速度，a为声速。
4、根据权利要求2所述的一种空天飞机高空爬升纵向控制方法，其特征在于所述步骤(3)中得到的攻角α的跟踪曲线方程为
其中，α为攻角，γ为空天飞机的爬升角。
5、根据权利要求1或2所述的一种空天飞机高空爬升纵向控制方法，其特征在于所述步骤(4)中俯仰角速率q的跟踪曲线方程为q＝0。
6、根据权利要求1或2所述的一种空天飞机高空爬升纵向控制方法，其特征在于所述步骤(5)中攻角α和俯仰角速率q的跟踪控制律采用基于特征模型的全系数自适应控制方法得到，步骤为
首先分别建立攻角α和俯仰角速率q的特征模型如下
α(k+2)＝f11(k)α(k+1)+f12(k)α(k)+g1(k)u(k)
q(k+2)＝f21(k)q(k+1)+f22(k)q(k)+g2(k)u(k)
式中，fi1(k)∈[1.4331，1.9974]，i＝1，2；fi2(k)∈[-0.9999，-0.5134]，i＝1，2；
然后确定gi(k)，i＝1，2的范围，使用最小二乘方法确定式中的6个参数，和
最后，利用得到的和设计全系数自适应控制律uj＝ugj+uij+udj，j＝1，2，
式中，
uij(k)＝uij(k-1)-kijej(k)，j＝1，2
udj(k)＝-kdj(ej(k)-ej(k-1))，j＝1，2
l1＝0.382，l2＝0.618
e1(k)＝α(k)-αr(k)，α(k)为状态变量，αr(k)表示攻角跟踪曲线；
e2(k)＝q(k)-qr(k)，q(k)为状态变量，qr(k)表示俯仰角速度跟踪曲线；
kijl，j，l＝1，2，cdj，ldj，j＝1，2为需调试参数。
全文摘要
本发明公开了一种空天飞机高空爬升纵向控制方法，步骤为(1)按照被控对象动力学方程和高空气动特点计算速度上升曲线；(2)沿着速度上升曲线计算爬升角的跟踪曲线；(3)按照爬升角的跟踪曲线设计攻角的跟踪曲线；(4)设计俯仰角速率跟踪曲线；(5)设计攻角和俯仰角速率跟踪控制律。本发明的控制方法为跟踪曲线设计奠定了理论依据，不仅实现了协调控制的目标，而且满足俯仰角和过载限制，克服了现有飞机控制中的调参问题，降低了设计的复杂性。
文档编号G05B13/04GK101393458SQ20081022517
公开日2009年3月25日 申请日期2008年10月30日 优先权日2008年10月30日
发明者果 李, 斌 孟, 吴宏鑫, 王大轶, 孙承启, 李智斌, 倪茂林, 杨俊春, 谈树萍, 欧阳高翔 申请人:北京控制工程研究所