一种月球天梯总体方案优化方法与流程

本发明提供一种月球天梯总体方案优化方法，它涉及一种基于改进重采样粒子群法的月球天梯总体方案优化方法，可以快速、有效地实现月球天梯总体方案优化设计，属于航空航天领域。

背景技术：

在人类探索太空的整个工程技术体系中，运载技术一直占据着及其重要的地位，运载能力的强弱直接决定了航天器的体量。从开始探索太空至今，运载火箭一直是主要的运载工具，而本发明中所针对的月球天梯是一种新概念运载技术。月球天梯可以理解为在月球上建立的太空电梯，其主体是一条缆绳，一端固定连接在月面平台，另一端同向太空与配种相连，整个缆绳保持紧绷。攀爬器可以在缆绳上运行，可以在月球表面与空间站之间运送货物和人员。月球天梯相对于传统的运载方式的优点有：运输成本极低，经过初步估算，使用天梯将1kg货物送入太空只需100美元，不到运载火箭的1％；运输效率高，基本可以实现不间断运行。月球天梯的建造对于月球基地的建造有着重要的意义，也可以为其他星球天梯的建造提供经验，它将进一步促进人类对太空的探索和利用。

在月球天梯的设计和建造过程中，其总体方案的评价和总体方案的优化是关键的步骤。对月球天梯总体方案的评价方法和系统已经另外申请专利，本专利提出一种月球天梯总体方案的优化设计方法。

所谓月球天梯总体方案的优化设计，就是利用计算机等辅助工具，采用数值计算和分析的方法寻找满足要求的最优设计方案。月球天梯是一个十分庞大而复杂的工程系统，可能有很多个局部最优解，同时其设计参数变化范围很大，所以传统的优化方法，如牛顿法、序列二次规划法等的求解效果一般。因此需要采用现代智能优化算法，现代智能优化算法是指通过计算机编程模拟自然现象，模仿动物乃至人类的社会行为和进化机制，从而实现对复杂优化问题求解的一大类算法的统称。其中重采样粒子群优化算法，是模拟鸟类种群觅食行为并融合了粒子滤波领域中重采样技术的一种现代智能优化算法。其优点是原理简单，实现容易，而且它在大型工程优化问题中体现出了良好的效果，但是在月球天梯总体方案的优化设计中，重采样粒子群优化方法并不能求解出满足要求的最优解，有时计算时间很长，有时陷入局部最优区域从而不能得到最优的设计方案。

为了解决上述问题，本发明在重采样粒子群优化方法中引入变异操作和震动操作，并定义了种群聚集度和粒子活跃值两个指标，提出了一种月球天梯总体方案优化方法。

技术实现要素：

1、目的

本发明提供一种月球天梯总体方案优化方法，以克服传统优化设计方法在一定程度上计算速度慢，以及容易陷入局部最优解的缺点，从而高效并可靠地完成月球天梯总体方案的优化设计。

2、技术方案

为了实现上述发明目的，本发明采用以下技术方案。

本发明提供一种月球天梯总体方案优化方法，即一种基于改进重采样粒子群法的月球天梯总体方案优化方法，主要包括以下几个步骤：

步骤一：建立月球天梯总体方案的优化模型

所述月球天梯总体方案包含了对月球天梯各个主要部分的形状、材料、性能等参数的选择和计算，并最终可以抽象表达为一组数学公式和多段软件程序的集合；

所述“优化模型”主要包含设计变量、约束条件和优化目标；

该“设计变量”是指月球天梯总体方案中各个主要部分的形状、材料、性能等参数，如缆绳的长度、缆绳的横截面积、缆绳的应力等。最优的设计变量是通过优化方法求解出来的，一旦得到最优的设计变量也就得到了最优的总体方案。

该“约束条件”是指月球天梯总体方案必须满足的条件，如成本不能超过某一个上限、运载能力必须高于某一个下限等。

该“优化目标”是指月球天梯总体方案优劣的评价分数，方案的评价分数越高，方案越好。方案评价分数的求解方法和系统已经另外申请专利，这里不详细说明，本发明主要是对上述优化模型的优化求解方法；

上述优化模型最终可以表示为：

maxf(x)

s.t.g(x)≤0

