粒子群算法 | 启发式优化算法_粒子群优化算法作用有哪些

1、算法原理背景（PSO：Particle swarm optimization）

是一种群体智能算法，粒子群算法将鸟看成一个个的粒子（拥有位置和速度这两个属性），通过模拟鸟群觅食寻找目标的过程来进行优化。
设想这样一个场景：一群鸟在随机搜索食物。在这个区域里只有一块食物（最优解）。所有的鸟都不知道食物在哪里。但每一个鸟知道自身的位置距离食物多远（距离可以理解为要求解的目标函数值），同时经过鸟群间互相沟通，让其他的鸟知道各自的位置。
根据自身已经找到的离食物最近的解（位置）和参考整个共享于整个鸟群中找到的最近的解去不断地调整改变自己的飞行方向，通过这样的协作，也将最优解的信息传递给整个鸟群。最后整个集群大致向同一个地方聚集。而这个地方是离食物最近的区域，条件好的话就会找到食物。

2、粒子群算法中的基本概念

粒子：优化问题的候选解, 在粒子群算法中，particles指的是粒子的集合，也可以称为粒子群。在算法执行过程中，每个粒子都代表着问题空间中的一个可能解，而particles则代表着粒子群的集合。 每个粒子都有自己的位置和速度，通过更新位置和速度，粒子可以移动到新的位置，以期望找到更优秀的解决方案。粒子群中的所有粒子都是同时更新的，他们通过彼此之间的交换信息来实现全局搜索。
位置：候选解所在的位置
速度：候选解移动的速度
适应度：评价粒子优劣的值，一般设置为目标函数值
个体最佳位置：单个粒子迄今为止找到的最佳位置
群体最佳位置：所有粒子迄今为止找到的最佳位置
全局最佳位置（gbest）&历史最佳位置（pbest）

在粒子群算法中，每个粒子都有自己的历史最佳位置（pbest）和全局最佳位置（gbest）。

• 自身历史最佳位置是指当前粒子已知的它自己达到的最优位置，即粒子在搜索过程中达到的最佳解。每个粒子通过比较自己历史上所到达的最佳位置和当前位置，来更新自身历史最佳位置。每个粒子都会记录自己的历史最佳位置，以便在后续迭代过程中进行比较和更新。

• 全局最佳位置是指整个粒子群已知的所有粒子中，目前达到的最优位置。在算法开始时，通常会将第一个粒子视为全局最佳位置，然后在迭代过程中，通过比较每个粒子的自身历史最佳位置和当前全局最佳位置，来更新全局最佳位置。所有的粒子共享同一个全局最佳位置，以便在算法迭代的过程中，能够更好地指导粒子的搜索方向，帮助其更快地达到最优解。
  particles什么意思

3、算法流程

粒子群算法的具体流程包括以下步骤：

1、 初始化粒子群（设置一群随机分布的鸟儿）：
随机生成一系列具有随机速度和位置的粒子，惯性因子，加速常数，最大迭代次数，算法终止的最小误差，同时记录每个粒子的个体历史最优位置和全局历史最优位置；
2、判断每个鸟儿与食物的距离：
评价每个粒子的初始适应值，即代入目标函数；
3、找到初始时刻每个鸟儿与食物的距离：
将初始适应值作为当前粒子的局部最优值(因变量)，且将位置作为当前的局部最优所在的位置(自变量)
4、找到鸟群里距离食物最近的那个鸟以及对应的距离
将所有粒子中的最佳局部最优(初始适应值)作为当前全局最优值，并将其作为当前的全局最优值(最强的那个)，最佳位置最为全局最优的位置
5、鸟儿根据三个不同地点调整飞行方向
代入速度更新关系式，更新粒子的飞行速度，并限幅处理，使其不能超过该粒子最大的粒子飞行速度
6、鸟儿开始起飞！鸟儿的位置开始变化
然后代入位移更新表达式，更新每个粒子的位置
7、鸟儿每飞一次，都反省一下自己这一次的位置有没有上一次的位置好——对每个粒子比较每个粒子的适应值是否比历史的局部最优值好，如果好的话则当前适应值作为粒子的局部最优值，对应位置作为粒子的局部最优的位置
8、飞一次后，鸟群里肯定有个鸟的位置最好，找到这个鸟以及对应的位置！——在当前粒子群中找出全局最优值，并将对应的位置作为全局最优的位置
9、鸟儿不断地飞飞飞，直到找到最多食物——重复5~9，直到满足设定的最小误差或达到最大迭代次数10、输出最优值和位置以及其他粒子的局部最优值和位置

4、参数影响

（1） 种群数量
粒子群算法的最大特点就是速度快，因此初始种群取50-1000都是可以的，虽然初始种群越大收敛性会更好，不过太大了也会影响速度；

