CN108958275A - 一种刚柔液耦合系统姿态控制器和机动路径联合优化方法 - Google Patents

Abstract

本发明提出了一种刚柔液耦合系统姿态控制器和机动路径联合优化方法，包括以下步骤：首先建立充液挠性航天器的动力学模型；再获取充液挠性航天器的角加速度曲线、角速度曲线、得到多段的角位置曲线，对充液挠性航天器姿态机动路径进行规划；再计算角位置多段曲线中每段曲线的表达式；采用PD控制系统对充液挠性航天器进行姿态控制；联合优化充液挠性航天器控制器和机动路径的参数：采用基于自适应网格的多目标粒子群优化算法对充液挠性航天器的控制器和机动路径参数进行联合优化。本发明的方法减少了姿态机动对挠性附件振动和液体晃动的激发，实现了充液挠性航天器姿态大角度快速机动快速稳定控制。

Description

一种刚柔液耦合系统姿态控制器和机动路径联合优化方法

技术领域

本发明属于航天器姿态控制领域，特别是一种刚柔液耦合系统姿态控制器和机动路径联合优化方法。

背景技术

为完成复杂任务并延长在轨运行寿命，现代航天器需要携带大量液体燃料和轻质附件，这种充液挠性航天器的姿控系统具有刚柔液耦合的复杂非线性动力学特性，姿态机动极易诱发液体晃动和附件振动，这对实现航天器姿态大角度快速机动、快速稳定控制提出了巨大挑战。

对充液挠性航天器动力学的研究表明，在姿态快速机动过程中，姿态角加速度的突变(即不连续性)常导致挠性帆板的强烈振动和液体的强烈晃动，同时，在机动过程中角加速度、角速度以及控制力矩易受到执行机构的限制。因此，需对航天器姿态机动的路径进行规划，寻找最优的姿态控制器参数及机动路径参数，以减少姿态机动对挠性附件振动和液体晃动的激发。

目前迫切需要一种刚柔液耦合系统姿态控制器和机动路径联合优化方法，使得充液挠性航天器能够快速完成大角度姿态机动并获得高的稳定度。但是现有技术中尚无相关描述。

发明内容

本发明所解决的技术问题在于提供一种刚柔液耦合系统姿态控制器和机动路径联合优化方法，以减少姿态机动对挠性附件振动和液体晃动的激发，实现充液挠性航天器姿态大角度快速机动快速稳定控制。

实现本发明目的的技术方案为：

一种刚柔液耦合系统姿态控制器和机动路径联合优化方法，包括以下步骤：

步骤1、建立充液挠性航天器的动力学模型：一类带有挠性太阳能帆板及液体燃料的三轴航天器的姿态动力学方程、挠性振动方程和液体燃料晃动方程；

步骤2、获取充液挠性航天器的角加速度曲线、角速度曲线、得到多段的角位置曲线，对充液挠性航天器姿态机动路径进行规划；

步骤3、计算角位置多段曲线中每段曲线的表达式；利用角加速度最大值、角加速度微分曲线中非零段时间和的一半或角加速度曲线中非零段时间和的一半、航天器角速度最大值、期望机动的角度，解算出航天器每段曲线中角位置的表达式；

步骤4、采用PD控制系统对充液挠性航天器进行姿态控制；

步骤5、联合优化充液挠性航天器控制器和机动路径的参数：采用基于自适应网格的多目标粒子群优化算法对充液挠性航天器的控制器和机动路径参数进行联合优化。

本发明与现有技术相比，其显著优点：

(1)对充液挠性航天器的大角度姿态机动路径进行规划，结合PD控制器，使得充液挠性航天器快速机动到稳定状态。

(2)对控制器参数及机动路径参数进行优化，寻找出较优的姿态控制器参数及机动路径参数，减少挠性附件振动和液体晃动的对充液挠性航天器姿态控制的影响，提高系统机动完成后的稳定性。

下面结合附图对本发明作进一步详细描述。

附图说明

图1为本发明方法的流程图。

图2为基于正弦型角加速度的分段路径图。

图3为充液挠性航天器姿态控制系统结构示意图。

图4为基于自适应网格的多目标粒子群优化算法流程图。

图5(a-c)分别为实施例中航天器滚动轴、俯仰轴、偏航轴姿态角位置仿真结果图。

图6(a-c)分别为实施例中航天器滚动轴、俯仰轴、偏航轴姿态角位置误差仿真结果图。

图7(a-c)分别为实施例中航天器滚动轴、俯仰轴、偏航轴姿态角速度误差仿真结果图。

具体实施方式

为了说明本发明的技术方案及技术目的，下面结合附图及具体实施例对本发明做进一步的介绍。

