CN110083170A - 一种使用固体微推力器进行轨道保持的优化控制方法 - Google Patents

Abstract

本发明公开了一种使用固体微推力器进行轨道保持的优化控制方法，属于最优化控制领域。首先初始化粒子群算法种群和设置参数，然后针对某微纳卫星上布置的固体微推力器阵列，根据固体微推力器不可重复使用的控制机制，以某个周期内所有固体微推力器的最小能耗为目标建立适应度函数。最后以最小化适应度函数为粒子群算法的优化目标，通过寻优模型得到最优加权阵Q和R，代入线性二次型控制算法，求解得到最优的轨道保持控制器，带入每一时刻的卫星轨道状态X(t)得到对应的最优反馈控制力U(t)，达到轨道保持的目的。本发明大大减少了Q、R参数选取的盲目性和不确定性，能自动给出与给定初始卫星轨道条件相匹配的最优Q、R参数值，效率更高，控制性能更优。

Description

一种使用固体微推力器进行轨道保持的优化控制方法

技术领域

本发明属于最优化控制领域，涉及航天器轨道保持控制，具体是一种使用固体微推力器进行轨道保持的优化控制方法。

背景技术

卫星推进系统的作用是为卫星轨道转移、位置保持提供所需要的推力，为姿态控制提供所需的力矩，是航天器最重要的分系统之一。不同的能量来源和转换方式，决定了不同的推进方式，主要有化学推进、冷气推进、电推进、核推进、动量转换推进和无工质推进等。推进系统的推进剂携带量和有效比冲制约着航天器的寿命，在深空探测中直接决定任务的成败。

微纳卫星具有一定的优点：生产和发射成本低，研发周期短，大量微纳卫星可以编队飞行，增加任务的可靠性和可行性，用来完成空间和地球科学任务，以及对其他卫星的检测和探测等。

但是在微纳卫星研制工作中推进系统成为制约其发展的一个重要因素，传统推进系统体积和质量比较大，不适合应用在微纳卫星上。现有微纳型卫星基本不配备推进系统，或者推进系统只提供极其有限的机动能力，控制能力和机动能力较差，无法完成复杂的航天任务。

目前航天推进领域已经运用MEMS技术发展了双组元、单组元、冷气推力器以及水蒸气推力器等，但是这些推力器存在技术难点，结构复杂。而近些年新出现的MEMS固体推力器阵列具有结构简单、质量小和高度高的优点，用其进行卫星轨道保持控制可以使微纳卫星具备较强的机动能力，从而完成更加复杂的航天任务；见文献[1]：杨灵芝，魏延明，刘旭辉在2016年2月发表的期刊论文《MEMS固体微推力器阵列发展研究》。

微纳卫星的轨道控制不同于现有的大型和中型卫星，其对控制精度和燃料消耗的优化都有很高的要求。因此，需要针对MEMS固体推力器阵列这种轨道控制执行机构设计专用的控制方法。传统的卫星轨道控制一般是对卫星的轨道六要素进行修正，精度较差，能量消耗多，不适应于MEMS固体推力器阵列这种微小推力执行机构。

为了达到微纳卫星的燃料消耗和控制精度最优，线性二次型最优控制算法(LQR)是一种可行的方法。目前的LQR控制理论中，加权矩阵Q和R的选择与控制性能之间尚未形成明确的理论关系，目前大都是通过试凑计算的方法来确定Q和R阵，非常依赖于设计人员的工作经验，当更换初始条件时往往需要重新寻找或试凑新的Q和R阵，操作起来十分繁琐和复杂。

针对LQR算法中Q和R阵的优化选取方法，目前已有的方法使用粒子群等智能搜索技术解决的都是微纳卫星轨道控制领域以外的问题，见文献[2]：张玉分，卢家暄等在2017年10月发表的期刊论文《基于粒子群优化的主动悬架最优控制研究》更没有考虑到MEMS固体为推力器这种特殊的脉冲式执行机构的特点，因此均不能用于以固体为推力器为执行机构的微纳卫星轨道控制设计。

发明内容

本发明是通过当前微纳卫星轨道控制系统的状态和预设计好的控制性能指标，自动地搜索出符合设计条件的最优Q和R阵，可以有效减少设计控制器时的试凑次数，降低对设计人员经验的依赖；具体是一种使用固体微推力器进行轨道保持的优化控制方法。

