CN114488788B - 概率自适应花授粉优化的高超声速飞行器再入段控制器 - Google Patents

Abstract

本发明公开了一种概率自适应花授粉优化的高超声速飞行器再入段控制器，该控制器由飞行器海拔高度传感器、飞行器飞行速度传感器、飞行器飞行航道倾角传感器、飞行器水平飞行距离传感器、飞行器MCU、飞行器攻角控制器构成。高超声速飞行器到达再入段空域后，各传感器开启获得飞行器海拔高度、速度、飞行航道倾角和飞行水平距离信息，MCU根据所获得的信息自动执行内部概率自适应花授粉优化算法，得到水平飞行距离最长的轨迹优化控制策略并转换为控制指令，发送给攻角控制器执行。本发明能够快速地得到轨迹优化控制策略，使高超声速飞行器获得更长的水平飞行距离。

Description

概率自适应花授粉优化的高超声速飞行器再入段控制器

技术领域

本发明涉及高超声速飞行器技术领域，尤其涉及一种概率自适应花授粉优化的高超声速 飞行器再入段轨迹优化控制器。

背景技术

高超声速飞行器MCU是实现远程快速精确打击以及全球快速到达的新型飞行器，在未 来的军事、政治和经济中具有十分重要的战略地位，已经成为世界航空航天领域一个极其重 要的发展方向，研究和发展高超声速飞行器在开发太空和国家安全方面具有非常重要的意义。在高超声速飞行器MCU的研究中，轨迹优化不仅有利于提高飞行器飞行品质以满足既定任 务要求，同时也是完成飞行任务的重要保证和实现机动飞行的必要条件。由于从大气从外缘 进入大气层，高度和速度的变化范围很大，高超声速飞行器面临各种严峻的再入环境，再入 段轨迹优化控制技术则是保证高超声速飞行器完成飞行任务的关键，对提高其打击范围和落点精度具有更重要的实用价值，其中如何获得更长的水平飞行距离，来提高打击范围，是当 前国内外研究的热点和难点。

发明内容

本发明的目的在于针对现有技术的不足，提供了一种概率自适应花授粉优化的高超声速 飞行器再入段控制器。本发明在高超声速飞行器到达再入段后能够给出高超声速飞行器轨迹 优化控制策略并转换为飞行器攻角控制指令，使高超声速飞行器获得更长的水平飞行距离，提高高超声速飞行器的打击范围。

本发明的目的是通过以下技术方案来实现的：一种概率自适应花授粉优化的高超声速飞 行器再入段轨迹优化控制器，根据高超声速飞行器再入段初始海拔高度、速度、飞行航道倾 角和飞行水平距离状态获取轨迹优化控制策略，通过控制飞行器攻角使高超声速飞行器获得更长的水平飞行距离。所述的轨迹优化控制器由飞行器海拔高度传感器、飞行器飞行速度传 感器、飞行器飞行航道倾角传感器、飞行器水平飞行距离传感器、飞行器MCU(Micro Control Unit，微控制单元)、飞行器攻角控制器构成，所述飞行器海拔高度传感器、飞行器飞行速度 传感器、飞行器飞行航道倾角传感器、飞行器水平飞行距离传感器、飞行器攻角控制器均通 过数据总线与飞行器MCU相连。所述控制器的运行过程如下：

步骤1)：在高超声速飞行器MCU中输入对应于该飞行器的气动系数模型、飞行器性能 约束条件、指定优化目标；

步骤2)：高超声速飞行器到达再入段后，开启飞行器海拔高度传感器、飞行器飞行速度 传感器、飞行器飞行航道倾角传感器和飞行器水平飞行距离传感器，得到高超声速飞行器当 前的海拔高度、速度、飞行航道倾角和飞行水平距离状态信息；

步骤3)：飞行器MCU根据设定的海拔高度、速度、飞行航道倾角要求自动执行内部概 率自适应花授粉优化算法，得到使高超声速飞行器水平飞行距离最长的轨迹优化控制策略；