其中，x是设计变量，g(x)是约束条件，f(x)是优化目标，也就是方案的评价分数；

步骤二：初始化粒子种群

所述初始化粒子种群包括：确定粒子种群中的粒子个数n、生成每个粒子的初始位置坐标x和初始速度矢量v、以及初始化每个粒子的历史最优位置pbest和群体的最优位置gbest；

该“粒子个数n”的取值与具体优化问题有关，一般取10～30；

该“初始位置坐标x”的生成方法为：

xid＝xmin(d)+rand1id·(xmax(d)-xmin(d))

其中，xid是第i个粒子第d维的坐标值，xmin(d)和xmax(d)分别是粒子第d维坐标值的下限和上限，rand1是一组0～1之间的随机数。

该“初始速度矢量v”的生成方法为：

vid＝xmin(d)+rand2id·(xmax(d)-xmin(d))-xid

其中，vid是第i个粒子第d维的速度值，xmin(d)和xmax(d)以及xid的含义同上，rand2是一组0～1之间的随机数；

该“每个粒子的历史最优位置pbest”的初始化方法为：记粒子的初始位置为粒子历史最优位置pbest的初始值，即pbest＝x。同时求出目标函数值，称之为历史最优值，记为pbest_f；

该“群体的最优位置gbest”的初始化方法为：比较上述每个粒子的历史最优值，其中历史最优值最小的粒子的位置为群体最优位置gbest的初始值，该历史最优值为当前的全局最优值，记为gbest_f；

步骤三：计算种群聚集度并对粒子种群进行重采样操作

所述“种群聚集度”是本发明中定义的用来描述种群中粒子多样性的宏观指标，种群聚集度越高说明粒子分布越集中，种群多样性越差；种群聚集度越低说明粒子分布越分散，种群多样性越好；种群聚集度的计算公式为：

式中，n为粒子个数，d为维数，disi为第i个粒子到中心位置的距离，是所有粒子到中心位置的平均距离，xij是第i个粒子第j维的坐标值，是中心位置第j维的坐标值，中心位置是在搜索空间内随机产生的；上式的含义为：种群聚集度是种群中所有粒子到一个随机产生的中心位置的距离的方差的倒数；

随着计算的推进，种群的多样性逐渐降低，粒子聚集度不断提高，当粒子聚集度大于临界值，即ad＞adth时，对粒子种群进行重采样操作；

该“重采样操作”是为了克服粒子群优化算法在处理月球天梯总体方案优化设计问题时所存在的收敛速度慢和易陷入局部最优的缺点而进行的特殊操作，其具体实施过程为：

首先对每个粒子按照高斯分布的规律进行权值分配，距离当前群体最优位置越近的粒子权值越大，越远的粒子权值越小，分配公式为：

其中qi为对第i个粒子赋予的权值，f(xi)为适应度函数，gbest为当前群体最优位置，σ为以f(xi)-gbest为样本计算所得的方差，qi为第i个粒子归一化后的权值；

然后对每个粒子的权值进行判断，当某个粒子的权值小于给定的阀值qt时，便以pr的概率随机产生新的粒子取代之：

当qi＜qt时，

其中，是新的粒子位置坐标，其确定方法为：

其中，t为当前迭代次数，xmin和xmax分别是粒子坐标值的下限和上限，rand3是一组0～1之间的随机数；

同时根据自适应速度修正公式修正该粒子的速度矢量：

其中t为最大迭代次数，t为当前迭代次数，为新引入的粒子速度；

其中，新引入的粒子速度的确定方法为：

其中，xmin和xmax分别是粒子坐标值的下限和上限，rand4是一组0～1之间的随机数，是上一步产生的新的粒子位置；

步骤四：更新粒子位置和速度

所述更新粒子位置和速度为：每个粒子在各自历史最优位置pbest和群体最优位置gbest的影响下根据给定规律生成新的位置坐标和速度矢量，速度矢量的更新方式为：

vi(t+1)＝χ{vi(t)+c1r1[pbesti-xi(t)]+c2r2[gbest-xi(t)]}

该“速度矢量更新公式”由三项构成，第一项vi(t)为原速度，第二项c1r1[pbesti-xi(t)]为个体历史最优位置对速度的影响，第三项c2r2[pg-xi(t)]为群体最优位置对速度的影响；