结合图1，本发明的一种刚柔液耦合系统姿态控制器和机动路径联合优化方法，包括以下步骤：

步骤1、建立充液挠性航天器的动力学模型：一类带有挠性太阳能帆板及液体燃料的三轴航天器的姿态动力学方程、挠性附件振动方程和液体燃料晃动方程，分别为：

其中：I∈R^3×3为航天器转动惯量；ω∈R^3×1为航天器角速度；ω^×∈R^3×3表示ω的斜对称矩阵；C₀∈R^3×n为挠性振动运动与航天器转动运动耦合系数矩阵； G₁∈R^2×m，C₁∈R^m×3，C₂∈R^m×3均为液体晃动运动与航天转动之间的耦合系数矩阵；τ∈R^3×1为控制力矩矢量；T_d∈R^3×1为环境干扰力矩；ζ∈R^n×n为挠性附件振动模态阻尼比；Λ∈R^n×n为挠性附件振动模态频率矩阵；Λ₂∈R^m×m为液体晃动模态频率平方矩阵；η∈R^n×1为挠性体振动模态矢量；q∈R^m×1为液体晃动模态矢量；α为航天器的姿态角矢量，

步骤2、获取充液挠性航天器的角加速度曲线、角速度曲线、得到多段的角位置曲线，对充液挠性航天器姿态机动路径进行规划；

优选的，对角加速度微分曲线划分为多段并对每段曲线进行积分，即得整个机动路径的角加速度多段曲线，再对角加速度曲线积分，即得整个机动路径的角速度多段曲线；再对角速度曲线积分，即得整个机动路径的角位置多段曲线，实现对充液挠性航天器姿态机动路径的规划。

优选的，结合图2，采用基于正弦型角加速度的分段路径，将角加速度微分曲线按波形分成7段，路径的角加速度微分曲线(图2-I)是由正弦函数组成的，I-①和I-⑦是正弦函数的正半周期，I-③和I-⑤是正弦函数的负半周期，I-②、I-④和I-⑥为0；对这七段角加速度微分曲线(图2-I)积分，即得整个机动路径的角加速度曲线(图2-Ⅱ)，对角加速度曲线(图2-Ⅱ)积分，即得整个机动路径的角速度曲线(图2-Ⅲ)，再对角速度曲线(图 2-Ⅲ)积分，即得整个机动路径的角位置曲线(图2-Ⅳ)；

在一些实施方式中，也可直接对角加速度曲线积分，即得整个机动路径的角速度曲线；再对角速度曲线积分，即得整个机动路径的角位置曲线；实现对充液挠性航天器姿态机动路径的规划。

步骤3、计算角位置多段曲线中每段曲线的表达式：利用角加速度最大值、角加速度微分曲线中非零段时间和的一半或角加速度曲线中非零段时间和的一半、航天器角速度最大值、期望机动的角度，解算出航天器每段曲线中角位置的表达式。

步骤3.1、定义j轴期望机动的角度其中，j＝x,y,z分别表示航天器的滚动轴、俯仰轴、偏航轴；

步骤3.2、定义j轴航天器角速度最大值V_{max_j}；

步骤3.3、获取j轴两个机动路径参数：角加速度最大值a_{max_j}、角加速度微分曲线中非零段时间和的一半T_j或角加速度曲线中非零段时间和的一半T_j；

步骤3.4、解算出航天器每段曲线中角位置的表达式。

结合步骤2优选的实施例，确定j轴(j＝x,y,z，分别表示航天器的滚动轴、俯仰轴、偏航轴，下文同理)机动到指定角度确定j轴航天器角速度最大值V_{max_j}；采用基于正弦型角加速度的分段路径，将角加速度微分曲线按波形分成7段，确定j轴角速度最大值a_{max_j}和图2-I中正弦函数的周期T_j，图中Ⅳ-②和Ⅳ-⑥时间相同，均为t_j1；Ⅳ-④时间为t_j2，可以根据a_{max_j}、T_j、V_{max_j}算出再根据a_{max_j}、T_j、V_{max_j}、t_j1算出从而算出角位置(图2-Ⅳ)的表达式为：

步骤4、结合图3，采用PD控制系统对充液挠性航天器进行姿态控制。PD控制律为：

其中，q_d为规划的欧拉角角位置输出，为规划的欧拉角角速度输出，为规划的欧拉角角加速度输出；q_p为实际的欧拉角角位置输出，为实际的欧拉角角速度输出； e＝q_d-q_p为欧拉角角位置误差向量；为欧拉角角速度误差向量。K_p为比例系数矩阵，K_d为微分系数矩阵，K_p＝diag(K_px,K_py,K_pz)，K_d＝diag(K_dx,K_dy,K_dz)，其中 K_px、K_py、K_pz分别表示航天器滚动轴、俯仰轴、偏航轴的控制器的比例系数，K_dx、K_dy、 K_dz分别表示航天器滚动轴、俯仰轴、偏航轴的控制器的微分系数。