步骤4)：飞行器MCU将所获得的轨迹优化控制策略发送给控制指令输出模块，并转换 为控制指令发送给飞行器攻角控制器执行。

所述的高超声速飞行器MCU，包括信息采集模块、初始化模块、系统状态求解模块、概 率自适应花授粉优化模块、收敛性判断模块、控制指令输出模块。其中，信息采集模块包括 飞行器海拔高度和速度采集子模块、飞行器飞行航道倾角和飞行水平距离采集子模块、飞行器海拔高度和速度设定采集子模块、飞行器飞行航道倾角设定采集子模块、飞行器的气动系 数模型和性能约束条件以及指定优化目标参数采集子模块。

所述的高超声速飞行器MCU自动产生攻角控制指令的概率自适应花授粉优化算法运行 步骤如下：

步骤1)：高超声速飞行器MCU到达再入段后，飞行器海拔高度传感器、飞行器飞行速 度传感器、飞行器飞行航道倾角传感器和飞行器水平飞行距离传感器开启，信息采集模块获 取高超声速飞行器当前的海拔高度、速度、飞行航道倾角和飞行水平距离状态信息；

步骤2)：初始化模块开始运行，对各项参数初始化，包括控制变量参数化参数和概率自 适应花授粉优化方法参数。控制变量参数化参数的初始化设置如下：设置时间段的分段数为 NE，设置控制变量参数化的插值拟合方式：分段常数或者分段线性；概率自适应花授粉优化方法参数的初始化设置如下：设置种群规模N，种群中第i个个体的初始位置向量，种 群初始历史最优位置，其中，J(·)表示要最大化的目标函数。 设置个体寻优最大迭代次数max_iters，花粉密度λ，初始迭代次数k＝0。

步骤3)：通过系统状态求解模块获得第k次迭代种群个体i的目标函数值，更新 第k次迭代种群历史最优位置g^(k),best，记录最小目标函数值

步骤4)：执行收敛性判断模块，如果当前迭代次数k达到最大迭代次数max_iters，则收 敛性条件得到满足，当前种群历史最优位置g^(k),best作为最优控制策略的离散解，经控制指令 输出模块输出；如果收敛性条件不满足，则迭代次数k加1。

步骤5)：执行概率自适应花授粉优化模块，依据种群分布状况计算表征个体聚集程度的 指标，并计算相应的转换概率，由随机数生成器产生均匀分布随机数，根据随机数与 转换概率p的相对大小更新个体的位置。该步骤执行完毕后跳转至步骤3，直至收敛性条件 得到满足为止。

其中，所述的系统状态求解模块，采用的是四步四阶Runge-Kutta方法，计算公式为：

其中，t_m表示四步四阶Runge-Kutta方法选择的积分时刻，t_m+1表示位于时刻t_m后的积分时刻， 并且t_m+1＝t_m+h，h为积分步长，F(·)是描述状态微分方程的函数，K1、K2、K3、K4分 别表示Runge-Kutta法积分过程中的4个节点的函数值，x(·)是描述系统状态变量的函数，且x(t_m)表示时刻t_m的系统状态变量，x(t_m+1)表示时刻t_m+1的系统状态变量。

进一步地，所述的概率自适应花授粉优化模块，采用如下步骤实现：

步骤1)：定义加入惩罚项后的增广性能指标函数，形式为：

其中，t表示积分时刻，t_f表示终值时刻，X(t)表示加入x_n+1(t)后的新状态向量，Φ₀[X(t_f)] 表示新状态向量在终值时刻处的终值函数，表示n_u维控制 向量，表示n_x维状态向量，ρ为惩罚因子，l表示光滑化 函数参数。H(t,x(t),u(t),ρ,l)为约束相关的惩罚项，具体形式为：