分以下步骤：

步骤一、初始化粒子群算法种群和设置参数。

参数包括粒子数，最大速度，学习因子和惯性权重；

步骤二、针对某微纳卫星上布置的固体微推力器阵列，根据固体微推力器不可重复使用的控制机制，以某个周期内所有固体微推力器的最小能耗为目标建立适应度函数。

适应度函数为：

ge(i)为第i次施加控制时所消耗的3轴方向的推力器个数。n代表整个周期内的总控制次数。

步骤三、以最小化适应度函数为粒子群算法的优化目标，通过寻优模型得到最优加权阵Q和R。

寻优模型包括三种：单变量R参数寻优模型，两变量Q、R寻优模型和九变量Q、R寻优模型。

单变量R参数寻优模型指的是：在LQR控制参数设计时，固定Q参数为1，即Q矩阵取单位阵I_6×6，R矩阵的对角线元素相等，即

寻优变量只有1个，即R矩阵对角线元素的指数。

两变量Q、R寻优模型为：

Q＝x₁×I_6×6，

寻优变量为2维：

x＝[x₁x₂]

九变量Q、R寻优模型为：

寻优变量为9维：

x＝[x₁ x₂ x₃ x₄ x₅ x₆ x₇ x₈ x₉]

将最小化适应度函数f_fitness作为标准粒子群算法的优化目标，选择上述三种寻优模型中的任一种，搜索得到向量x，求得最优加权阵Q和R。

步骤四、将最优加权阵Q和R代入线性二次型控制算法，求解得到最优的轨道保持控制器，带入每一时刻的卫星轨道状态X(t)得到对应的最优反馈控制力U(t)，达到轨道保持的目的。

具体步骤为：

步骤401、针对某时刻t，计算该时刻下线性化的卫星轨道动力学系统中状态X(t)的系统矩阵A以及反馈控制力U(t)的系数矩阵B；

状态Δx，Δy和Δz是卫星和标准轨道之间的3轴位置误差；

系数矩阵

其中，μ是地球常数；r是该时刻卫星的地心矩；

系数矩阵

步骤402、根据系数矩阵A和系数矩阵B，以及最优加权阵Q和R，求解黎卡提方程，得到矩阵P(t)；

步骤403、利用矩阵P(t)的数值解，得出最优反馈增益控制矩阵；

U(t)＝-R-¹(t)B^T(t)P(t)X(t)

u_x，u_y和u_z是卫星在当前时刻受到的3轴控制力。

本发明的优点在于：

一种使用固体微推力器进行轨道保持的优化控制方法，利用粒子群参数辨识算法，大大减少了卫星轨道LQR最优控制问题中Q、R参数选取的盲目性和不确定性，能够自动给出与给定初始卫星轨道条件相匹配的最优Q、R参数值，相比人工选择效率更高，控制性能更优。

附图说明

图1为本发明一种使用固体微推力器进行轨道保持的优化控制方法的步骤流程图；

图2为本发明实施例中卫星轨道为h＝700km选用单变量寻优模型的误差图；

图3为本发明实施例中一个周期12小时内推力器消耗情况图；

图4为本发明实施例中卫星轨道为h＝700km选用九变量寻优模型的误差图；

图5为本发明实施例中单变量寻优模型和九变量寻优模型的误差对比图。

具体实施方式

下面将结合附图和实例对本发明作进一步的详细说明。

本发明一种使用固体微推力器进行轨道保持的优化控制方法，以控制高集成脉冲式固体微推力器的最小能耗为目标，建立适应度函数，并以最小化适应度函数为粒子群算法的优化目标，对整个解空间进行搜索寻优，当MEMS固体微推力器消耗个数最少时，选取最优Q、R参数代入线性二次型控制算法，求解得到最优的轨道保持控制器，带入每一时刻的卫星轨道状态得到对应的最优反馈控制力，达到轨道保持的目的。具有效率高，控制效果好，燃料消耗少优势。

如图1所示，具体步骤如下：

步骤一、初始化粒子群算法种群和设置参数。

参数包括粒子数，最大速度V_max，学习因子和惯性权重；

粒子数：一般取20～40，对于比较难的问题，粒子数可以取到100或200。

最大速度V_max决定了粒子在一个循环中最大的移动距离，通常小于粒子的范围宽度。较大的V_max可以保证粒子种群的全局搜索能力，较小的V_max则保证粒子种群的局部搜索能力加强。