然后对输出进行限幅，控制力矩向量为u_limit为实际系统所允许的最大控制力矩输出向量；

步骤5、联合优化充液挠性航天器控制器和机动路径的参数：采用基于自适应网格的多目标粒子群优化算法对充液挠性航天器的控制器和机动路径参数进行联合优化；该算法通过外部粒子群以指导群体外其他粒子的飞行，并利用Pareto档案进化策略(PAES) 中的自适应网格法进行外部粒子群的维护，以平衡全局搜索和局部搜索能力。基于自适应网格的多目标粒子群优化算法本质上是一种带惯性权重的粒子群算法。

步骤5.1、建立多目标粒子群优化算法的表达式：

式中：d表示粒子的维数；i表示第i个粒子；和分别代表粒子在k时刻的运动速度向量和位置向量；和分别代表粒子在k+1时刻的运动速度向量和位置向量；为粒子在k时刻的自身最优位置；为群体在k时刻的全局最优位置；w∈(0,1)为惯性权重，代表着粒子先前速度项对粒子速度的影响；c₁和c₂为加速因子，c₁代表着粒子先前的本身最优位置对当前的影响，c₂是代表着群体先前的最优位置对当前的影响；r₁和 r₂为(0,1)之间的独立随机数。

步骤5.2、对充液挠性航天器的滚动轴、俯仰轴、偏航轴的姿态控制器和机动路径进行联合优化，需要优化的参数有12个，分别为K_px、K_py、K_pz、K_dx、K_dy、K_dz、a_{max_x}、 a_{max_y}、a_{max_z}、T_x、T_y、T_z。要求充液挠性航天器在规定的最大姿态机动时间t_max完成机动，在稳定时间段内达到指定的指向精度p_t，优化的性能指标是滚动轴、俯仰轴、偏航轴三轴角位置分别达到指定的指向精度p_t的时间f_1x、f_1y、f_1z以及在三轴分别在稳定时间段内的稳定度f_2x、f_2y、f_2z。优化的目标是f_1x、f_1y、f_1z、f_2x、f_2y、f_2z同时达到较小的值。

即，求取：ζ＝[a_{max_x},a_{max_y},a_{max_z},T_x,T_y,T_z,K_px,K_py,K_pz,K_dx,K_dy,K_dz]

s.t.F(ζ)＝min[f_1x,f_1y,f_1z,f_2x,f_2y,f_2z]

其中，ζ满足约束条件：

①K_px∈[K_{px_min},K_{px_max}]、K_py∈[K_{py_min},K_{py_max}]、K_pz∈[K_{pz_min},K_{pz_max}]，其中K_{px_min}、K_{px_max}分别为设置的滚动轴比例系数的最小值和最大值，K_{py_min}、K_{py_max}分别为为设置的俯仰轴比例系数的最小值和最大值，K_{pz_min}、K_{pz_max}分别为设置的偏航轴比例系数的最小值和最大值；

②K_dx∈[K_{dx_min},K_{dx_max}]、K_dy∈[K_{dy_min},K_{dy_max}]、K_dz∈[K_{dz_min},K_{dz_max}]，其中K_{dx_min}、K_{dx_max}分别为设置的滚动轴微分系数的最小值和最大值，K_{dy_min}、K_{dy_max}分别为设置的俯仰轴微分系数的最小值和最大值，K_{dz_min}、K_{dz_max}分别为设置的偏航轴微分系数的最小值和最大值；

③a_{max_x}∈(0,a_max]、a_{max_y}∈(0,a_max]、a_{max_z}∈(0,a_max]，其中a_max为设置的路径规划角加速度的最大值；

④T_x∈(0,t_max]、T_y∈(0,t_max]、T_z∈(0,t_max]，其中t_max为设置的最大姿态机动时间；角位置曲线中各段曲线的时间和不大于t_max。

⑤航天器角位置在稳定时达到指定机动角度

步骤5.3、结合图4，基于自适应网格的多目标粒子群优化算法，对充液挠性航天器的控制器和机动路径参数进行联合优化：

5.3.1、初始化：

①初始化粒子群的初始位置和初始速度：粒子的初始速度均为零，初始位置为一个 12维的向量：K_px、K_py、K_pz、K_dx、K_dy、K_dz、a_{max_x}、a_{max_y}、a_{max_z}、T_x、T_y、T_z；