其中，Q_i(t,x(t),u(t),ρ,l)为表征约束违背量的函数，θ(Q_i(·))为多阶段配置函数，m₁表示 标准等式约束的个数，m₂表示标准不等式约束的个数。ε为光滑化因子，其取值为一个很小 的整数。

γ(Q_i(·))为惩罚函数指数，按以下方式更新：

θ(Q_i(·))为多阶段配置函数，定义如下：

基于以上方式，惩罚项因子的值随约束违背量的大小自适应地变化。当约束违背量较大 时，采取较大的惩罚因子，将更多的候选解推至可行区域内；当约束为违背量逐渐减小时， 惩罚因子也随之减小，保证进化算法在可行最优解周围搜索，更利于发现可行域边界周围的全局最优值。

步骤2)：按如下公式获得个体i(i＝1,2,…,N)与其他个体的平均欧式距离d_i为：

其中，N表示种群规模，u_i、u_j表示种群中的两个个体。

步骤3)：定义密集度指数ρ表征种群密集度，按如下公式获得：

其中，d_g为当前迭代下全局最优个体与其他个体间的平均距离，d_min、d_max为d_i(i＝1,2,…,N) 的最大和最小值。

步骤4)：按如下公式获得转换概率p为：

步骤5)：由随机数生成器产生均匀分布随机数

步骤6)：根据与p的相对大小更新个体的位置。如果概率p大于随机数，个体依据 式(5)按照Lévy飞行方式更新；概率p小于随机数时，由下式计算更新后位置；

其中，gbest表示花粉中的最优个体，L(λ)是与花粉密度λ相关的Lévy飞行步长，可由Lévy 分布求得：

其中，Γ(λ)表示标准γ函数，s为步长，s₀为初始步长。

本发明的有益效果主要表现在：本发明获得使高超声速飞行器水平飞行距离更长的轨迹 优化攻角控制指令，提高高超声速飞行器打击范围。

附图说明

图1是本发明的结构示意图；

图2是本发明高超声速飞行器MCU内部模块结构图；

图3是实施例1的高超声速飞行器状态轨迹图，其中，(a)为飞行器的海拔高度随时间 变化图，(b)为飞行器的速度随时间变化图，(c)为飞行器的航向角随时间变化图，(d)为飞行器的飞行航程随时间变化图。

具体实施方式

高超声速飞行器再入段轨迹优化问题可以描述为：

其中，t为时间，x(t)表示高超声速飞行器的状态变量，x₁(t)表示飞行器海拔高度、x₂(t) 表示飞行器飞行速度、x₃(t)表示飞行器飞行航道倾角、x₄(t)表示飞行器水平飞行距离，u(t) 表示飞行器的攻角控制量，为本问题的控制变量；表示状态变量x(t)的一阶导数， F(x(t),u(t),t)是根据高超声速飞行器再入段三维空间运动方程建立的微分方程组数学模型； t₀表示再入段轨迹优化开始的时间点，h₀表示优化开始时刻飞行器的初始海拔高度，v₀表示 优化开始时刻飞行器的初始速度，γ₀表示优化开始时刻飞行器的初始飞行航道角，r₀表示优 化开始时刻飞行器的初始水平飞行距离，t_f表示再入段轨迹优化结束时间点，h_f表示优化结 束时刻飞行器的海拔高度，v_f表示优化结束时刻飞行器的速度，γ_f表示优化结束时刻飞行器 的飞行航道角；J[u(t)]表示高超声速飞行器轨迹优化的目标函数即优化结束时刻飞行器的水平飞行距离，G[u(t),x(t),t]是高超声速飞行器再入段过程的约束条件，u_min和u_max表示攻 角控制范围的下限值和上限值。

本发明在高超声速飞行器MCU中集成了概率自适应花授粉优化算法，在高超声速飞行器到达再入段后能够给出速飞行器攻角的控制指令，使高超声速飞行器获得更长的水平飞行 距离。

实施例1

