本文对LBAS算法的改进

嗅探前随机化 $d_0$ 的改进

在方向向量生成后、计算须位置前，引入随机化因子：

其中 $r \sim U(0,1)$ 是随机数。这一改进实际上将固定长度的方向向量变为随机长度的方向向量，使得天牛须的”嗅探”行为更加多样化。

二次进化距离的随机化改进

引入随机因子 $r \sim U(0,1)$，使二次进化的步长具有随机性：

随机因子 $r$ 的引入使得二次进化步骤具有了自适应步长调整的能力，平衡了探索与开发，针对本问题能够特别有效地搜索变量空间。

图5.2 展示了改进的LBAS算法的二次进化过程，其中灰色虚线为触角、红线为天牛个体主动移动，五角星为最优记录，蓝色点线为更优引导，橙色点划线为最优引导，棕线为最终二次进化路线。

模型求解

LBAS算法参数设置

针对本问题，LBAS算法的参数设置如下：

优化结果

使用改进的LBAS算法求解问题二，得到最优参数组合：

最优投放点坐标：$(17912.26, 9.81, 1800.00)\text{m}$
最优起爆点坐标：$(17929.22, 11.29, 1799.91)\text{m}$

六、问题三的模型建立与求解

优化模型建立

问题三要求在单无人机FY1上投放3枚烟幕干扰弹，实施对导弹M1的干扰，目标是最大化总有效遮蔽时间。这是一个8维优化问题，决策变量包括：

1. 无人机飞行方向 $\theta$（方位角）
2. 无人机飞行速度 $v$（m/s）
3. 第一枚烟幕弹投放时间 $T_{L1}$（s）c
4. 第一枚烟幕弹引信时间 $t_{F1}$（s）
5. 第二枚烟幕弹相对投放间隔 $\Delta T_{L2}$（s）
6. 第二枚烟幕弹引信时间 $t_{F2}$（s）
7. 第三枚烟幕弹相对投放间隔 $\Delta T_{L3}$（s）
8. 第三枚烟幕弹引信时间 $t_{F3}$（s）

其中 $T_{L2} = T_{L1} + 1 + \Delta T_{L2}$，$T_{L3} = T_{L2} + 1 + \Delta T_{L3}$。$I_{\text{shielding}}$为求区间函数。$\sum$为重叠区间求和。

基于对FY1位置和导弹轨迹的分析，我们提出以下建模策略：

1. 由于FY1位置相对靠后，导弹会快速越过无人机，因此烟雾弹必须在短时间内迅速投放。将投放时间和引信时间上限设为5秒，大幅缩减搜索空间。

2. 通过将投放间隔表示为”最小间隔+增量”的形式处理投放间隔约束：

   $$T_{L2} = T_{L1} + 1 + \Delta T_{L2}, \quad T_{L3} = T_{L2} + 1 + \Delta T_{L3}$$

   其中 $\Delta T_{L2}, \Delta T_{L3} \geq 0$，确保满足至少1秒间隔的要求。

3. 两阶段优化策略：

   • 第一阶段：添加无效烟雾弹惩罚，引导算法找到三个烟雾弹都有效的解区域
   • 第二阶段：取消惩罚，在可行解附近进行精细优化，找到全局最优解

模型求解

针对8维高维优化问题，对LBAS算法进行以下改进：

1. 增大种群规模：从30增加到300，增强全局探索能力
2. $d_{\eta} = \delta_{\eta} = 0.97$，减缓参数衰减速度
3. $P = 0.6$，增强学习能力
4. 方位角搜索范围扩大到 $[0, \frac{5\pi}{2}]$，避免最优空间在角度上割裂。

通过两阶段LBAS优化，得到的最优结果为：

总遮蔽时间：6.403160

七、问题四的模型建立与求解

模型分析与猜想

问题四要求利用FY1、FY2、FY3三架无人机各投放1枚烟幕干扰弹，实施对导弹M1的协同干扰。通过对三架无人机初始位置的空间分析：

• FY1位置：(17800, 0, 1800) m
• FY2位置：(12000, 1400, 1400) m
• FY3位置：(6000, -3000, 700) m
• M1位置：(20000, 0, 2000) m
• 假目标位置：(0, 0, 0) m
• 真目标位置：(0, 200, 0) m

基于三架无人机的初始位置分析，三架无人机分布在导弹飞行路径的不同区段，且彼此间距较大。基于此，可以提出以下猜想：

猜想7.1：由于三架无人机空间分布差异显著，各自的最优遮蔽策略相互独立，协同优化可分解为三个独立的子优化问题。

猜想7.2：三架无人机的最优投放时间存在先后顺序，与各自空间位置相关，形成时序不连续的遮蔽效果。预测三架无人机的遮蔽时间段不会重合，每架无人机负责导弹飞行路径的不同区段，形成分段遮蔽策略。

增强奖励函数

与问题二不同，该问题中F2,F3无人机的投放和引信时间变量空间更大，无法限制范围，解集更加稀疏。标准遮蔽时间函数 $T_{\text{shielding}}$ 在优化过程中存在明显缺陷：当烟雾弹完全未遮蔽导弹时，函数值为零，无法提供梯度信息，导致优化算法在平坦区域陷入停滞。为此可以设计一个增强的奖励函数，对于完全未遮蔽的情况下给予额外惩罚。为此，我们设计增强奖励函数，其数学表达式为：