其中c1、c2为加速度系数，代表着个体历史最优位置和群体最优位置对速度影响作用的大小，其取值和具体优化问题有关；r1、r2为随机因子，是0～1之间的随机数；χ是压缩因子，其确定方法为：

位置更新的方式为：

xi(t+1)＝xi(t)+vi(t+1)

上式由两项构成，第一项xi(t)为原位置，第二项vi(t+1)为根据上述速度更新公式得到的新的速度；

步骤五：更新每个粒子的历史最优位置和群体最优位置

经过步骤三和步骤四的操作，粒子的位置更新，导致每个粒子的历史最优位置pbest和群体最优位置gbest有所变化，所以要更新每个粒子的历史最优位置pbest和群体最优位置gbest；

所述“每个粒子的历史最优位置pbest”的更新过程为：

首先，求出每个粒子当前位置所对应的目标函数值f(xi(t+1))；

然后，比较f(xi(t+1))和该粒子的历史最优值pbest_fi；

iff(xi(t+1))＜pbest_fi

thenpbesti＝xi(t+1)

pbest_fi=f(xi(t+1))

上式中，脚标i代表第i个粒子；上式的含义为，当粒子当前位置对应的目标函数值小于该粒子历史最优值时，用当前位置作为该粒子的历史最优位置，用该目标函数值作为该粒子的历史最优值；否则，pbest和pbest_f保留原值。

所述“群体的最优位置gbest”的更新过程为：比较上述每个粒子的历史最优值，其中历史最优值最小的粒子的位置为群体最优位置gbest，该历史最优值为当前的全局最优值，记为gbest_f；

步骤六：计算粒子活跃值并进行变异操作

所述“粒子活跃值”是用来描述每个粒子活跃程度的微观指标，如果在某次迭代中粒子的目标函数值有所增加，那么该粒子为活跃粒子，其活跃值重置为最大值，记为ac＝acmax；如果在某次迭代中粒子的目标函数值没有增加，那么该粒子的活跃值减1；当某粒子的活跃值等于0时，该粒子为不活跃粒子，此时对该粒子进行变异操作同时重置其活跃值；

所述“变异操作”是指在粒子当前位置和当前群体最优位置的中点附近以正态分布形式随机产生新的粒子并取代之，其计算的公式为:

其中x^*表示目前为止所找到的最优解，norm(μ，σ)表示均值为μ、方差为σ的高斯随机数；

步骤七：进行震动操作

每当种群找到更优的位置，就进行小范围震动操作，其作用是通过在当前最优点的震动寻找它附近更优的位置，从而提高算法的局部探索能力；其具体过程为：先随机选择k维，k是1到预先设定的kmax上的随机数，在选出的每一维上产生一个震动点x’

x′(d)＝x(d)+norm(0，1)·βl(d)

其中β为震动系数；如果没有震动点比当前最优点更优，则保留当前最优点，否则用更优的震动点替换当前最优点作为新的最优点；

步骤八：如果不满足精度要求且尚未达到最大迭代次数，迭代次数加一，返回步骤三，否则记录并输出结果；

所述“结果”主要包括该月球天梯总体方案优化设计问题的最优设计点和目标函数最优值；所述“最优设计点”即当前的群体最优位置，它是月球天梯总体设计的最优方案，所述“目标函数最优值”即当前全局最优值，它是该最优方案的评价分数。

通过以上流程和步骤，可以高效并可靠地处理月球天梯总体方案优化设计问题，而且由于本发明引入了变易操作和震动操作，一方面加快求解速度，提高效率，另一方面改善了全局搜索的能力，提高了优化精度。

其中，在步骤一中所述的“建立月球天梯总体方案的优化模型”，其建立的过程如下：

首先根据月球天梯方案的总体要求，确定优化模型的设计变量，如缆绳的长度、缆绳的横截面积、缆绳的应力等；

然后根据月球天梯方案的总体需求，确定优化模型的优化目标，优化目标是一个或多个关键的指标，在这里优化目标就是根据月球天梯总体方案评价方法获得的方案的评价分数，该“月球天梯总体方案评价方法”已另外申请专利，在这里不再详细描述；

