CG伪谱法在欠驱动航天器姿态机动最优控制中的应用*

石青鑫 戈新生

(北京信息科技大学机电工程学院, 北京 100192)

引言

欠驱动航天器是指仅依靠两个执行机构控制实现姿态机动的航天器.Crouch[1]研究了欠驱动航天器的可控性,证明了在系统的动量矩等于零时,欠驱动航天器可控.Morin[2]设计了指数收敛的时变状态反馈稳定控制器,实现了欠驱动航天器的姿态稳定控制.戈新生等[3]利用Ritz理论和粒子群算法,研究了带有两动量轮的欠驱动航天器和太阳帆板展开过程的姿态控制问题.王冬霞等[4]利用分层滑模控制方法,为带两个控制执行机构的欠驱动刚体航天器的姿态控制系统设计了一种三轴稳定控制器.庄宇飞等[5]应用伪谱法解决了欠驱动刚性航天器的时间最优轨迹规划问题,表明伪谱法能够较好地满足各种约束条件,而且计算精度高、速度快,具有良好的实时性.近年来,伪谱法越来越成为解决最优控制问题(OCP)和非线性方程的重要方法,其根本原因是NLP问题的Karush-Kuhn-Tucker(KKT)条件与离散的哈密顿边值问题的一阶最优条件等价得到了证明[6].伪谱法主要包括Legendre伪谱法、Gauss伪谱法、Radau伪谱法、Chebyshev伪谱法,其中Legendre伪谱法和Gauss伪谱法使用较为普遍.唐小军等[7]提出以Chebyshev-Gauss(CG)伪谱法解决OCP的直接轨迹优化,严格推导了约束乘子估计与NLP问题的KKT条件的等价性.

航天器的姿态机动可以通过喷气推力器或动量轮实现.使用CG伪谱法研究欠驱动航天器姿态最优控制问题,采用CG配置点,利用 Clenshaw-Curtis积分近似得到性能指标函数中的积分项,并用快速傅里叶变换(FFT)以快速有效计算相应的积分权；为提高计算效率和数值稳定性,应用重心拉格朗日插值代替经典的拉格朗日插值去逼近状态变量和控制变量,并采用三角函数来减小计算配点数较大时的状态微分矩阵带来的误差,将连续OCP离散为具有代数约束的参数优化NLP问题,再通过SQP算法求解得到控制输入的规律.通过数值仿真得出,对两类欠驱动航天器的姿态最优控制均能达到设计的零边界控制要求,得到的姿态最优曲线与反向验证得到的曲线几乎完全重叠.

1 欠驱动航天器模型

设航天器系统由一个航天器主体B及两个动量轮Wi(i=1,2)组成,并且两轮的自转轴分别与航天器主体的前两个主轴重合(如图1所示),动量轮质心为Ci(i=1,2).以航天器系统质心建立惯性系(O-XYZ),以动量轮的质心建立连体坐标系分别为(Oi-xiyizi)(i=1,2),并设x1,y2分别为两动量轮的自转轴.

图1 带两动量轮的欠驱动航天器简化图Fig.1 Simplified representation of underactuated spacecraft with two momentum wheels

使用欧拉角ψ,θ,φ分别代表偏航向角,俯仰角,滚动角,以3-2-1次序旋转,得到方向矩阵[8]

(1)

式中sin(*)=s*,cos(*)=c*,航天器姿态运动学方程为[9,10]：

(2)

(3)

将式(2)展开可得以欧拉角描述姿态的运动学方程：

(4)

设航天器系统对质心O的总动量矩为H,可表示为：

(5)

其中JB=diag(JB1,JB2,JB3)和JWi分别为航天器主体、两动量轮对质心O转动惯量矩阵,JCi为动量轮对其质心的转动惯量矩阵,Ωi为动量轮相对航天器主体转动角速度矢量.

由于动量轮相对航天器主体作定轴转动,式(5)右侧第三项可简化为：

(6)

(7)

其中J为航天器系统对质心O的转动惯量矩阵,表示如下：

(8)

根据动量矩定理有：

(9)

式中M表示施加在系统质心的外力矩,将式(7)对时间求导：

(10)

(11)