其中：

• $t_m = T_L + t_F + \dfrac{T_s}{2}$ 为评估中点时刻
• $d_{\text{smoke-to-route}} = \dfrac{\|(\mathbf{S}(t_m) - \mathbf{C}_0) \times (\mathbf{M}(t_m) - \mathbf{C}_0)\|}{\|\mathbf{M}(t_m) - \mathbf{C}_0\|}$ 为烟雾云到导弹航线的距离
• $\mathbf{C}_0 = (0, 200, 5)$ 为目标中心点
• $\gamma = \begin{cases} 2 & \text{if } \|\mathbf{M}(t_m) - \mathbf{C}_0\| < \|\mathbf{S}(t_m) - \mathbf{C}_0\| \\ 1 & \text{otherwise} \end{cases}$ 为导弹不及时惩罚因子

增强奖励函数在未遮蔽区域仍能提供有意义的梯度信息，引导算法向有效解方向搜索。同时针对不及时造成的惩罚加倍，引导算法使无人机朝更早方向探索。

优化模型建立与求解

基于空间独立性猜想，我们将原问题分解为三个子问题：

其中 $R_i$ 为第 $i$ 架无人机的增强奖励函数。特殊地，在猜想之上，关于无人机FY1的结果可以直接沿用问题二的结果。

采用增强奖励函数和分解协调策略，得到的最优结果为：

总遮蔽时间：11.6792 s

猜想验证与讨论

验证猜想7.1：三架无人机的最优参数差异显著，且各自遮蔽时间接近独立优化时的最大值，证实了空间独立性猜想的正确性。无人机间的最优策略确实相互独立，协同效应主要体现在有效遮蔽时间段的分离上。

验证猜想7.2：三架无人机的投放时间呈现明显的时间序列：FY1(0.91s) → FY2(10.50s) → FY3(29.91s)，这与各自空间位置吻合：距离导弹起点最近的FY1最先投放，距离最远的FY3最后投放。

八、问题五的模型建立与求解

问题描述与建模思路

问题复杂性分析

问题五要求利用5架无人机（FY1-FY5）各投放至多3枚烟幕干扰弹，实施对3枚来袭导弹（M1-M3）的协同干扰。这是一个高度复杂的优化问题。每架无人机有8个决策变量（方位角、速度、3枚烟雾弹的投放时间和引信时间），5架无人机共40个决策变量。需要同时优化对3枚导弹的遮蔽效果，目标间存在竞争关系。且存在烟雾弹投放给不适宜的导弹导致无法达到最优解（如第3或第2、3枚烟雾弹屏蔽到第1或第1、2个导弹的情形）

图8.1中，M1M3为导弹位置，F1F5为无人机位置，T为真目标位置。红线表示烟雾弹需要遮蔽的导弹探测视线。

分层优化策略

基于问题复杂性，我们提出分层优化策略：

第一层：单无人机多导弹优化（问题五子问题）

• 扩展至单无人机对多导弹遮蔽
• 采用”一弹一目标”分配策略，避免“烟雾弹错位遮蔽”
• 使用问题四优化的奖励函数，惩罚烟雾弹未有效遮蔽的情形，辅助快速收敛到解集密集区。

第二层：多无人机多导弹优化

• 整合各无人机优化结果
• 全局协调微调

资源分配策略

基于各无人机位置特征，制定如下资源分配方案：

该分配方案基于以下几何考虑：

• FY1位置最靠前，难以在M2和M3之前再布置烟雾弹。故可专注于遮蔽路线最近的M1
• FY2、FY4位置偏北，优先遮蔽M2和M1
• FY3、FY5位置偏南，优先遮蔽M3和M1

优化模型建立

对于每架无人机，定义其目标函数为对分配导弹的总遮蔽奖励：

全局目标函数为各无人机目标函数的和：

其中：

• $M_{ij}$ 为第$i$架无人机第$j$枚烟雾弹分配的导弹目标
• $R$ 为增强的奖励函数。采用问题四中设计的增强奖励函数：

两阶段优化流程

第一阶段：增强奖励函数优化

• 使用增强奖励函数快速找到解集密集区
• 获得初始解 $\mathbf{x}^{(0)}$

第二阶段：标准目标函数优化

• 使用标准计算重复时间区间的方法进行精细优化
• 以 $\mathbf{x}^{(0)}$ 为初始点，寻找全局最优解

优化结果与分析

结论与启示

通过问题五的建模与求解，我们得到以下结论：

1. 改进的LBAS算法在解决此类高维、非线性、多约束的几何优化问题上展现出显著优势。
2. 如图8.2所示，本文提出的两阶段优化流程在实际应用中表现出色。针对生效时机相近的F2、F4、F5在M2导弹上的屏蔽展现了明显的优化偏移痕迹。
3. 通过三维运动学分析发现，F2和F4无法对第三枚导弹形成有效遮蔽并非算法局限，而是物理约束下的最优选择。以F2为例，如图8.3所示：