②初始化公式(1)中各参数：惯性权重w、加速因子c₁和c₂；

③设置相关参数：算法的最大迭代次数T_imax；外部档案NP的大小NP_size；自适应网格参数cubelen；粒子群规模popsize，同时令初始外部档案NP为空集，设置循环次数t_n＝1。

5.3.2、计算粒子的适应度值，即优化的6个性能指标，将产生的非劣解加入到外部档案NP中，若数量超出NP_size，则随机选出NP_size个粒子组成外部档案NP；同时初始化粒子的局部最优pbest。

5.3.3、将自适应网格空间划分为多个格子，并同时确定各粒子所在的格子；为所有包含至少一个NP成员的格子计算适应度值：使用轮盘赌算法选择一个格子，并从格子中随机选择一个粒子作为局部最优gbest。

5.3.4、利用公式(1)对粒子的速度和位置进行更新。

5.3.5、保证粒子不飞出搜索空间，当粒子将要飞出时，将粒子保持在搜索边界，同时反转粒子的飞行方向。

5.3.6、进行局部最优和全局最优的更新：

局部最优pbest的更新：当新粒子的6个性能指标均小于pbest时，则新粒子为新的pbest；如果pbest的6个性能指标均小于新粒子，则pbest不改变；否则从两个粒子中随机选择一个作为pbest。

全局最优gbest的更新：将目标空间分成很多格子，为每个至少包含一个外部粒子群的格子定义其适应度值，然后采用轮盘赌算法选择格子以及粒子的gbest。

5.3.7、如果t_n＝T_imax，则输出结果12个优化参数，结束算法；否则令t_n＝t_n+1，跳转至5.3.2，继续执行。

实施例

步骤1、建立充液挠性航天器的动力学模型：一类带有挠性太阳能帆板及液体燃料的三轴航天器的姿态动力学方程、挠性附件振动方程和液体燃料晃动方程分别可描述为：

其中：I∈R^3×3为航天器转动惯量；ω∈R^3×1为航天器角速度；ω^×∈R^3×3表示ω的斜对称矩阵；C₀∈R^3×n为挠性振动运动与航天器转动运动耦合系数矩阵； G₁∈R^2×m，C₁∈R^m×3，C₂∈R^m×3均为液体晃动运动与航天转动之间的耦合系数矩阵；τ∈R^3×1为控制力矩矢量；T_d∈R^3×1为环境干扰力矩；ζ∈R^n×n为挠性附件振动模态阻尼比；Λ∈R^n×n为挠性附件振动模态频率矩阵；Λ₂∈R^m×m为液体晃动模态频率平方矩阵；η∈R^n×1为挠性体振动模态矢量；q∈R^m×1为液体晃动模态矢量；α为航天器的姿态角矢量，采用李英波.挠性充液卫星动力学分析与姿态控制研究[D].上海交通大学, 2001.中的航天器，相关参数如下：

航天器转动惯量：

挠性振动运动与星体转动运动耦合系数矩阵：

挠性附件振动模态频率矩阵：

Λ＝diag{1.02464,1.23670,1.91610,2.85637,3.87904}

挠性附件振动模态阻尼比：

ζ＝diag{0.001,0.001,0.001,0.001,0.001}

液体晃动模态频率平方矩阵：

Λ₂＝10^-4diag{1.0589,1.0589}

液体晃动运动与星体转动之间的耦合系数矩阵：

环境干扰力矩仅考虑太阳辐射压力矩，取：

T_d＝10^-4[-0.185 0.731 0.0555]^T

步骤2、获取充液挠性航天器的角加速度曲线、角速度曲线、得到多段的角位置曲线，对充液挠性航天器姿态机动路径进行规划；

结合图2，采用基于正弦型角加速度的分段路径，将角加速度微分曲线按波形分成7段，路径的角加速度微分曲线(图2-I)是由正弦函数组成的，I-①和I-⑦是正弦函数的正半周期，I-③和I-⑤是正弦函数的负半周期，I-②、I-④和I-⑥为0；对这七段角加速度微分曲线(图2-I)积分，即得整个机动路径的角加速度曲线(图2-Ⅱ)，对角加速度曲线 (图2-Ⅱ)积分，即得整个机动路径的角速度曲线(图2-Ⅲ)，再对角速度曲线(图2-Ⅲ)积分，即得整个机动路径的角位置曲线(图2-Ⅳ)；

步骤3、计算角位置多段曲线中每段曲线的表达式：确定j轴(j＝x,y,z，分别表示航天器的滚动轴、俯仰轴、偏航轴)机动到指定角度定义滚动轴从-30°机动到30°，俯仰轴和偏航轴从0.5°机动到0°；定义j轴航天器角速度最大值V_{max_j}＝2.5°/s；

