从0开始机器学习--9.启发式算法（pulp，粒子群PSO、模拟退火SA、遗传GA算法剖析+对比，TOPSIS、NSGA-II组合算法介绍，含代码）

PuLP 是一个用于在 Python 中定义和解决线性规划 (LP) 和整数规划 (IP) 问题的库。它为定义线性规划问题提供了简单的接口，并可以与多个求解器（如 COIN-OR, CPLEX, Gurobi）配合使用。PuLP 的强大之处在于它能够使用 Python 编程的灵活性来定义线性规划模型，并将其转换为标准的数学规划格式，供底层求解器进行优化。

代码

import pulp
# 创建一个线性规划问题实例（最大化问题）
prob = pulp.LpProblem("Maximize_Profit", pulp.LpMaximize)
# 定义决策变量A和B，均为非负
A = pulp.LpVariable('A', lowBound=0, cat='Continuous')
B = pulp.LpVariable('B', lowBound=0, cat='Continuous')
# 定义目标函数
prob += 40 * A + 30 * B, "Total_Profit"
# 定义约束条件
prob += 2 * A + 1 * B <= 8, "Labor_Constraint"
prob += 3 * A + 2 * B <= 12, "Material_Constraint"
# 求解问题
prob.solve()
# 输出结果
print("Status:", pulp.LpStatus[prob.status])
# 输出解
print(f"Optimal number of product A to produce: {A.varValue}")
print(f"Optimal number of product B to produce: {B.varValue}")
# 输出最大利润
print(f"Maximum Profit: {pulp.value(prob.objective)}")

定义问题：通过 pulp.LpProblem("问题名称", pulp.LpMaximize) 定义一个线性规划的最大化问题。
决策变量：定义了两个决策变量。
𝐴和𝐵，代表生产的产品数量，它们都是非负的连续变量。
目标函数：我们通过 prob += 40 * A + 30 * B, "Total_Profit" 来定义目标函数，表示最大化利润。
约束条件：通过 prob += 添加约束条件，例如劳动时间和原材料的限制。
求解：使用 prob.solve() 调用默认的 COIN-OR 求解器来解决线性规划问题。
输出结果：通过 A.varValue 和 B.varValue 获取最优解，通过 pulp.value(prob.objective) 获取最大化的利润。

输出：

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

At line 2 NAME          MODEL
At line 3 ROWS
At line 7 COLUMNS
At line 14 RHS
At line 17 BOUNDS
At line 18 ENDATA
Problem MODEL has 2 rows, 2 columns and 4 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Presolve 2 (0) rows, 2 (0) columns and 4 (0) elements
0  Obj -0 Dual inf 85 (2)
0  Obj -0 Dual inf 85 (2)
1  Obj 180
Optimal - objective value 180
Optimal objective 180 - 1 iterations time 0.002
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.02

Status: Optimal
Optimal number of product A to produce: 0.0
Optimal number of product B to produce: 6.0
Maximum Profit: 180.0

粒子群算法（PSO）

粒子群算法（PSO）是群智能算法的一种，它通过模拟自然界中生物群体的行为来搜索最优解。是基于迭代的方法

原理

鸟群：假设一个区域，所有的鸟都不知道食物的位置，但是它们知道当前位置离食物还有多远
PSO：每个解看作一只鸟，称为粒子，所有粒子都有一个适应值（由目标函数决定），每个粒子都有一个速度决定它们的飞翔方向和距离，粒子们追随当前最优粒子在解空间中搜索。粒子们知道当前自己的位置和自己已发现的最好位置（pbest）以及所有pbest中最好的gbest。粒子通过自己的经验和同伴中最好的经验来决定下一步的运动。

参数（利用网格搜索法确定最佳）

粒子数：粒子群的大小，一般取20-40，对于较难或特定类别的问题可以取100-200（若过小：结果不准；过大：效率很慢）；
最大速度：决定粒子在一个循环中最大的移动距离，通常设定为粒子的范围宽度（若过大，可能结果会不在有效范围之内）；
终止条件：最大循环数以及最小错误要求；
惯性权重：w，控制粒子的搜索范围，有能力探索新的区域，较大的w有利于跳出局部极值，而较小的有利于算法收敛。因此随着迭代的进行，线性的减小w： $w=w_{max}-(w_{max}-w_{min})/i_{max}*i_{now}$ 其中， $w_{max}$ 、 $w_{min}$ :w的最大最小值； $i_{max}$ :最大迭代次数； $i_{now}$ :当前迭代次数；
c1、c2：加速因子，两个大于0的常数，属于0-1之间，将粒子推向统计加速项的权重，较大会导致粒子突然冲向或越过目标区域，较小则允许粒子再被拉回之前在目标区域外徘徊；
控制参数fac=c1+c2，fac越大，则粒子的位置变化越快；实验证明fac=4.1（通常c1=2.0,c2=2.0）具有很好的收敛效果；