图8.3中，红色为导弹位置与轨迹，蓝色为无人机位置与轨迹，黑色为烟雾弹位置及轨迹，绿色为导弹的探测视线。此图中，烟雾弹1和2均在导弹越过无人机之前被遮挡，而当烟雾弹3到达其遮挡的导弹之前、遮挡已经结束。F2的速度已达到140m/s的上限，但在此速度下仍然无法有效地在到达M1之前使烟雾弹达到M1视线附近。故抛弃对此导弹屏蔽在全局范围内是具有策略性的

九、模型优缺点

模型的优点

1. 本研究提出的改进LBAS算法在多个方面展现创新性：

• 种群机制改进：引入排序分组和角色划分，将单一天牛搜索扩展为群体协同搜索，全局搜索能力大幅提升。
• 二次进化策略：通过向更优个体和引领者学习，有效在控制迭代次数总是不高于300次的情况下找到全局最优。
• 自适应扰动机制：有效避免早熟收敛，在多种算法中收敛表现优异。

2. 模型采用两阶段分解协调策略，将40维问题降维处理，实时性能表现优异。

模型的缺点

1. 模型基于以下简化假设，在实际应用中可能存在偏差：

• 忽略空气阻力：烟幕弹运动采用理想抛体模型，实际弹道可能因空气阻力产生偏移，结合研究中常见的滞空时间（约5s），使用优化的烟幕弹弹身设计，预计最大偏差可能达到8-12%

图9.1展示了大型烟幕弹(m=10kg,C=0.02,v0=100.0m/s)的平抛运动下是否加入空气阻力的理论影响。其中理论为有空气阻力，理想为无空气阻力。

2. 图8.3展示了无人机的部分烟幕弹无法到达部分导弹，但未尝试使用这枚烟幕弹屏蔽其他导弹的情形。

十一、参考文献

[1] Jiang, X., & Li, S. (2018). BAS: beetle antennae search algorithm for optimization problems. International Journal of Robotics and Control.

[2] 李辉,殷文明.引领者扰动天牛须算法及其应用J.数学建模及其应用,2025,14(02)

[3] 廖列法,杨红. 天牛须搜索算法研究综述 [J]. 计算机工程与应用, 2021, 57 (12): 54-64.

附录

支撑材料列表

LBAS算法实现代码

"""
算法来自论文
李辉,殷文明.引领者扰动天牛须算法及其应用[J].数学建模及其应用,2025,14(02):20-27.
有修改
"""

import time
import kawaiitb.autoload

from matplotlib.animation import FuncAnimation
import numpy as np
import matplotlib.pyplot as plt
from tqdm import trange

plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False

feature_config = {
    # 是否禁用在特定迭代次数后开始扰动最优点
	'disable_turbulent': False,

    # 是否禁用在未二次进化时重新roll点
    # True时更偏向每个个体自己的探索。让每个个体自行在P概率下移动到最优。
	'disable_regenerate': False,  # 30维球测试条件可证明禁用更好，高维不敏感

    # 是否只通过引领者进行进化。
    # 待测试：使个体要么往最优方向，要么往上位方向。指向性最好强一点
	'only_forward_best': False,  # 似乎对结果不明显。
    
    # 手动复现论文的公式(6)的bug
	'recurrent_bug_from_raw_paper': False,  # 任何情况禁用都更好。

    # 二次进化的距离是否使用随机数。论文没提但意思应该是那样的
    # 使用固定值会导致随机性很低的向优趋势。
    'secondary_evolution_use_random': True,

    # 是否在“嗅探”前就立刻地随机化d0。论文里公式这里不这么做，但退火时这样是反直觉且梯度强假设的
    'randomly_tester_d0': True,
}