假设高超声速飞行器到达再入段空域，高超声速飞行器海拔高度传感器3、速度传感器4、 飞行航道倾角传感器5、水平飞行距离传感器6和MCU2均已开启。信息采集模块9立即采 集飞行器进入再入段时的初始海拔高度、速度、飞行航道倾角和水平飞行距离，设当前初始 时刻t₀＝0s，海拔高度传感器传入MCU的海拔高度为h₀＝80,000m，速度传感器传入MCU 的速度为v₀＝6400m/s，飞行航道倾角传感器传入MCU的飞行航道倾角为γ₀＝-0.052rad， 水平飞行距离传感器感器传入MCU的水平飞行距离为r₀＝0m；终值时刻t_f高超声速飞行器 需要满足的条件为海拔高度设定为h_f＝24000m，速度设定为v_f＝760m/s，飞行航道倾角设 定为γ_f＝-0.08rad；结合飞行器的三维空间运动方程、气动系数模型、飞行器性能约束条件 和指定优化目标，得到该问题的数学模型如下：

max J[u(t)]＝x₄(t_f)

C_L＝-0.15+3.44u(t)

C_D＝0.29-1.51u(t)+5.87u(t)²

x₁(0)＝80×10³,x₁(t_f)＝24×10³

x₂(0)＝6.4×10³,x₂(t_f)＝760

x₃(0)＝-0.052,x₃(t_f)＝-0.08

x₄(0)＝0

其中,L表示升力，D表示阻力，C_L表示升力系数，C_D表示阻力系数，h表示飞行海拔高度， v表示飞行速度，γ表示飞行航道倾角，u表示飞行器的水平飞行距离，x₁表示无量纲地心距， x₂表示无量纲速度，x₃表示航迹角，x₄表示经度。

为了便于表述，采用F(x(t),u(t),t)表示高超声速飞行器再入段三维空间运动方程建立 的微分方程组数学模型，即：

其中,L表示升力，D表示阻力，h表示飞行海拔高度，v表示飞行速度，γ表示飞行航道倾 角。

采用G[u(t),x(t),t]表示高超声速飞行器再入段过程的约束条件，为：

其中,L表示升力，D表示阻力，h表示飞行海拔高度，v表示飞行速度，γ表示飞行航道倾 角。

此外，J[u(t)]表示高超声速飞行器轨迹优化的目标函数即优化结束时刻飞行器的水平飞行距离。

高超声速飞行器MCU自动产生攻角控制指令的概率自适应花授粉优化算法如图2所示， 其运行步骤如下：

步骤1)：在高超声速飞行器MCU2中输入对应于该飞行器的气动系数模型、飞行器性能 约束条件、指定优化目标1；

步骤2)：高超声速飞行器到达再入段后，开启飞行器海拔高度传感器3、飞行器速度传 感器4、飞行器飞行航道倾角传感器5和飞行器水平飞行距离传感器6，得到高超声速飞行器 当前的海拔高度、速度、飞行航道倾角和飞行水平距离状态信息；

步骤3)：飞行器MCU 2根据设定的海拔高度、速度、飞行航道倾角要求自动执行内部 概率自适应花授粉优化算法13，得到使高超声速飞行器水平飞行距离最长的轨迹优化控制策 略；

步骤4)：高超声速飞行器MCU 2将获得的轨迹优化控制策略发送给控制指令输出模块 14，并转换为控制指令发送给飞行器攻角控制器8执行。

所述的高超声速飞行器MCU部分如图2所示，包括信息采集模块9、初始化模块10、系统状态求解模块11、收敛性判断模块12、概率自适应花授粉优化模块13、控制指令输出模块14。其中，信息采集模块9包括飞行器海拔高度和飞行速度采集、飞行器飞行航道倾角和飞行水平距离采集、飞行器海拔高度和速度设定采集、飞行器飞行航道倾角设定采集、飞行器的气动系数模型和性能约束条件以及指定优化目标参数采集五个子模块。