第i个粒子的速度

v(i+1) = w * v(i) + c1 * r1(i) * (pbest(i) - x(i)) + c2 * r2(i) * (gbest - x(i))，其中：

w、c1、c2是超参数，r1、r2是随机数；
x(i)是第i个粒子的位置，v(i-1)是第i-1个粒子的速度，pbest(i)是第i个粒子的最佳位置，gbest是所有pbest中最好的位置。
整个c1项为认知部分，即来源于自己的经验部分；c2项为社会部分，即来源于群体中其他优秀微粒的经验；c1越大，则粒子越容易接受认知信息，c2越大，则粒子越容易接受社会信息。

第i个粒子的位置

x(i+1) = x(i) + v(i)

算法步骤

初始化：随机生成一组粒子，并赋予其初始位置和速度。
评估：计算每个粒子的适应度值，即目标函数值。
选择：根据适应度值选择出最优的粒子。
交换：将最优的粒子的位置和速度信息传递给其他粒子。（设置最优的pbest为gbest）
更新：根据最优的粒子的位置信息更新粒子的位置和速度信息。
重复：重复步骤2-5，直到收敛（全局最优解）或达到最大迭代次数。

优点

需要调整的参数相对较少；
适用于连续和离散的优化问题，能处理多维、复杂约束的优化问题；
可以很容易地进行并行处理。

缺点

粒子群算法的收敛速度慢，在大型问题中，需要较多的迭代次数才能找到全局最优解；
粒子群算法的初始条件很敏感，依赖于这些参数的设置，初始条件的选择会影响最终的结果；
粒子群算法依赖于初始解，初始解的选择会影响最终的结果。

适应问题类型

连续优化问题：PSO 特别适合连续解空间的问题，如函数优化、参数估计等。它通过粒子的速度和位置更新，能够快速找到近似最优解。
多峰优化问题：PSO 可以处理具有多个峰值的复杂优化问题，通过全局和局部搜索的结合，较好地避免陷入局部最优。
动态优化问题：PSO 能够适应动态变化的环境，通过粒子位置和速度的实时更新来应对目标函数的变化。

典型应用场景

控制系统优化：在机器人控制、自动驾驶、无人机导航等领域，PSO 可以用于优化控制策略和路径规划。
信号处理和图像处理：用于滤波器设计、图像分割、特征提取等问题。
经济和金融领域：如投资组合优化、金融时间序列预测、市场模型优化等。

代码寻找Rastrigin函数的最小值位置

Rastrigin函数：具有多个局部最小值，适合测试优化算法的性能。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 定义目标函数（Rastrigin函数）
def objective_function(x):
    return x[0]**2 + x[1]**2 - 10 * np.cos(2 * np.pi * x[0]) - 10 * np.cos(2 * np.pi * x[1]) + 20

# 粒子群算法参数
num_particles = 30  # 粒子数量
num_iterations = 100  # 迭代次数
dim = 2  # 维度
w = 0.5  # 惯性权重
c1 = 1.5  # 个体学习因子
c2 = 1.5  # 社会学习因子

# 初始化粒子位置和速度
particles_position = np.random.uniform(-5.12, 5.12, (num_particles, dim))
particles_velocity = np.random.uniform(-1, 1, (num_particles, dim))

# 初始化个体最优位置和全局最优位置
personal_best_position = particles_position.copy()
personal_best_value = np.apply_along_axis(objective_function, 1, personal_best_position)
global_best_position = personal_best_position[np.argmin(personal_best_value)]
global_best_value = np.min(personal_best_value)

# PSO算法迭代
for _ in range(num_iterations):
    r1 = np.random.rand(num_particles, dim)
    r2 = np.random.rand(num_particles, dim)
    particles_velocity = (w * particles_velocity +
                          c1 * r1 * (personal_best_position - particles_position) +
                          c2 * r2 * (global_best_position - particles_position))
    particles_position += particles_velocity
    current_value = np.apply_along_axis(objective_function, 1, particles_position)
    better_mask = current_value < personal_best_value
    personal_best_position[better_mask] = particles_position[better_mask]
    personal_best_value[better_mask] = current_value[better_mask]
    if np.min(personal_best_value) < global_best_value:
        global_best_position = personal_best_position[np.argmin(personal_best_value)]
        global_best_value = np.min(personal_best_value)

# 输出最优解
print("Global Best Position:", global_best_position)
print("Global Best Value:", global_best_value)