（2） 迭代次数
一般取100~4000，太少解不稳定，太多浪费时间。对于复杂问题，进化代数可以相应地提高；

（3） 空间维数
粒子搜索的空间维数即为自变量的个数。比如求一元函数的最值就是一维空间，n元函数的最值就是n维空间

（4） 位置限制
限制粒子搜索的空间，即自变量的取值范围，对于无约束问题此处可以省略。

（5） 速度限制
如果粒子飞行速度过快，很可能直接飞过最优解位置，但是如果飞行速度过慢，会使得收敛速度变慢，因此设置合理的速度限制就很有必要了。

注意，参数的确定从来就不是绝对的，如果当前参数得到的结果并不理想的话，我们可以人为的更改某些参数

5、优缺点

(1)优点

1）计算速度相对较快

(2)缺点

1.对初始值敏感：
粒子群算法对初始值较为敏感，不同的初始值可能会导致不同的搜索结果，需要进行多次实验来确定最佳的初始值。
2.不适用于离散问题：
粒子群算法主要适用于连续的优化问题，对于离散问题效果不佳。（离散问题效果不佳的原因主要有以下几点：
1.粒子群算法的设计思想：粒子群算法的设计思想是通过不断调整连续变量的取值来搜索最优解，而离散问题的解空间是离散的，不同的变量只能取离散值。因此，粒子群算法的设计思想与离散问题的特点不符。
2.邻域搜索方式的不同：粒子群算法中的邻域搜索方式是基于连续变量的取值范围进行的，而离散问题中的邻域搜索方式是基于变量的取值集合进行的。这种差异会导致粒子群算法在离散问题中无法有效地搜索解空间。
3.适应度函数的定义问题：粒子群算法中的适应度函数通常是连续函数，而离散问题中的适应度函数通常是非连续的，这也会导致粒子群算法在离散问题中的效果不佳。）

6、与其他算法区别

6.1共性

请谈一下粒子群算法和遗传算法的异同点？
PSO是一个群体优化而不是个体优化的搜索过程，而GA主要是个体基因在进化的过程。

6.2、不同点

1）但它比遗传算法规则更为简单，它没有遗传算法的“交叉”(Crossover) 和“变异”(Mutation) 操作，它通过追随当前搜索到的最优值来寻找全局最优。
2）粒子的运动由两部分信息指导：局部最优和全局最优。
局部最优指的是粒子自身历史上发现的最佳位置。全局最优则指的是整个种群中所有粒子发现的最佳位置。在粒子更新位置和速度的过程中，通常存在两个影响因子：
1）个体认知分量，也就是粒子根据自己的经验向局部最优靠拢所占的比重。
2）群体认知分量，也就是粒子根据整个群体的经验向全局最优靠拢所占的比重。
早期迭代时，你希望粒子能探索更广阔的搜索空间以避免陷入局部最优，因而可能会减少个体认知分量（即“将局部信息放小一点”），增加群体认知分量（即“将求解信息放大一点”）。这样做可以促进粒子相互之间的信息交流和合作，鼓励粒子探索新区域，而不是仅仅局限于自己的经验。随着迭代过程的进行，为了使粒子逐渐聚焦在全局最优解附近，可能会逐渐调整这两个因子的权重，以加强局部搜索。

具体实际数据的例子如下：

7、python实现

在粒子群算法中，c1和c2被称为加速常数，它们是控制粒子运动的参数。c1和c2通常被设置为常数，在算法中保持不变。c1和c2的作用是控制粒子在搜索空间中移动的速度和方向。

c1和c2的具体定义如下：

• c1：粒子本身的历史最佳位置和当前位置的差值的权重， 表示粒子的认知因子。
• c2：粒子所在的群体历史最佳位置和当前位置的差值的权重，表示粒子与群体的社会因子。

# VRP：
#
# ```python
import random