# 引领者扰动天牛须算法 (LBAS) 实现
class LBAS:
    def __init__(self, obj_func, pop_size, bounds, max_iter, init_x=None,
                 seed=None,

                 # 超参数
                 P=None,  # 更大时有更大机会二次进化，朝更优者和引领者方向移动，否则重生。
                 d0=None, 
                 delta0=None, 
                 d0_init_ratio=0.5, 
                 delta0_init_ratio=0.5, 
                 d_eta=0.95, 
                 delta_eta=0.95, 
                 tol=1e-10, 

                 with_plot=False, 
                 tqdm=False, 
                 history_per_iter: int|None = 100):
        # 优化算法参数
        self.obj_func = obj_func  # 优化目标
        bounds = np.array(bounds)
        self.dim = len(bounds)  # 问题维度
        self.pop_size = pop_size  # 种群大小
        self.bounds = np.array(bounds)  # 搜索边界 [[min, max], ...]
        self.availd_domain = bounds[:, 1] - bounds[:, 0]
        self.max_iter = max_iter  # 最大迭代次数
        self.tol = tol  # 早停阈值
        self.best_history = []  # 记录每代最佳适应度
        self.seed = seed or int(time.time())
        self.rng = np.random.default_rng(self.seed)

        # 可视化选项和历史
        self.tqdm = tqdm
        self.history_per_iter = history_per_iter
        self.with_plot = with_plot
        if with_plot:
            if dim != 2:
                raise ValueError("2维问题才能绘制动画")
            self.with_plot = True
            self.positions_history = []
            self.onframe_lines = []
        
        # 天牛须算法参数
        self.d0 = d0 if d0 is not None else (
            (self.availd_domain) * d0_init_ratio
        )  # 初始两须距离
        self.delta0 = delta0 if delta0 is not None else (
            (self.availd_domain) * delta0_init_ratio
        )  # 初始步长
        self.d_eta = d_eta  # 距离衰减系数
        self.delta_eta = delta_eta  # 步长衰减系数
        
        # LBAS特定参数
        self.P = P or 0.5  # 二次进化概率
        self.t0 = 30  # 初始阶段迭代次数（之后开始扰动）
        self.reducing = 1  # 用于记录衰减系数
        
        # 初始化种群
        self.population = self.init_population(init_x)
        self.fitness = np.array([self.obj_func(ind) for ind in self.population])
        
        # 记录最佳解
        self.best_idx = np.argmin(self.fitness)
        self.best_solution = self.population[self.best_idx].copy()
        self.best_fitness = self.fitness[self.best_idx]
        
        # 记录收敛过程
        self.convergence = []
        self.iteration = 0
    
    def init_population(self, init_x=None):
        """初始化种群"""
        
        population = self.rng.uniform(
            low=self.bounds[:, 0],
            high=self.bounds[:, 1],
            size=(self.pop_size, self.dim)
        )
        
        if init_x != None:
            if len(init_x) != self.dim:
                raise ValueError("init_x的维度必须与问题维度一致")
            population[0] = self.clip(init_x)
        
        return population

    def clip(self, x):
        """确保解在有效域内"""
        return np.clip(x, self.bounds[:, 0], self.bounds[:, 1])

    def iterate(self):
        """执行一次迭代"""
        thisframe_lines = []
        self.iteration += 1
        if self.with_plot:
            self.positions_history.append(self.population.copy())
        # 对种群按适应度排序
        sorted_indices = np.argsort(self.fitness)
        sorted_population = self.population[sorted_indices]
        sorted_fitness = self.fitness[sorted_indices]
        
        # 更新全局最优
        if sorted_fitness[0] < self.best_fitness:
            self.best_fitness = sorted_fitness[0]
            self.best_solution = sorted_population[0].copy()
        
        self.best_history.append((self.best_solution.copy(), self.best_fitness))
        
        # 对每个个体进行更新
        for i in range(self.pop_size):
            # 当前个体
            x = self.population[i].copy()
            
            # 1. 生成随机方向向量
            b = self.rng.uniform(-1, 1, self.dim)
            b = b / np.linalg.norm(b)  # 单位化
            if feature_config.get('randomly_tester_d0', False):
                b = b * self.rng.random()  # 随机化
            
            # 2. 计算左右须位置
            x_left = x - self.d0 * b
            x_right = x + self.d0 * b
            
            # 确保在边界内
            x_left = self.clip(x_left)
            x_right = self.clip(x_right)
            
            # 3. 计算适应度
            f_left = self.obj_func(x_left)
            f_right = self.obj_func(x_right)
            
            # 4. 决定移动方向并应用惯性权重
            inertia = self.delta0 * b
            # 如果在计算时已经随机化步长，那现在就不随机化了
            if not feature_config.get('randomly_tester_d0', False):
                inertia = inertia * self.rng.random()
            if f_right > f_left:  # 向左移动
                x_new = x - inertia
            else:  # 向右移动
                x_new = x + inertia
            if self.with_plot:
                thisframe_lines.append((x_left, x_right, {  # 天牛须
                    'color': 'black',
                    'linestyle': '--',
                    'linewidth': 1
                }))
                thisframe_lines.append((x, x_new, {  # 天牛主动行动
                    'color': 'red',
                    'linestyle': '-',
                    'linewidth': 4
                }))
            
            # 确保在边界内
            x_new = self.clip(x_new)
            f_new = self.obj_func(x_new)
            
            # 5. 按概率P进行二次进化
            idx_in_sorted = np.where(sorted_indices == i)[0][0]
            if idx_in_sorted > 0:  # 如果不是最优个体
                # 找到排在当前个体前面的个体
                if self.rng.random() < self.P:
                    # 获取前一个个体和全局最优个体
                    x_prev = sorted_population[idx_in_sorted - 1]
                    x_best = self.best_solution
                    
                    # 二次进化：向更优个体和全局最优个体学习
                    # 使用正确的公式，基于上一次的位置
                    dx_secondary = 0.5 * (x_prev - x_new) + 0.5 * (x_best - x_new)
                    if feature_config.get('only_forward_best', False):
                        dx_secondary = 0 * (x_prev - x_new) + 1 * (x_best - x_new)
                    
                    use_rand = self.rng.random() if feature_config.get('secondary_evolution_use_random', False) else 1
                    x_secondary = x_new + dx_secondary * use_rand
                    if feature_config.get('recurrent_bug_from_raw_paper', False):
                        x_secondary = dx_secondary * use_rand
                    
                    if self.with_plot:
                        thisframe_lines.append((x_new, x_prev, {  # 二次进化，向更优个体学习
                            'color': 'cyan',
                            'linestyle': ':',
                            'linewidth': 2
                        }))
                        thisframe_lines.append((x_new, x_best, {  # 二次进化，向全局最优学习
                            'color': 'salmon',
                            'linestyle': '-.',
                            'linewidth': 2
                        }))
                        thisframe_lines.append((x_new, x_new + dx_secondary, {  # 二次进化总方向
                            'color': 'brown',
                            'linestyle': '-',
                            'linewidth': 2
                        }))
                        thisframe_lines.append((x_new, x_secondary, {  # 完整二次进化
                            'color': 'brown',
                            'linestyle': '-',
                            'linewidth': 4
                        }))
                    
                    x_secondary = self.clip(x_secondary)
                    f_secondary = self.obj_func(x_secondary)
                    
                    # 选择较好的解
                    if f_secondary < f_new:
                        x_new = x_secondary
                        f_new = f_secondary
                elif not feature_config.get('disable_regenerate', False):
                    # 未进行二次进化时，生成随机新个体
                    # 从最优个体开始往随机方向前进当前的两须长度再带个随机量
                    random_direction = self.rng.uniform(-1, 1, self.dim)
                    # 归一化方向向量
                    norm = np.linalg.norm(random_direction)
                    if norm > 0:
                        random_direction = random_direction / norm
                    
                    # 前进距离
                    step_length = self.delta0 * self.rng.random()  # 随机量设为0到d之间的值
                    
                    # 生成新个体
                    x_random = self.best_solution + random_direction * step_length
                    x_random = self.clip(x_random)
                    f_random = self.obj_func(x_random)
                    
                    # 如果新个体更优则接受
                    if f_random < f_new:
                        x_new = x_random
                        f_new = f_random
            

            # 6. 更新个体
            self.population[i] = x_new
            self.fitness[i] = f_new
        
        # 7. 对最优个体进行扰动（如果超过初始阶段）
        if (not feature_config.get('disable_turbulent', False)) and self.iteration > self.t0:
            # 对最优个体添加扰动
            perturbation = self.rng.uniform(-1, 1, self.dim) * self.availd_domain * np.exp(-4 * self.iteration / self.max_iter)
            x_perturbed = self.best_solution + perturbation
            x_perturbed = self.clip(x_perturbed)
            f_perturbed = self.obj_func(x_perturbed)
            
            # 如果扰动后的解更好，则替换最优解
            if f_perturbed < self.best_fitness:
                self.best_solution = x_perturbed
                self.best_fitness = f_perturbed
                
                # 也用扰动后的解替换种群中最差的个体
                worst_idx = np.argmax(self.fitness)
                self.population[worst_idx] = x_perturbed
                self.fitness[worst_idx] = f_perturbed
        
        # 早停
        # if self.best_fitness < self.tol:
        #     return True
        
        # 8. 更新步长和两须距离
        self.reducing = self.reducing * self.d_eta
        self.d0 = self.d_eta * self.d0
        self.delta0 = self.delta_eta * self.delta0
        
        # 记录当前最优适应度
        self.convergence.append(self.best_fitness)
        if self.with_plot:
            self.onframe_lines.append(thisframe_lines)
        return False
    
    def optimize(self):
        """执行优化过程"""
        if self.with_plot:
            self.positions_history = []
        self.convergence = []
        iters = trange if self.tqdm else range
        for self.iteration in iters(self.max_iter):
            if self.iterate():
                break
            if self.history_per_iter is not None and self.iteration % self.history_per_iter == 0:
                print(f"I{self.iteration:<6} | Fit={self.best_fitness:.4e}, rdc={self.reducing:.2e}")
        
        return self.best_solution, self.best_fitness, self.convergence
        
    def create_contour_plot(self):
        """创建等高线图"""
        x = np.linspace(self.bounds[0][0], self.bounds[0][1], 100)
        y = np.linspace(self.bounds[1][0], self.bounds[1][1], 100)
        X, Y = np.meshgrid(x, y)
        Z = np.zeros_like(X)
        
        for i in range(X.shape[0]):
            for j in range(X.shape[1]):
                Z[i, j] = self.obj_func(np.array([X[i, j], Y[i, j]]))
        
        plt.figure(figsize=(12, 10))
        contour = plt.contour(X, Y, Z, 50, cmap='viridis', alpha=0.6)
        plt.colorbar(contour, label='Function Value')
        plt.xlabel('X')
        plt.ylabel('Y')
        plt.title('LBAS Optimization Process')
        plt.grid(True, alpha=0.3)
        
        return plt.gca(), X, Y, Z
    
    def animate_optimization(self):
        """创建优化过程动画"""
        fig, ax = plt.subplots(figsize=(12, 10))
        ax: plt.Axes  # type: ignore
        
        # 创建等高线
        x = np.linspace(self.bounds[0][0], self.bounds[0][1], 100)
        y = np.linspace(self.bounds[1][0], self.bounds[1][1], 100)
        X, Y = np.meshgrid(x, y)
        Z = np.zeros_like(X)
        
        for i in range(X.shape[0]):
            for j in range(X.shape[1]):
                Z[i, j] = self.obj_func(np.array([X[i, j], Y[i, j]]))
        
        contour = ax.contour(X, Y, Z, 30, cmap='viridis', alpha=0.6)
        plt.colorbar(contour, ax=ax, label='Function Value')
        
        # 初始化图形元素
        beetles_scatter = ax.scatter([], [], c='blue', s=30, label='Beetles')
        best_beetle_scatter = ax.scatter([], [], c='red', s=100, marker='*', label='Best Beetle')
        antennae_lines = [ax.plot([], [], 'k-', alpha=0.5)[0] for _ in range(self.pop_size)]
        prev_lines = 100
        move_lines = [ax.plot([], [], 'r-', alpha=0.5)[0] for _ in range(prev_lines)]
        
        # 添加文本框显示迭代信息
        iter_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12,
                           bbox=dict(facecolor='white', alpha=0.7))
        best_text = ax.text(0.02, 0.88, '', transform=ax.transAxes, fontsize=10,
                           bbox=dict(facecolor='white', alpha=0.7))
        
        ax.set_xlim(self.bounds[0])
        ax.set_ylim(self.bounds[1])
        ax.set_xlabel('X')
        ax.set_ylabel('Y')
        ax.set_title('Leader Disturbance Beetle Antennae Search Algorithm')
        ax.legend()
        ax.grid(True, alpha=0.3)

        # 动画更新函数
        def update(frame):
            if frame < len(self.positions_history) - 1:  # 最后一帧有位置但没动作，索性全都不更新
                positions = self.positions_history[frame]
                best_pos, best_fit = self.best_history[frame]
                going_arrows: list[tuple[np.ndarray, np.ndarray, dict]] = self.onframe_lines[frame]
                
                # 更新天牛位置
                beetles_scatter.set_offsets(positions)
                
                # 更新最优天牛位置
                best_beetle_scatter.set_offsets([best_pos])
                
                # 更新移动箭头
                # for i, arrow in enumerate(going_arrows):
                #     start, end = arrow
                #     move_lines[i].set_data([start[0], end[0]], [start[1], end[1]])
                for line, arrow in zip(move_lines, going_arrows + [((-9999, -9999), (-9999, -9999), {})] * prev_lines):
                    start, end, style = arrow
                    line.set_data([start[0], end[0]], [start[1], end[1]])
                    line.set_color(style.get('color', 'black'))
                    line.set_linestyle(style.get('linestyle', '-'))
                    line.set_linewidth(style.get('linewidth', 1))

                # 更新文本信息
                iter_text.set_text(f'Iteration: {frame+1}/{self.max_iter}')
                best_text.set_text(f'Best Fitness: {best_fit:.6f}\nBest Position: ({best_pos[0]:.4f}, {best_pos[1]:.4f})')
            
            return beetles_scatter, best_beetle_scatter, *antennae_lines, iter_text, best_text
        
        ani = FuncAnimation(fig, update, frames=len(self.positions_history), 
                          interval=200, blit=True)
        return ani