采用基于正弦型角加速度的分段路径，将角加速度微分曲线按波形分成7段，确定j轴角速度最大值a_{max_j}和图2-I中正弦函数的周期T_j(即角加速度微分曲线中非零段时间和的一半)，图中Ⅳ-②和Ⅳ-⑥时间相同，均为t_j1；Ⅳ-④时间为t_j2，可以根据a_{max_j}、 T_j、V_{max_j}算出再根据a_{max_j}、T_j、V_{max_j}、t_j1算出从而算出角位置(图2-Ⅳ)的表达式为：

步骤4、结合图3，采用PD控制系统对充液挠性航天器进行姿态控制。PD控制律为：

其中，q_d为规划的欧拉角角位置输出，为规划的欧拉角角速度输出，为规划的欧拉角角加速度输出；q_p为实际的欧拉角角位置输出，为实际的欧拉角角速度输出； e＝q_d-q_p为欧拉角角位置误差向量；为欧拉角角速度误差向量。K_p为比例系数矩阵，K_d为微分系数矩阵，K_p＝diag(K_px,K_py,K_pz)，K_d＝diag(K_dx,K_dy,K_dz)，其中 K_px、K_py、K_pz分别表示航天器滚动轴、俯仰轴、偏航轴的控制器的比例系数，K_dx、K_dy、 K_dz分别表示航天器滚动轴、俯仰轴、偏航轴的控制器的微分系数。

然后对输出进行限幅，控制力矩向量为u_limit为实际系统所允许的最大控制力矩输出向量，设置u_limit为25Nm；

步骤5、联合优化充液挠性航天器控制器和机动路径的参数：采用基于自适应网格的多目标粒子群优化算法对充液挠性卫星的控制器和机动路径参数进行联合优化。基于自适应网格的多目标粒子群优化算法表达式如下：

步骤5.1、建立多目标粒子群优化算法的表达式：

式中：d表示粒子的维数；i表示第i个粒子；和分别代表粒子在k时刻的运动速度向量和位置向量；和分别代表粒子在k+1时刻的运动速度向量和位置向量；为粒子在k时刻的自身最优位置；为群体在k时刻的全局最优位置；w∈(0,1)为惯性权重，代表着粒子先前速度项对粒子速度的影响；c₁和c₂为加速因子，c₁代表着粒子先前的本身最优位置对当前的影响，c₂是代表着群体先前的最优位置对当前的影响；r₁和 r₂为(0,1)之间的独立随机数。

步骤5.2、对充液挠性航天器的滚动轴、俯仰轴、偏航轴的姿态控制器和机动路径进行联合优化，需要优化的参数有12个，分别为K_px、K_py、K_pz、K_dx、K_dy、K_dz、a_{max_x}、 a_{max_y}、a_{max_z}、T_x、T_y、T_z。本实施例要求充液挠性航天器在规定的最大姿态机动时间t_max＝70s完成机动，在稳定时间段80s～130s内达到指定的指向精度p_t＝0.005°，优化的性能指标是三轴角位置分别达到指定的指向精度p_t＝0.005°的时间f_1x、f_1y、f_1z以及三轴在稳定时间段80s～130s内的稳定度f_2x、f_2y、f_2z。优化的目标是f_1x、f_1y、f_1z、f_2x、 f_2y、f_2z同时达到较小的值。

即，求取：ζ＝[a_{max_x},a_{max_y},a_{max_z},T_x,T_y,T_z,K_px,K_py,K_pz,K_dx,K_dy,K_dz]

s.t.F(ζ)＝min[f_1x,f_1y,f_1z,f_2x,f_2y,f_2z]

其中，ζ满足约束条件：

①K_px∈[K_{px_min},K_{px_max}]、K_py∈[K_{py_min},K_{py_max}]、K_pz∈[K_{pz_min},K_{pz_max}]，其中K_{px_min}、K_{px_max}分别为设置的滚动轴比例系数的最小值和最大值，本实施例采用的K_{px_min}＝1、 K_{px_max}＝20，K_{py_min}、K_{py_max}分别为设置的俯仰轴比例系数的最小值和最大值，本实施例采用的K_{py_min}＝1、K_{py_max}＝20，K_{pz_min}、K_{pz_max}分别为设置的偏航轴比例系数的最小值和最大值，本实施例采用的K_{pz_min}＝1、K_{pz_max}＝20；

