第七章使用现代优化算法求解组合优化经典问题

在本章解决的问题都属于NP难问题，即无法在多项式时间复杂度内求解出来，只能用近似的算法或者我们即将使用的现代算法来进行求解，在以前的文章中我们也遇到了一个NP难问题，那就是旅行商问题TSP。在这里我们也会用现代算法来尝试求解TSP问题。

在开始前，请确保Python已经安装了scikit-opt。

pip install scikit-opt

本章我们不会深入的去探讨每个算法的实现原理，而是着重放在解决问题上。

scikit-opt

差分进化算法

R. Storn; K. Price. Differential Evolution - a Simple and Efficient Adaptive Scheme for Global Optimization over Continuous Spaces. Technical Report TR-95-012, ICSI, March 1995.

差分进化算法是一种后设启发式算法，不能够保证找到全局最优解。

差分进化算法之目的为求解优化问题，使用突变、交叉、选择计算以演化多个可能的解。首先，产生足量的随机变量，做为初始的可能解。接着，依序进行突变、交叉、选择计算，做完一轮后，检查某个终止条件。若终止条件尚未满足，则回到突变、交叉、选择计算，否则终止差分进化算法，输出最后一轮的最佳解。

例题简单的数学规划

既然这个算法能够求解优化问题，那么我们不妨让其从最简单的开始，这是一道非线性规划问题：

$\begin{align} &\max z = x_1^2+x_2^2+x_3^2\\ s.t.~~~&x_1x_2\ge1\\ &x_1x_2\le5\\ &x_1+x_2+x_3 = 1\\ &0\le x_1,x_2,x_3\le5 \end{align}$

import sko.DE as DE

# 定义目标函数
def func(p):
    x1, x2, x3 = p
    return x1 ** 2 + x2 ** 2 + x3 ** 2

# 定义约束条件，为等式约束
constraint_eq = [
    lambda x: 1 - x[0] - x[1] - x[2]
]

# 定义约束条件，为不等式约束
constraint_ueq = [
    lambda x: 1 - x[0] * x[1],
    lambda x: x[0] * x[1] - 5
]

# 调用DE算法，DE是差分进化算法
de = DE.DE(func=func, n_dim=3, size_pop=50, max_iter=2000, lb=[0, 0, 0], ub=[
           1, 1, 1], constraint_eq=constraint_eq, constraint_ueq=constraint_ueq)

# 运行算法
best_x, best_y = de.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)

我们得到的最终结果为：

best_x: [4.99999723e-01 5.00000277e-01 9.95851947e-14] 
 best_y: [0.5]

也就是说我们取到的优化结果为 $z=0.5$ ，其中 $x_1=0.499,x_2=0.501,x_3 = 9.95\times10^{-14}$

遗传算法

遗传算法是一种模拟自然进化过程的搜索优化方法，它借鉴了生物学中的遗传、变异、选择和杂交等现象，从一个初始的解集合（称为种群）出发，通过不断地产生新的解（称为个体），并根据适应度函数评价解的优劣，从而使种群向更好的方向进化，最终得到一个或多个优化解。遗传算法可以用于解决复杂的组合优化、函数优化、机器学习、人工智能等领域的问题。

多维背包问题MKP

多维背包问题是一种组合优化的问题，它涉及到一个背包和一些物品，每个物品有多个属性，比如重量、体积、价值等。背包也有多个限制条件，比如最大重量、最大体积等。问题的目标是选择一些物品放入背包，使得物品的总价值最大，同时满足背包的限制条件。

在之前学习动态规划，我们就涉及过背包问题：

数学建模系列(3)——数学规划模型 | LittleAO学习小站

当时只是简单的分析了一下，接下来我们会用差分进化算法来尝试求解多维背包问题，本次的问题约束分别为重量和体积。

我们给定背包的约束和所有的物品，如下面的代码所示：