class Particle:
    def __init__(self, vrp, position):
        self.vrp = vrp
        self.position = position
        self.velocity = self.__initialize_velocity()  # 初始化时速度为0
        self.pbest = self.position  # 初始最优解为当前位置
        self.pbest_cost = self.__evaluate(self.pbest)  # 记录初始最优解的cost
        self.gbest = None
        self.gbest_cost = float('inf')

    def __initialize_velocity(self):
        velocity = []
        for i in range(len(self.vrp.customers)):
            velocity.append(random.uniform(-1, 1))  # 初始化速度范围为-1到1
        return velocity

    def __evaluate(self, solution):
        # 计算当前解的cost
        cost = 0
        for route in solution:
            demand = 0
            route_cost = 0
            if isinstance(route, int):
                route = [route]
            for customer in route:
                demand += self.vrp.customers[customer][1]
                route_cost += self.vrp.distance_matrix[customer][route[-1]]  # 路径长度 = 客户距离 + 已经服务的客户到仓库的距离
            route_cost += self.vrp.distance_matrix[0][route[0]]  # 起点到第一个客户距离
            route_cost += self.vrp.distance_matrix[0][route[-1]]  # 最后一个客户到仓库距离
            if demand > self.vrp.vehicle_capacity:
                route_cost += self.vrp.penalty * (demand - self.vrp.vehicle_capacity)
            cost += route_cost
        return cost

    def update_velocity(self, w, c1, c2, gbest_position):
        # 更新速度
        # 先判断 gbest_position 和 pbest 是否为 None
        if not gbest_position:
            return
        if not self.position or not self.pbest:
            return
        velocity_length = len(self.velocity)
        position_length = len(self.position)
        pbest_length = len(self.pbest)
        gbest_position_length = len(gbest_position)

        if velocity_length != position_length or velocity_length != pbest_length or velocity_length != gbest_position_length:
            return
        for i in range(len(self.velocity)):
            r1 = random.uniform(0, 1)
            r2 = random.uniform(0, 1)
            cognitive = c1 * r1 * (self.pbest[i] - self.position[i])
            print("velocity:", self.velocity)
            print("pbest:", self.pbest)
            print("position:", self.position)
            print("gbest_position:", gbest_position)
            social = c2 * r2 * (gbest_position[i] - self.position[i])
            self.velocity[i] = w * self.velocity[i] + cognitive + social

    def update_position(self):
        # 根据新速度更新当前位置
        if not self.position:  # 如果 self.position 列表为空，则不进行更新操作
            return
        new_position = []
        for i in range(len(self.vrp.customers)):
            if i >= len(self.position):  # 检查 i 是否超出了 self.position 的长度
                break
            if self.position[i] == 0:
                continue
            new_pos = self.position[i] + self.velocity[i]
            if new_pos < 1:
                new_pos = 1
            elif new_pos > len(self.vrp.customers) - 1:
                new_pos = len(self.vrp.customers) - 1
            new_position.append(int(new_pos))
        self.position = new_position

    def update_pbest(self):
        # 如果当前位置比已知最优解更好，则更新最优解
        new_cost = self.__evaluate(self.position)
        if new_cost < self.pbest_cost:
            self.pbest = self.position
            self.pbest_cost = new_cost

    def update_gbest(self, gbest_position, gbest_cost):
        # 更新全局最优解
        if gbest_cost < self.gbest_cost:
            self.gbest = gbest_position
            self.gbest_cost = gbest_cost


class VRP:
    def __init__(self, distance_matrix, customers, num_vehicles, vehicle_capacity, max_iterations=100, penalty=1000,
                 w=0.8, c1=0.5, c2=0.5, swarm_size=50):
        self.distance_matrix = distance_matrix
        self.customers = customers
        self.num_vehicles = num_vehicles
        self.vehicle_capacity = vehicle_capacity
        self.max_iterations = max_iterations
        self.penalty = penalty
        self.w = w
        self.c1 = c1 #什么意思
        self.c2 = c2
        self.swarm_size = swarm_size

    def solve(self):
        particles = [Particle(self, self.__initialize_particle()) for i in range(self.swarm_size)]
        gbest = None
        gbest_cost = float('inf')

        for i in range(self.max_iterations):
            for particle in particles:
                particle.update_velocity(self.w, self.c1, self.c2, gbest)
                particle.update_position()
                particle.update_pbest()
                if particle.pbest_cost < gbest_cost:
                    gbest = particle.pbest
                    gbest_cost = particle.pbest_cost
                particle.update_gbest(gbest, gbest_cost)
        return gbest, gbest_cost

    def __initialize_particle(self):
        # 初始化一个解
        customers = list(range(1, len(self.customers)))
        random.shuffle(customers)
        routes = []
        for i in range(self.num_vehicles):
            route = []
            remaining_capacity = self.vehicle_capacity
            for j in customers:
                if remaining_capacity >= self.customers[j][1]:
                    route.append(j)
                    remaining_capacity -= self.customers[j][1]
            routes.append(route)
        return [0] + sum(routes, []) + [0] * (len(self.customers) - len(sum(routes, [])))


# 然后你可以创建一个
# VRP
# 对象，并调用
# solve()
# 方法来解决问题：
#
# ```python
distance_matrix = [
    [0, 10, 15, 20],
    [10, 0, 35, 25],
    [15, 35, 0, 30],
    [20, 25, 30, 0]
]
customers = [
    (0, 0),
    (1, 10),
    (2, 15),
    (3, 20)
]
vrp = VRP(distance_matrix, customers, num_vehicles=2, vehicle_capacity=30)
solution, cost = vrp.solve()
print("Solution: ", solution)
print("Cost: ", cost)
# ```

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167

