一种基于序列凸优化和模型预测控制的卫星集群重构控制方法

1.本发明属于卫星集群重构控制领域，具体涉及一种基于序列凸优化和模型预测控制的卫星集群重构控制方法。


背景技术：

2.卫星集群重构控制问题指的是卫星从一组轨道转移到另一组轨道的控制问题，其中卫星之间保持相对有界，即最远不能超过所规定的上界，最近不小于最小安全距离。卫星集群重构机动问题是一个复杂的多目标优化问题。通常是以推进剂消耗、轨道转移时间、卫星集群燃料消耗均衡性等作为目标，寻找最佳重构策略和机动策略，以实现特定任务。
3.卫星编队是卫星集群的一种特殊情况，编队所要求的不仅是有界性约束，并且包含严格的位置关系约束，因此，卫星集群在不追求最优性的情况下也可以使用编队飞行的一些控制方法。然而，以往的编队飞行工作通常涉及的航天器数量很少，最多的时候只有十几个，因此，将现有方法应用于大规模卫星集群轨道控制时算法的可扩展性较差，进行大规模卫星集群轨道控制时必须不同于以往的研究，需要同时考虑航天器数量多、单个微小卫星的能力有限以及复杂的动态环境干扰等问题。具体地说，大量的航天器使碰撞避免成为一大挑战。此外，每颗小卫星有限的计算能力要求制导律和控制算法非常简单，以便它们可以在卫星上实时运行。
4.大规模卫星集群的重构问题本质上也是一类求解具有非线性动态和非凸路径约束的多智能体系统最优控制问题。解决这类问题的方法有很多，包括间接法和直接法。由于航天器群的复杂的非线性动力学，间接方法需要推导出最优性的一阶必要条件，同时需要良好的初始猜测，因而使用起来非常困难。因此，许多最优控制问题都是用直接方法来解决的，这种方法将控制空间，有时甚至是状态空间参数化，将问题归结为一个非线性优化问题。伪谱法很适合处理非线性动力学，但随着卫星数量的增加，它们在计算上变得难以实现，而且还没有实时实现。到目前为止，非线性求解器的速度已经足够快，可以在机上实现适定问题，但不能很好地适合多智能体问题。混合整数线性规划(milp)可用于强制碰撞避免约束，并且已被实时实现以及用于预先规划轨迹。然而，随着航天器数量的增加，这些算法的可扩展性也很差，这是由于避碰约束数量的增加导致整数变量的增加。
5.在线求解最优控制问题依赖于当前不断进步的信息技术。最近，凸优化已被应用于多飞行器轨迹设计，凸规划问题可以在多项式时间内可靠地求解到全局最优。更重要的是，最近的进展表明，这些问题可以通过通用的二阶锥规划求解器和利用特定问题结构的定制求解器来实时解决。有研究用凸优化方法提出了一个凸问题的滚动时域控制器。此外，凸优化已被用于为队形重构和机器人运动规划寻找无碰撞轨迹。但是大量的研究文献的假设情况都是事先已经获得空间中的障碍物信息，且假设不存在执行机构误差、测量误差等不确定性误差，因此，控制方法的控制精度较低、鲁棒性较差。
6.在单一动态系统的各种控制方法中，模型预测控制(mpc)是一种应用广泛的控制
技术，由于其在处理多个控制目标和约束方面的灵活性和有效性而被广泛应用于多智能体系统。分布式模型预测控制是一种将模型预测控制(mpc)和分布式控制思想有机结合起来的控制方法，其具备显式处理约束、抑制干扰、滚动优化等mpc法的优势，且兼具分布式控制结构容错性高、结构灵活性强等特点。在dmpc设计过程中，可根据具体的任务需求，结合系统模型的特点，设计相应的目标函数、约束条件，通过滚动求解优化问题得到控制指令。与传统的控制器设计方法相比，基于mpc的方法具有更好的控制性能和更快、更平稳的响应。dmpc不需要数据处理中心，采用点对点通信，在大规模网络中具有可扩展性，对突发节点或链路故障具有鲁棒性。在输入约束、通信约束、通信延迟等不同条件下，基于dmpc的多智能体系统控制得到了广泛的研究。提出的dmpc框架中，每个本地控制器也考虑其相邻智能体的假设方案，并且所有控制器都按特定的顺序执行，而不是并行执行。有学者针对非线性连续时间多智能体系统，提出了一种自适应事件触发分布式模型预测控制算法。为了节省大量通信资源，提出了一种事件触发dmpc方法，保证了整个闭环系统的递归可行性和渐近稳定性。基于分布式模型预测控制(mpc)的联网车辆队列控制算法，重点研究了异常通信条件下的切换通信拓扑结构和控制策略。有学者提出了一种同步dmpc方案，利用相邻多智能体的估计信息来解决具有避碰性的跟踪和编队问题。然而，它不能将状态和输入约束强加到dmpc。
7.综上所述，现有方法仍然存在以下问题：
8.1)现有方法通常涉及的航天器数量很少，最多的时候只有十几个，方法的可扩展性差。
9.2)大规模卫星集群轨道控制会受到复杂的动态环境干扰，但是在实际任务中，每颗卫星的传感设备的传感范围和通信范围有限，卫星都只能得到周围一定范围的环境信息，无法预知全部信息，且没有考虑除了j2摄动和大气阻力的其他小扰动和执行机构误差，因此控制精度低。
10.3)现有方法在任务开始之前就已经直接给出了整个重构机动过程的控制序列，这就要求所有卫星对整个重构过程中的空间环境预先知道，方法的鲁棒性差。