②K_dx∈[K_{dx_min},K_{dx_max}]、K_dy∈[K_{dy_min},K_{dy_max}]、K_dz∈[K_{dz_min},K_{dz_max}]，其中K_{dx_min}、K_{dx_max}分别为设置的滚动轴比例系数的最小值和最大值，本实施例采用的K_{dx_min}＝1、 K_{dx_max}＝20，K_{dy_min}、K_{dy_max}分别为设置的俯仰轴比例系数的最小值和最大值，本实施例采用的K_{dy_min}＝1、K_{dy_max}＝20，K_{dz_min}、K_{dz_max}分别为设置的偏航轴比例系数的最小值和最大值，本实施例采用的K_{dz_min}＝1、K_{dz_max}＝20；

③a_{max_x}∈(0,a_max]、a_{max_y}∈(0,a_max]、a_{max_z}∈(0,a_max]，其中a_max为路径规划角加速度的最大值，本实施例采用的a_max＝0.4°/s²；

④T_x∈(0,t_max]、T_y∈(0,t_max]、T_z∈(0,t_max]，其中t_max为规定的最大姿态机动时间，设置t_max＝70s；2T_x+t_x1+t_x2∈(0,t_max]、2T_y+t_y1+t_y2∈(0,t_max]、2T_z+t_z1+t_z2∈(0,t_max]；

⑤航天器角位置在稳定时达到指定机动角度即滚动轴达到30°，俯仰轴和偏航轴达到0°。

步骤5.3、结合图4，基于自适应网格的多目标粒子群优化算法，对充液挠性卫星的控制器和机动路径参数进行联合优化：

5.3.1、初始化：

①初始化粒子群的初始位置和初始速度：粒子的初始速度均为零，初始位置为一个 12维的向量：K_px、K_py、K_pz、K_dx、K_dy、K_dz、a_{max_x}、a_{max_y}、a_{max_z}、T_x、T_y、T_z；

②初始化速度位置更新模型各参数：惯性权重ω＝0.6、加速因子c₁和c₂为 c₁＝c₂＝0.4；

③设置相关参数：算法的最大迭代次数T_imax＝10；外部档案NP的大小NP_size＝20；自适应网格参数cubelen＝30；粒子群规模popsize＝20，同时令初始外部档案NP为空集，设置循环次数t_n＝1。

5.3.2、计算粒子的适应度值，即优化的6个性能指标，将产生的非劣解加入到外部档案NP中，若数量超出NP_size，则随机选出NP_size个粒子组成外部档案NP；同时初始化粒子的局部最优pbest。

5.3.3、将自适应网格空间划分为多个格子，并同时确定各粒子所在的格子；为所有包含至少一个NP成员的格子计算适应度值：使用轮盘赌算法选择一个格子，并从格子中随机选择一个粒子作为gbest。

5.3.4、利用公式(1)对粒子的速度和位置进行更新。

5.3.5、保证粒子不飞出搜索空间，当粒子将要飞出时，将粒子保持在搜索边界，同时反转粒子的飞行方向。

5.3.6、进行局部最优和全局最优的更新：局部最优pbest的更新：当新解的6个性能指标均小于pbest时，则新解为新的pbest；如果pbest的6个性能指标均小于新解，则pbest 不改变；否则从两个解中随机选择一个作为pbest。全局最优gbest的更新：将目标空间分成很多格子，为每个至少包含一个外部粒子群的格子定义其适应度值，然后采用轮盘赌算法选择格子以及粒子的gbest。

5.3.7、如果t_n＝T_imax，则输出结果12个优化参数，结束算法；否则令t_n＝t_n+1，跳转至Step2，继续执行。

用Matlab进行仿真，执行优化算法，得到两组满足航天器性能指标的控制器参数和路径参数的Parato解，每组解包含12个优化的参数，将其应用于充液挠性航天器的姿态控制系统中。表1给出了基于自适应网格的多目标粒子群优化算法的两组Pareto解，表2 给出了参数优化结果对应的性能指标：

表1基于自适应网格的多目标粒子群优化算法联合优化参数优化结果

表2参数优化结果对应的性能指标

由以上两表可以看出航天器都能在49.25s以内达到稳定状态，并且稳定度都能达到 10^-6°/s数量级，姿态机动路径规划及参数优化效果较好，两组解各有优点。

将Pareto解1作为路径参数的充液挠性航天器的仿真结果如图5～7所示，分析仿真结果：由图5(a-c)可以看出，航天器三轴均能快速机动到指定角度；由图6(a-c)可以看出，航天器三轴角位置误差均能在47.05s内进入0.005°误差带，由图7(a-c)和表 2可以看出，航天器三轴稳定度均能够达到10^-6°/s数量级。