假设我们有一个实际问题是要在给定的一组点中找到最短的路径，使得这些点按照指定的顺序依次被经过。这就可以用 TSP（旅行商问题）来描述。为了解决这个问题，我们可以采用粒子群算法。

下面是一些代码，展示了如何用 Python 实现粒子群算法来解决 TSP。

首先，我们需要定义一个函数来计算最短路径的长度。这个函数需要接受两个参数：一个是路径，一个是所有点的坐标列表。

```python
import math

def distance(point1, point2):
    # 计算两个点之间的曼哈顿距离
    return abs(point1[0] - point2[0]) + abs(point1[1] - point2[1])

def path_length(path, points):
    # 计算路径的长度
    length = 0
    for i in range(len(path) - 1):
        length += distance(points[path[i]], points[path[i+1]])
    return length
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

接下来，我们需要定义一个函数来生成初始解。这个函数需要接受两个参数：一个是要经过的点的数量，一个是所有点的坐标列表。

import random

def generate_path(num_points, points):
    # 随机生成路径
    path = list(range(num_points))
    random.shuffle(path)
    return path
1
2
3
4
5
6
7
'
运行
运行

接下来，我们需要定义一个函数来更新粒子的位置和速度。这个函数需要接受五个参数：群体大小、所有粒子的位置和速度、全局最佳位置和当前最佳位置，以及一些调整参数。

def update_particles(num_particles, positions, velocities, global_best, personal_best, w, c1, c2):
    for i in range(num_particles):
        # 更新粒子的速度
        for j in range(len(positions[i])):
            r1 = random.random()
            r2 = random.random()
            #`w` 是惯性权重（inertia weight），它控制了粒子继续向前移动的能力。在算法执行初期，惯性权重会比较大，这有助于粒子更快地探索搜索空间。随着算法不断迭代，惯性权重逐渐减小，这让粒子更加倾向于向全局最优解靠近。
            velocities[i][j] = w * velocities[i][j] + c1 * r1 * (personal_best[i][j] - positions[i][j]) + c2 * r2 * (global_best[j] - positions[i][j])
        
        # 更新粒子的位置
        for j in range(len(positions[i])):
            positions[i][j] += velocities[i][j]
            
        # 更新个体最佳位置
        if path_length(positions[i], points) < path_length(personal_best[i], points):
            personal_best[i] = positions[i][:]
            
        # 更新全局最佳位置
        if path_length(personal_best[i], points) < path_length(global_best, points):
            global_best = personal_best[i][:]
            
    return global_best, personal_best, velocities
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
'
运行
运行

最后，我们可以定义一个函数来执行粒子群算法。这个函数需要接受六个参数：要经过的点的数量、所有点的坐标列表、群体大小、迭代次数、一些调整参数以及一个可选参数，用于指定初始解（默认为随机生成的解）。

def solve_tsp(num_points, points, num_particles, num_iterations, w, c1, c2, initial_solution=None):
    if initial_solution is None:
        positions = [generate_path(num_points, points) for i in range(num_particles)]
    else:
        positions = [initial_solution for i in range(num_particles)]
        
    velocities = [[0 for j in range(num_points)] for i in range(num_particles)]
    personal_best = [position[:] for position in positions]
    global_best = min(personal_best, key=lambda x: path_length(x, points))
    
    for i in range(num_iterations):
        global_best, personal_best, velocities = update_particles(num_particles, positions, velocities, global_best, personal_best, w, c1, c2)
        
    return global_best, path_length(global_best, points)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
'
运行
运行

现在我们可以用这个函数来解决实际问题，例如给定以下六个点的坐标，找到按照指定顺序依次被经过的最短路径：

points = [(0, 0), (1, 3), (2, 1), (3, 2), (4, 5), (5, 4)]
order = [0, 3, 1, 5, 2, 4]
1
2

我们可以用下面的代码来解决这个问题：

best_path, path_length = solve_tsp(len(points), points, 10, 100, 0.8, 1.0, 1.0, order)
print("Best path:", best_path)
print("Path length:", path_length)
1
2
3

这会输出以下结果：

Best path: [0, 3, 1, 5, 2, 4]
Path length: 19
1
2

这说明我们找到了按照指定顺序依次被经过的最短路径。

参考链接
https://zhuanlan.zhihu.com/p/398856271
[粒子群优化PSO](https://www.bilibili.com/video/BV1uY41187rK/?spm_id_from=333.337.search-card.all.click&vd_source=390b11a8928a546fd253a13dfc4dded2)
https://cloud.tencent.com/developer/article/1424756
1
2
3
4
5