技术实现要素：

11.本发明的目的是为解决现有控制方法的控制精度低、可扩展性差以及鲁棒性差的问题，而提出的一种基于序列凸优化和模型预测控制的卫星集群重构控制方法。
12.本发明为解决上述技术问题所采取的技术方案是：
13.一种基于序列凸优化和模型预测控制的卫星集群重构最优控制方法，所述方法具体包括以下步骤：
14.步骤s1、基于罚函数建立惩罚项代价函数；
15.步骤s2、建立凸子问题模型的各项约束；
16.步骤s3、基于步骤s1的惩罚项代价函数和步骤s2的各项约束建立基于模型预测的凸子问题；
17.步骤s4、利用分布式序列凸优化方法求解步骤s3的凸子问题，获得卫星集群重构控制序列。
18.本发明的有益效果是：
19.本发明方法可以应用于大规模卫星集群轨道控制，避免了现有方法的可扩展性差的问题，一旦时间更新到下一个预测时域，卫星将更新环境和邻居信息，考虑了空间环境的误差和卫星自身机构的误差，减少了任务执行过程中的误差，提高了控制精度，而且不需要所有卫星对整个重构过程中的空间环境预先知道，进而提高了鲁棒性。
20.同时，将分布式模型预测控制与序列凸优化相结合，只在预测时域内考虑避碰和避障约束，减少了约束数量，减少了凸子问题求解时间，提高了求解效率。
附图说明
21.图1是本发明的一种基于序列凸优化和模型预测控制的卫星集群重构最优控制方法的流程图；
22.图中，x
j,current
代表在每一次滚动优化过程中第j颗卫星当前的状态信息，x
j,f
代表第j颗卫星的终端状态信息，tf代表卫星重构机动总时间；
23.图2是凸子问题时域划分图。
具体实施方式
24.具体实施方式一、结合图1说明本实施方式。本实施方式所述的一种基于序列凸优化和模型预测控制的卫星集群重构控制方法，所述方法具体包括以下步骤：
25.步骤s1、基于罚函数建立惩罚项代价函数；
26.步骤s2、建立凸子问题模型的各项约束；
27.步骤s3、基于步骤s1的惩罚项代价函数和步骤s2的各项约束建立基于模型预测的凸子问题；
28.步骤s4、利用分布式序列凸优化方法求解步骤s3的凸子问题，获得卫星集群重构控制序列。
29.首先，由于每颗卫星对环境的局部认知，对其应用预测控制策略来实时更新其检测环境和局部有效轨迹。其次，基于当前预测范围内检测到的环境信息，卫星构建自己的解耦优化子问题，其中可能包含一些依赖于其他机器人未来状态的参数。为了独立构建自己的问题，机器人通过与邻居通信来预测邻居的运动，以检索他们的当前和/或历史状态，而不是等待邻居的预测信息，有限时域说明该卫星仅考虑在最近一段时间内与其他邻居卫星或障碍物的碰撞避免约束，极大的减少了凸子问题的大小。第三，利用分布式序列凸规划方法对子问题进行求解。由于凸子问题的建立和求解不依赖于其他机器人的信息，因此每颗卫星都可以分布式地求解子问题并执行其动作策略。第四，一旦时间更新到下一个预测时域，卫星将更新环境和邻居信息，考虑了空间环境的误差和卫星自身机构的误差，提高了鲁棒性。
30.具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤s1的具体过程为：
[0031][0032]
其中，j
penalty
为基于l1罚函数方法得到的惩罚项，ω
eq
为等式约束的罚系数，ω
ineq
为不等式约束的罚系数，为等式约束，n
eq
为等式约束的个数，为不等式约束，n
ineq
为不等式约束的个数，
[0033]
x代表所有卫星在离散点处的状态量和控制量的集合，将x表示为：
[0034]
x＝(x1,x2,...,xn)
t
,xj＝(xj[0],...,xj[k],uj[0],...,uj[k]),j＝1,...,n
ꢀꢀ
(2)
[0035]
其中，xj代表第j颗卫星在所有离散点处的状态量和控制量的集合，j＝1,...,n，n代表集群卫星的总个数，xj[k]代表第k步处第j颗卫星在lvlh坐标系三个坐标轴方向的位置速度信息，uj[k]代表第k步处第j颗卫星在lvlh坐标系三个坐标轴方向施加的推力加速度，k是离散点的总个数。
[0036]
将等式约束和不等式约束乘以一个非常大的罚系数，可以提高问题求解的可能性。
[0037]
其它步骤及参数与具体实施方式一相同。
[0038]
具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述步骤s2的具体过程为：
[0039]
步骤s21、建立动力学约束、推力幅值约束、速度幅值约束和避碰避障约束，并将建立的约束凸化；
[0040]
建立动力学约束并凸化：
[0041][0042]
其中，k0是在每一次滚动优化时的初始时刻离散点；
[0043]
上式为经过线性化和离散化处理以后得到的动力学约束，式中：
[0044][0045]
其中，uj(tk)代表第j颗卫星在第k步所对应的时间点tk处三个坐标轴方向施加的推力加速度，o(tk)为主星在tk处的轨道要素，为已知量，δt是离散步长，e是自然对数的底数，b＝[03×3,i3×3]
t
，xj[k]代表第j颗卫星在第k步时的基准轨迹，ξ是积分变量；aj[k]表示离散的状态转移矩阵，bj[k]表示离散的控制矩阵，sj[k]表示分段常值向量，和都是与基准轨迹有关的常量，得到基准轨迹后便可求得；
[0046]
假设主星轨道要素已知，且不对主星施加推力作用，则主星运动状态完全已知；
[0047]
凸化后的推力幅值约束：
[0048]
||uj[k]||
∞
≤u
max
,k＝k0,...,k-1,j＝1,...,n
ꢀꢀ
(5)
[0049]
式中，||
·
||
∞
为无穷范数，u
max
为推力器幅值上限；
[0050]
凸化后的速度幅值约束：
[0051]
||hxj[k]||2≤v
max
,k＝k0,...,kh,j＝1,...,n
ꢀꢀ
(6)其中，h＝[03×3,i3×3]，v
max
为速度最大值，||
·
||2为2范数，kh代表预测步长；
[0052]
凸化后的离散点处卫星间避碰约束：
[0053][0054]
式中，上角标t代表转置，代表第j颗卫星在第k步时的基准轨迹，代表第i颗卫星在第k步时的基准轨迹，g＝[i3×3,03×3]
t
，i3×3为单位矩阵，r
col
为卫星之间需满足的最小碰撞约束距离，代表前一次迭代过程中(每次滚动重新求解凸子问题时均需要进行若干次迭代才能获得最优控制序列)卫星j的邻居集合；
[0055]
只需要对集合内的卫星进行碰撞约束即可，可大大减少避碰约束；
[0056]
凸化后的离散点处避免外在障碍物约束：
[0057][0058]
其中，oi′
代表第i
′
个障碍物的位置状态信息，o代表所有障碍物位置状态信息集合，r
obs,i
′
代表第i
′
个障碍物的最小安全距离，kh代表预测步长；
[0059]
凸化后的任意两离散点之间的避障约束，保证重构全过程都满足避障约束：
[0060][0061]
其中，oi′
代表第i
′
个障碍物的位置信息，o代表所有障碍物状态信息集合，代表第j颗卫星在第k-1步时的基准轨迹，r
obs,i
′
代表第i
′
个障碍物的最小安全距离，kh代表预测步长。
[0062]
其它步骤及参数与具体实施方式一或二相同。
[0063]
具体实施方式四：结合图2说明本实施方式。本实施方式与具体实施方式一至三之一不同的是，所述步骤s3的具体过程为：
[0064]
步骤s31、当预测时域不包含终点时，建立的凸子问题p1如下：
[0065][0066][0067]
式(11)为信赖域约束，iter代表迭代次数，x
j,iter
[k]代表第iter次迭代时第j颗卫星在第k步时刻的状态信息，代表第iter次迭代时第j颗卫星在第k步时刻的基准轨迹，为信赖域半径，信赖域半径随着迭代次数增加逐渐减小；||
·
||
p
代表p范数，p
＝1，||
·
||
∞
代表无穷范数，δt1和δt2分别是预测时域内和预测时域外的离散时间步长。
[0068]
为等式约束，包含动力学约束，为式(3)，为不等式约束，包含推力幅值约束(式(5))、速度幅值约束(式(6))、离散点处卫星间避碰约束(式(7))和避障约束(式(8))以及任意两离散点之间避障约束(式(9))。基准轨迹xj是对实际轨迹xj的初始猜测，并且用于凸化碰撞规避约束，k代表离散化过程中离散区间数。
[0069]
目标代价函数分为三项，第一项为预测时域内累计的目标代价，注意在预测时域内考虑全部等式和不等式约束，确保卫星间和与障碍物之间的碰撞避免。第二项为预测时域以后累计的代价函数，不考虑卫星间碰撞规避以及和障碍物碰撞规避约束。第三项为包含不等式约束和等式约束的惩罚项，用以提升可行解的概率。
[0070]
步骤s32、当预测时域包含终点时，建立的凸子问题p2如下：
[0071][0072]
为等式约束，包含动力学约束，为式(3)，为不等式约束，包含推力幅值约束(式(5))、速度幅值约束(式(6))、离散点处卫星间避碰约束(式(7))和避障约束(式(8))以及任意两离散点之间避障约束(式(9))。其中，卫星间避碰约束和避障约束需要修正，其余约束均保持不变。
[0073]
对离散点处卫星间避碰约束和避障约束进行修正如下：
[0074][0075][0076][0077]
其中，为考虑不确定性干扰时卫星i与卫星j之间需满足的最小碰撞约束距离。
[0078]
当预测时域终点到达任务终点时，即凸子问题p2所示，目标代价函数分为两项，与p1相比减少的是预测时域外的代价函数项。
[0079]
其它步骤及参数与具体实施方式一至三之一相同。
[0080]
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述步骤s4的具体过程为：
[0081]
步骤s41、判断预测时域是否包含终点，即判断是否满足：
zj]
t
，η，ηj，ξ，ξj，rj和r
jz
是为了简化势能项；
[0099][0100]
式中，μ代表地球引力常数，μ＝398600.4418，单位是km3/s2，k
j2
＝2.633
×
10
10
，单位是km5/s2。
[0101]
通过公式(17)可以解出xj,yj,zj和也就是更新的卫星集群最新状态(就是位置速度)信息。
[0102]
其它步骤及参数与具体实施方式一至五之一相同。
[0103]
本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。