最后根据月球天梯方案的总体要求，确定优化模型的约束条件，它是我们在设计中必须满足的条件，其具体指标由用户给出，如成本不能超过某一个上限、运载能力必须高于某一个下限等。

其中，在步骤二中所述的“初始化粒子种群”，其初始化的作法如下：

先确定粒子个数n，一般取10～30；

然后生成初始位置坐标x，其方法为：

xid＝xmin(d)+rand1id·(xmax(d)-xmin(d))

其中，xid是第i个粒子第d维的坐标值，xmin(d)和xmax(d)分别是粒子第d维坐标值的下限和上限，rand1是一组0～1之间的随机数。

再生成初始速度矢量v，其方法为：

vid＝xmin(d)+rand2id·(xmax(d)-xmin(d))-xid

其中，vid是第i个粒子第d维的速度值，xmin(d)和xmax(d)以及xid的含义同上，rand2是一组0～1之间的随机数；

然后初始化每个粒子的历史最优位置pbest，其方法为：记粒子的初始位置为粒子历史最优位置pbest的初始值，即pbest＝x。同时求出目标函数值，称之为历史最优值，记为pbest_f；

最后初始化群体的最优位置gbest，其方法为：比较上述每个粒子的历史最优值，其中历史最优值最小的粒子的位置为群体最优位置gbest的初始值，该历史最优值为当前的全局最优值，记为gbest_f；

其中，在步骤六中所述的“计算粒子活跃值并进行变异操作”，其作法如下：

首先，计算每个粒子的活跃值aci，

然后，判断aci是否等于0，如果aci不为0，则不进行任何操作，如果aci等于0，则根据以下公式求出该粒子经过变异操作后的粒子。

。

3、优点及功效

本发明提供的一种月球天梯总体方案优化方法，它涉及一种基于改进重采样粒子群法的月球天梯总体方案优化方法，其主要优点是：可以实时监测种群的多样性和每个粒子的活跃程度并做出相应处理，从而合理分配计算能力，提高计算效率；同时通过变易操作提高了全局搜索能力，通过震动操作提高了局部探索能力；最终实现了高效并可靠地完成月球天梯总体方案的优化设计的目的。

附图说明

图1本发明所述方法流程图。

图2重采样操作流程图。

图3变异操作流程图。

图4震动操作流程图。

图5月球天梯方案综合优化流程图。

具体实施方式

以下结合附图和一个月球天梯总体优化设计实例对本方法作进一步描述，但本实例并不用于限制本方法，绝大部分月球天梯优化设计问题都可以用本方法优化求解。

由图1可以看出，本发明提供的基于改进重采样粒子群优化算法的月球天梯总体方案优化方法，主要包括以下八个步骤：

步骤一：建立月球天梯总体方案的优化模型

为了阐明方法的逻辑过程，这里结合一个月球天梯总体优化设计实例对本方法作详细描述，该实例的综合优化流程如图5所示。

本实例的设计变量有三个：天梯缆绳的长度、天梯缆绳的初始横截面积、天梯缆绳的应力。

本实例的约束条件有五个：天梯缆绳长度在1000000km到3000000km之间，天梯缆绳的初始横截面积在20mm²到50mm²之间，天梯缆绳的应力在20gpa到43.3gpa之间，总成本不超过2000亿美元，运载能力不小于5000吨。

本实例的优化目标是方案的评价分数。方案评价分数的求解方法和系统已经另外申请专利，这里不详细说明，本发明主要是对上述优化模型的优化求解方法。

上述优化模型最终可以表示为：

其中，r是方案的评价分数，l是缆绳长度，a0是缆绳的初始横截面积，σ是缆绳应力，c是总成本，m是运载能力。lmin和lmax是缆绳长度的上下限，amin和amax是缆绳初始横截面积的上下限，[σ]min和[σ]max是缆绳应力的上下限，cmax是最大成本，mmin是最小运载能力。且lmin＝1000000km，lmax＝3000000km，amin＝20mm²，amax＝50mm²，[σ]min＝20gpa，[σ]max＝43.3gpa，cmax＝2000亿美元，mmin＝5000t。

步骤二：初始化粒子种群

所述初始化粒子种群包括：确定粒子种群中的粒子个数n、生成每个粒子的初始位置坐标x和初始速度矢量v、以及初始化每个粒子的历史最优位置pbest和群体的最优位置gbest。

