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

7
import time
8
import kawaiitb.autoload
9

10
from matplotlib.animation import FuncAnimation
11
import numpy as np
12
import matplotlib.pyplot as plt
13
from tqdm import trange
14

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

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

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

26
    # 是否只通过引领者进行进化。
27
    # 待测试：使个体要么往最优方向，要么往上位方向。指向性最好强一点
28
  'only_forward_best': False,  # 似乎对结果不明显。
29

30
    # 手动复现论文的公式(6)的bug
31
  'recurrent_bug_from_raw_paper': False,  # 任何情况禁用都更好。
32

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

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

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

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

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

72
        # 可视化选项和历史
73
        self.tqdm = tqdm
74
        self.history_per_iter = history_per_iter
75
        self.with_plot = with_plot
76
        if with_plot:
77
            if dim != 2:
78
                raise ValueError("2维问题才能绘制动画")
79
            self.with_plot = True
80
            self.positions_history = []
81
            self.onframe_lines = []
82

83
        # 天牛须算法参数
84
        self.d0 = d0 if d0 is not None else (
85
            (self.availd_domain) * d0_init_ratio
86
        )  # 初始两须距离
87
        self.delta0 = delta0 if delta0 is not None else (
88
            (self.availd_domain) * delta0_init_ratio
89
        )  # 初始步长
90
        self.d_eta = d_eta  # 距离衰减系数
91
        self.delta_eta = delta_eta  # 步长衰减系数
92

93
        # LBAS特定参数
94
        self.P = P or 0.5  # 二次进化概率
95
        self.t0 = 30  # 初始阶段迭代次数（之后开始扰动）
96
        self.reducing = 1  # 用于记录衰减系数
97

98
        # 初始化种群
99
        self.population = self.init_population(init_x)
100
        self.fitness = np.array([self.obj_func(ind) for ind in self.population])
101

102
        # 记录最佳解
103
        self.best_idx = np.argmin(self.fitness)
104
        self.best_solution = self.population[self.best_idx].copy()
105
        self.best_fitness = self.fitness[self.best_idx]
106

107
        # 记录收敛过程
108
        self.convergence = []
109
        self.iteration = 0
110

111
    def init_population(self, init_x=None):
112
        """初始化种群"""
113

114
        population = self.rng.uniform(
115
            low=self.bounds[:, 0],
116
            high=self.bounds[:, 1],
117
            size=(self.pop_size, self.dim)
118
        )
119

120
        if init_x != None:
121
            if len(init_x) != self.dim:
122
                raise ValueError("init_x的维度必须与问题维度一致")
123
            population[0] = self.clip(init_x)
124

125
        return population
126

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

131
    def iterate(self):
132
        """执行一次迭代"""
133
        thisframe_lines = []
134
        self.iteration += 1
135
        if self.with_plot:
136
            self.positions_history.append(self.population.copy())
137
        # 对种群按适应度排序
138
        sorted_indices = np.argsort(self.fitness)
139
        sorted_population = self.population[sorted_indices]
140
        sorted_fitness = self.fitness[sorted_indices]
141

142
        # 更新全局最优
143
        if sorted_fitness[0] < self.best_fitness:
144
            self.best_fitness = sorted_fitness[0]
145
            self.best_solution = sorted_population[0].copy()
146

147
        self.best_history.append((self.best_solution.copy(), self.best_fitness))
148

149
        # 对每个个体进行更新
150
        for i in range(self.pop_size):
151
            # 当前个体
152
            x = self.population[i].copy()
153

154
            # 1. 生成随机方向向量
155
            b = self.rng.uniform(-1, 1, self.dim)
156
            b = b / np.linalg.norm(b)  # 单位化
157
            if feature_config.get('randomly_tester_d0', False):
158
                b = b * self.rng.random()  # 随机化
159

160
            # 2. 计算左右须位置
161
            x_left = x - self.d0 * b
162
            x_right = x + self.d0 * b
163

164
            # 确保在边界内
165
            x_left = self.clip(x_left)
166
            x_right = self.clip(x_right)
167

168
            # 3. 计算适应度
169
            f_left = self.obj_func(x_left)
170
            f_right = self.obj_func(x_right)
171

172
            # 4. 决定移动方向并应用惯性权重
173
            inertia = self.delta0 * b
