月球最优软着陆两点边值问题的数值解法.doc

月球最优软着陆两点边值问题的数值解法 3王大轶李铁寿马兴瑞中国航天科技集团公司 , 北京 100830北京控制工程研究所 , 北京 100080摘 要 对于月球软着陆 , 燃耗最优是制导过程的基本要求 。文中首先应用极大值原理设计了最优着陆制导控制律 , 此时求解最优轨迹变成一个两点边值问题 ( TPBVP) 。本文利用一种基于初值猜测技术的打靶法求解这个两点边值问题 , 得 到软着陆最优轨迹 。结果表明该方法可有效改善迭代计算 , 具有一定的优越性 。 主题词 月球软着陆 最优轨迹 两点边值问题Optimal L unar Soft Landing 3Numerical Sol ution of TPBVP inWang Dayi Li TieshouBeijing Institute of Control Engineering , 100080Ma XingruiChina Aerospace Science Corporation , 100830Abstract Minimal f uel guidance is the primary demand f or lunar sof t landing . First , theMaximum principle is used to generate an optimal guidance law f or lunar landing , and the TPBVP is to be solved as a result of the numerical solution f or optimal trajectory . In this pa2 per , shooting methods based on an initial variable guess technique are proposed to solve the TPBVP, and the optimal landing trajectory is obtained . A simulation result is given todemonstrate the f easibility of the improved method f or iteration calculation .Subject termsL unar sof t landingOptimal trajectory TPBVP引言1月球表面软着陆是实现月球探测任务的一项关键技术1 。软着陆过程一般从一条环月停泊轨道开始 。一旦选定着陆点 , 则通过V 变轨转入一条椭圆轨道 。当到达近月点时 ,指令着陆制动发动机点火 , 开始动力下降着陆 。燃耗最优是月球软着陆过程的基本要求2 ,3。文中首先应用极大值原理设计了最优着陆制导控制律 , 通过解软着陆方程得到最优着陆轨迹 。最优着陆轨迹的求解是一个关于两点边值问题的求解过程 , 一般需要经过大3 国家自然科学基金资助项目 (19782004)收稿日期1999 年 12 月 44量的迭代计算才能得出。解两点边值问题的方法有很多 , 主要有4 : 插值法 (L1Fox) 、变量法 (R1S1Varga) 、配置法 (L1Collatz ) 、伪 线 性 化 法 ( R1Bellman ) 和 打 靶 法 ( P1B1Bailey , H1B1Keller , T1R1Goodman etc) 等。其中 , 打靶法由于简单、易于编程且能保证局部收敛性 , 因此应用较多。但是 , 应用打靶法解两点边值问题时需要首先猜测未知状态变量的初值 , 当 猜测值与真值相差较大时 , 计算过程往往会陷入局部极值点 , 或者使计算过程发散。本文在一种初值猜测技术的基础上利用打靶法求解这种关于月球软着陆最优控制的两 点边值问题 。具体做法是 : 将对共轭变量初值的猜测问题转变为对有物理意义的量 (如推 力和推力方向角初值) 的猜测问题 。通过解方程组得到一组初值 , 用打靶法经迭代计算可 得到初始共轭变量的真值 , 此时将状态方程进行积分便得最优着陆轨迹 。仿真结果表明 , 这种解最优着陆过程两点边值问题的方法具有一定的优越性 。燃耗最优制导律2着陆器近月运行轨道为远月点高度 100km 、近月点高度 15km 的椭圆轨道 , 建立图 1 所示的着陆坐标系 。 假设 , 着陆轨迹在纵向面内 , 忽略月球自转 。取月心 O为坐标原点 , Oy1 指向近月点的星下点 , Ox1 指向着陆器运动方向 。着陆器质心运动方程为 :r = vv = ( F/ m) sin - / r2 + r2= = - ( ( F/ m) co s + 2 v) / r m = - F/ Isp( 1)图 1纵向面软着陆极坐标系其中 , v 是着陆器距月心矢径 r 方向上的速度 ; 是着陆器方位角的角速度 ; m 为着陆器质量 ; 为月球引力常数 ; F 为常值制动推力器的推力 , F 取 0 或 Fmax ; Isp 为制动推力器的 比冲 , 是一个常值 ; 为推力方向角 。将式 (1) 表示成状态方程 x = f ( x , u) 的形式 , 其中 x =r , v , m T , u =F , T 。最优软着陆问题设计的目的 , 就是要寻找一个制导律 u , 来调整推力的大小和方向 ,使着陆器在月面实现燃耗最优着陆 。下面应用极大值原理设计这个最优制导律 。定义如下 性能指标 :tJ = f m( t) d t = - ( m ( t ) - m ( t ) )(2)0ft0软着陆的初始条件由霍曼变轨椭圆轨道的近月点确定 , 终端条件为着陆器降落到月面时实现软着陆 。定义初始时刻 t0 = 0 , 终端时刻 tf 不定 。r (0) = r0 , v (0) = 0 , (0) = 0 , (0)r ( tf ) = rf , v ( tf ) = vf , ( tf ) = 0= 0 ,m (0)(3)(4)=m0哈密顿函数为H ( x , u) = Tf ( x , u) -(5)F/ Isprvm T 是下列共轭方程的解 :其中 , 共轭变量 = 9H ( x , u)= -( 6)9xr ( tf ) -v ( tf ) -定义终端约束条件rfvf0G ( x ( tf ) ) = 0( 7)( tf )-横截条件T 9S ( xf )9G ( xf )(8)=+f9x9x式中 , 是拉格朗日乘子 , S ( xf ) = 0 。将方程式 (1) 代入式(5) 得 F co s - 2 v F F F2H = rv + v m sin - v r2+ v r+ - m r- m I -rIspspu ( t) = u 3 ( t) , t t0 , tf 在控制输入达到最优时 , 哈密顿函数 H 取绝对极大 , 即 :H ( x 3 ( t) ,( t) , u 3 ( t) ) = max H ( x 3 ( t) ,( t) , u ( t) )u ( t) Ux 3 为最优状态轨迹 , u 为控制域 。由此可得如下最优制导律 :v r = - arctanL ( t ) L ( t )L ( t)Fmax0不定 0 0= 0(9)v sinco sm1L ( t) =F =-mm rIIsp sp可见 , 只要在区间 t0 , tf 内不存在 L ( t) = 0 的子区间 , 则 F ( t) 是开关函数 。制动推力器要么以最大推力工作 , 要么以最小推力工作 。这样 , 将最优制导律式 ( 9) 代入方 程式 (1) , 利用式 (3) 、(4) 和 (8) , 对状态方程和共轭方程进行积分 , 便可得到最优着陆轨迹 。此时 , 求最优轨迹即变成求解一个关于两点边值的问题 。下面利用一种基于初值猜测技术的打靶法来解两点边值问题 。两点边值问题的解3对于实际的两点边值问题 , 其解析解通常不易求得 , 因此 , 寻找一种好的数值解法就显得异常重要5 。下面给出的这种数值解法基于初值猜测技术和打靶法 。311初值猜测技术用打靶法解两点边值问题 , 需要首先对未知的变量初值进行猜测 。在解关于最优控制 的两点边值问题时 , 一般需要猜测的未知变量初值是共轭变量初值 。而共轭变量一般没有 物理意义 。本文采用通过控制变量来对共轭变量初值进行猜测的方法 。不失一般性 , 考虑系统 (1) 为自治非线性系统 x= f ( x , u) , 终端条件 G ( xf ) = 0 ,t初始条件 x ( t0) = x0 。选择控制变量 u ( t) 使性能指标 J = S ( xf ) +Q ( x , u) d t 达到极ft0值 。下面应用极大值原理解这个问题 。定义哈密顿函数 H = Tf ( x , u)+ Q ( x , u) , 其取绝对极值的控制量 u 3 就是最优控制 。相应的共轭方程= - 9H/ 9x , 横截条件为 f = 469S ( xf ) / 9x + ( 9G ( xf ) / 9x) T,是拉格朗日乘子 。在初始时刻 t0 的某些个邻域内 , 将共轭变量分别进行一阶泰勒展开 , 得如下一组方程n= 0 + | t t n ,( n = 1 , 2 , N )(10)0这里 , t n 是 t0 领域内的某一时刻 , N 不小于未知初值的共轭变量个数 。当非线性系统状态方程不显含时间变量时 , 哈密顿函数在最优控制下是常值H ( t )将哈密顿函数和共轭方程代入式 (10)0 -t t0 , tf (11)= c ,中 , 可得Tf = c9Ht(12)n9xnt0如果控制变量 u ( t) 在初始时刻的值可大致估算出 , 则根据初始条件 , f n 就可通过对系统方程在初始时刻的邻域内进行积分而获得 。当在 t0 的邻域内选定时间点 t1 , t2 , tN 后 ,就能够得到一个 N 维的线性方程组 , 如式 (12) 所示 。通过解这个方程组就可得出未知共轭变量初值的一组估值 。312打靶法 考虑非线性方程组 :y i= gi ( y1 , y2 , y n , t) ,i = 1 , 2 ,( 13), n初始条件为yi ( t0)yi ( tf ) =i = 1 , 2 ,m = 1 , 2 ,( 14)(15)=ci , r, n -终端条件为如果 yi ( t) ,cirmm, n ; t0 t tf, 是方程 (13) 的解 , 假设它的修正解为 yi ( t )i = 1 , 2 ,+yi ( t ) ,这里变分yi ( t ) 是对 yi ( t ) 的一阶修正项 。这样做的目的是将i = 1 , 2 , n ,yi ( t) 看成是相对已知初始条件式 (14) 和未知初始条件估值的方程解 , 将yi ( t) 看作是对于 yi ( t ) 的必要修正以得到边值问题式 (13) (15) 的真实解 。 将修正解代入方程式 (13) 中 , 将其右端进行泰勒展开 , 保留一阶项得n 9gy i = 9y yj ,ii = 1 , 2 , n(16)jj = 1式 (16) 是线性时变微分方程 , 称为变分方程 。根据解线性两点边值问题的结论4变分方程式 (16) 的辅助方程, 构造n 9gpi = - 9y pj ,ji = 1 , 2 , n( 17)ij = 1它也是一个线性时变微分方程 。建立关系式 , 将yi ( t) 与 pi ( t ) 的初值和终值相关联nnpi ( tf )yi ( tf )- pi ( t0)yi ( t0)(18), n 的修正量yi ( t0) ,= 0i = 1i = 1利用关系式 (18) 可求出对于未知初始条件 yi ( t0) , i =r + 1 , n 。选择辅助变量的终端条件i = r + 1 ,1 ,0 ,i = imi imp ( m)( tf )(19)=, m = 1 , 2 , n - ri由此 , 对辅助方程 (17) 进行反向积分 , 可得一组初值 , 它们是线性无关的 。根据式 (18) 及对辅助方程(17) 的 n - r 次反向积分 , 得到下面一组关于 n -, n 的 n - r 维线性方程组 :r 个未知修正量y ( k) ( t ) , i=r + 1 ,i 0y i ( tf )1( k)p (1) ( t )p (1)(1)y r +1 ( t0)( k)r +2 ( t0)p n ( t0)r +1 0y ( k) ( tp (2)(2)(2)( k)y r +2 ( t0)f )r +1 ( t0)p r +2 ( t0)p n ( t0)i2(20)=p ( n - r)( n - r)( n - r)( k)p n ( t0) y n ( t0)y ( k) ( t )( t0)( t0)p r +2r +1i n - rf其中 , 上标 k 代表第 k 次迭代过程 。方程式 (20) 的右端可由下式给出 :y ( k)( k)( tf ) = yi ( tf ) - y i ( tf ) ,mmm = 1 , 2 ,(21)(20) 左端系, n - rim右端第一项由式 (15) 给出 , 第二项为第 k 次迭代过程的终端计算值 。若式数矩阵可逆 , 则可得出未知的初值修正量yi ( t0) ,对于第 k + 1 次迭代过程 , 新的初始条件为, n 。i = r + 1 ,y ( k +1)( t0) = yi ( t0)i = 1 , 2 ,= ci , r ;iy ( k +1)( k)( k)= y i ( t0) +y i ( t0) , i = r + 1 ,( t0)( 22), ni313计算过程由打靶法解两点边值问题时 , 采用构造辅助方程的方法解非线性微分方程的步骤如下 :11 确定偏微分项gi /yj , i , j = 1 , 2 , n 的解析表达式 。21 初始化迭代过程计数器 , 设定 k = 0 。猜测未知的初值 , y (0) ( t ) ,31 对于 k = 0 时 ,i = r + 1 , r + 2 , , n 。i 041 由初值 y ( k) ( t )= c , i = 1 , 2 , , r ; y ( k) ( t ) , n , 对式 (13) 进行积i = r + 1 ,i 0i 0i分 , 保存状态变量 y ( k) ( t ) , n , t t t 。i = 1 , 2 ,i f0f51 设定辅助方程积分计数器 , m = 1 。61 确定式 (20) 左端系数矩阵第 m 行的值 。根据保存的状态变量轨迹 , 计算偏微分( m)项的值 , 由终端条件式 (19) , 从 tf 到 t0 , 对辅助方程式 (17) 进行积分 , 保存 p i ( t0) , n 。i = r + 1 ,71 根据式81 如果 m91 构造 n -101 根据式111 设 k =(21) 确定式 (20)右端第 m 行的值 。= m + 1 , 返回步骤 6 。 n - r , 则令 mr 维线性代数方程 (20) , 解得y ( k) ( t ) , n 。= r + 1 ,ii0(22) 得到下次迭代的状态变量初值 。k + 1 , 并返回到步骤 4 。当 Maxy ( k) ( t ) , i121 终止计算条件 :, n 小于某一设定值 , 或迭代计数= r + 1 ,i f器 k 大于某一最大迭代设置值 。仿真结果4对于月球最优着陆制导 , 如式 (1) (9) 所述的两点边值问题采用基于初值猜测技术的打靶法求解。由系统方程 (1) 和共轭方程 (6) 共同组成一个 10 阶规范方程。状态变量初始和终 端约束式 (3) 和 (4) 的取值分别为 r0 = 1753km; v0 = 0 ;0 = 0 ;0 = 9165 10- 4rad/ s ; m0 = 48600kg ; rf = 1738km; vf = 0 ;f = 0 。月球半径 Rm = 1738km ,月球引力常数 = 4902175km3/ s2 。制动发动机采用常推力开关式控制 , Fmax = 1500N ,比冲 Isp = 300 918m/ s 。首先根据初值猜测技术解共轭变量初值 。由于终端时间 tf 自由 , 哈密顿函数式 (11)等于零 。由横截条件式 (8) 及共轭方程有 = 0 。令m ( 0) = 1 , 每 5s 取一个时间点 , N= 3 。着陆时间 tf 是根据将制动推力用于抵消水平速度的假设进行估算的 , 故有vh ( tf ) - vh ( t)Cm ( t )1 - exptf = t +CFmax其中 , vh 表示着陆器的水平速度 。根据式(12) 可解得初始共轭变量的一组估值 。然后按照打靶法的解题步骤 , 经过迭代运算可得出初始共轭变量的迭代真值及相应的最优着陆曲线 。表 1 为月球最优着陆初始共轭变量的解 。图 2 为最优着陆曲线 。图 2 最优着陆曲线表 1 月球最优着陆共轭变量的初值下转第 55 页rvmHtf猜测值- 41088e - 4- 01094- 31603e + 5110051713 ( s)迭代真值- 51539e - 4010058- 31603e + 511001003853912 ( s)结论共进行八次轨道控制 , 在每次轨控前后 , 调用了 CB ERS 卫星 AOCS 飞控方案研制的七个软件包 , 成功地完成了 CB ERS 的飞行任务 。姿轨控系统的轨控飞行模式运行正常 , 轨 道控制策略和实施方案及软件正确 。轨道维持阶段轨控大系统相对误差 (其中包括测轨与发动机推力模型和卫星动力学模型误差在内) 在 2 %左右 , 标定系数为 019485 至 019665 之间 。这表明短脉冲轨控已经达 到了很高的水平 。参考文献1 杨嘉樨 1 航天器轨道动力学与控制 (上) 1 导弹与航天丛书卫星工程系列 1 宇航出版社 1 1995 年 12月2 刘林 1 人造卫星轨道力学 1 高等教育出版社 , 1993 年3 杨维廉 1 ZY - 1 卫星轨道设计 1 中国空间技术研究院飞行器总体设计部 1 1988 年 9 月4 Jiachi Yang. Development of Intelligent Autonomous Control Technology for the Chinese Space Program. IFAC IACA95 Conference , Beijing5 Jiachi Yang , Zhang Guofu , Sun Chengqi . “Three - Axis Stabilized Attitude Control System for Chinese Near EarthOrbit Satellite”. VIII IFAC Symposium on Automatic Control in Space UK J uly 2 - 6 1979上接第 4