# 绘制函数的3D图像
X = np.linspace(-5.12, 5.12, 400)
Y = np.linspace(-5.12, 5.12, 400)
X, Y = np.meshgrid(X, Y)
Z = X**2 + Y**2 - 10 * np.cos(2 * np.pi * X) - 10 * np.cos(2 * np.pi * Y) + 20

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')

# 标记出全局最优解
ax.scatter(global_best_position[0], global_best_position[1], global_best_value, color='r', s=100)
ax.set_title('Rastrigin Function')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

输出：

Global Best Position: [-1.36850095e-10 -1.20921201e-10]
Global Best Value: 0.0

对于梯度的介绍，我在《6.1神经网络基础》中具体介绍过。

模拟退火算法（SA）

用了蒙特卡洛优化（在《8.强化学习》中具体介绍过）方法，是对爬山算法的一种改进；

爬山算法

简单的贪心搜索算法，从当前黑的临近空间解中选择一个最优解作为当前解直到局部最优，但不一定能达到全局最优。

原理

模拟退火算法来源于固体退火原理，引入温度这一变量初始状态设计一个高温，升温内能增大，非最优解也会愿意尝试；随着温度下降，粒子趋于有序，尝试的活性下降。相对于爬山算法，增大了找到全局最优解的概率

是一个双重迭代的算法：

外循环：退火过程（模拟温度下降的过程）；
内循环：Metropolis算法（模拟某一温度下粒子的跳动过程）。

内循环中随机移动后的解如果更好(增量大于0)就接受他，更新位置；如果结果更差(增量小于0)，就以p的概率（ $exp(-df/T)$ ）接受他(产生随机概率值(0<r<1)，如果p>r则接受新解，否则接受旧解)。内循环的最后再更新温度。

优点

有概率跳出局部最优解并趋于全局最优；
可以灵活调整参数以平衡搜索质量和速度；
能处理复杂的连续或离散优化问题；
解和初始状态无关。

缺点

不同的问题可能需要不同的设置，需要合适的初始解和参数调整；
搜索空间的复杂性：搜索空间越大，算法的收敛速度越慢，可能碰到的局部最优解最多；
运行结果受随机概率影响可能有一定变化，随机性强；
SA与其他优化算法不同，不需要初始化种群操作，收敛速度慢。

适应的问题类型

复杂的组合优化问题：SA 适用于需要全局搜索能力的复杂组合优化问题，尤其是目标函数具有多个局部最优解时，如大规模的旅行商问题。
非线性和非凸优化问题：SA 在处理非凸函数优化时表现良好，适用于那些目标函数可能具有多个局部最优的非线性问题。
目标函数不光滑或不可导的优化问题：SA 不依赖于目标函数的导数信息，适用于优化不可导函数或目标函数不光滑的问题。

典型应用场景

资源配置和调度问题：如生产调度、物流运输优化、任务调度等，需要在多个可能解中找到最优配置。
电路设计：用于优化复杂电路的参数设置和布局设计。
生物信息学和物理学问题：如蛋白质折叠、晶体结构预测等问题，这些问题通常具有复杂的能量表面。

double T = 2000;  // 初始温度
double dT = 0.99; // 降温速率，越接近1降温越慢
double eps = 1e-14; // 终止温度，越小越好
while (T > eps) {
    // 这里是每次退火的操作，例如：
    // 1. 生成新解
    // 2. 计算新解的代价
    // 3. 根据代价和当前温度决定是否接受新解,有一定概率接受比当前解较差的解
    T *= dT; // 降温 T减小，接受较差解的概率减小
}

代码寻找Rastrigin函数的最小值位置

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 定义Rastrigin函数作为目标函数
def objective_function(x):
    return x[0]**2 + x[1]**2 - 10 * np.cos(2 * np.pi * x[0]) - 10 * np.cos(2 * np.pi * x[1]) + 20

# 邻域函数，产生一个新的二维解
def neighbor(x):
    return x + np.random.normal(0, 1, size=2)  # 生成一个二维向量

# 模拟退火算法
def simulated_annealing(objective_function, initial_solution, temp, cooling_rate, max_iterations, inner_iterations):
    current_solution = initial_solution
    current_value = objective_function(current_solution)
    best_solution = current_solution
    best_value = current_value

    for iteration in range(max_iterations):
        for _ in range(inner_iterations):
            # 生成新解
            new_solution = neighbor(current_solution)
            new_value = objective_function(new_solution)
            # 根据概率判断是否接受新解
            if new_value < current_value or np.random.rand() < np.exp(-(new_value - current_value) / temp):
                current_solution = new_solution
                current_value = new_value
            # 更新最优解
            if current_value < best_value:
                best_solution = current_solution
                best_value = current_value
        # 降低温度
        temp *= cooling_rate
        print(f'Iteration {iteration + 1}, Temp: {temp:.4f}, Best Value: {best_value:.4f}')
    
    return best_solution, best_value