学习因子包括为局部学习因子c₁和全局学习因子c₂；一般c₂＞c₁。

惯性权重：一个大的惯性权值有利于展开全局寻优，而一个小的惯性权值有利于局部寻优。当粒子的最大速度V_max很小时，使用接近于1的惯性权重。当V_max不是很小时，使用权重0.8较好，或者使用时变权重。如果在迭代过程中采用线性递减惯性权值，则粒子群算法在开始时具有良好的全局搜索性能，能够迅速定位到接近全局最优点的区域，而在后期具有良好的局部搜索性能，能够精确地得到全局最优解。经验表明，惯性权重采用从0.90线性递减到0.10的策略，会获得比较好的算法性能。

步骤二、针对某微纳卫星上布置的固体微推力器阵列，根据固体微推力器不可重复使用的控制机制，以某个周期内所有固体微推力器的最小能耗为目标建立适应度函数。

在寻优过程中，有一个重要问题是如何选取适应度函数。适应度函数用于计算每次迭代得到的解是否达到最优。在普通的粒子群算法中，每次迭代都是向着使适应度函数值变大的方向发展，最终优化的结果是找到使适应度函数值最大的解。需要注意的是，由于在工程应用中需要解决的最优问题一般都是求最小值，比如最小能耗、最短距离等等，所以在MATLAB软件中对适应度函数做了一点改进，使得最终优化的结果是适应度函数值达到最小，称这个改进后的适应度函数为目标函数。这样做的好处是可以把需要求某些最小值的工程问题抽象化后直接作为优化算法的目标函数，具有一定的实际意义。因此，需要注意对粒子群算法中的适应度函数和实际应用MATLAB中的目标函数加以区别。

适应度函数的选取遵循如下规则：

利用固体微推力器进行轨道控制时，由于推力器是一种一次性点火使用的控制机构，因此需要尽量节省使用，从而能使卫星在轨道上运行足够长的时间。所以，任务目标是在满足误差精度要求的情况下，尽量使推力器消耗最少，延长卫星寿命。因此，以最优能耗为优化目标。在仿真过程中发现，卫星的位置误差随着卫星轨道周期呈现出一定的周期性变化，推力器的消耗量也有同样的规律。所以任务目标是需要尽量减少一个轨道周期内的推力器消耗量，令一个轨道周期内的推力器消耗量作为目标函数：

ge(i)为第i次施加控制时所消耗的3轴方向的推力器个数。n代表整个周期内的总控制次数。

步骤三、以最小化适应度函数为粒子群算法的优化目标，通过寻优模型得到最优加权阵Q和R。

需要优化的Q矩阵为6×6维对角半正定矩阵，R矩阵为3×3维对角正定矩阵。由于固体微推力器是一种不可再生控制力，应该引入较大的控制量权值矩阵R。根据先验知识，一般的R矩阵对角线元素取值范围为10⁷～10¹¹。因此，在进行系统辨识寻优时，可以选取Q阵的对角线元素(简称Q参数)和R阵对角线元素作为寻优变量x。但是，由于R矩阵的对角线元素绝对值很大，若直接将其作为寻优变量，会导致寻优时的区间范围过大而影响寻优速度及精度。因此选取R阵对角线元素的指数(简称R参数)作为寻优变量。即：

x＝[x₁ x₂ x₃ x₄ x₅ x₆ x₇ x₈ x₉]

设计以单个卫星周期内推力器消耗最少为目标的自适应LQR控制，可以实现轨道高度为500～1000km内太阳同步轨道的LQR参数辨识，根据轨道初始信息，可以自动辨识出当前轨道一个周期内控制消耗最少时的Q、R全局最优解，大大提高了LQR控制理论的通用性，不再需要盲目的试凑Q和R矩阵，而且能够保证Q、R参数的可用性。

根据Q和R矩阵的复杂程度，建立了3种不同的寻优模型：单变量R参数寻优模型，两变量Q、R寻优模型和九变量Q、R寻优模型。

单变量R参数寻优模型指的是：在LQR控制参数设计时，固定Q参数为1，即Q矩阵取单位阵I_6×6，R矩阵的对角线元素相等，即

寻优变量只有1个，即R矩阵对角线元素的指数。这样做的目的是为了简化寻优模型，减少运算量。