# 背包约束
knapsack_const = {'weight': 600, 'volume': 600}
# 商品信息, 商品有不同的属性，如重量weight、体积volume，且有各自的价值value
items_info = [{'name': 'item_0', 'value': 1898, 'weight': 45, 'volume': 30},
 {'name': 'item_1', 'value': 440, 'weight': 5, 'volume': 20},
 {'name': 'item_2', 'value': 22507, 'weight': 85, 'volume': 125},
 {'name': 'item_3', 'value': 270, 'weight': 150, 'volume': 5},
 {'name': 'item_4', 'value': 14148, 'weight': 65, 'volume': 80},
 {'name': 'item_5', 'value': 3100, 'weight': 95, 'volume': 25},
 {'name': 'item_6', 'value': 4650, 'weight': 30, 'volume': 35},
 {'name': 'item_7', 'value': 30800, 'weight': 12, 'volume': 73},
 {'name': 'item_8', 'value': 615, 'weight': 170, 'volume': 12},
 {'name': 'item_9', 'value': 4975, 'weight': 20, 'volume': 15},
 {'name': 'item_10', 'value': 1160, 'weight': 40, 'volume': 15},
 {'name': 'item_11', 'value': 4225, 'weight': 25, 'volume': 40},
 {'name': 'item_12', 'value': 510, 'weight': 20, 'volume': 5},
 {'name': 'item_13', 'value': 11880, 'weight': 3, 'volume': 10},
 {'name': 'item_14', 'value': 479, 'weight': 7, 'volume': 10},
 {'name': 'item_15', 'value': 440, 'weight': 25, 'volume': 12},
 {'name': 'item_16', 'value': 490, 'weight': 12, 'volume': 10},
 {'name': 'item_17', 'value': 330, 'weight': 22, 'volume': 9},
 {'name': 'item_18', 'value': 110, 'weight': 25, 'volume': 10},
 {'name': 'item_19', 'value': 560, 'weight': 9, 'volume': 20},
 {'name': 'item_20', 'value': 24355, 'weight': 165, 'volume': 60},
 {'name': 'item_21', 'value': 2885, 'weight': 2, 'volume': 40},
 {'name': 'item_22', 'value': 11748, 'weight': 85, 'volume': 50},
 {'name': 'item_23', 'value': 4550, 'weight': 15, 'volume': 36},
 {'name': 'item_24', 'value': 750, 'weight': 9, 'volume': 49},
 {'name': 'item_25', 'value': 3720, 'weight': 2, 'volume': 40},
 {'name': 'item_26', 'value': 1950, 'weight': 4, 'volume': 19},
 {'name': 'item_27', 'value': 10500, 'weight': 100, 'volume': 150}]

我们可以写出这类问题的数学表达式：

$\begin{aligned} &\max z = \sum\limits^n_{j=1}c_jx_j \\ s.t.& \sum\limits^n_{j=1}a_{ij}x_j\le b_i,~~&i\in M=\{1,2,...,m\}\\ & x_j\in\{0,1\}, &j\in N=\{1,2,...,n\} \end{aligned}$

符号	含义
$m$	背包约束的个数
$n$	物品的个数
$b_i$	背包约束 $i$
$a_{ij}$	物品 $j$ 对背包的约束 $i$ 的消耗
$c_i$	物品的价值
$x_j$	决策变量， $x_j=1$ 表示物品 $j$ 放入背包； $x_j=0$ 表示物品 $j$ 未放入背包。

代码实现如下：

# 差分进化算法求解MKP
from sko.DE import DE
from sko.GA import GA
import matplotlib.pyplot as plt
import time

# 读取数据
# 背包约束
knapsack_const = {'weight': 600, 'volume': 600}
# 商品信息, 商品有不同的属性，如重量weight、体积volume，且有各自的价值value
items_info = [{'name': 'item_0', 'value': 1898, 'weight': 45, 'volume': 30},
              {'name': 'item_1', 'value': 440, 'weight': 5, 'volume': 20},
              {'name': 'item_2', 'value': 22507, 'weight': 85, 'volume': 125},
              {'name': 'item_3', 'value': 270, 'weight': 150, 'volume': 5},
              {'name': 'item_4', 'value': 14148, 'weight': 65, 'volume': 80},
              {'name': 'item_5', 'value': 3100, 'weight': 95, 'volume': 25},
              {'name': 'item_6', 'value': 4650, 'weight': 30, 'volume': 35},
              {'name': 'item_7', 'value': 30800, 'weight': 12, 'volume': 73},
              {'name': 'item_8', 'value': 615, 'weight': 170, 'volume': 12},
              {'name': 'item_9', 'value': 4975, 'weight': 20, 'volume': 15},
              {'name': 'item_10', 'value': 1160, 'weight': 40, 'volume': 15},
              {'name': 'item_11', 'value': 4225, 'weight': 25, 'volume': 40},
              {'name': 'item_12', 'value': 510, 'weight': 20, 'volume': 5},
              {'name': 'item_13', 'value': 11880, 'weight': 3, 'volume': 10},
              {'name': 'item_14', 'value': 479, 'weight': 7, 'volume': 10},
              {'name': 'item_15', 'value': 440, 'weight': 25, 'volume': 12},
              {'name': 'item_16', 'value': 490, 'weight': 12, 'volume': 10},
              {'name': 'item_17', 'value': 330, 'weight': 22, 'volume': 9},
              {'name': 'item_18', 'value': 110, 'weight': 25, 'volume': 10},
              {'name': 'item_19', 'value': 560, 'weight': 9, 'volume': 20},
              {'name': 'item_20', 'value': 24355, 'weight': 165, 'volume': 60},
              {'name': 'item_21', 'value': 2885, 'weight': 2, 'volume': 40},
              {'name': 'item_22', 'value': 11748, 'weight': 85, 'volume': 50},
              {'name': 'item_23', 'value': 4550, 'weight': 15, 'volume': 36},
              {'name': 'item_24', 'value': 750, 'weight': 9, 'volume': 49},
              {'name': 'item_25', 'value': 3720, 'weight': 2, 'volume': 40},
              {'name': 'item_26', 'value': 1950, 'weight': 4, 'volume': 19},
              {'name': 'item_27', 'value': 10500, 'weight': 100, 'volume': 150}]