# 参数设置
initial_solution = np.random.uniform(-5.12, 5.12, size=2)  # 初始解，二维
initial_temp = 100  # 初始温度
cooling_rate = 0.9  # 降温速率
max_iterations = 50  # 最大外循环次数
inner_iterations = 10  # 每次降温时的内循环次数

# 运行模拟退火算法
best_solution, best_value = simulated_annealing(objective_function, initial_solution, initial_temp, cooling_rate, max_iterations, inner_iterations)

# 输出最优解
print("Best Solution:", best_solution)
print("Best Value:", best_value)

# 绘制Rastrigin函数的三维图像
X = np.linspace(-5.12, 5.12, 400)
Y = np.linspace(-5.12, 5.12, 400)
X, Y = np.meshgrid(X, Y)
Z = X**2 + Y**2 - 10 * np.cos(2 * np.pi * X) - 10 * np.cos(2 * np.pi * Y) + 20

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')

# 标记出模拟退火算法找到的最优解
ax.scatter(best_solution[0], best_solution[1], best_value, color='r', s=100)
ax.set_title('Rastrigin Function with Simulated Annealing')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

输出：

Iteration 1, Temp: 90.0000, Best Value: 16.0347
Iteration 2, Temp: 81.0000, Best Value: 13.4949
Iteration 3, Temp: 72.9000, Best Value: 13.4949
Iteration 4, Temp: 65.6100, Best Value: 13.4949
Iteration 5, Temp: 59.0490, Best Value: 13.4949
Iteration 6, Temp: 53.1441, Best Value: 13.4949
Iteration 7, Temp: 47.8297, Best Value: 13.4949
Iteration 8, Temp: 43.0467, Best Value: 13.4949
Iteration 9, Temp: 38.7420, Best Value: 13.4949
Iteration 10, Temp: 34.8678, Best Value: 13.4949
Iteration 11, Temp: 31.3811, Best Value: 13.4949
Iteration 12, Temp: 28.2430, Best Value: 13.4949
Iteration 13, Temp: 25.4187, Best Value: 8.5710
Iteration 14, Temp: 22.8768, Best Value: 5.9105
Iteration 15, Temp: 20.5891, Best Value: 5.9105
Iteration 16, Temp: 18.5302, Best Value: 5.9105
Iteration 17, Temp: 16.6772, Best Value: 5.9105
Iteration 18, Temp: 15.0095, Best Value: 5.9105
Iteration 19, Temp: 13.5085, Best Value: 5.9105
Iteration 20, Temp: 12.1577, Best Value: 5.9105
Iteration 21, Temp: 10.9419, Best Value: 5.9105
Iteration 22, Temp: 9.8477, Best Value: 5.9105
Iteration 23, Temp: 8.8629, Best Value: 5.9105
Iteration 24, Temp: 7.9766, Best Value: 5.9105
Iteration 25, Temp: 7.1790, Best Value: 2.2807
Iteration 26, Temp: 6.4611, Best Value: 1.8803
Iteration 27, Temp: 5.8150, Best Value: 1.8803
Iteration 28, Temp: 5.2335, Best Value: 1.8803
Iteration 29, Temp: 4.7101, Best Value: 1.8803
Iteration 30, Temp: 4.2391, Best Value: 1.8803
Iteration 31, Temp: 3.8152, Best Value: 1.8803
Iteration 32, Temp: 3.4337, Best Value: 0.4254
Iteration 33, Temp: 3.0903, Best Value: 0.4254
Iteration 34, Temp: 2.7813, Best Value: 0.4254
Iteration 35, Temp: 2.5032, Best Value: 0.4254
Iteration 36, Temp: 2.2528, Best Value: 0.4254
Iteration 37, Temp: 2.0276, Best Value: 0.4254
Iteration 38, Temp: 1.8248, Best Value: 0.4254
Iteration 39, Temp: 1.6423, Best Value: 0.4254
Iteration 40, Temp: 1.4781, Best Value: 0.4254
Iteration 41, Temp: 1.3303, Best Value: 0.4254
Iteration 42, Temp: 1.1973, Best Value: 0.4254
Iteration 43, Temp: 1.0775, Best Value: 0.4254
Iteration 44, Temp: 0.9698, Best Value: 0.4254
Iteration 45, Temp: 0.8728, Best Value: 0.4254
Iteration 46, Temp: 0.7855, Best Value: 0.4254
Iteration 47, Temp: 0.7070, Best Value: 0.4254
Iteration 48, Temp: 0.6363, Best Value: 0.4254
Iteration 49, Temp: 0.5726, Best Value: 0.4254
Iteration 50, Temp: 0.5154, Best Value: 0.4254
Best Solution: [0.01949567 0.04213167]
Best Value: 0.42543128700373956