在单变量R参数寻优模型中，固定了Q参数的取值为1，这样做的好处是寻优速度较快，但是在一定程度上也具有局限性，无法确定二次型积分指标式中Q矩阵对固体微推力器消耗量是否有影响。因此研究Q、R参数两变量寻优。

两变量Q、R寻优模型为：

Q＝x₁×I_6×6，

寻优变量为2维：

x＝[x₁ x₂]

从理论上分析来看，Q参数的取值并不会影响最终结果，假设Q参数取a(a≠1)。即有

代入二次型指标函数中，得

即当Q参数不取1时，只需在二次型指标函数中提取出来一个常数，内部又转化成Q为单位阵的二次型指标函数，R矩阵对角线元素对应的缩小a倍。因此，经过验证结果也表明，若Q矩阵和R矩阵均对应放大或缩小相同倍数，最终计算出来反馈矩阵K是完全相同的。

上述单变量和两变量寻优都是对原寻优模型的一种简化，虽然大大缩减了寻优时间，但是寻优结果与原模型相比仍存在一定的差距，为了更进一步贴合实际，增加模型的准确性，采用九变量同时寻优。

九变量Q、R寻优模型为：

寻优变量为9维：

x＝[x₁ x₂ x₃ x₄ x₅ x₆ x₇ x₈ x₉]

即不再假设加权矩阵的对角线元素相同。

将最小化适应度函数f_fitness作为标准粒子群算法的优化目标，选择上述三种寻优模型中的任一种，搜索得到向量x，求得最优解加权阵Q和R。

步骤四、将最优加权阵Q和R代入线性二次型控制算法，求解得到最优的轨道保持控制器，带入每一时刻的卫星轨道状态X(t)得到对应的最优反馈控制力U(t)，达到轨道保持的目的。

具体步骤为：

步骤401、针对某时刻t，计算该时刻下线性化的卫星轨道动力学系统中状态X(t)的系统矩阵A以及反馈控制力U(t)的系数矩阵B；

首先，已知的惯性系下的线性轨道动力学状态方程为：

状态Δx，Δy和Δz是卫星和标准轨道之间的3轴位置误差；

系数矩阵

其中，μ是地球常数；r是该时刻卫星的地心矩；

系数矩阵u_x，u_y和u_z是卫星在当前时刻受到的3轴控制力。

步骤402、根据系数矩阵A和系数矩阵B，以及最优加权阵Q和R，求解黎卡提方程，得到矩阵P(t)；

步骤403、利用矩阵P(t)的数值解，得出最优反馈增益控制矩阵；

U(t)＝-R-¹(t)B^T(t)P(t)X(t)

变更最优反馈增益控制矩阵为：

K(t)＝R-¹(t)B^T(t)P(t)

则有

U(t)＝-K(t)X(t)。

选取太阳同步卫星轨道，为了验证自适应LQR控制算法的通用性，选取轨道高度分别为500km、600km、700km、800km、900km、1000km的6条太阳同步轨道进行仿真验证，卫星控制周期为60秒，卫星轨道初始误差设为0.3m，卫星参数见表1：

表1

卫星质量	1kg
		卫星尺寸	0.2m×0.2m×0.2m
卫星面质比	0.02m<sup>2</sup>/kg
		单个推力器冲量	1*10-<sup>4</sup>N·s
推力器点火延迟	0.0017s
		推力器点火持续时间	0.015s

6条太阳同步轨道的初始轨道根数如表2所示，考虑到太阳同步轨道的冻结性，偏心率和轨道倾角都是根据轨道高度计算得到的。

表2

采用单变量R参数寻优模型，固定Q参数为1，R参数寻优范围为8～11；以最优能耗即单个周期内推力器消耗最少为目标，辨识结果如表3所示：

表3

以700km轨道高度为例，用粒子群算法计算出的R参数进行12小时的长时间仿真，误差变化情况如图2所示，可以看出三轴误差基本能控制在02m以内，x方向和z方向误差随卫星周期呈周期性变化。其中x方向和z方向的误差大小相近，在相位上z方向的误差变化和x方向相差1/4个周期。y方向的误差模值较小。仿真过程中后期稳定时误差最大值为0.2144m，符合0.3m的精度要求范围。

一个周期12小时内的推力器消耗图如图3所示，可以看出，在12小时内推力器总共消耗了69个，平均每个周期消耗了9个，其中绝大部分的推力器都消耗在了卫星地心轨道系中的z轴方向。