所述的高超声速飞行器MCU自动产生攻角控制指令的概率自适应优化算法运行步骤如 下：

步骤1)：高超声速飞行器MCU到达再入段后，飞行器海拔高度传感器、飞行器飞行速 度传感器、飞行器飞行航道倾角传感器和飞行器水平飞行距离传感器开启，信息采集模块9 获取初始时刻t₀＝0s时高超声速飞行器海拔高度h₀＝80,000m，飞行速度为v₀＝6400m/s， 飞行航道倾角为γ₀＝-0.052rad，水平飞行距离传感器感器水平飞行距离设置为r₀＝0m；终 值时刻t_f高超声速飞行器海拔高度要求设定为h_f＝24000m，飞行速度要求设定为 v_f＝760m/s，飞行航道倾角要求设定为γ_f＝-0.008rad；

步骤2)：初始化模块10开始运行，对各项参数初始化，设置时间段的分段数为NE，设 置控制变量参数化的插值拟合方式：分段常数或者分段线性；概率自适应花授粉优化方法参 数初始化如下：设置种群规模N，种群中第i个个体的初始位置向量，种群初始历史最优 位置，其中，J(·)表示要最大化的目标函数。设置个体寻优最 大迭代次数max_iters，花粉密度λ，初始迭代次数k＝0。

步骤3)：通过系统状态求解模块11，获得第k次迭代个体i的目标函数值，更新 第k次迭代种群历史最优位置g^(k),best，记录最大目标函数值和最小目标函数值

步骤4)：执行收敛性判断模块12，如果当前迭代次数k达到最大迭代次数max_iters，则 收敛性条件得到满足，当前种群历史最优位置g^(k),best作为最优控制策略的离散解，经控制指 令输出模块输出；如果收敛性条件不满足，则迭代次数k加1。

步骤5)：执行概率自适应花授粉优化模块13，依据种群分布状况计算表征个体聚集程度 的指标，并计算相应的转换概率，由随机数生成器产生均匀分布随机数，根据与p的 相对大小更新个体的位置。该步骤执行完毕后跳转至步骤3，直至收敛性条件得到满足为止。

所述的概率自适应花授粉优化模块13，采用如下步骤实现：

步骤1)：定义加入惩罚项后的增广性能指标函数，形式为：

其中，t表示积分时刻，t_f表示终值时刻，X(t)表示加入x_n+1(t)后的新状态向量，Φ₀[X(t_f)] 表示新状态向量在终值时刻处的终值函数，表示n_u维控制 向量，表示n_x维状态向量，ρ为惩罚因子，l表示光滑化 函数参数。H(t,x(t),u(t),ρ,l,ε)为约束相关的惩罚项，具体形式为：

其中，Q_i(t,x(t),u(t),ρ,l,ε)为表征约束违背量的函数。θ(Q_i(·))为多阶段配置函数，m₁表 示标准等式约束的个数，m₂表示标准不等式约束的个数。ε为光滑化因子，其取值为一个很 小的整数。

γ(Q_i(·))为惩罚函数指数，按以下方式更新：

θ(Q_i(·))为多阶段配置函数，定义如下：

基于以上方式，惩罚项因子的值随约束违背量的大小自适应地变化。当约束违背量较大 时，采取较大的惩罚因子，将更多的候选解推至可行区域内；当约束为违背量逐渐减小时， 惩罚因子也随之减小，保证进化算法在可行最优解周围搜索，更利于发现可行域边界周围的全局最优值。

步骤2)：按如下公式获得个体i与其他个体的平均欧式距离为：

其中，N表示种群规模，u_i、u_j表示种群中的两个个体。

步骤3)：定义密集度指数ρ表征种群密集度，按如下公式获得：

其中，d_g为当前迭代下全局最优个体与其他个体间的平均距离，d_min、d_max为d_i(i＝1,2,…,N) 的最大和最小值。