遗传算法（GA）

同样用了蒙特卡洛优化（在《8.强化学习》中具体介绍过）方法。

概念

染色体：问题中个体的某种字符串形式的编码表示；
基因：染色体（字符串）中的单个字符。

算法步骤

初始化种群：设定种群规模，编码染色体，产生厨师种群。随机生成N个染色体，每个染色体由N个基因组成。
计算适应度：定义适应度函数，对于每个染色体，计算其对应的适应度值（代入到目标函数中的值），适应度值越高，染色体的优劣程度越高。
判断是否终止（迭代次数或误差最小的那一项满足条件）。
选择-复制：对于一个规模为N的种群S，按每个染色体的选择概率P=f(xi)/求和f(xj)[j从1到N]所决定的选择机会，分N次从S中随机选定1个染色体（赌轮选择法），并进行复制，组成群体S1。其中，f:目标函数，即适应度值。
交叉：按交叉率从S1中随机确定c个染色体，配对进行交叉操作，产生新的染色体代替原染色体，得到S2。
变异：按变异率从S2中随机确定m个染色体，分别进行变异操作，产生新的染色体代替原染色体，得到S3。
生成新一代种群，即S3。迭代次数加一。
迭代：重复步骤2-7，直到满足终止条件。

控制参数

种群大小N：种群的大小决定了算法的复杂度。
最大换代数（可理解为迭代次数）。
交叉率（crossover rate）：参加交叉运算的染色体个数占全体染色体总数的比例，记为Pc，一般取0.4-0.99。
变异率（mutation rate）：发生变异的基因位数所占全体染色体的基因总位数的比例，记为Pm，一般取0.0001-0.1。
如计算可以变异的基因位数：4（个染色体）*5（位基因）*Pm=0.0001*4*5=0.002，不足一位，所以本轮遗传操作不变异。

优点

全局搜索能力强：GA 通过模拟自然选择和遗传变异过程，具有较强的全局搜索能力，能够有效探索解空间。
适应性强：GA 能够处理复杂的非线性、多峰优化问题，适用于广泛的应用场景。
多样性维护：GA 通过种群多样性保持（如交叉、变异等操作），能够避免过早收敛到局部最优解。
无需梯度信息：GA 不需要目标函数的梯度信息，适合处理不可导或不连续的优化问题。

缺点

收敛速度较慢：GA 的进化过程往往需要较多的迭代次数，导致整体收敛速度较慢，尤其是在优化问题较复杂时。
参数设置复杂：GA 需要设置较多参数，如种群大小、交叉概率、变异概率等，这些参数的选择对算法性能有较大影响。
易陷入局部最优：尽管 GA 具有全局搜索能力，但在复杂问题中，种群可能会逐渐趋同，导致陷入局部最优解。
计算代价较高：由于需要维护和操作一个种群，GA 的计算开销较大，尤其是在大规模问题中。

适应问题类型

组合优化问题：如旅行商问题（TSP）、背包问题、排课问题、工厂调度问题等，这类问题通常有离散的解空间，GA 的种群多样性和进化机制使其能够有效探索这些复杂的解空间。
非线性优化问题：遗传算法适用于处理复杂的非线性优化问题，尤其是在目标函数具有多个局部最优解的情况下。
多目标优化问题：GA 可以同时处理多个目标，通过进化机制生成一组解，满足不同目标的平衡。
规则发现和符号回归：GA 可以用来发现数据中的规则，或进行符号回归，从而在数据中找到潜在的函数关系。

典型应用场景

人工智能和机器学习：用于特征选择、超参数优化、神经网络结构优化等。
工程设计优化：如结构优化、机械设计、电子电路设计等需要考虑多目标的复杂工程问题。
生物信息学：基因序列比对、基因组拼接、蛋白质结构预测等。

代码寻找Rastrigin函数的最小值位置

适应度函数：将 Rastrigin 函数定义为优化目标（注意在遗传算法中为了选择最大值，这里用负值做适应度）
种群初始化：生成二维的解（x 和 y）。
交叉操作和变异操作：这些操作都在二维向量上进行。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 定义Rastrigin函数作为目标函数
def fitness_function(x):
    return -(x[0]**2 + x[1]**2 - 10 * np.cos(2 * np.pi * x[0]) - 10 * np.cos(2 * np.pi * x[1]) + 20)

# 初始化种群，二维
def initialize_population(size, x_min, x_max):
    return np.random.uniform(x_min, x_max, (size, 2))

# 选择父代（轮盘赌选择法）
def select_parents(population, fitness_values):
    total_fitness = np.sum(fitness_values)
    selection_prob = fitness_values / total_fitness
    selected_idx = np.random.choice(len(population), size=len(population), p=selection_prob)
    return population[selected_idx]