174
            # 如果在计算时已经随机化步长，那现在就不随机化了
175
            if not feature_config.get('randomly_tester_d0', False):
176
                inertia = inertia * self.rng.random()
177
            if f_right > f_left:  # 向左移动
178
                x_new = x - inertia
179
            else:  # 向右移动
180
                x_new = x + inertia
181
            if self.with_plot:
182
                thisframe_lines.append((x_left, x_right, {  # 天牛须
183
                    'color': 'black',
184
                    'linestyle': '--',
185
                    'linewidth': 1
186
                }))
187
                thisframe_lines.append((x, x_new, {  # 天牛主动行动
188
                    'color': 'red',
189
                    'linestyle': '-',
190
                    'linewidth': 4
191
                }))
192

193
            # 确保在边界内
194
            x_new = self.clip(x_new)
195
            f_new = self.obj_func(x_new)
196

197
            # 5. 按概率P进行二次进化
198
            idx_in_sorted = np.where(sorted_indices == i)[0][0]
199
            if idx_in_sorted > 0:  # 如果不是最优个体
200
                # 找到排在当前个体前面的个体
201
                if self.rng.random() < self.P:
202
                    # 获取前一个个体和全局最优个体
203
                    x_prev = sorted_population[idx_in_sorted - 1]
204
                    x_best = self.best_solution
205

206
                    # 二次进化：向更优个体和全局最优个体学习
207
                    # 使用正确的公式，基于上一次的位置
208
                    dx_secondary = 0.5 * (x_prev - x_new) + 0.5 * (x_best - x_new)
209
                    if feature_config.get('only_forward_best', False):
210
                        dx_secondary = 0 * (x_prev - x_new) + 1 * (x_best - x_new)
211

212
                    use_rand = self.rng.random() if feature_config.get('secondary_evolution_use_random', False) else 1
213
                    x_secondary = x_new + dx_secondary * use_rand
214
                    if feature_config.get('recurrent_bug_from_raw_paper', False):
215
                        x_secondary = dx_secondary * use_rand
216

217
                    if self.with_plot:
218
                        thisframe_lines.append((x_new, x_prev, {  # 二次进化，向更优个体学习
219
                            'color': 'cyan',
220
                            'linestyle': ':',
221
                            'linewidth': 2
222
                        }))
223
                        thisframe_lines.append((x_new, x_best, {  # 二次进化，向全局最优学习
224
                            'color': 'salmon',
225
                            'linestyle': '-.',
226
                            'linewidth': 2
227
                        }))
228
                        thisframe_lines.append((x_new, x_new + dx_secondary, {  # 二次进化总方向
229
                            'color': 'brown',
230
                            'linestyle': '-',
231
                            'linewidth': 2
232
                        }))
233
                        thisframe_lines.append((x_new, x_secondary, {  # 完整二次进化
234
                            'color': 'brown',
235
                            'linestyle': '-',
236
                            'linewidth': 4
237
                        }))
238

239
                    x_secondary = self.clip(x_secondary)
240
                    f_secondary = self.obj_func(x_secondary)
241

242
                    # 选择较好的解
243
                    if f_secondary < f_new:
244
                        x_new = x_secondary
245
                        f_new = f_secondary
246
                elif not feature_config.get('disable_regenerate', False):
247
                    # 未进行二次进化时，生成随机新个体
248
                    # 从最优个体开始往随机方向前进当前的两须长度再带个随机量
249
                    random_direction = self.rng.uniform(-1, 1, self.dim)
250
                    # 归一化方向向量
251
                    norm = np.linalg.norm(random_direction)
252
                    if norm > 0:
253
                        random_direction = random_direction / norm
254

255
                    # 前进距离
256
                    step_length = self.delta0 * self.rng.random()  # 随机量设为0到d之间的值
257

258
                    # 生成新个体
259
                    x_random = self.best_solution + random_direction * step_length
260
                    x_random = self.clip(x_random)
261
                    f_random = self.obj_func(x_random)
262

263
                    # 如果新个体更优则接受
264
                    if f_random < f_new:
265
                        x_new = x_random
266
                        f_new = f_random
267

268

269
            # 6. 更新个体
270
            self.population[i] = x_new
271
            self.fitness[i] = f_new
272

273
        # 7. 对最优个体进行扰动（如果超过初始阶段）