当航天器系统以动量轮作为执行机构时,式(11)中M=0,此时令动量轮的角加速度作为系统的控制变量,即：

(12)

则带动量轮的欠驱动航天器系统动力学模型为：

(13)

其中：

(14)

该系统对应的状态变量为:

(15)

该系统对应的状态变量为:

2 欠驱动航天器姿态最优控制问题

为使欠驱动航天器既实现期望的姿态位置,又满足执行工作任务的其它要求(譬如运行时间最短、消耗燃料最少等),最优控制问题由此而生,在此选取Bolza形式的最优控制：

(16)

上述方程依次代表：性能指标函数、状态约束方程、边界约束方程和路径约束方程.其中x∈n为状态变量,用u∈m代表上文提到的控制变量和M,t为实际任意时间,t0为实际初始时间,tf为实际终端时间.函数Φ和g是标量,f∈n为n维向量函数,φ∈q为q维向量函数,C∈c为c维向量函数.

根据可控性秩条件(Chow定理)[11]可知：控制系统存在一个满秩可控性李代数,原则上能求解系统的姿态机动问题.设控制系统存在最优解u*∈L2([0,T]),其中L2([0,T])为可测向量函数u(t),t∈[0,T]构成的Hilbert空间.根据最优控制原理,选取优化指标为：

(17)

对于带动量轮的航天器,给定系统初末位形、角速度和动量轮的角速度,即x0,xf∈8；而对于带推力器的航天器,给定初末位形和角速度,即x0,xf∈6.通过目标函数寻求控制输入u(t)∈2,t∈[0,T]T,从而可确定系统从x0机动到xf的姿态运动轨迹.控制航天器的执行机构在启动和关闭瞬间的控制值应均为零,即u0=uf=0.另外,控制值不可能无限大,应设置最大值,因此还需考虑

3 Chebyshev-Gauss伪谱法

CG伪谱法通过一系列变换将连续最优控制问题离散为具有代数约束的NLP问题.

首先进行时域变换,需要引入新的时间变量τ∈[-1,1],将时间区间[t0,tf]转换为[-1,1],定义时间变量t为：

(18)

其次对状态变量与控制变量作如下处理.取K个CG点以及τ0=-1作为CG伪谱法的节点,其中CG点[12]的定义为：

(19)

CG点在区间(-1,1)上的分布特点为两端稠密,中间稀疏,能有效抑制拉格朗日插值时两端容易出现的Runge现象.用此K+1个节点构造拉格朗日插值多项式,并以此为基函数来逼近状态变量有：

(20)

为提高计算效率和数值稳定性,采用重心拉格朗日插值[13]：

(21)

其中重心权ξi的定义为：

(22)

对于式(22),在数值计算τi-τj时会产生浮点误差[13],为避免误差作如下处理,设：

(23)

考虑到τK+1=1,由式(22)和(23)可得：

ξi=(τi-τK+1)ξi′

(24)

(25)

基于时域变换、状态变量与控制变量的近似变换,式(16)写作：

(26)

然后将运动学与动力学微分方程进行离散,在CG伪谱法中,由全局正交多项式近似状态变量,可通过对式(20)求导来近似得到：

(27)

式中τk为CG点.文中采用重心拉格朗日插值,因此微分矩阵D∈RK×(K+1)表示如下,该矩阵具有较好的数值稳定性[14]：

(28)

式中k=1,…,K,i=0,…,K.联立式(27)、(28)以及(26)中的状态约束方程,得到离散方程：

(29)

再将终端状态约束离散化,由于式(26)中包含了终端状态约束,而式(20)的时间区间为[-1,1),因而状态变量的离散表达式中未包含末端时间节点X(τf),对(16)中的状态约束方程在区间[-1,1]上积分：

(30)

将式(30)左边展开得：

(31)

对式(31)中的积分项采用Clenshaw-Curtis积分[15]求其近似值：

(32)

(33)

(34)

(35)

同样利用Clenshaw-Curtis积分方法将指标函数离散化,近似式(26)中的性能指标函数,可得：

J= Φ(X0,t0,Xf,tf)+

(36)

考虑到系统模型有两个控制变量,即U1K,U2K∈K,再根据式(17),式(36)变为：

(37)