# 交叉操作
def crossover(parent1, parent2):
    alpha = np.random.rand()
    return alpha * parent1 + (1 - alpha) * parent2

# 变异操作
def mutate(child, mutation_rate, x_min, x_max):
    if np.random.rand() < mutation_rate:
        return np.random.uniform(x_min, x_max, size=2)
    return child

# 遗传算法主循环
def genetic_algorithm(pop_size, x_min, x_max, generations, mutation_rate):
    population = initialize_population(pop_size, x_min, x_max)
    for generation in range(generations):
        fitness_values = np.array([fitness_function(x) for x in population])
        parents = select_parents(population, fitness_values)

        next_generation = []
        for i in range(0, pop_size, 2):
            parent1, parent2 = parents[i], parents[i+1]
            child1, child2 = crossover(parent1, parent2), crossover(parent2, parent1)
            child1, child2 = mutate(child1, mutation_rate, x_min, x_max), mutate(child2, mutation_rate, x_min, x_max)
            next_generation.extend([child1, child2])
        population = np.array(next_generation)
    
    # 返回最优解
    best_individual = population[np.argmax([fitness_function(x) for x in population])]
    best_fitness = fitness_function(best_individual)
    return best_individual, best_fitness

# 参数设置
pop_size = 20
x_min, x_max = -5.12, 5.12
generations = 100
mutation_rate = 0.1

# 运行遗传算法
best_solution, best_fitness = genetic_algorithm(pop_size, x_min, x_max, generations, mutation_rate)
print(f"最优解: x = {best_solution}, 适应度值 = {best_fitness}")

# 绘制Rastrigin函数的三维图像
X = np.linspace(-5.12, 5.12, 400)
Y = np.linspace(-5.12, 5.12, 400)
X, Y = np.meshgrid(X, Y)
Z = X**2 + Y**2 - 10 * np.cos(2 * np.pi * X) - 10 * np.cos(2 * np.pi * Y) + 20

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')

# 标记出遗传算法找到的最优解
ax.scatter(best_solution[0], best_solution[1], -best_fitness, color='r', s=100)  # -best_fitness，因为fitness是负值
ax.set_title('Rastrigin Function with Genetic Algorithm')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

输出：

最优解: x = [-0.92258478 -1.0867697 ], 适应度值 = -4.641796195714294

PSO、SA、GA对比总结

不同

并行：PSO可（粒子的独立性），SA不可（串行降温），GA部分可（按种群分为子种群）；
搜索机制：PSO通过群体合作和个体经验共享搜索最优解，SA利用概率性接受劣解的机制寻找最优解，GA模拟自然选择和遗传变异，通过交叉和变异等操作产生新解；
全局最优：PSO相对表现最差；
初始条件依赖性：PSO对初始解和参数敏感,SA初始温度和冷却速率的选择对算法结果有较大影响,GA对初始种群的多样性要求较高；
收敛速度：PSO通常在早期收敛速度较快，而GA和SA的收敛速度相对较慢。SA在搜索后期的收敛速度尤其慢，因为温度逐渐降低，导致搜索步长减小。

相同

适应性：三者都有很强的适应性；
都是从任一解出发，按照某种机制，以一定概率在整个求解空间中探索最优解。

TOPSIS，NSGA-II组合算法

TOPSIS（Technique for Order Preference by Similarity to Ideal Solution）和NSGA-II（Non-dominated Sorting Genetic Algorithm II）的组合算法是一种将多标准决策分析和进化算法结合的优化方法，用于解决复杂的多目标优化问题。它们的结合能够兼顾多个冲突目标，得到接近最优解的 Pareto 前沿，同时提供具体方案的优劣排序。

TOPSIS

TOPSIS 是一种经典的多标准决策分析方法，其基本思想是通过计算各备选方案与“理想解”和“负理想解”的距离，来对备选方案进行排序。理想解是各指标的最好值，负理想解是各指标的最差值。TOPSIS 的步骤如下：

标准化决策矩阵：将所有指标标准化，使得它们具有可比性。
加权标准化矩阵：考虑各指标的重要性，乘以权重。
计算理想解和负理想解：找到各指标的最大值和最小值，分别构成理想解和负理想解。
- 理想解：在每个评价指标上，选择最优值（例如最小的成本或最大的效益）。
- 负理想解：在每个评价指标上，选择最差值（例如最大的成本或最小的效益）。
计算与理想解和负理想解的欧几里得距离。
计算相对接近度：根据与理想解和负理想解的距离计算每个方案的相对接近度，距离越小方案越优。
排序：根据相对接近度对所有方案进行排序。