# 测试算法
if __name__ == "__main__":
    
    from pprint import pprint

    pprint(feature_config)
    def sphere_function(x):
        def function(x):
            return np.sum(x**2)
        return function, (-100, 100)
    # 参数设置
    center = 0
    test_function = sphere_function
    obj, func_bound = test_function(center)
    dim = 2  # 问题维度
    bounds = np.array([func_bound] * dim)  # 求解范围
    pop_size = 3  # 种群大小
    max_iter = 50  # 最大迭代次数
    
    # 创建并运行算法
    eta = 0.9
    lbas = LBAS(obj, pop_size, bounds, max_iter,
                d_eta=eta, delta_eta=eta, tol=1e-5, with_plot=True)
    best_solution, best_fitness, convergence = lbas.optimize()
    
    print()
    print(f"Test params: {dim=},{center=},{test_function.__name__=},{bounds=},")
    print(f"Algorithm params: {pop_size=},{max_iter=}")
    print(f"solution: {best_solution}")
    print(f"fitness: {best_fitness:.6e}")
    print(f"error: {np.linalg.norm(best_solution - center)}")
    print(f"iter: {lbas.iteration} / {lbas.max_iter}")
    
    # 绘制收敛曲线
    # plt.figure(figsize=(10, 6))
    # plt.semilogy(convergence)
    # plt.title('Convergence Curve')
    # plt.xlabel('Iteration')
    # plt.ylabel('Best Fitness (log scale)')
    # plt.grid(True)
    # plt.savefig(f'fig/fig{time.time()}.png')
    # plt.show()

    ani = lbas.animate_optimization()
    ani.save('lbas_optimization.gif', writer='pillow', fps=1)
    print("动画已保存为 'lbas_optimization.gif'")
    ani.event_source.stop()
    del ani
    plt.close()