本实施例采用基于自适应网格的多目标粒子群优化算法，输出12个优化参数，即航天器的控制器参数和路径参数，使得航天器在47.05s以内达到稳定状态，并且稳定度都能达到10^-6°/s数量级，减小了液体晃动和挠性附件振动对航天器控制性能的影响，使航天器系统快速机动快速稳定，验证了本发明的有效性。

Claims (9)

1.一种刚柔液耦合系统姿态控制器和机动路径联合优化方法，其特征在于，包括以下步骤：

步骤1、建立充液挠性航天器的动力学模型：一类带有挠性太阳能帆板及液体燃料的三轴航天器的姿态动力学方程、挠性附件振动方程和液体燃料晃动方程；

步骤2、获取充液挠性航天器的角加速度曲线、角速度曲线、得到多段的角位置曲线，对充液挠性航天器姿态机动路径进行规划；

步骤3、计算角位置多段曲线中每段曲线的表达式；利用角加速度最大值、角加速度微分曲线中非零段时间和的一半或角加速度曲线中非零段时间和的一半、航天器角速度最大值、期望机动的角度，解算出航天器每段曲线中角位置的表达式；

步骤4、采用PD控制系统对充液挠性航天器进行姿态控制；

步骤5、联合优化充液挠性航天器控制器参数和机动路径的参数：采用基于自适应网格的多目标粒子群优化算法对航天器的控制器参数和机动路径参数进行联合优化。

2.根据权利要求1所述的刚柔液耦合系统姿态控制器参数和机动路径参数联合优化方法，其特征在于，步骤1建立的姿态动力学方程、挠性附件振动方程和液体燃料晃动方程，分别为：

其中：I为航天器转动惯量；ω∈R^3×1为航天器角速度；ω^×∈R^3×3表示ω的斜对称矩阵；C₀∈R^3×n为挠性振动运动与航天器转动运动耦合系数矩阵；G₁∈R^2×m，C₁∈R^m×3，C₂∈R^m×3均为液体晃动运动与航天转动之间的耦合系数矩阵；τ∈R^3×1为控制力矩矢量；T_d∈R^3×1为环境干扰力矩；ζ∈R^n×n为挠性附件振动模态阻尼比；Λ∈R^n×n为挠性附件振动模态频率矩阵；Λ₂∈R^m×m为液体晃动模态频率平方矩阵；η∈R^n×1为挠性体振动模态矢量；q∈R^m×1为液体晃动模态矢量；α为航天器的姿态角矢量，

3.根据权利要求2所述的刚柔液耦合系统姿态控制器和机动路径联合优化方法，其特征在于，步骤2对充液挠性航天器姿态机动路径进行规划，采用将角加速度微分曲线划分为多段并对每段曲线进行积分，即得整个机动路径的角加速度多段曲线，再对角加速度曲线积分，即得整个机动路径的角速度多段曲线；再对角速度曲线积分，即得整个机动路径的角位置多段曲线，实现对充液挠性航天器姿态机动路径的规划。

4.根据权利要求3所述的刚柔液耦合系统姿态控制器和机动路径联合优化方法，其特征在于，步骤3计算角位置多段曲线中每段曲线的表达式，具体包括以下步骤：

步骤3.1、定义j轴期望机动的角度其中，j＝x,y,z分别表示航天器的滚动轴、俯仰轴、偏航轴；

步骤3.2、定义j轴航天器角速度最大值V_{max_j}；

步骤3.3、获取j轴两个机动路径参数：角加速度最大值a_{max_j}、角加速度微分曲线中非零段时间和的一半T_j或角加速度曲线中非零段时间和的一半T_j；

步骤3.4、解算出航天器每段曲线中角位置的表达式。

5.根据权利要求4所述的刚柔液耦合系统姿态控制器和机动路径联合优化方法，其特征在于，PD控制律为：

其中，q_d为规划的欧拉角角位置输出，为规划的欧拉角角速度输出，为规划的欧拉角角加速度输出；q_p为实际的欧拉角角位置输出，为实际的欧拉角角速度输出；e＝q_d-q_p为欧拉角角位置误差向量；为欧拉角角速度误差向量；K_p为比例系数矩阵，K_d为微分系数矩阵；

然后对输出进行限幅，控制力矩向量为u_limit为实际系统所允许的最大控制力矩输出向量。

6.根据权利要求5所述的刚柔液耦合系统姿态控制器和机动路径联合优化方法，其特征在于，步骤5联合优化充液挠性航天器控制器和机动路径的参数，具体包括以下步骤：

步骤5.1、建立多目标粒子群优化算法的表达式：