该“粒子个数n”的取值与具体优化问题有关，一般取10～30，在本例中取n＝20。本例的维数为d＝3。

该“初始位置坐标x”的生成方法为：

xid＝xmin(d)+rand1id·(xmax(d)-xmin(d))

其中，xid是第i个粒子第d维的坐标值，xmin(d)和xmax(d)分别是粒子第d维坐标值的下限和上限，rand1是一组0～1之间的随机数。

该“初始速度矢量v”的生成方法为：

vid＝xmin(d)+rand2id·(xmax(d)-xmin(d))-xid

其中，vid是第i个粒子第d维的速度值，xmin(d)和xmax(d)以及xid的含义同上，rand2是一组0～1之间的随机数。

该“每个粒子的历史最优位置pbest”的初始化方法为：记粒子的初始位置为粒子历史最优位置pbest的初始值，即pbest＝x。同时求出目标函数值，称之为历史最优值，记为pbest_f。

该“群体的最优位置gbest”的初始化方法为：比较上述每个粒子的历史最优值，其中历史最优值最小的粒子的位置为群体最优位置gbest的初始值，该历史最优值为当前的全局最优值，记为gbest_f。

步骤三：计算种群聚集度并对粒子种群进行重采样操作

所述种群聚集度是本发明中定义的用来描述种群中粒子多样性的宏观指标，种群聚集度越高说明粒子分布越集中，种群多样性越差；种群聚集度越低说明粒子分布越分散，种群多样性越好。种群聚集度的计算公式为：

式中，n为粒子个数，d为维数，disi为第i个粒子到中心位置的距离，是所有粒子到中心位置的平均距离，xij是第i个粒子第j维的坐标值，是中心位置第j维的坐标值。中心位置是在搜索空间内随机产生的。上式的含义为：种群聚集度是种群中所有粒子到一个随机产生的中心位置的距离的方差的倒数；

随着计算的推进，种群的多样性逐渐降低，粒子聚集度不断提高，当粒子聚集度大于临界值，即ad＞adth时，对粒子种群进行重采样操作；

该“重采样操作”是为了克服粒子群优化算法在处理月球天梯总体方案优化设计问题时所存在的收敛速度慢和易陷入局部最优的缺点而进行的特殊操作，其流程如图2所示，其具体实施过程为：

首先对每个粒子按照高斯分布的规律进行权值分配，距离当前群体最优位置越近的粒子权值越大，越远的粒子权值越小，分配公式为：

其中qi为对第i个粒子赋予的权值，f(xi)为适应度函数，gbest为当前群体最优位置，σ为以f(xi)-gbest为样本计算所得的方差。qi为第i个粒子归一化后的权值。

然后对每个粒子的权值进行判断，当某个粒子的权值小于给定的阀值qt时，便以pr的概率随机产生新的粒子取代之：

当qi＜qt时，

其中，是新的粒子位置坐标，其确定方法为：

其中，t为当前迭代次数，xmin和xmax分别是粒子坐标值的下限和上限，rand3是一组0～1之间的随机数。

同时根据自适应速度修正公式修正该粒子的速度矢量：

其中t为最大迭代次数，t为当前迭代次数，为新引入的粒子速度。

其中，新引入的粒子速度的确定方法为：

其中，xmin和xmax分别是粒子坐标值的下限和上限，rand4是一组0～1之间的随机数。是上一步产生的新的粒子位置

步骤四：更新粒子位置和速度

所述更新粒子位置和速度为：每个粒子在各自历史最优位置pbest和群体最优位置gbest的影响下根据给定规律生成新的位置坐标和速度矢量，速度矢量的更新方式为：

vi(t+1)＝χ{vi(t)+c1r1[pbesti-xi(t)]+c2r2[gbest-xi(t)]}

该“速度矢量更新公式”由三项构成，第一项vi(t)为原速度，第二项c1r1[pbesti-xi(t)]为个体历史最优位置对速度的影响，第三项c2r2[pg-xi(t)]为群体最优位置对速度的影响；

其中c1、c2为加速度系数，代表着个体历史最优位置和群体最优位置对速度影响作用的大小，其取值和具体优化问题有关；r1、r2为随机因子，是0～1之间的随机数；χ是压缩因子，其确定方法为：