问题一的模型代码

from functools import lru_cache
from math import sin, cos
import numpy as np
from warnings import warn

from utils import bool_func_singularity, calc_total_true_length, total_interval_sum

g = 9.8  # 重力加速度(m/s^2)
TARGET_POSITION = np.array([0, 200, 0])  # 真目标位置
GRENADE_MIN_INTERVAL = 1  # 烟雾弹最小间隔(s)

# 真目标检查点
TARGET_CHECKPOINTS = [
    TARGET_POSITION + np.array([7,  0,  0]),
    TARGET_POSITION + np.array([0,  7,  0]),
    TARGET_POSITION + np.array([-7, 0,  0]),
    TARGET_POSITION + np.array([0, -7,  0]),
    TARGET_POSITION + np.array([7,  0, 10]),
    TARGET_POSITION + np.array([0,  7, 10]),
    TARGET_POSITION + np.array([-7, 0, 10]),
    TARGET_POSITION + np.array([0, -7, 10]),
]

DECOY_POSITION = np.array([0, 0, 0])  # 假目标位置
MISSILE_VELOCITY = 300  # 导弹飞行速度(m/s)
SMOKE_SHIELDING_TIME = 20  # 烟雾云屏蔽时间(s)
SMOKE_VELOCITY = [0, 0, -3]  # 烟雾云速度(m/s)
SMOKE_RADIUS = 10  # 烟雾云半径(m)