274
        if (not feature_config.get('disable_turbulent', False)) and self.iteration > self.t0:
275
            # 对最优个体添加扰动
276
            perturbation = self.rng.uniform(-1, 1, self.dim) * self.availd_domain * np.exp(-4 * self.iteration / self.max_iter)
277
            x_perturbed = self.best_solution + perturbation
278
            x_perturbed = self.clip(x_perturbed)
279
            f_perturbed = self.obj_func(x_perturbed)
280

281
            # 如果扰动后的解更好，则替换最优解
282
            if f_perturbed < self.best_fitness:
283
                self.best_solution = x_perturbed
284
                self.best_fitness = f_perturbed
285

286
                # 也用扰动后的解替换种群中最差的个体
287
                worst_idx = np.argmax(self.fitness)
288
                self.population[worst_idx] = x_perturbed
289
                self.fitness[worst_idx] = f_perturbed
290

291
        # 早停
292
        # if self.best_fitness < self.tol:
293
        #     return True
294

295
        # 8. 更新步长和两须距离
296
        self.reducing = self.reducing * self.d_eta
297
        self.d0 = self.d_eta * self.d0
298
        self.delta0 = self.delta_eta * self.delta0
299

300
        # 记录当前最优适应度
301
        self.convergence.append(self.best_fitness)
302
        if self.with_plot:
303
            self.onframe_lines.append(thisframe_lines)
304
        return False
305

306
    def optimize(self):
307
        """执行优化过程"""
308
        if self.with_plot:
309
            self.positions_history = []
310
        self.convergence = []
311
        iters = trange if self.tqdm else range
312
        for self.iteration in iters(self.max_iter):
313
            if self.iterate():
314
                break
315
            if self.history_per_iter is not None and self.iteration % self.history_per_iter == 0:
316
                print(f"I{self.iteration:<6} | Fit={self.best_fitness:.4e}, rdc={self.reducing:.2e}")
317

318
        return self.best_solution, self.best_fitness, self.convergence
319

320
    def create_contour_plot(self):
321
        """创建等高线图"""
322
        x = np.linspace(self.bounds[0][0], self.bounds[0][1], 100)
323
        y = np.linspace(self.bounds[1][0], self.bounds[1][1], 100)
324
        X, Y = np.meshgrid(x, y)
325
        Z = np.zeros_like(X)
326

327
        for i in range(X.shape[0]):
328
            for j in range(X.shape[1]):
329
                Z[i, j] = self.obj_func(np.array([X[i, j], Y[i, j]]))
330

331
        plt.figure(figsize=(12, 10))
332
        contour = plt.contour(X, Y, Z, 50, cmap='viridis', alpha=0.6)
333
        plt.colorbar(contour, label='Function Value')
334
        plt.xlabel('X')
335
        plt.ylabel('Y')
336
        plt.title('LBAS Optimization Process')
337
        plt.grid(True, alpha=0.3)
338

339
        return plt.gca(), X, Y, Z
340

341
    def animate_optimization(self):
342
        """创建优化过程动画"""
343
        fig, ax = plt.subplots(figsize=(12, 10))
344
        ax: plt.Axes  # type: ignore
345

346
        # 创建等高线
347
        x = np.linspace(self.bounds[0][0], self.bounds[0][1], 100)
348
        y = np.linspace(self.bounds[1][0], self.bounds[1][1], 100)
349
        X, Y = np.meshgrid(x, y)
350
        Z = np.zeros_like(X)
351

352
        for i in range(X.shape[0]):
353
            for j in range(X.shape[1]):
354
                Z[i, j] = self.obj_func(np.array([X[i, j], Y[i, j]]))
355

356
        contour = ax.contour(X, Y, Z, 30, cmap='viridis', alpha=0.6)
357
        plt.colorbar(contour, ax=ax, label='Function Value')
358

359
        # 初始化图形元素
360
        beetles_scatter = ax.scatter([], [], c='blue', s=30, label='Beetles')
361
        best_beetle_scatter = ax.scatter([], [], c='red', s=100, marker='*', label='Best Beetle')
362
        antennae_lines = [ax.plot([], [], 'k-', alpha=0.5)[0] for _ in range(self.pop_size)]
363
        prev_lines = 100