# 目标函数

def func(x):
    # x为决策变量，是一个一维数组，长度为商品数量，每个元素的取值为0或1，表示是否选择该商品
    # 该函数返回值为目标函数值
    value = 0
    weight = 0
    volume = 0
    for i in range(28):
        if x[i] == 1:
            value += items_info[i]['value']
            weight += items_info[i]['weight']
            volume += items_info[i]['volume']
    if weight > knapsack_const['weight'] or volume > knapsack_const['volume']:
        value = 0
    return -value

# 不等式约束

def con(x):
    sum = 0
    for i in range(28):
        sum += x[i] * items_info[i]['weight']
    return sum - knapsack_const['weight']

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")

# 连续运行10次，取最优解
best_x = []
best_y = 0
start = time.time()
print(f'开始计算')
for i in range(15):
    ga = GA(func=func, n_dim=28, size_pop=80, max_iter=2000, lb=[
        0]*28, ub=[1]*28, constraint_eq=[con], precision=1, prob_mut=0.02)
    x, y = ga.run()
    print(f'第{i+1}次计算结果：{-y}, 共耗时{time.time()-start:.2f}')
    if y < best_y:
        best_x = x
        best_y = y
print('best_x:', best_x, '\n', 'best_y:', -best_y)

# 输出方案
weight = 0
volume = 0
value = 0
for i in range(28):
    if best_x[i] == 1:
        weight += items_info[i]['weight']
        volume += items_info[i]['volume']
        value += items_info[i]['value']
        print(items_info[i]['name'])
print('weight:', weight, 'volume:', volume, 'value:', value)

运行结果：

开始计算
第1次计算结果：[87470], 共耗时12.48
第2次计算结果：[110503], 共耗时24.77
第3次计算结果：[104905], 共耗时37.27
第4次计算结果：[127724], 共耗时50.11
第5次计算结果：[110288], 共耗时63.04
第6次计算结果：[87432], 共耗时78.55
第7次计算结果：[125978], 共耗时92.91
第8次计算结果：[114195], 共耗时107.31
第9次计算结果：[90532], 共耗时120.17
第10次计算结果：[134945], 共耗时132.67
第11次计算结果：[115270], 共耗时146.43
第12次计算结果：[126319], 共耗时160.79
第13次计算结果：[117359], 共耗时174.22
第14次计算结果：[112465], 共耗时186.81
第15次计算结果：[88326], 共耗时198.96
best_x: [1. 0. 1. 0. 1. 0. 1. 1. 0. 1. 0. 1. 1. 1. 1. 0. 1. 1. 0. 0. 1. 0. 1. 0.
 0. 0. 1. 0.]
 best_y: [134945]
item_0
item_2
item_4
item_6
item_7
item_9
item_11
item_12
item_13
item_14
item_16
item_17
item_20
item_22
item_26
weight: 600 volume: 571 value: 134945

请注意，我们得到的结果不一定是全局最优解，只能无限逼近于最优解。

模拟退火算法

模拟退火算法是一种模拟自然界中固体物质的退火过程的优化算法，它利用温度参数和概率因素来控制搜索过程，从而能够在解空间中跳出局部最优解，寻找全局最优解。它的基本思想是：

从一个较高的初始温度开始，随机生成一个初始解，并计算其目标函数值。
在每个温度下，重复多次以下步骤：
- 对当前解施加一个小的随机扰动，得到一个新解，并计算其目标函数值。
- 如果新解比当前解更优，或者满足一定的概率条件（Metropolis准则），则接受新解作为当前解，否则保持当前解不变。
降低温度参数，直到达到预设的终止条件（如最小温度或最大迭代次数）。