def missile_coord(M, t):
    """
    计算导弹在t时刻的位置
    M: 导弹初始位置
    t: 时间(s)
    """
    M = np.array(M)
    assert M.shape == (3,), f"M的形状必须为(3,), 实际为{M.shape}"
    # 导弹飞向假目标(0,0,0)
    direction = DECOY_POSITION - M
    distance = np.linalg.norm(direction)
    unit_vector = direction / distance
    
    # 导弹以300m/s的速度飞行
    flight_distance = MISSILE_VELOCITY * t
    if flight_distance >= distance:
        warn("导弹到达假目标")
        return DECOY_POSITION
    else:
        return M + unit_vector * flight_distance

def uav_coord(FA, theta, v, t):
    """
    计算无人机在t时刻的位置
    FA: 无人机初始位置
    theta: 朝向
    v: 飞行速度(m/s)
    t: 时间(s)
    """
    # 无人机朝向azimuth飞行
    vector = np.array([cos(theta), sin(theta), 0])  # 模长为1
    
    # 无人机以恒定速度飞行
    return FA + vector * v * t

def grenade_coord(FA, theta, v, t_launch, t_fuse, t):
    """
    计算烟幕弹在t时刻的位置
    FA: 无人机初始位置
    azimuth: 朝向假目标方向
    speed: 飞行速度(m/s)
    t_launch: 投放时间(s)
    t_fuse: 引信时间(s)(从投放后算起)
    t: 时间(s)
    """
    # 投放时刻无人机的位置和速度
    launch_pos = uav_coord(FA, theta, v, t_launch)
    
    # 无人机速度向量
    uav_velocity = np.array([np.cos(theta), np.sin(theta), 0]) * v
    
    # 爆炸时刻
    t_burst = t_launch + t_fuse
    
    if t > t_burst:
        warn("计算烟幕弹位置时发现烟幕弹已爆炸")
        t = t_burst
    if t < t_launch:
        warn("计算烟幕弹位置时发现烟幕弹未投放")
        t = t_launch
    
    # 平抛运动
    delta_t = t - t_launch
    dH = uav_velocity * delta_t
    dZ = np.array([0, 0, -0.5 * g * delta_t**2])  # 负号在这
    return launch_pos + dH + dZ

def smoke_coord(FA, theta, v, t_launch, t_fuse, t):
    """
    计算烟雾云在t时刻的位置
    FA: 无人机初始位置
    azimuth: 朝向假目标方向
    speed: 飞行速度(m/s)
    t_launch: 投放时间(s)
    t_fuse: 引信时间(s)(从投放后算起)
    t: 时间(s)
    """
    t_burst = t_launch + t_fuse
    if t < t_burst:
        warn("计算烟雾云位置时发现烟雾弹未爆炸")
        t = t_burst
    if t > t_burst + SMOKE_SHIELDING_TIME:
        warn("计算烟雾云位置时发现烟雾弹已爆炸")
        t = t_burst + SMOKE_SHIELDING_TIME
    grenade_pos = grenade_coord(FA, theta, v, t_launch, t_fuse, t_launch + t_fuse)
    return grenade_pos + np.array(SMOKE_VELOCITY) * (t - t_burst)  # 烟雾云高度100m

def is_shielding(smoke_coord, missile_coord):
    """
    判断导弹是否被烟雾云遮蔽
    1. 判断导弹是否在烟雾云内部，返回True
    2. 判断导弹是否比烟雾云更快，返回False
    3. 判断导弹与目标的直线与烟雾云的距离是否大于烟雾云半径，返回False
    返回True
    smoke_coord: 烟雾云位置
    missile_coord: 导弹位置
    """
    # 1. 判断导弹在烟雾云内部
    if np.linalg.norm(missile_coord - smoke_coord) <= SMOKE_RADIUS:
        return True
    
    for checkpoint in TARGET_CHECKPOINTS:
        ms2cp = missile_coord - checkpoint
        cp2sc = checkpoint - smoke_coord
        # 2. 判断导弹是否比烟雾云更快
        distance_to_checkpoint = np.linalg.norm(ms2cp)
        distance_to_smoke = np.linalg.norm(cp2sc)
        if distance_to_checkpoint < distance_to_smoke:
            return False
        # 3. 判断导弹与目标的直线与烟雾云的距离是否大于烟雾云半径
        d_smoke_to_route = np.linalg.norm(np.cross(cp2sc, ms2cp)) / distance_to_checkpoint
        if d_smoke_to_route > SMOKE_RADIUS:
            return False
    
    return True

def shielding_para(FA, M, azimuth, speed, launch_time, fuse_time):
    """
    计算有效遮蔽时间
    FA: 无人机初始位置
    M: 导弹初始位置
    azimuth: 朝向假目标方向
    speed: 飞行速度(m/s)
    launch_time: 投放时间(s)
    fuse_time: 引信时间(s)
    """
    # 计算有效遮蔽时间
    t_burst = launch_time + fuse_time
    interval = (t_burst, t_burst + SMOKE_SHIELDING_TIME)
    shielding_paras = bool_func_singularity(
        lambda t: is_shielding(
            smoke_coord(FA, azimuth, speed, launch_time, fuse_time, t), 
            missile_coord(M, t)
        ),
        interval, tol=1e-6, max_depth=6
    )
    return shielding_paras