其中μcc为Clenshaw-Curtis积分权向量.设置状态变量和控制变量初末值及控制约束条件：

x0-X0=0,xf-Xf=0;

(38)

式中X0,Xf分别代表近似状态变量的初末值,相应的控制变量初末值分别为U0,Uf.Umax为控制变量的最大值.

至此,基于CG伪谱法的最优控制问题转化为NLP问题,即在满足式(29)、(32)和式(38)的约束条件下寻求最优控制输入以使式(37)取得最小值.

4 数值仿真

取系统质量几何参数为[17]：

JB=diag(86.215,85.07,113.565)kg·m2

JW1=diag(0.5,0.45,0.45)kg·m2

JW2=diag(0.45,0.5,0.45)kg·m2

J1=J2=0.5kg·m2

假设欠驱动航天器系统进行rest-to-rest的姿态机动,设置机动时间为20s,给定航天器姿态角ψ,θ,φ的初始值依次为0,0,0,末端值为0,π/6,0,控制输入初始和终端值均为零.优化算法采用从可行解到最优解的串行优化策略,可行解选取CG点个数为6个,最优解选取CG点个数为60个,根据上文的理论公式,在内存为4.0G,CPU频率为2.1GHz的Windows7操作系统中使用MATLAB R2015b编程仿真.

带动量轮的欠驱动航天器姿态最优控制结果如图2～图5所示.

图2 带动量轮欠驱动航天器姿态角的机动曲线Fig.2 Optimal trajectory of the underactuated spacecraft with momentum wheels

图中CGPM代表CG伪谱法.为验证文中姿态优化曲线的合理性及有效性,将求解得到的图5中的控制变量代入到状态方程式(4)和式(13)中使用ode45进行积分,积分后得出的状态变量结果如图2～图4中的蓝色虚线所示,可以看出得到的积分曲线与用CG伪谱法优化得到的曲线几乎是重叠的,以图2姿态角ψ曲线图为例,观察其局部放大图可以发现两曲线误差在10-2数量级,可以看出该算法求解得到的姿态运动以及控制输入都是合理且有效的.

图3 带动量轮欠驱动航天器姿态角速度的机动曲线Fig.3 Optimal angular velocity of the underactuated spacecraft with momentum wheels

为直观体现CG伪谱法的优势,文中同时给出采用Gauss伪谱法计算的结果.Gauss伪谱法的配置点为Legendre-Gauss(LG)点；采用Gauss积分来离散性能指标函数中的积分项；应用经典拉格朗日插值去逼近状态和控制变量.Gauss伪谱法同样采用串行优化策略,可行解选取LG点个数为6个,最优解选取LG点个数为60个,状态和控制变量的参数设置均与CG伪谱法的相同.CG伪谱法和Gauss伪谱法可行解的初始猜测值均采用MATLAB的随机函数生成,结果如表1所示.表中的参数Itf和Jf分别表示可行解求解阶段的迭代次数以及目标函数值；参数Ito和Jo则分别代表最优解求解阶段的迭代次数以及目标函数值；而参数Ae和tCPU则分别表示最终的末端姿态最大误差以及CPU总的运行时间.表1直观地表明,CG伪谱法的目标函数的可行解和最优解值均比Gauss伪谱法的小,求解时间比Gauss伪谱法明显缩短.

图4 两动量轮转动角速度曲线Fig.4 Optimal angular velocity of two momentum wheels

图5 最优控制输入(两动量轮角加速度)曲线Fig.5 Optimal control input for two momentum wheel accelerations

表1 CG伪谱法与Gauss伪谱法的比较Table 1 Comparisonof CG Pseudospectral Method and Gauss Pseudospectral Method

带推力器的欠驱动航天器姿态最优控制结果如图6～图8所示.

图6 带推力器欠驱动航天器姿态角的机动曲线Fig.6 Optimal trajectory of the underactuated spacecraft with thrusters

图7 带推力器欠驱动航天器姿态角速度的机动曲线Fig.7 Optimal angular velocity of the underactuated spacecraft with thrusters

从图6～图8中的曲线可看到,姿态角和姿态角速度均可精确机动到预定的末端值,且相应的轨迹曲线是光滑的.

5 结论