式中：d表示粒子的维数；i表示第i个粒子；和分别代表粒子在k时刻的运动速度向量和位置向量；和分别代表粒子在k+1时刻的运动速度向量和位置向量；为粒子在k时刻的自身最优位置；为群体在k时刻的全局最优位置；w∈(0,1)为惯性权重，代表着粒子先前速度项对粒子速度的影响；c₁和c₂为加速因子；

步骤5.2、对充液挠性航天器的滚动轴、俯仰轴、偏航轴的姿态控制器和机动路径进行联合优化，需要优化的参数有12个，分别为K_px、K_py、K_pz、K_dx、K_dy、K_dz、a_{max_x}、a_{max_y}、a_{max_z}、T_x、T_y、T_z。要求充液挠性航天器在规定的最大姿态机动时间t_max完成机动，在稳定时间段内达到指定的指向精度p_t，优化的性能指标是滚动轴、俯仰轴、偏航轴三轴角位置分别达到指定的指向精度p_t的时间f_1x、f_1y、f_1z以及在三轴分别在稳定时间段内的稳定度f_2x、f_2y、f_2z。优化的目标是f_1x、f_1y、f_1z、f_2x、f_2y、f_2z同时达到较小的值。

即，求取：ζ＝[a_{max_x},a_{max_y},a_{max_z},T_x,T_y,T_z,K_px,K_py,K_pz,K_dx,K_dy,K_dz]

s.t.F(ζ)＝min[f_1x,f_1y,f_1z,f_2x,f_2y,f_2z]

其中，ζ满足约束条件：

①K_px∈[K_{px_min},K_{px_max}]、K_py∈[K_{py_min},K_{py_max}]、K_pz∈[K_{pz_min},K_{pz_max}]，其中K_{px_min}、K_{px_max}分别为设置的滚动轴比例系数的最小值和最大值，K_{py_min}、K_{py_max}分别为为设置的俯仰轴比例系数的最小值和最大值，K_{pz_min}、K_{pz_max}分别为设置的偏航轴比例系数的最小值和最大值；

②K_dx∈[K_{dx_min},K_{dx_max}]、K_dy∈[K_{dy_min},K_{dy_max}]、K_dz∈[K_{dz_min},K_{dz_max}]，其中K_{dx_min}、K_{dx_max}分别为设置的滚动轴微分系数的最小值和最大值，K_{dy_min}、K_{dy_max}分别为设置的俯仰轴微分系数的最小值和最大值，K_{dz_min}、K_{dz_max}分别为设置的偏航轴微分系数的最小值和最大值；

③a_{max_x}∈(0,a_max]、a_{max_y}∈(0,a_max]、a_{max_z}∈(0,a_max]，其中a_max为设置的路径规划角加速度的最大值；

④T_x∈(0,t_max]、T_y∈(0,t_max]、T_z∈(0,t_max]，其中t_max为设置的最大姿态机动时间；角位置曲线中各段曲线的时间和不大于t_max；

⑤航天器角位置在稳定时达到指定机动角度

步骤5.3、基于自适应网格的多目标粒子群优化算法，对充液挠性航天器的控制器和机动路径参数进行联合优化。

7.根据权利要求6所述的刚柔液耦合系统姿态控制器和机动路径联合优化方法，其特征在于，步骤5.3对充液挠性航天器的控制器和机动路径参数进行联合优化，具体包括以下步骤：

5.3.1、初始化：初始化粒子群的初始位置和初始速度、惯性权重、加速因子、算法的最大迭代次数T_imax；外部档案NP的大小NP_size；自适应网格参数cubelen；粒子群规模popsize，同时令初始外部档案NP为空集，设置循环次数t_n＝1；

5.3.2、计算粒子的适应度值；

5.3.3、选择一个粒子作为局部最优gbest；

5.3.4、利用多目标粒子群优化算法对粒子的速度和位置进行更新；

5.3.5、将粒子保持在搜索边界；

5.3.6、进行局部最优和全局最优的更新；

5.3.7、如果t_n＝T_imax，则输出结果12个优化参数，结束算法；否则令t_n＝t_n+1，跳转至5.3.2，继续执行。

8.根据权利要求4所述的刚柔液耦合系统姿态控制器和机动路径联合优化方法，其特征在于，采用基于正弦型角加速度的分段路径。

9.根据权利要求8所述的刚柔液耦合系统姿态控制器和机动路径联合优化方法，其特征在于，将角加速度微分曲线按波形分成7段，角位置的表达式为：

Ⅳ-①：

Ⅳ-②：

Ⅳ-③：

Ⅳ-④：

Ⅳ-⑤：

Ⅳ-⑥：

Ⅳ-⑦：

Ⅳ-⑧：

其中，