技术领域
[0001] 本
发明属于电磁编队卫星姿态和轨道相对控制技术领域，具体涉及一种电磁编队卫星姿态和轨道相对控制方法。
背景技术
[0002] 卫星电磁编队飞行(EMFF)具有无需推进剂，而且没有光学
传感器污染等多种优点，在近十年得到了广泛关注。电磁编队飞行这一概念能够应用于多种未来的任务，例如分布式
航天器，天文的干涉量度分析法，以及相机装配等。
[0003] Miller教授的团队首先提出了电磁编队飞行这一概念，随后，MIT的空间系统实验室开始对其进行研究，主要研究了磁偶极子模型并进行了动态分析、
控制器的设计以及地面模拟试验，这些研究为电磁编队飞行的研究打下了坚实的
基础。随后，Sedwick教授的团队提出了谐振感应近场发电系统，这是电磁编队飞行在微重
力环境下的首次验证。近年来，在中国也进行了很多关于电磁编队飞行的研究，这些研究主要是关于电磁编队飞行的相对平衡状态、动力学性能和对其控制方面。研究提出了多种控制器，例如Zeng et al.于2012年提出的有限时间控制器，Cai et al.于2013年提出的误差反馈控制器，Zhang et al.于2014年提出的自适应控制器和Huang et al.于2015年提出的基于线性矩阵不等式的控制器等。
[0004] 电磁编队卫星的姿态能够通过给分布式卫星配备三个
正交的反作用
飞轮来实现，分布式卫星的相对
位置可以通过电磁相互作用来控制。以往的研究主要是假定能够对其姿态进行很好的控制，来研究运动动力学模型和对电磁编队飞行的控制。针对卫星姿态控制系统，主要都是用反作用飞轮来实现，反作用飞轮通过低速摩擦补偿和鲁棒控制器的设计来进行控制。通过以上设计方法得到的卫星姿态和轨道相对控制系统，仍然具有较大的局限性，难以更为高效的有针对性的对卫星姿态和轨道进行相对控制。
发明内容
[0005] 针对
现有技术存在的
缺陷，本发明提供一种电磁编队卫星姿态和轨道相对控制方法，可有效提高电磁编队卫星姿态和轨道相对控制
精度。
[0006] 本发明采用的技术方案如下：
[0007] 本发明提供一种电磁编队卫星姿态和轨道相对控制方法，包括以下步骤：
[0008] 步骤1：建立姿态动力学方程；姿态动力学方程包括：
[0009] (一)用四元数描述的线性形式的姿态动力学方程：
[0010]
[0011] 其中：
[0012] —描述航天器姿态四元数的一阶导数；
[0013] ωxb、ωyb、ωzb分别代表体
坐标系三个轴的转动
角速度；
[0014] q1＝cos(φ/2)；
[0015] q2＝exsin(φ/2)；
[0016] q3＝eysin(φ/2)；
[0017] q4＝ezsin(φ/2)；
[0018] φ—欧拉轴的旋转角；
[0019] [ex ey ez]—欧拉轴的单位矢量；
[0020] (二)体坐标系下的姿态动力学方程：
[0021]
[0022] 其中：
[0023] Ixb、Iyb、Izb分别为体坐标系下三个轴的
转动惯量；
[0024] Lxb、Lyb、Lzb分别为体坐标系下三个轴的力矩；
[0025] (三)刚体卫星的动力学方程：
[0026]
[0027] 其中：
[0028] IT—卫星及飞轮整体的转动惯量；
[0031] Iw—反作用飞轮的转动惯量；
[0032] ωB—卫星的角速度；
[0033] ωw—体坐标系下飞轮的角速度；
[0034] LB—体坐标系下的力矩；
[0035] 步骤2：根据步骤1建立的四元数描述的线性形式的姿态动力学方程、体坐标系下的姿态动力学方程和刚体卫星的动力学方程，设计得到卫星姿控系统反馈控制律；
[0036] 本步骤具体包括：
[0037] 步骤2.1，假定当前四元数为q，期望四元数为qc，则误差四元数qe的定义为：
[0038]
[0039] 假定当前角速度为ω，期望角速度为ωc；则误差角速度ωe的定义为：
[0040] ωe＝ω-ωc (2)
[0041] 步骤2.2，基于误差四元数和误差角速度的姿态动力学模型可被描述为[0042]
[0043] 其中：
[0044] I—整星的惯量矩阵；
[0045] —误差角加速度；
[0046] —期望角加速度；
[0047] H—整星的
角动量，H＝IBωB+Iwωw；
[0048] IB—卫星的转动惯量；
[0049] u—反作用飞轮的控制力矩，u＝Iwωw；
[0050] Md—外部干扰力矩；
[0052] —误差四元数的一阶导数；
[0053] 步骤2.3，将地球磁场近似为一个以地球为中心的偶极子；非线性控制系统的非线性部分能够利用反
馈线性化方法被抵消，从而提出反馈控制律：
[0054]
[0055] 其中：
[0056] P，D—正定增益矩阵，通常被设计为对角矩阵；
[0057] qe123—误差四元数中的三个分量qe1、qe2、qe3组成的列向量，qe123＝[qe1 qe2 qe3]T；
[0058] 步骤3：根据步骤2提出的反馈控制律，对电磁编队卫星姿态和轨道进行相对控制。
[0059] 优选的，步骤3具体为：
[0060] 在根据步骤2提出的反馈控制律对电磁编队卫星姿态和轨道进行相对控制时，在控制过程中，实时判断是否发生磁饱和现象，如果发生，采用飞轮角动量磁卸载
算法卸载磁力矩。
[0061] 优选的，在控制过程中，实时判断是否发生磁饱和现象，如果发生，采用飞轮角动量磁卸载算法卸载磁力矩，具体为：
[0062] 步骤3.1，假设反作用飞轮的冗余转动惯量为ΔH；在控制过程中，实时测量飞轮所处磁场B，判断是否满足以下的磁卸载条件：
[0063] 45°＜(B,ΔH)＜135° (30)
[0064] 如果满足，进一步判断磁场B是否垂直于ΔH，如果不垂直，则按式(37)计算得到磁卸载力矩值Tunloading，并将该磁卸载力矩值Tunloading施加于飞轮中，使其反向转动，进行磁卸载：
[0065]
[0066] |B|—飞轮所处磁场B的大小；
[0067] K—增益系数；
[0068] 如果B垂直于ΔH，则按下式(38)计算得到磁卸载力矩值Tunloading，并将该磁卸载力矩值Tunloading施加于飞轮中，使其反向转动，进行磁卸载：
[0069] Tunloading＝-KΔH(K＞0) (38)
[0070] 优选的，还包括：
[0071] 步骤4，设计大角度姿态机动方法，包括：
[0072] 步骤4.1，假设机动前，即初始时刻t0时刻的理论姿态为q(t0)，理论角速度为ω(t0)＝[0 0 0]T；机动后，即tf时刻目标姿态为q(tf)，目标角速度为ω(tf)＝[0 0 0]T；
[0073] 步骤4.2，根据下式计算得到控制力矩Lc：
[0074]
[0075] 其中：
[0076] IT—卫星及飞轮整体的转动惯量；
[0077] q123—q1、q2、q3组成的列向量，q123＝[q1 q2 q3]T；
[0078] Kd—大角度姿态机动控制律的微分系数；
[0079] Kp—大角度姿态机动控制律的比例系数；
[0080] 步骤4.3，将控制力矩Lc作用于飞轮，使其以最快速度按下式规划的路径进行姿态机动：
[0081]
[0082] 其中：
[0083] ωc(t)—期望角速度关于时间的函数；
[0084] —期望角加速度关于时间的函数；
[0085] qc(t)—期望四元数关于时间的函数；
[0086] —瞬时欧拉轴；
[0087] —绕欧拉轴旋转角加速度关于时间的函数；
[0088] —绕欧拉轴旋转角速度关于时间的函数；
[0089] θ(t)—绕欧拉轴旋转角度大小关于时间的函数。
[0090] 优选的，还包括：
[0091] 步骤5，建立相对运动动力学方程；
[0092] 相对运动动力学方程为：
[0093]
[0094] 其中：
[0095] ρ—环绕星相对参考星的相对位置矢量；
[0096] —环绕星相对参考星的相对速度矢量；
[0097] —环绕星相对参考星的相对加速度矢量；
[0098] ωn—参考卫星的轨道角速度；
[0099] r—参考星在地惯系中的位置矢量；
[0100] um—电磁控制相对加速度；
[0101] R—干扰项和未建模的动力学部分；
[0102] C(ωn)∈R3×3为
[0103]
[0104] R3×3—三维实数集；
[0105] 非线性矩阵N(·)∈R3×3定义为如下等式：
[0106]
[0107] μE—地球引力常数；
[0108] x、y、z—环绕星与参考星的相对位置矢量在Hill坐标系下的分量；
[0109] 步骤6，根据步骤5建立的相对运动动力学方程，设计得到电磁编队飞行卫星运动控制器；
[0110] 所述电磁编队飞行卫星运动控制器是一种终端滑模控制器，用于实现轨道
跟踪误差保持在零附近。
[0111] 优选的，轨道跟踪误差可被定义为
[0112]
[0113] 其中：
[0114] e—轨道跟踪位置误差；
[0115] —轨道跟踪速度误差；
[0116] ρd—目标相对位置矢量；
[0117] —目标相对速度矢量；
[0118] 因此，设计的终端滑模控制器描述为：
[0119]
[0120] 其中：β＞0，ki＞0，i＝1,2,3，0＜γ＜1；
[0121] us—终端滑模控制相对加速度；
[0122] —目标相对加速度矢量；
[0123] β—控制系数β＞0；
[0124] γ—控制系数0＜γ＜1；
[0125] k—控制系数ki＞0，i＝1,2,3；
[0126] sat(·)—饱和函数，
[0127] φ0—一严格正常数，表示所设计的
边界层的厚度；
[0128] S—滑模面，滑模面可描述为
[0129]
[0130] |S|—滑模面的模。
[0131] 优选的，为了消除振动，用连续的饱和函数代替了公式(51)中的饱和函数，连续的饱和函数可描述为
[0132]
[0133] 其中，设计系数τ满足0＜τ＜1；公式(51)可被写为
[0134]
[0135] 由此得到终端滑模控制器的另一种描述。
[0136] 优选的，还包括：
[0137] 步骤7，建立电磁力和电磁力矩的模型，包括：
[0138] 步骤7.1，卫星j对卫星i的电磁力Fij和电磁力矩Tij可被描述为
[0139]
[0140]
[0141] 其中：
[0142] μ0—
真空磁导率，μ0＝4π×10-7H/m；
[0143] rij—由卫星i到卫星j的相对位置矢量；
[0144] rij—由卫星i到卫星j的相对位置矢量的模；
[0145] μi—卫星i的磁偶极子的强度，μi＝NiIiSi；
[0146] μj—卫星j的磁偶极子的强度，μj＝NjIjSj；
[0147] Ni、Nj分别为偶极子A、B的绕线
匝数；
[0148] Ii、Ij分别为偶极子A、B的线圈
电流；
[0149] Si、Sj分别为偶极子A、B的线圈面积；
[0150] 由公式(55)和公式(56)，可得到如下三
自由度的远场电磁模型：
[0151]
[0152] 其中：
[0153] —偶极子A受到的电磁力在三个轴方向的分量；
[0154] —偶极子A受到的电磁力矩在三个轴方向的分量；
[0155] μA—卫星A的磁偶极子强度；
[0156] μB—卫星B的磁偶极子强度；
[0157] α、β分别为偶极子A、B与ox轴之间的夹角；
[0158] χ、δ分别为偶极子A、B绕ox轴的转角；
[0159] 步骤8，在步骤6设计的电磁编队飞行卫星运动控制器的基础上，基于步骤7的电磁力和电磁力矩的模型，为获取卫星间相对位置的控制过程中每颗卫星的磁偶极子，还设计两种计算磁偶极子的方法：
[0160] 第一种是根据
牛顿第三定律，将电磁编队飞行卫星的某一偶极子固定，计算可得其他偶极子的电磁力；
[0161] 第二种是将偶极子的计算问题转换为一个优化问题，利用一种粒子群算法和序列二次规划方法混合算法以实现自由偶极子的最优化，最终实现
电磁干扰力矩MT和
地磁场干扰力矩ME最小的优化目标。
[0162] 优选的，步骤7中，第一种计算磁偶极子的方法包括：
[0163] 根据牛顿第三定律，电磁编队飞行卫星运动控制器所对应的公式(54)代表针对n个卫星的3(n-1)个标量公式，因此，系统在此约束下方程会有无穷多解，如果将电磁编队飞行卫星的某一偶极子固定，其他偶极子就可以通过如下公式进行计算
[0164]
[0165] 其中：
[0166] μ0—真空磁导率，μ0＝4π×10-7H/m；
[0167] rij—由卫星i到卫星j的相对位置矢量；
[0168] rij—由卫星i到卫星j的相对位置矢量的模；
[0169] μi—卫星i的磁偶极子强度；
[0170] μj—卫星j的磁偶极子强度；
[0171] 步骤7中，第二种计算磁偶极子的方法为：
[0172] 建立以下目标函数：
[0173]
[0174] μoptimal—最优磁矩；
[0175] arg(·)—幅值；
[0176] n—卫星个数；
[0177] 采用粒子群算法和序列二次规划方法混合算法以实现自由偶极子的最优化，其中，首先采用粒子群算法设定初始参数值，然后采用粒子群算法进行渐进估算，当收敛值满足预先设定的临界值时，将收敛值提供给序列二次规划方法，采用序列二次规划方法进行非线性约束下问题的非线性解算。
[0178] 本发明提供的电磁编队卫星姿态和轨道相对控制具有以下优点：
[0179] 本发明提供的电磁编队卫星姿态和轨道相对控制方法，充分考虑了电磁编队卫星姿态和轨道相对控制方法设计每个阶段需要考虑的因素，建立了姿态动力学方程、相对运动动力学方程及电磁力和电磁力矩的模型，由此设计的卫星姿控系统反馈控制律、飞轮角动量磁卸载算法、大角度姿态机动方法、电磁编队飞行卫星运动控制器和磁偶极子的控制分配方法具有较好的工作性能，能够实现高效的电磁编队卫星姿态和轨道相对控制。
附图说明
[0180] 图1为本发明提供的电磁编队卫星姿态和轨道相对控制方法的流程示意图；
[0181] 图2为远场电磁模型；
[0183] 图4为角加速度随时间变化的曲线；
具体实施方式
[0185] 为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及
实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0186] 结合图1，本发明提供一种电磁编队卫星姿态和轨道相对控制方法，在建立姿态动力学方程、建立相对运动动力学方程和建立电磁力和电磁力矩的模型的基础上，提出以下创新点：1、设计一种卫星姿控系统反馈控制律；2、设计了一种飞轮角动量磁卸载算法；3、设计了一种大角度姿态机动方法；4、设计一种电磁编队飞行卫星运动控制器；5、设计一种磁偶极子的控制分配。按照本发明设计方法，能够根据电磁编队卫星的姿态控制方法，实现电磁编队卫星的对日定向；同时还能够实现磁卸载，电磁编队卫星
相位位置保持的非线性控制，以及磁矩分配。下面对上述五大部分分别详细介绍：
[0187] (一)卫星姿控系统反馈控制律
[0188] 步骤1：建立姿态动力学方程；姿态动力学方程包括：用四元数描述的线性形式的姿态动力学方程、体坐标系下的姿态动力学方程和刚体卫星的动力学方程，其建立过程包括以下步骤：
[0189] S1.1，用四元数描述航天器的姿态，将四元数描述的姿态动力学方程描述为线性形式；
[0190] 用四元数描述航天器的姿态为：
[0191] q＝(q1 q2 q3 q4)T (1)
[0192] 其中：
[0193] q—航天器姿态四元数；
[0194] q1＝cos(φ/2)；
[0195] q2＝exsin(φ/2)；
[0196] q3＝eysin(φ/2)；
[0197] q4＝ezsin(φ/2)；
[0198] φ—欧拉轴的旋转角；
[0199] [ex ey ez]—欧拉轴的单位矢量；
[0200] T—矩阵的转置；
[0201] 四元数满足如下约束：
[0202]
[0203] 因此，四元数描述下的姿态动力学可以描述为公式(6)所示的线性形式：
[0204]
[0205] 其中：
[0206] —描述航天器姿态的四元数的一阶导数；
[0207] ωxb、ωyb、ωzb分别代表体坐标系的三个轴上的转动角速度；
[0208] S1.2，根据轨道坐标系到体坐标系的转换矩阵，得到体坐标系下的卫星角动量；
[0209] 由轨道坐标系到体坐标系的转换矩阵BG为：
[0210]
[0211] 其中：
[0213] θ—航天器的俯仰角；
[0215] 且三个角均为欧拉角；
[0216] 由于欧拉角很小，因此上式可以近似为：
[0217]
[0218] 惯性坐标系下的卫星角动量为HI＝II·ωI，
[0219] II—惯性坐标系下的转动惯量；
[0220] ωI—惯性坐标系下的角速度；
[0221] 下标I—惯性坐标系下的向量；
[0222] S1.3，根据角动量定理，得到体坐标系下滚
转轴的动力学方程；
[0223] 根据刚体的角动量定理，得到如下等式：
[0224]
[0225] LI—惯性坐标系下的力矩；
[0226] 定义惯性坐标系到体坐标系的转换矩阵为BI，那么体坐标系下的转动惯量IB、体坐标系下的角速度ωB和体坐标系下的对
应力矩LB分别为：
[0227] IB＝BIIIBIT，ωB＝BI·ωI，LB＝BI·LI (7)
[0228] 因此，体坐标系下的卫星角动量HB为：
[0229] HB＝IB·ωB＝BIIIBITBI·ωI＝BI·HI (8)
[0230] 根据角动量定理，能够得到：
[0231]
[0232] 因此：
[0233]
[0234] Δθ、Δψ—航天器的滚转角、俯仰角、偏航角的角度变化量；则体坐标系下滚转轴的动力学方程为等式(14)，其中变量都定义在体坐标系下：
[0235]
[0236] 其中：
[0237] 体坐标系下的卫星角动量
[0238] 体坐标系下的力矩
[0239] 体坐标系下的卫星转动角速度
[0240] 其中：Ixb、Iyb、Izb分别为体坐标系下三个轴的转动惯量；
[0241] Lxb、Lyb、Lzb分别为体坐标系下三个轴的力矩；
[0242] Hxb、Hyb、HZb分别为体坐标系下三个轴的角动量；
[0243]
[0244] 因此，在体坐标系下，姿态动力学方程可被写作：
[0245]
[0246] 其中：
[0247] IT—卫星及飞轮整体的转动惯量；
[0248] —卫星的角加速度；
[0249] —飞轮的角加速度；
[0250] Iw—反作用飞轮的转动惯量；
[0251] ωB—卫星的角速度；
[0252] ωw—体坐标系下飞轮的角速度；
[0253] LB—体坐标系下的力矩；
[0254] S1.4，由于电磁编队飞行卫星的姿态能够通过三个正交的反作用飞轮控制，最后推导了刚性卫星的动力学公式；
[0255] 电磁编队飞行卫星的姿态能够通过三个正交的反作用飞轮控制，这三个反作用飞轮可以沿体坐标系的三个坐标轴安装；设反作用飞轮的转动惯量为Iw，三个轴的转动惯量分别为Iwx，Iwy和Iwz，体坐标系下飞轮的角速度为ωw，卫星的角速度为ωB，卫星及飞轮整体的转动惯量为IT；因此，
[0256] HB＝ITωB+Iwωw (13)
[0257] 那么，刚性卫星的动力学方程可被写作：
[0258]
[0259] 其中：
[0260] IT—卫星及飞轮整体的转动惯量；
[0261] —卫星的角加速度；
[0262] —飞轮的角加速度；
[0263] Iw—反作用飞轮的转动惯量；
[0264] ωB—卫星的角速度；
[0265] ωw—体坐标系下飞轮的角速度；
[0266] LB—体坐标系下的力矩；
[0267] 步骤2：根据步骤1建立的四元数描述的线性形式的姿态动力学方程、体坐标系下的姿态动力学方程和刚体卫星的动力学方程，设计得到卫星姿控系统反馈控制律；
[0268] 本步骤具体包括：
[0269] 步骤2.1：假定当前四元数为q，期望四元数为qc。那么，误差四元数qe的定义为：
[0270]
[0271] 假定当前角速度为ω，期望角速度为ωc；那么，误差角速度ωe的定义为：
[0272] ωe＝ω-ωc (2)
[0273] 步骤2.2，基于误差四元数和误差角速度的姿态动力学模型可被描述为[0274]
[0275] I—整星的惯量矩阵；
[0276] —误差角加速度；
[0277] —期望角加速度；
[0278] H—整星的角动量，H＝IBωB+Iwωw；
[0279] IB—卫星的转动惯量；
[0280] u—反作用飞轮的控制力矩，u＝Iwωw；
[0281] Md—外部干扰力矩；
[0282] ME—地磁场引起的干扰力矩；
[0283] —误差四元数的一阶导数；
[0284] 为便于计算ME，将地球磁场近似为一个以地球为中心的偶极子；非线性控制系统的非线性部分能够利用反馈线性方法被抵消，从而提出反馈控制律：
[0285]
[0286] 其中：
[0287] P，D—正定增益矩阵，通常被设计为对角矩阵；
[0288] qe123—误差四元数中的三个分量qe1、qe2、qe3组成的列向量，qe123＝[qe1 qe2 qe3]T；
[0289] 姿态控制系统中，在反馈线性化后的闭环控制是一个典型的二阶线性系统，因此能够确定增益矩阵P和D；卫星姿态控制回路中包括动力学模型、感应器、控制器和反作用飞轮模型。以俯仰控制回路为例，该闭环回路的框架图如图3所示；
[0290] 反作用飞轮和陀螺的模型可以视为传递函数，反作用飞轮的传递函数为惯性环节：
[0291]
[0293] s—拉氏算子；
[0294] 陀螺的传递函数Wgyro为振荡环节：
[0295]
[0296] ωg—陀螺的振荡角频率；
[0297] ζg—陀螺的阻尼系数；
[0298] 在图3中，姿态控制系统开环传递函数由动力学模型的传递函数、陀螺的传递函数、反作用飞轮的传递函数和PD控制环节组成，系统的开环传递函数Wopen可表示为[0299]
[0300] p—PD控制的环节比例系数；
[0301] d—PD控制的环节微分系数；
[0302] Iy—干扰力矩产生的转动惯量；
[0303] 系统的闭环传递函数Wclose为
[0304]
[0305] 外部干扰力矩Ty(s)到姿态输出θ(s)的闭环传递函数Wdisturb为
[0306]
[0307] 从而能够得到输入θc与外部干扰力矩Ty、输出θout的简化关系式
[0308]
[0309] 如公式(24)所示，p的增加能够抑制干扰对姿态的影响；与之相反，p的增加导致高频段的增益增大，这会影响弹性振动的
稳定性；因此，参数p需要合理设计。
[0310] (二)飞轮角动量磁卸载算法
[0311] 步骤3：根据步骤2提出的反馈控制律，对电磁编队卫星姿态和轨道进行相对控制。具体为：在根据步骤2提出的反馈控制律对电磁编队卫星姿态和轨道进行相对控制时，在控制过程中，实时判断是否发生磁饱和现象，如果发生，采用飞轮角动量磁卸载算法卸载磁力矩。
[0312] 电磁力矩是由卫星磁偶极子的相互作用和地球磁场产生的，对于低轨
小卫星，可以利用反作用飞轮对磁力矩器进行磁卸载；假定反作用飞轮的冗余转动惯量为ΔH，则其磁动量为
[0313]
[0314] K—增益系数；
[0315] B—飞轮所处磁场；
[0316] |B|—飞轮所处磁场B的大小；
[0317] 飞轮的磁力矩Tw为：
[0318]
[0319] 如果B垂直于ΔH，那么
[0320] Tw＝-KΔH(K＞0) (27)
[0321] 如果B不垂直于ΔH，B可以被分解为沿着ΔH的方向和垂直于ΔH的方向两部分，分别记为BΔH和BN。那么对卸载有效的力矩可写为
[0322]
[0323] 对卸载无效的力矩为
[0324]
[0325] 由等式(29)可得，如果B和ΔH之间的角度很小，那么不适合磁卸载；磁卸载通常在|b·Δh|＜0.707时进行，其中b和Δh分别是B和ΔH的单位向量；因此，得到磁卸载条件：
[0326] 45°＜(B,ΔH)＜135° (30)
[0327] 如果将比例反馈控制律用于磁卸载，那么会得到陀螺速度ωb和体坐标系下飞轮的角速度ωw；因此，整体的角动量可被表示为
[0328] H＝I·ωb+Iw·ωw (31)
[0329] 轨道的角速度为ω0，则理论角动量可被写为
[0330] H0＝I·ω0T (32)
[0331] 因此，多余的角动量为
[0332] ΔH＝H-H0≈Iw·ωw (33)
[0333] 通过如下比例反馈控制律，可以得到磁卸载的指令：
[0334] Tw＝-KΔH (34)
[0335] 由于Tw＝M×B，因此
[0336] B×T＝B×(M×B)＝M(B·B)-B(B·M) (35)
[0337] 如果M与B平行，那么不能够产生力矩，即B·M＝0；因此，
[0338] M＝B×Tw/B2 (36)
[0339] 磁卸载的施加的力矩为
[0340]
[0341] 如果B垂直于ΔH，则按下式(38)计算得到磁卸载力矩值Tunloading，并将该磁卸载力矩值Tunloading施加于飞轮中，使其反向转动，进行磁卸载：
[0342] Tunloading＝-KΔH(K＞0) (38)
[0343] (三)设计大角度姿态机动方法
[0344] 为使
太阳能帆板产生足够电流，调整电磁编队飞行卫星的姿态使其保持对日定向是最好的方法；传统的过程是调整姿态，使其由三轴稳定到对日定向；反作用飞轮用于控制姿态，因此，最大力矩和最大转动速度会导致大角度机动产生新问题；在飞轮的最大力矩和转动速度的限制下，基于欧拉旋转理论，能够设计出最理想的机动算法；此外，根据误差四元数和角速度的反馈控制律能够实现跟踪控制；
[0345] 步骤4.1，假设机动前，即初始时刻t0时刻的理论姿态为q(t0)，理论角速度为ω(t0)＝[0 0 0]T；机动后，即tf时刻目标姿态为q(tf)，目标角速度为ω(tf)＝[0 0 0]T；
[0346] 步骤4.2，根据欧拉定理，通过将瞬时欧拉轴旋转θ角，能够实现由q(t0)到q(tf)的姿态机动，且机动路径是最短的；为了满足最大力矩和转动速度的约束，机动路径可以规划为
[0347]
[0348] Lc—飞轮的转动力矩；
[0349] —飞轮的转动速度；
[0350] —飞轮的转动速度；
[0351] —绕欧拉轴旋转角加速度；
[0352] —瞬时欧拉轴；
[0353] 假定飞轮的最大输出力矩为Lcmax，那么绕欧拉轴的最大角加速度为
[0354]
[0355] 为以最快速度实现姿态机动，机动路径可以规划为
[0356]
[0357] 其中：
[0358] ωc(t)—期望角速度关于时间的函数；
[0359] —期望角加速度关于时间的函数；
[0360] qc(t)—期望四元数关于时间的函数；
[0361] —瞬时欧拉轴；
[0362] —绕欧拉轴旋转角加速度关于时间的函数；
[0363] —绕欧拉轴旋转角速度关于时间的函数；
[0364] θ(t)—绕欧拉轴旋转角度大小关于时间的函数。
[0365] 为满足飞轮的磁饱和约束，机动过程要满足如图4所示变化律；
[0366] 如图4所示，确定t1和tf以确定机动路径；通过忽略机动期间的干扰，可以保持角动量；因此，
[0367]
[0368] θf—机动的角度；
[0369] 为了确保飞轮在机动过程中不会达到磁饱和，应确定机动时间tf；假定最大角速度为 那么
[0370]
[0371]
[0372] t1—从初始时刻开始以最大角加速度 机动的结束时刻；
[0373] tf—第二次以最大角加速度 机动的结束时刻，机动角加速度与t0～t1相反；
[0374] 其中，t1和tf的计算过程流程图如图5所示；
[0375] 基于动量守恒定律能够得到ωwt的计算方程
[0376]
[0377] ωwt—t时刻飞轮的角速度；
[0378] ωBt—体坐标系下卫星t时刻的角速度；
[0379] ωw0—初始时刻飞轮的角速度；
[0380] ωB0—体坐标系下卫星初始时刻的角速度；
[0381] θt—t时刻机动的角度；
[0382] —t时刻机动的角速度；
[0383] 那么，控制律可写为
[0384]
[0385] 其中：
[0386] IT—卫星及飞轮整体的转动惯量；
[0387] q123—q1、q2、q3组成的列向量，q123＝[q1 q2 q3]T；
[0388] Kd—大角度姿态机动控制律的微分系数；
[0389] Kp—大角度姿态机动控制律的比例系数；
[0390] (四)电磁编队飞行卫星运动控制器
[0391] 步骤5，建立相对运动动力学方程；
[0392] 相对运动动力学方程为：
[0393]
[0394] 其中：
[0395] ρ—环绕星相对参考星的相对位置矢量；
[0396] —环绕卫星相对参考卫星的相对速度矢量；
[0397] —环绕卫星相对参考卫星的相对加速度矢量；
[0398] ωn—参考卫星的轨道角速度；
[0399] r—参考星在地惯系中的位置矢量；
[0400] um—电磁控制相对加速度；
[0401] R—干扰项和未建模的动力学部分；
[0402] C(ωn)∈R3×3为
[0403]
[0404] R3×3—三维实数集；
[0405] 非线性矩阵N(·)∈R3×3定义为如下等式：
[0406]
[0407] μE—地球引力常数；
[0408] x、y、z—环绕星与参考星的相对位置矢量在Hill坐标系下的分量；
[0409] 步骤6，根据步骤5建立的相对运动动力学方程，设计得到电磁编队飞行卫星运动控制器；
[0410] 所述电磁编队飞行卫星运动控制器是一种终端滑模控制器，用于实现轨道跟踪误差保持在零附近。轨道跟踪误差可被定义为
[0411]
[0412] 其中：
[0413] e—轨道跟踪位置误差；
[0414] —轨道跟踪速度误差；
[0415] ρd—目标相对位置矢量；
[0416] —目标相对速度矢量；
[0417] 因此，设计的终端滑模控制器描述为：
[0418]
[0419] 其中β＞0，ki＞0，i＝1,2,3，0＜γ＜1；
[0420] us—终端滑模控制相对加速度；
[0421] —目标相对加速度矢量；
[0422] β—控制系数β＞0；
[0423] γ—控制系数0＜γ＜1；
[0424] k—控制系数ki＞0，i＝1,2,3；
[0425] sat(·)—饱和函数，
[0426] φ0—一严格正常数，表示所设计的边界层的厚度；
[0427] S—滑模面，滑模面可描述为
[0428]
[0429] |S|—滑模面的模。
[0430] 注：为了消除振动，用连续的饱和函数代替了公式(51)中的饱和函数；
[0431]
[0432] 其中，设计系数τ满足0＜τ＜1；公式(51)可被写为
[0433]
[0434] 由此得到终端滑模控制器的另一种描述。
[0435] (五)磁偶极子的控制分配
[0436] 步骤7，建立电磁力和电磁力矩的模型；
[0437] 本发明提出了电磁力和电磁力矩的远场模型，如图5所示。远场模型在卫星距离大于线圈半径的6.76倍时是有效的，此时远场模型的误差小于10％；
[0438] 卫星j对卫星i的电磁力Fij和电磁力矩Tij可被描述为
[0439]
[0440]
[0441] μ0—真空磁导率，μ0＝4π×10-7H/m；
[0442] rij—由卫星i到卫星j的相对位置矢量；
[0443] rij—由卫星i到卫星j的相对位置矢量的模；
[0444] μi—卫星i的磁偶极子的强度，μi＝NiIiSi；
[0445] μj—卫星j的磁偶极子的强度，μj＝NjIjSj；
[0446] Ni、Nj分别为偶极子A、B的绕线匝数；
[0447] Ii、Ij分别为偶极子A、B的线圈电流；
[0448] Si、Sj分别为偶极子A、B的线圈面积；
[0449] 由公式(55)和公式(56)，可得到如下三自由度的远场电磁模型：
[0450]
[0451] —偶极子A受到的电磁力在三个轴方向的分量；
[0452] —偶极子A受到的电磁力矩在三个轴方向的分量；
[0453] μA—卫星A的磁偶极子强度；
[0454] μB—卫星B的磁偶极子强度；
[0455] α、β分别为偶极子A、B与ox轴之间的夹角；
[0456] χ、δ分别为偶极子A、B绕ox轴的转角；
[0457] 在公式(57)中，电磁力和电磁力矩由两颗卫星的磁偶极子乘积决定，并且随两卫星的相对距离和角度变化；电磁力与相对距离的四次方成反比，而电磁力矩与相对距离的三次方成反比；因此，电磁力和电磁力矩是强非线性化和强耦合的。
[0458] 步骤8，在步骤6设计的电磁编队飞行卫星运动控制器的基础上，基于步骤7的电磁力和电磁力矩的模型，为获取卫星间相对位置的控制过程中每颗卫星的磁偶极子，还设计两种计算磁偶极子的方法：
[0459] 第一种是根据牛顿第三定律，将电磁编队飞行卫星的某一偶极子固定在，计算可得其他偶极子的电磁力；
[0460] 根据牛顿第三定律，电磁编队飞行卫星运动控制器所对应的公式(54)代表针对n个卫星的3(n-1)个标量公式，因此，系统在此约束下方程会有无穷多解，如果将电磁编队飞行卫星的某一偶极子固定，其他偶极子就可以通过如下公式进行计算
[0461]
[0462] 其中，μ0—真空磁导率，μ0＝4π×10-7H/m；
[0463] rij—由卫星i到卫星j的相对位置矢量；
[0464] rij—由卫星i到卫星j的相对位置矢量的模；
[0465] μi—卫星i的磁偶极子强度；
[0466] μj—卫星j的磁偶极子强度；
[0467] 第二种是将偶极子的计算问题转换为一个优化问题，利用一种粒子群算法和序列二次规划方法(PSO-SQP)混合算法以实现自由偶极子的最优化，最终实现电磁干扰力矩MT和地磁场干扰力矩ME最小的优化目标。
[0468] 建立以下目标函数：
[0469]
[0470] μoptimal—最优磁矩；
[0471] arg(·)—幅值；
[0472] n—卫星个数；
[0473] 本发明提出的粒子群算法和序列二次规划方法(PSO-SQP)混合算法以实现该自由偶极子的最优化，其中粒子群算法是一种渐进估算算法，而序列二次规划方法是一种能够高度适用于非线性约束下问题的非线性方法；这一混合算法利用粒子群算法设定初始参数值，将其提供给序列二次规划方法；其优点在于序列二次规划方法能够减少计算时间并且提高收敛性能；当收敛值满足预先设定的临界值时，应用序列二次规划方法。
[0474] 本发明提供的电磁编队卫星姿态和轨道相对控制方法，充分考虑了电磁编队卫星姿态和轨道相对控制方法设计每个阶段需要考虑的因素，建立了姿态动力学方程、相对运动动力学方程及电磁力和电磁力矩的模型，由此设计的卫星姿控系统反馈控制律、飞轮角动量磁卸载算法、大角度姿态机动方法、电磁编队飞行卫星运动控制器和磁偶极子的控制分配方法具有较好的工作性能，能够实现高效的电磁编队卫星姿态和轨道相对控制。按照本发明提出的电磁编队卫星姿态和轨道相对控制方法，能够实现电磁编队卫星的对日定向；同时还能够实现磁卸载，电磁编队卫星相位位置保持的非线性控制，以及磁矩分配。具体具有以下优点：
[0475] (1)在飞轮的最大力矩和转动速度的限制下，基于欧拉旋转理论，能够设计出最理想的大角度姿态。此外，根据误差四元数和角速度的反馈控制律能够实现跟踪控制。将反作用飞轮用于控制姿态，能够实现最短机动路径的姿态机动，并满足飞轮的磁饱和约束，最终确定机动过程的控制律，实现姿态控制对日定向和大角度机动。
[0476] (2)提出了一种终端滑模控制器来实现轨道跟踪误差保持在零附近。设计的终端滑模控制器能够使给出的相对运动系统是大范围渐进稳定的，具有鲁棒性。
[0477] (3)提出两种计算磁偶极子的方法：一种是根据牛顿第三定律，将某一偶极子固定在电磁编队飞行卫星上，计算可得其他偶极子的电磁力；另一种是将偶极子的计算问题转换为一个优化问题，利用一种粒子群算法和序列二次规划方法(PSO-SQP)混合算法以实现该自由偶极子的最优化，最终实现飞轮力矩和地球磁场干扰力矩最小化的优化目标；这两种方法均能够实现磁矩的最优控制分配。
[0478] 以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。