本文利用CG伪谱法解决了带动量轮、推力器两类欠驱动航天器姿态机动最优控制问题.基于CG伪谱法把系统的连续姿态运动规划转化成NLP问题,再用求解非线性约束优化问题的SQP算法运算.通过对这两类模型的数值仿真,得到了系统从初始姿态机动到终端姿态的运动轨迹,终端角速度也达到了预定值,而且最优控制输入的初末值均等于零,可以说明该方法的可行性.将求得的最优控制解代回系统的姿态运动方程中,利用数值积分求解状态变量的结果,可以发现两种方法求解得到的曲线几乎完全重叠,进而反向验证了CG伪谱法对两类系统的姿态最优控制的合理性和有效性.另外,使用 Gauss伪谱法进行仿真对比,看出CG伪谱法的计算效率明显更高.

图8 推力器的最优控制输入(力矩)曲线Fig.8 Optimal control input for the thrusters

1Crouch P E. Application of geometric control theory to rigid body models.IEEETransactionsonAutomaticControl, 1984,29(4):87～95

2Morin P,Samson C. Time-varying exponential stabilization of a rigid spacecraft with two control torques.IEEETransactionsonAutomaticControl, 1997,42(4):528～534

3戈新生,陈立群,刘延柱. 带有两个动量飞轮刚体航天器的姿态非完整运动规划问题. 控制理论与应用, 2004,21(5):781～784 (Ge X S, Chen L Q, Liu Y Z. Nonholonomic motion planning for the attitude of rigid spacecraft with two momentum wheel actuator.ControlTheory&Applications, 2004,21(5):781～784 (in Chinese))

4王冬霞,贾英宏,金磊. 欠驱动航天器姿态稳定的分层滑模控制器设计. 宇航学报, 2013,34(1):17～24 (Wang D X, Jia Y H, Jin L. Hierarchical sliding-mode control for attitude stabilization of an underactuated spacecraft.JournalofAstronautics, 2013,34(1):17～24 (in Chinese))

5庄宇飞,马广富,黄海滨. 欠驱动刚性航天器时间最优轨迹规划设计. 控制与决策, 2010,25(10):1469～1473 (Zhuang Y F, Ma G F, Huang H B. Time-optimal motion planning of an underactuated rigid spacecraft.ControlandDecision, 2010,25(10):1469～1473 (in Chinese))

6Benson D A. A gauss pseudospectral transcription for optimal control[PhD Thesis]. Cambridge, Massachusetts：Massachusetts Institute of Technology, 2005

7Tang X J. A Chebyshev-gauss pseudospectral method for solving optimal control problems.ActaAutomaticaSinica, 2015,41(10):1778～1787

8Peterson C. Advances in underactuated spacecraft control[PhD Thesis]. Michigan: University of Michigan, 2016

9Hughes P C. Spacecraft attitude dynamics. New Jersey: Wiley, Hoboken, 1994:17～31

10 Sidi M J. Spacecraft dynamics and control. Cambridge: Cambridge University Press, 1997:25～28

11 Isidori A. Nonlinear control systems II. Berlin:Springer-Verlag, 1999

12 Ge X S, Yi Z G. Optimal control of attitude for coupled-rigid-body spacecraft via Chebyshev-Gauss pseudospectral method.AppliedMathematicsandMechanics, 2017,38(9):1257～1272

13 Berrut J P, Trefethen L N. Barycentric lagrange interpolation.SIAMReview, 2004,46(3):501～517

14 Costa B, Don W S. On the computation of high order pseudo-spectral derivatives.AppliedNumericalMathematics, 2000,33(1):151～159

15 Trefethen L N. Is Gauss quadrature better than Clenshaw-Curtis?SiamReview, 2008,50(1):67～87

16 Waldvogel J. Fast construction of the Fejér and Clenshaw-Curtis quadrature rules.BITNumericalMathematics, 2006,46(1):195～202

17 Krishnan H, Mc clamroch N H, Reyhanoglu M. Attitude stabilization of a rigid spacecraft using two momentum wheel actuators.JournalofGuidanceControl&Dynamics, 1995,18(2):256～263

18 Gong Q, Ross I M, Fahroo F. Costate computation by a chebyshev pseudospectral method.JournalofGuidanceControl&Dynamics, 2010,33(2):623～628