364
        move_lines = [ax.plot([], [], 'r-', alpha=0.5)[0] for _ in range(prev_lines)]
365

366
        # 添加文本框显示迭代信息
367
        iter_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12,
368
                           bbox=dict(facecolor='white', alpha=0.7))
369
        best_text = ax.text(0.02, 0.88, '', transform=ax.transAxes, fontsize=10,
370
                           bbox=dict(facecolor='white', alpha=0.7))
371

372
        ax.set_xlim(self.bounds[0])
373
        ax.set_ylim(self.bounds[1])
374
        ax.set_xlabel('X')
375
        ax.set_ylabel('Y')
376
        ax.set_title('Leader Disturbance Beetle Antennae Search Algorithm')
377
        ax.legend()
378
        ax.grid(True, alpha=0.3)
379

380
        # 动画更新函数
381
        def update(frame):
382
            if frame < len(self.positions_history) - 1:  # 最后一帧有位置但没动作，索性全都不更新
383
                positions = self.positions_history[frame]
384
                best_pos, best_fit = self.best_history[frame]
385
                going_arrows: list[tuple[np.ndarray, np.ndarray, dict]] = self.onframe_lines[frame]
386

387
                # 更新天牛位置
388
                beetles_scatter.set_offsets(positions)
389

390
                # 更新最优天牛位置
391
                best_beetle_scatter.set_offsets([best_pos])
392

393
                # 更新移动箭头
394
                # for i, arrow in enumerate(going_arrows):
395
                #     start, end = arrow
396
                #     move_lines[i].set_data([start[0], end[0]], [start[1], end[1]])
397
                for line, arrow in zip(move_lines, going_arrows + [((-9999, -9999), (-9999, -9999), {})] * prev_lines):
398
                    start, end, style = arrow
399
                    line.set_data([start[0], end[0]], [start[1], end[1]])
400
                    line.set_color(style.get('color', 'black'))
401
                    line.set_linestyle(style.get('linestyle', '-'))
402
                    line.set_linewidth(style.get('linewidth', 1))
403

404
                # 更新文本信息
405
                iter_text.set_text(f'Iteration: {frame+1}/{self.max_iter}')
406
                best_text.set_text(f'Best Fitness: {best_fit:.6f}\nBest Position: ({best_pos[0]:.4f}, {best_pos[1]:.4f})')
407

408
            return beetles_scatter, best_beetle_scatter, *antennae_lines, iter_text, best_text
409

410
        ani = FuncAnimation(fig, update, frames=len(self.positions_history),
411
                          interval=200, blit=True)
412
        return ani
413
# 测试算法
414
if __name__ == "__main__":
415

416
    from pprint import pprint
417

418
    pprint(feature_config)
419
    def sphere_function(x):
420
        def function(x):
421
            return np.sum(x**2)
422
        return function, (-100, 100)
423
    # 参数设置
424
    center = 0
425
    test_function = sphere_function
426
    obj, func_bound = test_function(center)
427
    dim = 2  # 问题维度
428
    bounds = np.array([func_bound] * dim)  # 求解范围
429
    pop_size = 3  # 种群大小
430
    max_iter = 50  # 最大迭代次数
431

432
    # 创建并运行算法
433
    eta = 0.9
434
    lbas = LBAS(obj, pop_size, bounds, max_iter,
435
                d_eta=eta, delta_eta=eta, tol=1e-5, with_plot=True)
436
    best_solution, best_fitness, convergence = lbas.optimize()
437

438
    print()
439
    print(f"Test params: {dim=},{center=},{test_function.__name__=},{bounds=},")
440
    print(f"Algorithm params: {pop_size=},{max_iter=}")
441
    print(f"solution: {best_solution}")
442
    print(f"fitness: {best_fitness:.6e}")
443
    print(f"error: {np.linalg.norm(best_solution - center)}")
444
    print(f"iter: {lbas.iteration} / {lbas.max_iter}")
445

446
    # 绘制收敛曲线
447
    # plt.figure(figsize=(10, 6))
448
    # plt.semilogy(convergence)
449
    # plt.title('Convergence Curve')
450
    # plt.xlabel('Iteration')
451
    # plt.ylabel('Best Fitness (log scale)')
452
    # plt.grid(True)
453
    # plt.savefig(f'fig/fig{time.time()}.png')
454
    # plt.show()
455

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