TOPSIS 可以帮助决策者在多个方案中找到离最优解最近的方案，但它仅适合处理多个目标指标，并不具备全局优化的能力。

适用于以下情况：

多属性决策：当一个决策问题涉及多个评价标准时（例如价格、质量、服务等），而决策者需要在这些标准之间进行权衡时。
确定最优方案：当我们有多个备选方案（例如多个产品、多个供应商等）时，TOPSIS 帮助选择最符合要求的方案。

应用场景

产品选择：比如选择最适合的智能手机，综合考虑价格、性能、用户评价、续航等因素。
供应商选择：例如选择一个供应商，综合考虑供应能力、交货期、价格、质量等多个标准。
投资决策：在多个投资选项中，评估每个选项在多个标准（如回报率、风险、流动性等）上的表现，选择最合适的投资方案。

代码

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler

def topsis(data, weights, impacts):
    # 归一化决策矩阵
    scaler = MinMaxScaler()
    norm_data = scaler.fit_transform(data)
    
    # 权重化
    weighted_data = norm_data * weights
    
    # 理想解和负理想解
    ideal_solution = np.max(weighted_data, axis=0) if impacts == 'max' else np.min(weighted_data, axis=0)
    negative_ideal_solution = np.min(weighted_data, axis=0) if impacts == 'max' else np.max(weighted_data, axis=0)
    
    # 计算距离
    distance_ideal = np.sqrt(np.sum((weighted_data - ideal_solution) ** 2, axis=1))
    distance_negative_ideal = np.sqrt(np.sum((weighted_data - negative_ideal_solution) ** 2, axis=1))
    
    # 计算相对接近度
    relative_closeness = distance_negative_ideal / (distance_ideal + distance_negative_ideal)
    
    return relative_closeness

# 示例数据：每行代表一个方案，每列代表一个属性
data = np.array([[250, 7.5, 300], 
                 [220, 8.2, 250], 
                 [240, 6.8, 280]])

# 属性的权重
weights = np.array([0.5, 0.3, 0.2])

# 属性影响类型（'max' 表示越大越好，'min' 表示越小越好）
impacts = ['max', 'max', 'min']

# TOPSIS计算
scores = topsis(data, weights, impacts)

# 排序方案
sorted_scores = np.argsort(scores)[::-1]

# 可视化TOPSIS结果
plt.bar(range(1, len(scores)+1), scores[sorted_scores])
plt.xlabel('方案')
plt.ylabel('相对接近度')
plt.title('TOPSIS评分')
plt.xticks(range(1, len(scores)+1), sorted_scores + 1)
plt.show()

数据：data 是多个备选方案的决策矩阵，每行表示一个备选方案，每列表示一个属性。
权重：weights 表示每个属性的相对重要性。
影响方向：impacts 表示每个属性是越大越好（'max'）还是越小越好（'min'）。
输出：代码通过计算每个方案的相对接近度，输出排序后的方案。

NSGA-II

NSGA-II 是一种常用于多目标优化的进化算法。它通过模拟自然选择和遗传操作来逐步进化出一组解决方案，称为 Pareto 前沿解，即一组在所有目标上都互相不被严格支配的解。NSGA-II 的步骤包括：

初始化：随机生成一组种群。
非支配排序：根据 Pareto 优势关系对种群进行分层，生成非支配排序层级。
拥挤度排序：在同一层内使用拥挤度距离来保持解的多样性。
选择、交叉、变异：使用选择、交叉、变异操作生成下一代种群。
迭代更新：重复非支配排序和拥挤度排序，最终得到分布均匀的 Pareto 前沿。

NSGA-II 能有效处理多个目标之间的冲突，得到一个多样化的解集，但它不提供对每个解的具体排序。

解决的问题：

NSGA-II 适用于以下类型的多目标优化问题：

多目标决策：当面临多个目标且目标之间存在冲突时，如何在这些目标之间找到最佳的平衡。
无单一最优解：多目标优化问题通常没有一个单一的最优解，而是存在一组互相之间没有被另一个解“支配”的解，称为“帕累托最优解”。
约束优化：当存在约束条件时，NSGA-II能够处理约束问题，并找到约束条件下的最优解。

应用场景

工程设计：在工程设计中，往往需要平衡多个目标，如结构强度、重量、成本等。
生产调度问题：如生产线调度优化，既要考虑时间的最小化，又要考虑成本的最小化。
资源分配：在多个资源分配问题中，优化目标可能包括时间、成本、效益等多个标准。
机器人路径规划：在机器人路径规划问题中，通常需要考虑路径的长度、能量消耗、时间等多个目标。

代码

import numpy as np
import matplotlib.pyplot as plt
from deap import base, creator, tools, algorithms
import random