步骤4)：按如下公式获得转换概率p为：

步骤5)：由随机数生成器产生均匀分布随机数

步骤6)：根据与p的相对大小更新个体的位置。如果概率p大于随机数，个体依据 式(28)按照Lévy飞行方式更新；概率p小于随机数时，由下式计算更新后位置；

其中，gbest表示花粉中的最优个体，L(λ)是与花粉密度λ相关的Lévy飞行步长，可由Lévy 分布求得：

其中，Γ(λ)表示标准γ函数，s为步长。

最后，飞行器MCU 2将获得的优化轨迹作为指令输出到控制指令输出模块14，转换为 控制指令发送给攻角控制器8，完成轨迹优化的执行。所得到的飞行器状态轨迹如图3的(a) -(d)所示

其中，实施例1所得到的攻角控制策略曲线如图3所示。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明 的具体实施只限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离发明 构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (1)

1.一种概率自适应花授粉优化的高超声速飞行器再入段轨迹优化控制器，根据高超声速飞行器再入段初始海拔高度、速度、飞行航道倾角和飞行水平距离状态获取轨迹优化控制策略，通过控制飞行器攻角使高超声速飞行器获得更长的水平飞行距离；其特征在于：所述的轨迹优化控制器由飞行器海拔高度传感器、飞行器飞行速度传感器、飞行器飞行航道倾角传感器、飞行器水平飞行距离传感器、飞行器MCU、飞行器攻角控制器构成，所述飞行器海拔高度传感器、飞行器飞行速度传感器、飞行器飞行航道倾角传感器、飞行器水平飞行距离传感器、飞行器攻角控制器均通过数据总线与飞行器MCU相连；所述控制器的运行过程如下：

步骤1)：在高超声速飞行器MCU中输入对应于该飞行器的气动系数模型、飞行器性能约束条件、指定优化目标；

步骤2)：高超声速飞行器到达再入段后，开启飞行器海拔高度传感器、飞行器飞行速度传感器、飞行器飞行航道倾角传感器和飞行器水平飞行距离传感器，得到高超声速飞行器当前的海拔高度、速度、飞行航道倾角和飞行水平距离状态信息；

步骤3)：飞行器MCU根据设定的海拔高度、速度、飞行航道倾角要求自动执行内部概率自适应花授粉优化算法，得到使高超声速飞行器水平飞行距离最长的轨迹优化控制策略；

步骤4)：飞行器MCU将所获得的轨迹优化控制策略发送给控制指令输出模块，并转换为控制指令发送给飞行器攻角控制器执行；

所述的高超声速飞行器MCU，包括信息采集模块、初始化模块、系统状态求解模块、概率自适应花授粉优化模块、收敛性判断模块、控制指令输出模块；其中，信息采集模块包括飞行器海拔高度和速度采集子模块、飞行器飞行航道倾角和飞行水平距离采集子模块、飞行器海拔高度和速度设定采集子模块、飞行器飞行航道倾角设定采集子模块、飞行器的气动系数模型和性能约束条件以及指定优化目标参数采集子模块；

所述的高超声速飞行器MCU自动产生攻角控制指令的概率自适应花授粉优化算法运行步骤如下：

步骤1)：高超声速飞行器MCU到达再入段后，飞行器海拔高度传感器、飞行器飞行速度传感器、飞行器飞行航道倾角传感器和飞行器水平飞行距离传感器开启，信息采集模块获取高超声速飞行器当前的海拔高度、速度、飞行航道倾角和飞行水平距离状态信息；

步骤2)：初始化模块开始运行，对各项参数初始化，包括控制变量参数化参数和概率自适应花授粉优化方法参数；控制变量参数化参数的初始化设置如下：设置时间段的分段数为NE，设置控制变量参数化的插值拟合方式：分段常数或者分段线性；概率自适应花授粉优化方法参数的初始化设置如下：设置种群规模N，种群中第i个个体的初始位置向量种群初始历史最优位置其中，J(·)表示要最大化的目标函数；设置个体寻优最大迭代次数max_iters，花粉密度λ，初始迭代次数k＝0；