位置更新的方式为：

xi(t+1)＝xi(t)+vi(t+1)

上式由两项构成，第一项xi(t)为原位置，第二项vi(t+1)为根据上述速度更新公式得到的新的速度；

步骤五：更新每个粒子的历史最优位置和群体最优位置

经过步骤三和步骤四的操作，粒子的位置更新，导致每个粒子的历史最优位置pbest和群体最优位置gbest有所变化，所以要更新每个粒子的历史最优位置pbest和群体最优位置gbest；

所述“每个粒子的历史最优位置pbest”的更新过程为：

首先，求出每个粒子当前位置所对应的目标函数值f(xi(t+1))；

然后，比较f(xi(t+1))和该粒子的历史最优值pbest_fi；

iff(xi(t+1))＜pbest_fi

thenpbesti＝xi(t+1)

pbest_fi＝f(xi(t+1))

上式中，脚标i代表第i个粒子。上式的含义为，当粒子当前位置对应的目标函数值小于该粒子历史最优值时，用当前位置作为该粒子的历史最优位置，用该目标函数值作为该粒子的历史最优值。否则，pbest和pbest_f保留原值；

所述“群体的最优位置gbest”的更新过程为：比较上述每个粒子的历史最优值，其中历史最优值最小的粒子的位置为群体最优位置gbest，该历史最优值为当前的全局最优值，记为gbest_f。

步骤六：计算粒子活跃值并进行变异操作

所述“粒子活跃值”是用来描述每个粒子活跃程度的微观指标，如果在某次迭代中粒子的目标函数值有所增加，那么该粒子为活跃粒子，其活跃值重置为最大值，记为ac＝acmax；如果在某次迭代中粒子的目标函数值没有增加，那么该粒子的活跃值减1；当某粒子的活跃值等于0时，该粒子为不活跃粒子，此时对该粒子进行变异操作同时重置其活跃值。

所述“变异操作”是指在粒子当前位置和当前群体最优位置的中点附近以正态分布形式随机产生新的粒子并取代之，其流程如图3所示，变异公式为:

其中x^*表示目前为止所找到的最优解，norm(μ，σ)表示均值为μ、方差为σ的高斯随机数。

步骤七：进行震动操作

每当种群找到更优的位置，就进行小范围震动操作，其作用是通过在当前最优点的震动寻找它附近更优的位置，从而提高算法的局部探索能力。其具体过程为：先随机选择k维，k是1到预先设定的kmax上的随机数，在选出的每一维上产生一个震动点x’，其流程如图4所示，震动公式为:

x′(d)＝x(d)+norm(0，1)·βl(d)

其中β为震动系数。如果没有震动点比当前最优点更优，则保留当前最优点，否则用更优的震动点替换当前最优点作为新的最优点。

步骤八：如果不满足精度要求且尚未达到最大迭代次数，迭代次数加一，返回步骤三，否则记录并输出结果。

所述“结果”主要包括该月球天梯总体方案优化设计问题的最优设计点和目标函数最优值。所述“最优设计点”即当前的群体最优位置，它是月球天梯总体设计的最优方案，所述“目标函数最优值”即当前全局最优值，它是该最优方案的评价分数。

通过数值计算和分析，得到本例的最有设计方案为：l＝204590km，a0＝44.75mm²，σ＝43.06gpa。

该方案的质量规模为19.6万吨，l1点年运载能力为12180吨，建造时间为19.9年，建造成本为2000亿美元。

通过以上流程和步骤，可以有效处理月球天梯总体方案优化设计的实际问题，而且由于本发明引入了两个指标：种群聚集度和粒子活跃值，可以实时监测种群的多样性和粒子的活跃程度，从而合理分配计算资源；同时本发明引入的两项操作：变异操作和震动操作，一方面改善了全局搜索能力，避免陷入局部最优解，一方面提高了局部探索能力，提高了收敛速度。最终实现了高效并可靠地完成月球天梯总体方案的优化设计的目的。

上面对本专利的实施方式作了详细说明，但是本专利并不限于上述实施方式，在本领域的普通技术人员所具备的知识范围内，还可以在不脱离本专利宗旨的前提下做出各种变化。