# 定义目标函数
def objective1(individual):
    return sum(x**2 for x in individual)  # 直接返回数值，目标越小越好

def objective2(individual):
    return sum((x - 2)**2 for x in individual)  # 直接返回数值，目标越小越好

# 设置遗传算法环境
creator.create("FitnessMulti", base.Fitness, weights=(-1.0, -1.0))  # 两个目标都要最小化
creator.create("Individual", list, fitness=creator.FitnessMulti)

# 生成初始个体
def create_individual():
    return [random.uniform(-5, 5) for _ in range(2)]  # 每个个体有两个基因

# 定义评估函数
def evaluate(individual):
    return objective1(individual), objective2(individual)  # 返回的是一个元组，但每个目标值是数值

# 创建工具箱
toolbox = base.Toolbox()
toolbox.register("individual", tools.initIterate, creator.Individual, create_individual)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.2)
toolbox.register("select", tools.selNSGA2)
toolbox.register("evaluate", evaluate)

# 遗传算法参数
population = toolbox.population(n=100)
probab_crossover = 0.7
probab_mutation = 0.2
number_of_generations = 50

# 运行遗传算法
result = algorithms.eaMuPlusLambda(population, toolbox, mu=100, lambda_=200,
                                   cxpb=probab_crossover, mutpb=probab_mutation,
                                   ngen=number_of_generations, stats=None, halloffame=None)

# 获取结果
front = tools.sortNondominated(population, len(population), first_front_only=True)[0]

# 可视化结果
front_fitness = [ind.fitness.values for ind in front]
front_fitness = np.array(front_fitness)

plt.scatter(front_fitness[:, 0], front_fitness[:, 1], c='r')
plt.xlabel('Objective 1')
plt.ylabel('Objective 2')
plt.title('NSGA-II Pareto Front')
plt.show()

目标函数：objective1 和 objective2 分别表示两个目标函数，目标是最小化这两个目标。
初始化种群：通过 create_individual 函数随机生成一组个体（解），每个个体表示一个潜在解。
选择、交叉、变异：NSGA-II 通过选择、交叉和变异操作生成新的个体，并通过非支配排序（non-dominated sorting）和拥挤度距离（crowding distance）来引导进化过程。
非支配排序：将个体按照支配关系进行排序，从而将其分为不同的等级（Pareto Front）。
进化过程：通过模拟自然选择过程，种群逐步优化，直到达到停止条件（如最大代数）。
输出：最终会输出一组帕累托最优解，表示在两个目标之间的最优折衷解。

gen     nevals
0       100
1       177   
2       181   
3       185   
4       177   
5       177   
6       171   
7       180   
8       183   
9       183   
10      185   
11      178   
12      186   
13      184   
14      181   
15      186   
16      178   
17      175   
18      177   
19      179   
20      179   
21      188   
22      177   
23      171   
24      175   
25      181   
26      181   
27      185   
28      185   
29      180   
30      181   
31      181   
32      179   
33      181   
34      181   
35      180   
36      179   
37      184   
38      175   
39      182   
40      183   
41      183   
42      187   
43      176   
44      168   
45      181   
46      175   
47      174   
48      177   
49      178   
50      174

TOPSIS 和 NSGA-II 的组合算法

将 TOPSIS 和 NSGA-II 结合，可以将 NSGA-II 生成的 Pareto 前沿解进一步细化和排序，帮助决策者在 Pareto 解集内做出具体的选择。具体步骤如下：

使用 NSGA-II 生成 Pareto 前沿：先使用 NSGA-II 生成一组非支配的 Pareto 前沿解。
对 Pareto 前沿解进行 TOPSIS 处理：将 Pareto 前沿解输入 TOPSIS 算法，根据距离理想解的接近度，对解进行排序。
综合考虑权重和距离：TOPSIS 能根据权重分配和指标优化目标，帮助决策者在 Pareto 前沿解集中做出更加合理的选择。

这种组合算法的优势在于 NSGA-II 提供了解的多样性，而 TOPSIS 提供了从多目标中优中选优的能力。

应用场景

这种算法组合在实际中常用于复杂的多目标优化问题，比如供应链优化、投资组合优化、工程设计等，需要在多个冲突目标之间进行权衡时使用。

总结

pulp求解器是一种基于线性规划的优化工具，能够在确定性条件下为线性或整数规划问题找到全局最优解。与此相比，启发式算法如粒子群优化（PSO）、模拟退火（SA）、遗传算法（GA）等则是基于生物或物理学启发的非确定性方法，通常用于解决复杂的非线性或高维度问题。这些算法通过随机性和局部搜索探索解空间，逐步逼近最优解，虽然不一定保证全局最优，但在处理高复杂度问题时通常能提供有效的近似解。