采用两变量Q、R寻优模型进行仿真，每种轨道高度计算出来的对应最少推力器消耗的Q、R参数如表4所示。

表4

由表中可以看出，两变量寻优模型的得到的R参数与单变量寻优模型得到的R参数不尽相同，而且对同一轨道进行寻优时，有时得到的结果不是每次都相同，但是控制效果是相同的；通过进行12小时仿真，结果同粒子群算法单变量寻优模型的结果一致。

由于利用九变量寻优模型仿真时间太长，因此只针对轨道高度700km的太阳同步轨道进行了粒子群优化算法仿真验证。寻优结果为：

x＝[1.7147 2.9062 1 8.3282 1 5.9246 11 11 11]

以该组参数进行12小时仿真，误差如图4所示。单变量寻优模型和九变量寻优模型的误差对比图如图5所示，可以看出，九变量寻优模型在控制精度上不如单变量寻优模型。从误差模值对比可以看出，九变量寻优模型的结果在前两个小时内有误差超过0.3m的情况，但是之后就稳定在0.3m以下震荡，而单变量寻优模型的结果，误差模值在0.2m附近震荡。

仿真验证结果表明，基于粒子群算法的Q、R参数识别轨道保持算法能够有效解决非线性系统的参数辨识问题，大大减少了LQR最优控制问题中Q、R参数选取的盲目性和不确定性，能够自动给出与给定初始条件相匹配的Q、R参数值，相比人工选择更加有效，更具有通用性。

Claims (3)

1.一种使用固体微推力器进行轨道保持的优化控制方法，其特征在于，包括以下步骤：

步骤一、初始化粒子群算法种群和设置参数；

步骤二、针对某微纳卫星上布置的固体微推力器阵列，根据固体微推力器不可重复使用的控制机制，以某个周期内所有固体微推力器的最小能耗为目标建立适应度函数；

适应度函数为：

ge(i)为第i次施加控制时所消耗的3轴方向的推力器个数；n代表整个周期内的总控制次数；

步骤三、以最小化适应度函数为粒子群算法的优化目标，通过寻优模型得到最优加权阵Q和R；

步骤四、将最优加权阵Q和R代入线性二次型控制算法，求解得到最优的轨道保持控制器，带入每一时刻的卫星轨道状态X(t)得到对应的最优反馈控制力U(t)，达到轨道保持的目的；

具体步骤为：

步骤401、针对某时刻t，计算该时刻下线性化的卫星轨道动力学系统中状态X(t)的系统矩阵A以及反馈控制力U(t)的系数矩阵B；

状态Δx，Δy和Δz是卫星和标准轨道之间的3轴位置误差；

系数矩阵

其中，μ是地球常数；r是该时刻卫星的地心矩；

系数矩阵

步骤402、根据系数矩阵A和系数矩阵B，以及最优加权阵Q和R，求解黎卡提方程，得到矩阵P(t)；

步骤403、利用矩阵P(t)的数值解，得出最优反馈增益控制矩阵；

U(t)＝-R^-1(t)B^T(t)P(t)X(t)

u_x，u_y和u_z是卫星在当前时刻受到的3轴控制力。

2.如权利要求1所述的一种使用固体微推力器进行轨道保持的优化控制方法，其特征在于，步骤一中所述的参数包括粒子数，最大速度，学习因子和惯性权重。

3.如权利要求1所述的一种使用固体微推力器进行轨道保持的优化控制方法，其特征在于，步骤三中所述的寻优模型包括三种：单变量R参数寻优模型，两变量Q、R寻优模型和九变量Q、R寻优模型；

单变量R参数寻优模型指的是：在LQR控制参数设计时，固定Q参数为1，即Q矩阵取单位阵I_6×6，R矩阵的对角线元素相等，即

寻优变量只有1个，即R矩阵对角线元素的指数；

两变量Q、R寻优模型为：

寻优变量为2维：

x＝[x₁ x₂]

九变量Q、R寻优模型为：

寻优变量为9维：

x＝[x₁ x₂ x₃ x₄ x₅ x₆ x₇ x₈ x₉]

将最小化适应度函数f_fitness作为标准粒子群算法的优化目标，选择上述三种寻优模型中的任一种，搜索得到向量x，求得最优加权阵Q和R。