摘要
针对复杂约束下航天器大角度姿态机动问题,提出一种基于三维特殊正交群的姿态高效机动路径规划和时域调节方法。针对航天器大角度机动过程中的姿态约束问题,设计了基于梯度的避障策略,得到了虚拟时域中的姿态路径和期望角速度。为了满足执行机构最大输出力矩约束,提出了迭代非线性时域调节方法,通过时间缩放得到了满足控制力矩约束的高效机动角速度/控制力矩。仿真结果表明,该方法在满足航天器姿态机动过程中姿态约束和控制力矩约束的同时,相较现有同类方法,显著缩短了机动时间,为受约论文拓展束的航天器高效机动规划与控制提供了新的思路。
Abstract
For the large-angle attitude maneuver of spacecraft under complex constraints, a efficient maneuver path planning for attitude and time scaling method on 3-dimensional special orthogonal group was proposed. Aiming at handling the attitude constraints during maneuvering, a gradient-based obstacle avoidance methodology was designed and the attitude routes and desired angular velocity trajectory was obtained on the virtual time domain. Considering the maximum output torque of the actuator, an iterative nonlinear time scaling method was proposed to adjust the angular velocity/control torque. Simulation results show that the proposed method not only satisfies the attitude constraints and the control torque constraints during the spacecraft maneuver process, but also significantly shortens the maneuver time comparing with existing methods. Novel results provide a new insight for efficient constrained attitude planning and controller design.
Keywords
航天器在轨运行期间经常需要执行大角度姿态机动来完成观测、通信等科学任务以延长其工作弧段,提高航天器效能。然而,航天器在姿态机动的过程中,某些光学敏感器性能易受杂光影响或损伤;同时为保证工作区间内的正常通信,天线指向也被要求在一定范围内,这类约束常被称为航天器姿态机动过程中的姿态约束。此外,航天器还受到其最大工作力矩、有效载荷姿态稳定性要求的角速度限制等约束,这些约束统称为动力学约束。同时考虑姿态约束和动力学约束情况下的航天器高效机动路径规划和控制问题研究具有重要的理论和工程应用价值,受到了极大的关注[1]。
针对复杂约束下航天器姿态机动规划问题,解决方案大致上可以分为以下四类:几何法、路径规划法、轨迹优化法和人工势场法。
几何法[2-3]期望通过几何计算得到航天器机动路径。该方法几何意义明确,简单直观,但其随着姿态约束的增加,计算量急剧增大,且不利于扩展到考虑动力学约束的情形。路径规划法通过将姿态参数空间离散化,采用基于网格的或者基于采样的算法搜索从初始姿态到目标姿态的可行路径。基于网格的路径规划算法是分辨率完备的,但在设定分辨率时需在计算量和解的质量之间进行权衡[4-5];基于采样的路径规划算法可以处理由姿态约束和动力学约束形成的复杂约束,但这种算法只能保证概率完备,并且不是一种确定性方法[6-8]。轨迹优化法将原问题视为最优控制问题,通过离散化方法转化为参数优化问题,使用数值优化算法进行求解。为了克服计算量较大的问题,近些年轨迹优化中的凸优化方法发展迅速[9],凸优化方法将非凸姿态约束凸化处理为凸约束,采用数值优化方法求解凸优化问题[10-12]。人工势场法则是通过建立姿态约束依赖的势函数,并将其代入控制器设计,从而达到障碍规避的效果[13-14]。本质上讲,前三类方法属于路径规划的范畴,而最后一类方法属于控制器设计的范畴。
上述文献使用不同方法部分解决了姿态规划问题,但其所用姿态参数几何直观性差,且不同时具有全局性和唯一性。三维特殊正交群(3-dimensional special orthogonal group,SO(3))上的姿态旋转矩阵可以全局地、唯一地表征姿态,相比其他姿态参数具有明显的优势[15-16]。目前直接在SO(3)上进行姿态规划的研究很少。Biggs等[17]基于SO(3)提出一种两步走的半分析方法。该方法预定义姿态轨迹的解析形式,并用试错法针对禁止区域迭代更新参数。Celani等[18]基于SO(3)使用基于梯度的优化方法和线性时域调节(linear time scaling,LTS),实现了指向约束和控制力矩约束下的“rest-to-rest”姿态运动规划,但是其得到的姿态机动时间较长。在此基础上,Celani等[19]进一步使用目标函数为显式最小化机动时间的无梯度优化算法求解姿态机动问题。但该方法计算量较大,且未充分利用执行机构能力,实际机动时间可进一步优化。
鉴于以上观察,为实现复杂约束下的航天器高效机动,同时避免现存文献中计算量大、机动时间长等问题,本文将复杂约束下的高效姿态机动规划问题转化为姿态约束下的路径规划和控制力矩约束下的时域调节问题。本文提出了梯度优化路径规划和迭代非线性时域调节(iterative nonlinear time scaling,INTS)相结合的高效姿态机动运动规划方法,在满足指向约束、控制力矩约束和边界条件约束的情况下,充分利用执行机构能力,相比现有同类方法大幅缩短了机动时间。
1 复杂约束下姿态机动问题描述
1.1 SO(3)上姿态运动学与动力学建模
定义惯性参考系为FI,原点位于航天器质心的体坐标系为FB。旋转矩阵R∈SO(3)定义了一种变换[20],它将向量z在本体系中的投影zB转换为其在惯性系中的投影zI,即zI=RzB。使用R表示航天器姿态,其姿态运动学方程[20]可表示为:
(1)
式中,,ω1、ω2、ω3是航天器角速度在本体系下的分量,矩阵A1、A2、A3是对应于李群SO(3)的李代数so(3)的基:
角速度和控制力矩的关系由姿态动力学方程给出,见式(2)。
(2)
式中,J为航天器在本体系内的转动惯量。
1.2 约束描述
用单位矢量r表示航天器敏感仪器在本体系下的指向。假设对于敏感仪器r有p个受限的指向,用惯性系下的单位矢量wi(i=1,2,···,p)表示。航天器的指向约束可描述为:
(3)
即航天器敏感元件的光轴方向与第i个受限指向之间的夹角大于规定值θi。
通常情况下,执行机构提供的控制力矩大小有限,对应的控制力矩约束为:
(4)
为了满足某些任务的需要,姿态机动中初始角速度和目标角速度必须满足等于零的边界条件约束,产生rest-to-rest机动。
综上所述,研究的问题归纳为在已知初始条件R(0)=R0、ω(0)=03×1和终端条件R(tf)=Rf、ω(tf)=03×1约束的情况下,求出有限时间段t=[0,tf]内的角速度ω(t)和力矩输入量T(t),使姿态机动同时满足指向约束,执行机构控制力矩约束,并且尽可能缩短姿态机动时间。
2 基于梯度优化的路径规划
路径规划阶段仅考虑航天器姿态运动学,将角速度视为控制量,同时只考虑初始条件、终端条件和指向约束。
在虚拟时域τ∈[0,1](对应实际时域时间τ∈[0,tf]),角速度表示为m个基函数vk(τ)的叠加。
(5)
式中,αjk为各基函数的权重系数。
针对式(5),有如下命题成立。
命题1 基函数必须满足vk(0)=vk(1)=0,使得角速度满足边界条件约束ω(0)=ω(1)=03×1。
证明:由式(5)可证明该命题。
命题2 在满足vk(0)=vk(1)=0的基础上,如果基函数还满足(dvk/dτ)(0)=(dvk/dτ)(1)=0,则控制力矩满足T(0)=T(1)=03×1。
证明:由式(5)和式(2)可证明该命题。
考虑到rest-to-rest机动和实际应用的要求,基函数被选为扁长椭球波函数[21]。这些函数较为光滑,并且可以通过参数设置使其在边界处自然趋向于零,符合rest-to-rest机动要求。同时,也可以通过参数设置使基函数的导数在端点处自然趋向于零,便于实现平滑控制。
于是,路径规划问题被建模为优化变量为权重系数αjk的优化问题。其中,需要最小化的目标函数为:
最小化该函数有助于规划得到合理的角速度曲线和姿态机动路径,从而减少航天器的实际机动时间[22]。
已知姿态机动的初始条件为:
(6)
为了保证航天器机动到期望的最终姿态,终端条件描述为如式(7)所示等式约束。
(7)
综上,受约束的航天器姿态路径规划问题建模为如下优化问题:
(8)
为了数值求解上述优化问题,需要将该优化问题进行离散化处理。
为满足机动过程中不等式约束式(3)表示的姿态约束,将虚拟时域τ∈[0,1]等分成n段,定义Δτ=1/n,τl=(l-1)Δτ(l=1,2,···,n+1),使离散时间点τl处的姿态R(τl)满足如下约束:
(9)
离散时间点τl处姿态R(τl)可由姿态运动学方程的Lie-Trotter乘积公式[23]近似得到。
(10)
式中,P(τl)表示无穷小旋转:
使用零阶保持近似ω(τ)时,该近似解表示姿态运动学的连续精确解。式(10)也被用来计算等式约束式(7)中的R(1)。
为了使用基于梯度的优化算法,需要得到约束相对于权重系数αjk的梯度解析表达式。
定义:
对各个离散时间点τl处的姿态求导,得:
(11)
进一步推导得到:
(12)
在进一步推导P(τq)/ ωj(τq)的表达式之前,先证明引理1。
引理1 设A和B为l×l斜埃尔米特矩阵(即AH=-A,BH=-B),λi,分别表示A的特征值和对应的单位特征向量。设,则对于任何m,n∈1,2,···,l,下式成立。
证明:
当λm=λn时,易得
当λm≠λn时,
综上所述,引理1得证。
P(τq)相对于ωj(τq)偏导数由命题3描述。
命题3 Ω(τq)是三维反对称矩阵,它的特征值为零和一对共轭纯虚数。设(i=1,2,3)为Ω(τq)的特征值和对应的单位特征向量。定义复矩阵 [u1 u2 u3],易知UHU=UUH=I3×3。于是
(13)
式中,E是一个3×3复矩阵,Emn是E矩阵中的元素。
证明:
式(13)等价于
进一步等价于
注意到
于是
进一步使用引理1的结论,该命题得证。
至此,利用式(10)~(12)和命题3可以得到约束式(7)和式(9)关于αjk的梯度的解析表达式。
3 基于迭代非线性时域调节的运动规划
在通过路径规划得到角速度后,可以直接通过逆动力学求解力矩输入量。但一般情况下虚拟时域控制力矩较大且分布不合理(过于集中)。通过迭代非线性时域调节,可对虚拟时域内的控制力矩曲线进行调整,减小其峰值(最大绝对值),而后将其映射到实际时域,以兼顾姿态约束和动力学约束,缩短航天器姿态机动时间。
下面推导时域调节前后角速度之间的关系和时域调节前后控制力矩之间的关系,然后给出运动规划的具体步骤。
3.1 角速度关系
设时域间的映射关系为τ=F(t),相应的姿态关系为R(t)=R*(τ),即时域调节不改变姿态。
已知姿态运动学方程:
(14)
(15)
其中,(·)*指时域τ中的物理量,否则为时域t中的物理量,下同。
下面推导Ω*(τ)和Ω(t)的关系,进一步得到ω*(τ)和ω(t)关系。首先根据姿态运动学关系得到:
(16)
结合式(14)~(16),整理得到:
(17)
利用式(17),结合以下两式:
不难得到:
(18)
从式(18)中看出时域调节前后角速度间的关系为简单的比例关系。
3.2 控制力矩关系
根据姿态动力学方程,分别建立两个时域内的姿态动力学方程:
(19)
(20)
将式(18)代入式(20),得到:
(21)
式(21)与式(19)对比得:
(22)
从式(22)中看出时域调节前后两个时域内控制力矩间的并不是简单的比例关系。
对于式(22),有以下命题成立。
命题4 当τ处于虚拟时域,即τ∈[0,1]时 式中,绝对值符号表示按向量的每个元素分别取绝对值,下标i指向量中第i个分量,下同。
证明:设,则τmax∈[0,1]。
已知初始条件ω*(0)=03×1,根据拉格朗日中值定理,存在ξ∈(0,τmax),使得:
于是
进一步
于是
命题得证。
命题5
式中,λmax(·)(λmin(·))表示矩阵的最大(最小)特征值。
证明:参见文献[24],命题得证。
一般情况下,λmax(·)与λmin(·)相差不会很大,由命题5可知的值较小,进而可以推断较小。一般比大一个数量级以上,于是有:
实际上,一般比大一个数量级以上。于是根据式(22),假设时域调节前后均处于虚拟时域,即τ∈[0,1]且t∈[0,1],则对的大小起主导作用的是项。于是可以通过设计使时域调节减小三轴控制力矩峰值,即使时域调节前后:
反复进行这个过程直至收敛。
3.3 迭代非线性时域调节
通常,直接在虚拟时域上由逆动力学得到的控制力矩峰值较大,并不满足执行机构控制力矩约束。因此需要通过时域调节调整控制力矩的值,使得实际时域中控制力矩满足约束。
式中,τ为虚拟时域,t为实际时域。
进一步得到:
将该映射关系分别代入式(18)和式(22),得到如下关系:
(23)
为了使控制力矩满足执行机构控制力矩约束式(4),tf必须满足:
一般选定
(24)
然后由式(23)得到实际时域t∈[0,tf]中的物理量。
线性时间调节采用将时间等比例延长的方式进行时域调节来减小控制力矩峰值。该方法使得所有时间段内的控制力矩以相同百分比减小,并未考虑控制力矩的实际分布特性。较大幅值的控制力矩一般仅集中分布于某一时间段,等比例缩放并没有对这种集中现象进行优化。所以线性时间调节在使实际时域控制力矩满足动力学约束的同时,实际机动时间变得很长,留下进一步优化的空间。
为了解决线性时间调节方法的不足,提出迭代非线性时域调节方法,先在虚拟时域内对控制力矩和角速度曲线进行修整,减小控制力矩曲线峰值,而后采用线性时间调节将它们映射到实际时域,如图1所示。整个运动规划分为两个阶段。
图1运动规划示意图
Fig.1Sketch map of the motion planning
阶段一:采用迭代非线性时域调节,在虚拟时域τi∈[0,1](i=1,···,ns)内迭代进行时域调节调整控制力矩形状,减少三轴控制力矩的峰值。
阶段二:采用线性时间调节,得到实际时域t∈[0,tf]的角速度和控制力矩。
阶段一中迭代非线性时域调节的时域间映射关系设计为:
(25)
式(25)为第i次迭代时,时域调节后τi+1与时域调节前τi之间的关系,其中a为常数,可达到调节控制力矩形状的目的。
式(25)中映射具有如下性质:
1)τi=0时,τi+1=0。
2)τi=1时,τi+1=0。
3) 。
4) 。
其中,性质1、性质2保证了虚拟时域映射前后在端点处满足时间长度不变的边界条件;性质3保证了该映射对虚拟时域内控制力矩曲线的塑形能力,当,,由3.2节可知,该时域调节将减小三轴控制力矩的峰值;性质4保证在一般迭代时域调节中,控制力矩关系式(22)中为主导项,的值通常比小一个数量级以上。
在第i次迭代时,设
(26)
(27)
为了充分发挥式(25)减小三轴控制力矩峰值的能力,应当调整a的值使
(28)
这样将使得的值尽可能减少,进而使第i次迭代时域调节后:。
令式(27)中k=1,然后将式(27)代入式(28),整理后得到:
(29)
所以a的值通过式(29)求解,式中τimax可在第i次迭代的时域调节时确定。
为进一步说明式(25)的功能,下面举例说明。将a=0时的图像绘制如图2所示。
图2a=0时取值
Fig.2 when a=0
a=0时,由式(25)得到τi+1=0.5对应τi=0.5,而在τi+1=0.5时取得最小值0.840 8。根据式(18),经过该时域调节后,ω(0.5)幅值减小。此外,经过该时域调节后,由于式(22)中主导项减小,T(0.5)幅值减小。
经过ns次非线性时域调节后,最大角速度和最大控制力矩幅值将不再发生明显变化,此时终止迭代。然后通过线性时间调节方法,得到姿态机动所需的真实角速度和控制力矩。
综上所述,运动规划阶段算法如算法1所示。在该算法中,运动规划阶段一中算法使用非线性时域调节方法使得角速度和控制力矩曲线更为平滑,期望起到角速度和控制力矩在时间域上平均的结果。从操作层面看,总可以选取时域调节后控制力矩峰值不大于时域调节前的结果,如果调节后控制力矩峰值大于调节前,则可以舍弃相应的时域调节结果。从实际效果看,运动规划阶段一中算法选取了连续nf次调节后,控制力矩峰值未减小作为算法终止条件,也可认为是算法的实际收敛条件。
算法1 运动规划算法
Alg.1 Motion planning algorithm
4 数值仿真
4.1 算法有效性验证
如文献[19]所述,侧摆机动是近地轨道航天器执行对地观测的一种典型机动模式。初始姿态设置为R0=I3×3,目标姿态设置为Rf=exp(3/4πA1),侧摆角度为135°。航天器配备了一个星敏感器,假定其在机动期间必须以规定的偏移角度避开太阳和月亮的方向。本体系中星敏感器的指向为:r=[0 1 0]T。两个禁止锥的参数是:
1)太阳锥:惯性系中方向w1=[-0.18 0.56 0.81]T,最小偏角θ1=47°。
2)月球锥:惯性系中方向w2=[0.95 0 0.31]T,最小偏角θ2=36°。
航天器的转动惯量是J=diag[3 000 4 500 6 000]kg·m2,沿每个本体轴的最大控制力矩设置为Tmax=0.25 N·m。
在计算路径规划中的优化问题时,vk(τl)选择扁长球面波函数的离散形式,称为Slepian序列[25]。这些序列通过长度N和半带宽参数W定义。本次仿真中vk(τl)使用N=500,W=0.015的前M=4个Slepian序列,如图3所示。
图3基函数vk(τl)(k=1,2,3,4;l=1,2,···,500)
Fig.3Basis functions vk (τl) (k=1,2,3,4;l=1,2,···,500)
由图3可知,选取的基函数有如下性质:
1) 。
2)基函数较为光滑。
根据命题1,性质1保证得到的姿态机动为rest-to-rest形式;根据命题2,性质1使得边界处控制力矩接近于零。性质2保证了由式(5)得到的角速度也较为光滑,进而由逆动力学得到的控制力矩也较为光滑,便于实际使用。
通过基于局部梯度的内点法求解优化问题式(8)得到决策变量αjk的值如表1所示,进而使用式(5)得到虚拟时域的角速度ω(τ),如图4所示。根据逆动力学式(20)得到虚拟时域的控制力矩T(τ),如图5所示。
表1αjk的解
Tab.1 Solution of αjk
由图4和图5可知,虚拟时域的角速度值和控制力矩值均较大,不满足航天器所受到的动力学约束,所需控制力矩超过航天器执行机构所能提供的最大力矩。
为实现在有限最大控制力矩情形下的姿态快速机动,采用迭代非线性时域调节方法进行时域调节,最后得到实际姿态角速度、控制力矩曲线,分别如图6和图7所示。迭代过程中姿态机动时间变化如图8所示。
图4虚拟时域角速度曲线
Fig.4Angular velocity curve on the virtual time domain
图5虚拟时域控制力矩曲线
Fig.5Control torque curve on the virtual time domain
图6实际时域角速度曲线
Fig.6Angular velocity curve on the real time domain
图7实际时域控制力矩曲线
Fig.7Control torque curve on the real time domain
图8姿态机动时间随迭代次数变化
Fig.8Attitude maneuver time after each iteration optimization
由图6和图7可知,经过运动规划算法,航天器姿态机动过程中的角速度和控制力矩分布更加均匀合理,且满足航天器最大控制力矩限制。整个机动过程中控制力矩连续可导。
边界处的值接近于零,相比其他方法(如时间最优规划)得到的不连续控制力矩(包括在边界处控制力矩不为零的情况),更适合实际应用[26]。
因为控制力矩的突变会激发航天器的柔性模态和溢出效应(姿态机动终止后的振动),从而降低规划的精度并可能损坏航天器上的精密器件[26-28]。从图8中可以看到,姿态机动所需要的机动时间随着迭代次数的增加逐渐减少,且收敛到一个稳态值。
该场景也被文献[18]和文献[19]中的方法使用。为了验证本文方法的有效性,将本文方法规划得到的实际机动时间与现有同类方法[18-19]得到的实际机动时间进行对比,如表2所示。从表2中看出,本文方法规划得到的机动时间较现有方法明显缩短。
表2规划结果对比
Tab.2 Comparison of planning results
惯性系中敏感仪器指向在单位球上的路径如图9所示。由图9可知,本文规划得到的航天器姿态机动在避开了禁止锥的同时选择了明显较短的路径。
图9敏感仪器指向在单位球上的路径
Fig.9Path of sensitive direction on the unit sphere
4.2 算法随机性测试
本节中部分场景的设定来源于文献[13]。假设航天器携带一个光学敏感仪器,其视轴方向与航天器本体系的z向固连,即r=[0 0 1]T。航天器沿每个本体轴的最大控制力矩设置为Tmax=0.1 N·m。航天器的转动惯量为J=diag[350 180 290] kg·m2。
在本场景中,航天器要进行姿态机动使其上敏感仪器再定向,在此过程中要避开四个天体。设置互不重叠的四个姿态禁止区域,详细信息如表3所示,其中指向相应天体的单位向量在惯性系下表示。
随机选取初始姿态和目标姿态并与上述环境参数组成一个完整的随机场景。仿真中,本文方法使用的基函数vk(τl)为N=500,W=0.015的前M=4个Slepian序列,与4.1节仿真相同;文献[18]方法的参数设置与文献[18]中应对较复杂的案例1时所采用的参数相同;文献[19]方法的参数设置与文献[19]中应对较复杂的案例2时所选取的参数相同。选取10个随机场景,使用本文方法规划得到的惯性系下敏感仪器路径如图10所示。
表3仿真参数
Tab.3 Simulation parameters
图10中红色标记表示初始点,绿色标记表示目标点。从图10中看出,所有路径均满足指向约束。
图10随机场景中敏感指向路径
Fig.10Path of sensitive direction in random scene
表4列出了10个随机场景中不同方法规划得到的实际机动时间tf。
从表4中看出,在各种随机场景中,本文方法表现均最好,大部分情况下文献[18]方法表现最差。以文献[18]方法作为基线方法,得到各规划方法的百分比如表5所示。10个随机场景在文献[18]、文献[19]和本文情况下的平均值分别为100%、92.6%和66.6%,标准差分别为0%、23.9%和1.8%。
以文献[18]方法作为基线方法,文献[19]方法规划得到的实际机动时间的百分比平均值和标准差较本文方法高,说明文献[19]方法结果较差、分布较为分散且稳定性差。注意到场景3中文献[19]方法规划所得结果比基线方法差,而本文方法并没有出现过这种情况。
表4实际机动时间对比
Tab.4 Comparison of actual maneuver time
表5规划表现对比
Tab.5 Comparison of planning performance
5 结论
本文研究了复杂约束下SO(3)流形上航天器高效机动路径规划和时域调节问题:将典型的rest-to-rest机动分为路径规划和运动规划两个步骤,实现航天器姿态约束和动力学约束解耦。在路径规划过程中,采用梯度优化法得到了航天器姿态受限情形下的机动路径;在运动规划过程中,利用所提出的迭代非线性时域调节将虚拟时域内角速度和控制力矩整形,使航天器姿态机动角速度和控制力矩分布更均匀合理,然后利用线性时间调节法,使航天器满足最大输出力矩约束。仿真结果表明,本文方法实现了航天器在姿态约束和最大力矩受限等复杂约束下的高效机动,较现有同类研究,大幅缩短了姿态机动所需时间。与现有运动约束下时间最优规划算法相比,其优点为:采用分层求解的思路,降低了问题的复杂度;规划得到的控制力矩在初始与终止时刻为零且较为平滑,适合于工程中使用;运动规划阶段算法输入角速度和控制力矩,输出相同姿态路径下更短时间的角速度和控制力矩,所以该算法可以在其他规划算法结果的基础上做进一步优化,具有广阔的应用前景。该方法为受约束的航天器高效机动规划与控制提供了新的思路。未来的工作将探索该方法在时间最优规划方面的应用。