步骤3)：通过系统状态求解模块获得第k次迭代种群个体i的目标函数值更新第k次迭代种群历史最优位置g^(k),^best，记录最小目标函数值

步骤4)：执行收敛性判断模块，如果当前迭代次数k达到最大迭代次数max_iters，则收敛性条件得到满足，当前种群历史最优位置g^(k),^best作为最优控制策略的离散解，经控制指令输出模块输出；如果收敛性条件不满足，则迭代次数k加1；

步骤5)：执行概率自适应花授粉优化模块，依据种群分布状况计算表征个体聚集程度的指标，并计算相应的转换概率，由随机数生成器产生均匀分布随机数根据随机数与转换概率p的相对大小更新个体的位置；该步骤执行完毕后跳转至步骤3，直至收敛性条件得到满足为止；

其中，所述的系统状态求解模块，采用的是四步四阶Runge-Kutta方法，计算公式为：

其中，t_m表示四步四阶Runge-Kutta方法选择的积分时刻，t_m+1表示位于时刻t_m后的积分时刻，并且t_m+1＝t_m+h，h为积分步长，F(·)是描述状态微分方程的函数，K1、K2、K3、K4分别表示Runge-Kutta法积分过程中的4个节点的函数值，x(·)是描述系统状态变量的函数，且x(t_m)表示时刻t_m的系统状态变量，x(t_m+1)表示时刻t_m+1的系统状态变量；

所述的概率自适应花授粉优化模块，采用如下步骤实现：

步骤1)：定义加入惩罚项后的增广性能指标函数形式为：

其中，t表示积分时刻，t_f表示终值时刻，X(t)表示加入x_n+1(t)后的新状态向量，Φ₀[X(t_f)]表示新状态向量在终值时刻处的终值函数，表示n_u维控制向量，表示n_x维状态向量，ρ为惩罚因子，l表示光滑化函数参数；H(t,x(t),u(t),ρ,l)为约束相关的惩罚项，具体形式为：

其中，Q_i(t,x(t),u(t),ρ,l)为表征约束违背量的函数，θ(Q_i(·))为多阶段配置函数，m₁表示标准等式约束的个数，m₂表示标准不等式约束的个数；ε为光滑化因子，其取值为一个很小的整数；

γ(Q_i(·))为惩罚函数指数，按以下方式更新：

θ(Q_i(·))为多阶段配置函数，定义如下：

基于以上方式，惩罚项因子的值随约束违背量的大小自适应地变化；当约束违背量较大时，采取较大的惩罚因子，将更多的候选解推至可行区域内；当约束为违背量逐渐减小时，惩罚因子也随之减小，保证进化算法在可行最优解周围搜索，更利于发现可行域边界周围的全局最优值；

步骤2)：按如下公式获得个体i(i＝1,2,…,N)与其他个体的平均欧式距离d_i为：

其中，N表示种群规模，u_i、u_j表示种群中的两个个体；

步骤3)：定义密集度指数ρ表征种群密集度，按如下公式获得：

其中，d_g为当前迭代下全局最优个体与其他个体间的平均距离，d_min、d_max为d_i(i＝1,2,…,N)的最大和最小值；

步骤4)：按如下公式获得转换概率p为：

步骤5)：由随机数生成器产生均匀分布随机数

步骤6)：根据与p的相对大小更新个体的位置；如果概率p大于随机数个体依据式(5)按照Lévy飞行方式更新；概率p小于随机数时，由下式计算更新后位置；

其中，gbest表示花粉中的最优个体，L(λ)是与花粉密度λ相关的Lévy飞行步长，可由Lévy分布求得：

其中，Γ(λ)表示标准γ函数，s为步长，s₀为初始步长。