模拟退火算法的优点是能够有效地避免陷入局部最优解，适用于各种复杂的组合优化问题。

免疫优化算法

免疫优化算法是一种模仿生物免疫系统的智能优化方法，它可以解决一些复杂的优化问题，如旅行商问题、多目标优化问题等。免疫优化算法的基本思想是利用抗体和抗原之间的亲和力和排斥力来选择和更新解，从而达到全局或局部最优。

旅行商问题TSP

想要学会车辆路径问题VRP，那么一定要先了解TSP。VRP是TSP的另一种类型。

在之前我们已经学习过TSP问题，使用了动态规划的思想，还初步接触了遗传算法。接下来我们要使用模拟退火算法来实现TSP问题：

数学建模系列(5)——图论和网格优化 | LittleAO学习小站

from sko.IA import IA_TSP
from matplotlib.ticker import FormatStrFormatter
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt

num_points = 30

points_coordinate = np.random.rand(
    num_points, 2)  # generate coordinate of points
distance_matrix = spatial.distance.cdist(
    points_coordinate, points_coordinate, metric='euclidean')

def cal_total_distance(routine):
    '''The objective function. input routine, return total distance.
    cal_total_distance(np.arange(num_points))
    '''
    num_points, = routine.shape
    return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])

ia_tsp = IA_TSP(func=cal_total_distance, n_dim=num_points,
                size_pop=50, max_iter=50, prob_mut=0.2)

# 每迭代一次就会记录绘图一次
for i in range(500):
    best_points, best_distance = ia_tsp.run()
    best_points_ = np.concatenate([best_points, [best_points[0]]])
    best_points_coordinate = points_coordinate[best_points_, :]
    plt.subplot(1, 2, 1)
    plt.plot(ia_tsp.generation_best_Y)
    plt.subplot(1, 2, 2)
    plt.plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1],
             marker='o', markerfacecolor='b', color='c', linestyle='-')
    plt.pause(0.01)
    plt.clf()

Untitled

车辆路径问题VRP

车辆路径问题（Vehicle Routing Problem，VRP）是一类经典的组合优化问题，涉及如何有效地分配一组车辆去访问多个客户点，并在满足约束条件的情况下最小化总行驶距离或成本。

给定一组客户点、车辆容量、车辆数量、起始点和终点，目标是找到使得所有客户点都被访问一次的最短路径方案。车辆路径问题的数学模型可以使用整数线性规划来表示。

$\begin{align} &\min z = \sum\limits^m_{k=1}\sum\limits^n_{i=0}\sum\limits^n_{j=0}c_{ij}\cdot x_{ij}^k \\ s.t.& \sum\limits^n_{j=1}x^k_{0j} = 1,~~k\in\{1,2,...,m\}\\ & \sum\limits^n_{j=1}x^k_{i0} = 1,~~k\in\{1,2,...,m\}\\ & \sum\limits^m_{k=1}\sum\limits^n_{j=0}x^k_{ij}=1,~~i\in\{1,2,...,n\}\\ &x^k_{ij}\in\{0,1\} \end{align}$

符号及式子解释：

符号	含义
$m$	车辆的数量
$n$	客户点的数量
$c_{ij}$	客户点 $i$ 到 $j$ 的距离
$x^k_{ij}$	决策变量，车辆 $k$ 在访问客户点 $i$ 后移动到 $j$ 。 $x^k_{ij}=1$ 为确认移动，否则为0。

式子	含义
(6)	总路径最小
(7)	车辆要经过起始点
(8)	车辆要经过终点
(9)	每个点只能被访问1次
(10)	0-1规划

目录CONTENT

数学建模系列(7)——现代优化算法及组合优化经典问题

第七章使用现代优化算法求解组合优化经典问题

差分进化算法

例题简单的数学规划

遗传算法

多维背包问题MKP

模拟退火算法

免疫优化算法

旅行商问题TSP

车辆路径问题VRP

评论区

数学建模系列(7)——现代优化算法及组合优化经典问题

第七章 使用现代优化算法求解组合优化经典问题

差分进化算法

例题 简单的数学规划

遗传算法

多维背包问题MKP

模拟退火算法

免疫优化算法

旅行商问题TSP

车辆路径问题VRP

评论区

第七章使用现代优化算法求解组合优化经典问题

例题简单的数学规划