def para_start_end(para):
    """
    计算有效遮蔽时间
    """
    leftval, n_point, points, interval = para
    if leftval and n_point == 0:
        return interval[0], interval[1]
    if not leftval and n_point == 0:
        return interval[0], interval[0]
    elif leftval and n_point == 1:
        return interval[0], points[0]
    elif not leftval and n_point == 1:
        return points[0], interval[1]
    elif leftval and n_point == 2:
        raise ValueError("不预期的两次遮挡")
    elif not leftval and n_point == 2:
        return points[0], points[1]
    else:
        raise ValueError(f"不预期的参数 {para}")

def shielding_time(FA, M, azimuth, speed, launch_time, fuse_time):
    """
    计算有效遮蔽时间
    FA: 无人机初始位置
    M: 导弹初始位置
    azimuth: 朝向假目标方向
    speed: 飞行速度(m/s)
    launch_time: 投放时间(s)
    fuse_time: 引信时间(s)
    """
    t_shielding = calc_total_true_length(*shielding_para(FA, M, azimuth, speed, launch_time, fuse_time))
    return t_shielding

def sum_of_shielding_times(*para_groups):
    """
    计算总有效遮蔽时间
    """
    intervals = [para_start_end(para) for para in para_groups]
    return total_interval_sum(intervals)

def shielding_reward(FA, M, azimuth, speed, launch_time, fuse_time, assert_no_shielding=False):
    """
    计算导弹被烟雾云遮蔽的奖励。**此函数为最大化函数**
    1. 判断导弹是否实际被烟雾云遮蔽，返回原始屏蔽时长作为奖励
    获取烟雾弹爆炸开始至消散的平均位置与导弹当时的平均位置
    2. 计算导弹航线与烟雾云的距离，返回基于此距离的惩罚。
    3. 计算是否是导弹更靠近检查点，如果是则遮蔽不到，惩罚加倍
    """
    # 1. 判断是否被烟雾云遮蔽
    if not assert_no_shielding:
        T_shielding = shielding_time(FA, M, azimuth, speed, launch_time, fuse_time)
        if T_shielding > 0:
            return T_shielding
    
    smoke_pos = smoke_coord(FA, azimuth, speed, launch_time, fuse_time, launch_time + fuse_time + SMOKE_SHIELDING_TIME / 2)
    missile_pos = missile_coord(M, launch_time + fuse_time + SMOKE_SHIELDING_TIME / 2)

    target_center = TARGET_POSITION + np.array([0,  0,  5])
    cp2ms = missile_pos - target_center
    sk2cp = target_center - smoke_pos
    # 2. 计算导弹航线与烟雾云的距离
    d_missile_target = np.linalg.norm(cp2ms)
    d_target_smoke = np.linalg.norm(sk2cp)
    d_smoke_to_route = np.linalg.norm(np.cross(sk2cp, cp2ms)) / d_missile_target
    penalty = -d_smoke_to_route
    # 3. 判断导弹是否比烟雾云更快
    if d_missile_target < d_target_smoke:
        return penalty * 2
    return penalty

# 问题1参数
FA = np.array([17800, 0, 1800])  # FY1初始位置
M1 = np.array([20000, 0, 2000])  # M1初始位置
azimuth = np.arctan2(DECOY_POSITION[1] - FA[1], DECOY_POSITION[0] - FA[0])
speed = 120  # 飞行速度(m/s)
launch_time = 1.5  # 投放时间(s)
fuse_time = 3.6  # 引信时间(s)

if __name__ == "__main__":
    # 解决问题1
    print("问题一计算结果：")
    print(f"飞行方向：{np.degrees(azimuth)}°")
    print(f"飞行速度：{speed} m/s")
    print(f"投放时间：{launch_time} s")
    print(f"引信时间：{fuse_time} s")
    effective_time = shielding_time(FA, M1, azimuth, speed, launch_time, fuse_time)
    print(f"遮蔽时间: {effective_time:.6f}秒")



if __name__ == "__main__1":  # 测试函数评估速度，即可优化性
    import timeit
    repeats = 5
    number_per_repeat = 20
    def test_shielding_time():
        return shielding_time(FA, M1, azimuth, speed, launch_time, fuse_time)
    
    times = timeit.repeat(test_shielding_time, number=number_per_repeat, repeat=repeats)
    
    print(f"重复测试 {repeats} 次，每次 {number_per_repeat} 轮")
    for i, time_taken in enumerate(times, 1):
        average_time = time_taken / number_per_repeat
        print(f"第 {i} 次测试 - 平均时间: {average_time:.6f} 秒")
    
    best_time = min(times) / number_per_repeat
    worst_time = max(times) / number_per_repeat
    avg_time = sum(times) / (repeats * number_per_repeat)
    
    print(f"\n最佳执行时间: {best_time:.6f} 秒")
    print(f"最差执行时间: {worst_time:.6f} 秒")
    print(f"平均执行时间: {avg_time:.6f} 秒")
    
    # 运行一次获取结果
    result = test_shielding_time()
    print(f"有效遮蔽时长: {result:.6f} 秒")