今年突然变成分赛区 (10%) 推国，
国奖结果还没出，感觉一等有点悬，论文写的太一般了我没时间去修。
结果出啦，我们是一等奖（2%）！！！很感谢评委们对 idea 的认可，也很庆幸这次没有被论文拖住。
碎碎念：4 月又被拉着不务正业打了次比赛，刚好这几天有闲暇，传一下之前写的解题思路，不过时间太久了，就不重新花时间整理完整代码了，只分享核心代码。这题的解题难点不在代码部分，理清思路使用简单的启发式算法应该就可以获得靠前的奖项。
以下部分是我写给队友的思路的修改版本，在看的时候可以将自己代入进行理解，希望对你们有所帮助。
代码部分的模拟退火算法是有冗余的，没编写复用代码，基本直接复制过去修改候选解，这里指出来方便对比差异。
原题目看采用的图片是出自于 18 年的一篇论文：PCI Planning Based on Binary Quadratic Programming in LTE/LTE-A Networks，发现这一点的时候我也很惊讶，这篇论文应该也提供了一种解决思路（不过因为当时有想法就没有去细看，如果出于学习目的的话可以去通读）
初步梳理
- 冲突MR数 - 如果两个小区被分配了相同的 PCI 且是邻区关系，会产生冲突。
- 干扰MR数 - 类似于冲突，如果两个同频小区的 PCI 模 3 值相同，会产生干扰。（不同频默认 0）
- 混淆MR数 - 如果两个小区被分配相同的 PCI，并且它们都是第三个小区的邻区，则会混淆。
符号定义
- a i , j a_{i,j} ai,j: 小区 i i i 和小区 j j j 之间的冲突MR数。
- b i , j b_{i,j} bi,j: 小区 i i i 和小区 j j j 之间的混淆MR数。
- c i , j c_{i,j} ci,j: 小区 i i i 和小区 j j j 之间的模3干扰MR数。
- x i , k x_{i,k} xi,k: 二进制变量，若小区 i i i 被分配了PCI k k k 则为1，否则为0。
- y i , j y_{i,j} yi,j: 二进制变量，表示小区 i i i 和小区 j j j 是否被分配了相同的PCI。
- z i , j z_{i,j} zi,j: 二进制变量，表示小区 i i i 和小区 j j j 的PCI值在模3的结果是否相同。
题目1
给这2067个小区重新分配 PCI，使得这 2067 个小区之间的冲突MR 数、混淆 MR 数和模 3 干扰 MR 数的总和最少。
建模
目标：使得这 2067 个小区之间的冲突 MR 数、混淆 MR 数和模 3 干扰 MR 数的总和最少。
计算描述：
- 若小区 i i i 和 j j j 分配相同的 PCI 值，则冲突数增加 a i , j + a j , i a_{i, j} + a_{j, i} ai,j+aj,i。
- 若小区 i i i 和 j j j 为第三方小区的共同邻区且被分配相同的 PCI，则混淆数增加 b i j + b j i b_{ij} + b_{ji} bij+bji。
- 若小区 i i i 和 j j j 分配的 PCI 值在模 3 上相同，则模 3 干扰数增加 c i j + c j i c_{ij} + c_{ji} cij+cji。
二进制决策变量
建模之前，定义一个二进制决策变量
x
i
,
k
x_{i,k}
xi,k，其中：
x
i
,
k
=
{
1
如果小区
i
被分配了PCI值
k
0
否则
x_{i,k} = \begin{cases} 1 & \text{如果小区 } i \text{ 被分配了PCI值 } k \\ 0 & \text{否则} \end{cases}
xi,k={10如果小区 i 被分配了PCI值 k否则
这里
i
i
i 索引小区，
k
k
k 索引可能的 PCI 值（0 到 1007）。
可以写出目标函数：
Minimize
F
o
b
j
=
∑
i
=
0
N
−
1
∑
j
=
i
+
1
N
−
1
∑
k
=
0
1007
∑
l
=
0
1007
(
(
a
i
,
j
+
a
j
,
i
)
⋅
x
i
,
k
⋅
x
j
,
k
+
(
b
i
,
j
+
b
j
,
i
)
⋅
x
i
,
k
⋅
x
j
,
k
+
(
c
i
,
j
+
c
j
,
i
)
⋅
1
(
k
m
o
d
3
=
l
m
o
d
3
)
⋅
x
i
,
k
⋅
x
j
,
l
)
\text{Minimize} \quad F_{obj} = \sum_{i=0}^{N-1} \sum_{j=i+1}^{N-1} \sum_{k=0}^{1007} \sum_{l=0}^{1007} \left( (a_{i,j} + a_{j,i}) \cdot x_{i,k} \cdot x_{j,k} + (b_{i,j} + b_{j,i}) \cdot x_{i,k} \cdot x_{j,k} + (c_{i,j} + c_{j,i}) \cdot \mathbf{1}_{(k \mod 3 = l \mod 3)} \cdot x_{i,k} \cdot x_{j,l} \right)
MinimizeFobj=i=0∑N−1j=i+1∑N−1k=0∑1007l=0∑1007((ai,j+aj,i)⋅xi,k⋅xj,k+(bi,j+bj,i)⋅xi,k⋅xj,k+(ci,j+cj,i)⋅1(kmod3=lmod3)⋅xi,k⋅xj,l)
公式 1 ( k m o d 3 = l m o d 3 ) \mathbf{1}_{(k \mod 3 = l \mod 3)} 1(kmod3=lmod3)代表一个指示函数。
指示函数定义：
1
(
k
m
o
d
3
=
l
m
o
d
3
)
=
{
1
如果
k
m
o
d
3
=
l
m
o
d
3
0
否则
\mathbf{1}_{(k \mod 3 = l \mod 3)} = \begin{cases} 1 & \text{如果 } k \mod 3 = l \mod 3 \\ 0 & \text{否则} \end{cases}
1(kmod3=lmod3)={10如果 kmod3=lmod3否则
约束条件
-
每个小区分配一个PCI：
- ∑ k = 0 1007 x i , k = 1 , ∀ i \sum_{k=0}^{1007} x_{i,k} = 1, \quad \forall i k=0∑1007xi,k=1,∀i
- 每个小区必须分配且只分配一个 PCI 值。
-
二进制决策变量约束：
- x i , k ∈ { 0 , 1 } , ∀ i , k x_{i,k} \in \{0, 1\}, \quad \forall i, k xi,k∈{0,1},∀i,k，整数规划约束，确保决策变量是二进制的。
引入 y i , j y_{i,j} yi,j 和 z i , j z_{i,j} zi,j
此时有：
Minimize
F
o
b
j
=
∑
i
=
0
N
−
1
∑
j
=
i
+
1
N
−
1
(
(
a
i
,
j
+
a
j
,
i
)
⋅
y
i
,
j
+
(
b
i
,
j
+
b
j
,
i
)
⋅
y
i
,
j
+
(
c
i
,
j
+
c
j
,
i
)
⋅
z
i
,
j
)
\text{Minimize} \quad F_{obj} = \sum_{i=0}^{N-1} \sum_{j=i+1}^{N-1} \left( (a_{i,j} + a_{j,i}) \cdot y_{i,j} + (b_{i,j} + b_{j,i}) \cdot y_{i,j} + (c_{i,j} + c_{j,i}) \cdot z_{i,j} \right)
MinimizeFobj=i=0∑N−1j=i+1∑N−1((ai,j+aj,i)⋅yi,j+(bi,j+bj,i)⋅yi,j+(ci,j+cj,i)⋅zi,j)
其中：
y
i
,
j
=
{
1
如果小区
i
和小区
j
分配了相同的 PCI
0
否则
y_{i,j} = \begin{cases} 1 & \text{如果小区 } i \text{ 和小区 } j \text{ 分配了相同的 PCI} \\ 0 & \text{否则} \end{cases}
yi,j={10如果小区 i 和小区 j 分配了相同的 PCI否则
z
i
,
j
=
{
1
如果小区
i
和小区
j
的 PCI 模3结果相同
0
否则
z_{i,j} = \begin{cases} 1 & \text{如果小区 } i \text{ 和小区 } j \text{ 的 PCI 模3结果相同} \\ 0 & \text{否则} \end{cases}
zi,j={10如果小区 i 和小区 j 的 PCI 模3结果相同否则
-
y
i
,
j
y_{i,j}
yi,j 定义：
- y i , j = ∑ k = 0 1007 x i , k ⋅ x j , k y_{i,j} = \sum_{k=0}^{1007}x_{i,k} \cdot x_{j,k} yi,j=k=0∑1007xi,k⋅xj,k
- 这表示小区 i i i 和 j j j 是否被分配了相同的PCI。
-
z
i
,
j
z_{i,j}
zi,j 定义：
- z i , j = ∑ k = 0 1007 ∑ l = 0 1007 ( 1 ( k m o d 3 = l m o d 3 ) ⋅ x i , k ⋅ x j , l ) z_{i,j} = \sum_{k=0}^{1007} \sum_{l=0}^{1007} \left( \mathbf{1}_{(k \mod 3 = l \mod 3)} \cdot x_{i,k} \cdot x_{j,l} \right) zi,j=∑k=01007∑l=01007(1(kmod3=lmod3)⋅xi,k⋅xj,l)，这表示小区 i i i 和 j j j 被分配的PCI模3后的值是否相同
使用辅助变量 u i , j , k u_{i,j,k} ui,j,k 和 v i , j , k , l v_{i,j,k,l} vi,j,k,l处理非线性情况
这部分的处理和2023MathorCup A题的QUBO模型三次转二次的过程基本一致
为了将非线性的乘积条件转化为线性规划模型能够处理的形式，我们引入辅助变量。
-
y
i
,
j
y_{i,j}
yi,j 的重新定义：
- y i , j = ∑ k = 0 1007 u i , j , k y_{i,j} = \sum_{k=0}^{1007} u_{i,j,k} yi,j=k=0∑1007ui,j,k
-
z
i
,
j
z_{i,j}
zi,j 的重新定义：
- z i , j = ∑ k = 0 1007 ∑ l = 0 1007 v i , j , k , l （ i f k m o d 3 = l m o d 3 ） z_{i,j} = \sum_{k=0}^{1007} \sum_{l=0}^{1007} v_{i,j,k,l}（if \ k \mod 3 = l \mod 3） zi,j=k=0∑1007l=0∑1007vi,j,k,l（if kmod3=lmod3）
约束条件
-
对于每对小区 i , j i, j i,j 和每个可能的PCI值 k k k，定义以下 u i , j , k u_{i,j,k} ui,j,k 约束：
- u i , j , k ≤ x i , k u_{i,j,k} \leq x_{i,k} ui,j,k≤xi,k
- u i , j , k ≤ x j , k u_{i,j,k} \leq x_{j,k} ui,j,k≤xj,k
- u i , j , k ≥ x i , k + x j , k − 1 u_{i,j,k} \geq x_{i,k} + x_{j,k} - 1 ui,j,k≥xi,k+xj,k−1
-
对于每对小区 i , j i, j i,j和每对可能的PCI值 k , l k, l k,l（如果 k m o d 3 = l m o d 3 k \mod 3 = l \mod 3 kmod3=lmod3），定义以下 v i , j , k , l v_{i,j,k,l} vi,j,k,l 约束：
- v i , j , k , l ≤ x i , k v_{i,j,k,l} \leq x_{i,k} vi,j,k,l≤xi,k
- v i , j , k , l ≤ x j , l v_{i,j,k,l} \leq x_{j,l} vi,j,k,l≤xj,l
- v i , j , k , l ≥ x i , k + x j , l − 1 v_{i,j,k,l} \geq x_{i,k} + x_{j,l} - 1 vi,j,k,l≥xi,k+xj,l−1（如果 k m o d 3 = l m o d 3 k \mod 3 = l \mod 3 kmod3=lmod3）
线性规划思路
（这里我没有放进论文，只是解题期间闲暇写给其他人理解的初始思路，可以结合代码中的拓展进行理解）
至此，可以写出线性规划的伪代码如下：
算法 1: PCI 分配问题的线性规划求解
输入:
- N N N: 小区的数量
- K K K: 可用PCI值的数量
- A , B , C A, B, C A,B,C: 分别为冲突、混淆和模3干扰的MR数矩阵
输出:
- 最优的PCI分配
- 目标函数的最小值
步骤:
-
初始化问题:
- 创建一个最小化目标的线性规划模型
-
定义决策变量:
- 对每个小区 i i i 和每个PCI k k k，定义二进制决策变量 x i , k x_{i,k} xi,k
- 对每对小区 i , j i, j i,j，定义二进制变量 y i , j y_{i,j} yi,j 和 z i , j z_{i,j} zi,j
- 定义辅助二进制变量 u i , j , k u_{i,j,k} ui,j,k 和 v i , j , k , l v_{i,j,k,l} vi,j,k,l
-
构建目标函数:
- 添加目标函数到模型： min ∑ i = 0 N − 1 ∑ j = i + 1 N − 1 ( ( a i , j + a j , i ) y i , j + ( b i , j + b j , i ) y i , j + ( c i , j + c j , i ) z i , j ) \min \sum_{i=0}^{N-1} \sum_{j=i+1}^{N-1} ((a_{i,j} + a_{j,i}) y_{i,j} + (b_{i,j} + b_{j,i}) y_{i,j} + (c_{i,j} + c_{j,i}) z_{i,j}) min∑i=0N−1∑j=i+1N−1((ai,j+aj,i)yi,j+(bi,j+bj,i)yi,j+(ci,j+cj,i)zi,j)
-
添加约束:
- 确保每个小区只被分配一个PCI： ∑ k = 0 K − 1 x i , k = 1 , ∀ i \sum_{k=0}^{K-1} x_{i,k} = 1, \quad \forall i ∑k=0K−1xi,k=1,∀i
- 为 y i , j y_{i,j} yi,j 和 z i , j z_{i,j} zi,j 添加线性化辅助变量的约束
-
求解模型:
- 使用线性规划求解器解决模型
-
输出结果:
- 如果找到最优解，输出PCI分配和目标函数的值
- 否则，报告问题无解
线性规划算法因为计算资源问题无法求得解的结果，所以转用启发式算法，线性规划后续用于启发式算法正确性的验证。
矩阵形式
注意到
y
i
,
j
y_{i,j}
yi,j 等价于
y
j
,
i
y_{j,i}
yj,i，
z
i
,
j
z_{i,j}
zi,j 等价于
z
j
,
i
z_{j,i}
zj,i，我们可以将目标函数先写成以下形式：
Minimize
F
o
b
j
=
∑
i
=
0
N
−
1
∑
j
=
0
N
−
1
(
a
i
,
j
⋅
y
i
,
j
+
b
i
,
j
⋅
y
i
,
j
+
c
i
,
j
⋅
z
i
,
j
)
\text{Minimize} \quad F_{obj} = \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} \left( a_{i,j} \cdot y_{i,j} + b_{i,j} \cdot y_{i,j} + c_{i,j} \cdot z_{i,j} \right)
MinimizeFobj=i=0∑N−1j=0∑N−1(ai,j⋅yi,j+bi,j⋅yi,j+ci,j⋅zi,j)
更进一步的，将其写成矩阵的形式。
- 计算相同PCI的情况（
Y
Y
Y）：
Y = X × X T Y = X \times X^T Y=X×XT
这里 Y [ i , j ] = y i , j Y[i, j]=y_{i,j} Y[i,j]=yi,j 表示小区 i i i 和小区 j j j 是否分配了相同的PCI，1是0否。即 Y Y Y 的每个元素由 X X X 的行向量的点乘得到，反映了两个小区是否被分配了相同的PCI。 - 计算模3相同的情况（
Z
Z
Z）：
首先定义一个矩阵 M M M，其中 M [ k , l ] = 1 M[k, l] = 1 M[k,l]=1 当且仅当 k m o d 3 = l m o d 3 k \mod 3 = l \mod 3 kmod3=lmod3。然后：
Z = ( X ⊙ M ) × ( X ⊙ M ) T Z = (X \odot M) \times (X \odot M)^T Z=(X⊙M)×(X⊙M)T
这里 ⊙ \odot ⊙ 表示哈达玛积（元素乘），即将 X X X 中的每个元素与 M M M 中对应的模3结果进行乘法操作，然后进行矩阵乘法。 - 目标函数的矩阵形式：
Minimize F o b j = sum ( ( A + B ) ⊙ Y ) + sum ( C ⊙ Z ) \text{Minimize} \quad F_{obj} = \text{sum}((A + B) \odot Y) + \text{sum}(C \odot Z) MinimizeFobj=sum((A+B)⊙Y)+sum(C⊙Z)
矩阵形式的目标函数在计算上可以利用 numpy 的广播形式进行迅速计算，采用合适的矩阵形式进行计算能够将时间缩短到原来的 1/70 。
数据观察后的思考
注意到所有频数为 156510 的小区与其他所有小区之间的所有 MR 数都为 0，故可以固定这个小区的 PCI 值（比如直接按顺序随便分配，意图在于防止只是数据漏报，如果一味的分配到同一个 PCI 可能会导致实际冲突 MR 剧增），并将它们排除 pcis 的分配行列，以减少运算消耗和加快启发式算法的工作速率。
核心代码
建议使用 Jupyter notebook 复现。
数据处理部分
from pathlib import Path
import pandas as pd
def read_excel_data(file_path, header_row=1):
return pd.read_excel(file_path, header=header_row)
# 路径名字
folder = Path('../附件')
file_names = ['附件1：小区基本信息.xlsx', '附件2：冲突及干扰矩阵数据.xlsx', '附件3：混淆矩阵数据.xlsx']
basic_info_data = read_excel_data(folder / file_names[0])
conflict_interference_data = read_excel_data(folder / file_names[1])
confusion_matrix_data = read_excel_data(folder / file_names[2])
# 筛选出频点156510的小区
specific_freq_cells = basic_info_data[basic_info_data['频点'] == 156510].copy()
specific_freq_cells['PCI'] = 1007 # 为这些小区分配PCI值1007
# 筛选优化小区相关数据，排除频点为156610的小区
optimized_cells = basic_info_data[(basic_info_data['是否优化小区'] == 1) & (basic_info_data['频点'] != 156510)]
optimized_conflict_interference = conflict_interference_data[
(conflict_interference_data['小区编号'].isin(optimized_cells['小区编号'])) &
(conflict_interference_data['邻小区编号'].isin(optimized_cells['小区编号']))
]
optimized_confusion = confusion_matrix_data[
(confusion_matrix_data['小区0编号'].isin(optimized_cells['小区编号'])) &
(confusion_matrix_data['小区1编号'].isin(optimized_cells['小区编号']))
]
# 检查数据大小
print("优化小区数据大小:", optimized_cells.shape)
print("优化冲突及干扰数据大小:", optimized_conflict_interference.shape)
print("优化混淆矩阵数据大小:", optimized_confusion.shape)
# 获取小区编号
cell_ids = optimized_cells['小区编号'].tolist()
# 打印出来看一下cell_ids的前几个元素
print("Cell IDs (first 10):", cell_ids[:10])
构造矩阵
# 先转成矩阵形式
import numpy as np
import pandas as pd
N = len(cell_ids)
A = np.zeros((N, N), dtype=int)
B = np.zeros((N, N), dtype=int)
C = np.zeros((N, N), dtype=int)
# A
for index, row in optimized_conflict_interference.iterrows():
i = cell_ids.index(row['小区编号'])
j = cell_ids.index(row['邻小区编号'])
A[i][j] += row['冲突MR数']
# B
for index, row in optimized_confusion.iterrows():
i = cell_ids.index(row['小区0编号'])
j = cell_ids.index(row['小区1编号'])
B[i][j] += row['混淆MR数']
# B[j][i] = B[i][j]
# C
for index, row in optimized_conflict_interference.iterrows():
i = cell_ids.index(row['小区编号'])
j = cell_ids.index(row['邻小区编号'])
C[i][j] += row['干扰MR数']
# 提前处理减少运算
A = A+A.T
B = B+B.T
C = C+C.T
AB = A+B
目标函数
初始学习版
def objective_function(pcis, A, B, C):
total_conflicts = 0
total_interferences = 0
total_confusions = 0
num_cells = len(pcis)
# 遍历所有小区组合，考虑题目描述，因此需要计算a_{ij} + a_{ji}等
for i in range(num_cells):
for j in range(i + 1, num_cells):
if pcis[i] == pcis[j]:
total_conflicts += (A[i][j] + A[j][i])
total_confusions += (B[i][j] + B[j][i])
if pcis[i] % 3 == pcis[j] % 3:
total_interferences += (C[i][j] + C[j][i])
return total_conflicts + total_interferences + total_confusions
矩阵计算版本
def objective_function(pcis, AB, C):
# 将pcis数组转换为NumPy数组
pci_array = np.array(pcis)
# 计算PCI相同的掩码
pci_conflict_mask = pci_array[:, None] == pci_array
# 计算模3相同的掩码
mod3_conflict_mask = pci_array[:, None] % 3 == pci_array % 3
# 计算总冲突、混淆和干扰MR数
total_conflicts_confusions = np.sum((AB) * pci_conflict_mask)
total_interferences = np.sum(C * mod3_conflict_mask)
return int((total_conflicts_confusions + total_interferences)/2)
查看当前PCI配置下的得分
# optimized_cells[optimized_cells['小区编号'] == id] 返回id对应的行
pcis = [optimized_cells[optimized_cells['小区编号'] == id]['现网PCI'].values[0] for id in cell_ids]
# 不用管顺序，
current_score = objective_function(pcis, AB,C)
# 打印当前得分
print("当前PCI配置下的目标函数得分:", current_score)
相关变量理解
pcis 就是按小区顺序对应的 PCI
display(optimized_cells[optimized_cells['小区编号'] == cell_ids[0]])
display(optimized_cells['现网PCI'][:10])
display(pcis[:10])
模拟退火
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
# plt.rcParams['font.sans-serif'] = ['PingFang HK']
# plt.rcParams['axes.unicode_minus'] = False
def simulated_annealing(AB, C, num_cells=1613, temp=1e3, initial_temp=1e3, cooling_rate=0.99, max_steps=10000, current_solution=None): # temp=1000, initial_temp=1000, cooling_rate=0.999, max_steps=1000
# 随机初始化PCI分配
current_cost = objective_function(current_solution, AB, C)
best_solution = current_solution.copy()
best_cost = current_cost
costs = [] # 用于存储每一步的成本
temperatures = [] # 用于存储每一步的温度
acceptance_rate = [] # 用于记录接受率
accepted_changes = 0 # 用于计算接受率
last_improvement_step = 0
for step in tqdm(range(max_steps), desc="Running Simulated Annealing Algorithm"):
base_change_rate = 0.01 # 初始变化率
min_change_rate = 0.001 # 最小变化率
change_rate = max(min_change_rate, base_change_rate * (temp / initial_temp))
num_changes = int(np.ceil(num_cells * change_rate))
change_indices = np.random.choice(num_cells, num_changes, replace=False)
# 产生新的候选解
candidate_solution = current_solution.copy()
candidate_solution[change_indices] = np.random.randint(0, 1008, size=num_changes)
# 计算新解的成本
candidate_cost = objective_function(candidate_solution, AB, C)
# 接受概率，以便于跳出局部最小值
if candidate_cost < current_cost or np.random.rand() < np.exp((current_cost - candidate_cost) / temp):
current_solution, current_cost = candidate_solution, candidate_cost
accepted_changes += 1
# 更新最优解和最低成本
if candidate_cost < best_cost:
best_solution = candidate_solution.copy()
best_cost = candidate_cost
last_improvement_step = step
costs.append(current_cost)
acceptance_rate.append(accepted_changes / (step + 1))
temperatures.append(temp)
# 自适应调整冷却率（简易版，可以自行修改）
# max_temp = initial_temp
if step - last_improvement_step > 100:
temp *= 0.999 # 减缓冷却
else:
temp *= cooling_rate # 正常冷却
return best_solution, best_cost, costs, temperatures, acceptance_rate
def plot_simulation_results(costs, temperatures, acceptance_rate):
fig, axs = plt.subplots(3, 1, figsize=(10, 12))
axs[0].plot(costs, label='Cost')
axs[0].set_title('Cost Over Time')
axs[0].set_xlabel('Iteration')
axs[0].set_ylabel('Cost')
axs[0].legend()
axs[1].plot(temperatures, label='Temperature', color='red')
axs[1].set_title('Temperature Over Time')
axs[1].set_xlabel('Iteration')
axs[1].set_ylabel('Temperature')
axs[1].legend()
axs[2].plot(acceptance_rate, label='Acceptance Rate', color='green')
axs[2].set_title('Acceptance Rate Over Time')
axs[2].set_xlabel('Iteration')
axs[2].set_ylabel('Acceptance Rate')
axs[2].legend()
plt.tight_layout()
plt.show()
多线程并行求解
from multiprocessing import Pool
def parallel_simulated_annealing(args):
AB, C, temp, initial_temp, cooling_rate, max_steps, initial_solution = args
return simulated_annealing(AB, C, temp=temp, initial_temp=initial_temp, cooling_rate=cooling_rate, max_steps=max_steps,current_solution=initial_solution)
# 修改参数前问下GPT这些会有什么影响
temp = initial_temp = 1e5
cooling_rate = 0.99
max_steps = 100000
n_cores = 6 # 注意cpu是几核
n = 6 # 这里的意思是总共跑6个，最好设置为n_cores的倍数，进度条跑完一次等于n_cores个都跑完
args = [(AB, C, temp, initial_temp, cooling_rate, max_steps, np.random.randint(0, 1008, size=N)) for _ in range(n)]
with Pool(processes=6) as pool:
results = pool.map(parallel_simulated_annealing, args)
拓展：线性规划
线性规划是无法求解整个问题的，一开始我是写来当作练手目的（第一天时间充足），顺便用于校验小型范围下模拟退火的正确性，放在这里仅供学习，建议看懂后再决定是否使用。
# 用于校验小型范围结果
N=5
A = AA[:N, :N]
B = BB[:N, :N]
C = CC[:N, :N]
K=3
import pulp as pl
import numpy as np
# 创建问题实例
model = pl.LpProblem("PCI_Planning", pl.LpMinimize)
# 创建决策变量
x = pl.LpVariable.dicts("x", (range(N), range(K)), cat=pl.LpBinary)
y = pl.LpVariable.dicts("y", (range(N), range(N)), cat=pl.LpBinary)
z = pl.LpVariable.dicts("z", (range(N), range(N)), cat=pl.LpBinary)
u = pl.LpVariable.dicts("u", (range(N), range(N), range(K)), 0, 1, pl.LpBinary) # 辅助变量 for y
v = pl.LpVariable.dicts("v", (range(N), range(N), range(K), range(K)), 0, 1, pl.LpBinary) # 辅助变量 for z
# 目标函数
model += pl.lpSum([
(A[i][j] + A[j][i]) * y[i][j] +
(B[i][j] + B[j][i]) * y[i][j] +
(C[i][j] + C[j][i]) * z[i][j]
for i in range(N) for j in range(i+1, N)
])
# 约束条件
for i in range(N):
model += pl.lpSum(x[i][k] for k in range(K)) == 1
for i in range(N):
for j in range(i+1, N):
model += y[i][j] == pl.lpSum(u[i][j][k] for k in range(K))
model += z[i][j] == pl.lpSum(v[i][j][k][l] for k in range(K) for l in range(K) if k % 3 == l % 3)
for k in range(K):
model += u[i][j][k] <= x[i][k]
model += u[i][j][k] <= x[j][k]
model += u[i][j][k] >= x[i][k] + x[j][k] - 1
for l in range(K):
if k % 3 == l % 3:
model += v[i][j][k][l] <= x[i][k]
model += v[i][j][k][l] <= x[j][l]
model += v[i][j][k][l] >= x[i][k] + x[j][l] - 1
# 求解模型
# model.solve()
model.solve(pl.PULP_CBC_CMD(msg=True, options=['-strategy', '1', '-maxIt', '10000']))
# 输出结果
if pl.LpStatus[model.status] == 'Optimal':
print(f"找到最优解！{pl.value(model.objective)}")
pci_assignments = {i: [k for k in range(K) if pl.value(x[i][k]) > 0.5][0] for i in range(N)}
print("PCI 分配结果：", pci_assignments)
else:
print("没有找到最优解。")
print("Status:", pl.LpStatus[model.status])
# 输出结果
for v in model.variables():
if v.varValue > 0:
print(f"{v.name} = {v.varValue}")
题目2
考虑冲突、混淆和干扰的不同优先级，给这2067个小区重新分配PCI。
首先保证冲突的 MR 数降到最低，在此基础上保证混淆的 MR 数降到最低，最后尽量降低模3干扰的MR 数。
建模
这里要先看一下视频理解概念，只看优先级法的部分就行，很简单
【【数学建模】15分钟速通多目标规划模型：理论+算法+例题+代码！】 【精准空降到 06:52】
这是一个分层优化问题，需要逐步解决，每个步骤都建立在前一个阶段的基础上。
首先，定义三个优化目标：
- 最小化冲突MR数
- 在冲突MR数最小化的基础上，最小化混淆MR数
- 在冲突和混淆MR数都最小化的基础上，最小化模3干扰MR数
定义目标函数和约束
将每个优化目标定义为一个单独的线性规划问题，在每个阶段，前一个阶段的解将用作当前阶段的约束。(这里依旧维持原来函数定义，后面再转矩阵)
第一阶段：最小化冲突MR数
实际上减少冲突 MR 数主要取决于小区间的相对位置及其 PCI 分配，模型假设不需要考虑小区之间的相对位置，所以通过重新分配 PCI，确保附表数据相关联的小区不共享相同的 PCI 值，是直接减少冲突 MR 数的最有效方法。
Minimize
F
conflict
=
∑
i
=
0
N
−
1
∑
j
=
i
+
1
N
−
1
(
a
i
,
j
+
a
j
,
i
)
y
i
,
j
\text{Minimize} \quad F_{\text{conflict}} = \sum_{i=0}^{N-1} \sum_{j=i+1}^{N-1} (a_{i,j}+a_{j,i})\ y_{i,j} \\
MinimizeFconflict=i=0∑N−1j=i+1∑N−1(ai,j+aj,i) yi,j
s.t.
∑
k
=
0
1007
x
i
,
k
=
1
,
∀
i
∈
N
y
i
,
j
=
∑
k
=
0
1007
x
i
,
k
x
j
,
k
,
∀
i
,
j
∈
N
x
i
,
k
∈
{
0
,
1
}
,
∀
i
∈
N
,
∀
k
∈
{
0
,
…
,
1007
}
y
i
,
j
∈
{
0
,
1
}
,
∀
i
,
j
∈
N
\begin{aligned} \text{s.t.} \quad \sum_{k=0}^{1007} x_{i,k} &= 1, && \forall i \in N \\ y_{i,j} &= \sum_{k=0}^{1007} x_{i,k} x_{j,k}, && \forall i, j \in N \\ x_{i,k} &\in \{0, 1\}, && \forall i \in N, \forall k \in \{0, \ldots, 1007\} \\ y_{i,j} &\in \{0, 1\}, && \forall i, j \in N \\ \end{aligned}
s.t.k=0∑1007xi,kyi,jxi,kyi,j=1,=k=0∑1007xi,kxj,k,∈{0,1},∈{0,1},∀i∈N∀i,j∈N∀i∈N,∀k∈{0,…,1007}∀i,j∈N
记录最终的 y y y 为 y i , j ∗ y_{i,j}^* yi,j∗。
第二阶段：最小化混淆MR数
- 维持第一阶段的解，确保之前设定的
y
i
,
j
y_{i,j}
yi,j 不被改变:
y i , j ≥ y i , j ∗ ∀ i , j y_{i,j} \geq y_{i,j}^* \quad \forall i, j yi,j≥yi,j∗∀i,j - 每个小区只分配一个PCI值（同1）：
∑ k = 0 1007 x i , k = 1 , ∀ i \sum_{k=0}^{1007} x_{i,k} = 1, \quad \forall i ∑k=01007xi,k=1,∀i - （定义）小区i和j是否分配了相同的PCI（同1）：
y i , j = ∑ k = 0 1007 x i , k ⋅ x j , k , ∀ i , j y_{i,j} = \sum_{k=0}^{1007} x_{i,k} \cdot x_{j,k}, \quad \forall i, j yi,j=∑k=01007xi,k⋅xj,k,∀i,j
记录最终的 y i , j y_{i,j} yi,j 为 y i , j ∗ ∗ y_{i,j}^{**} yi,j∗∗
在第一阶段的解决方案基础上，增加优化目标：
Minimize
F
confusion
=
∑
i
=
0
N
−
1
∑
j
=
i
+
1
N
−
1
b
i
j
y
i
j
\text{Minimize} \quad F_{\text{confusion}} = \sum_{i=0}^{N-1} \sum_{j=i+1}^{N-1} b_{ij} y_{ij} \\
MinimizeFconfusion=i=0∑N−1j=i+1∑N−1bijyij
s.t.
y
i
j
≥
y
i
j
∗
,
∀
i
,
j
∈
N
其余约束
同第一阶段
\begin{aligned} \text{s.t.} \quad y_{ij} &\geq y_{ij}^*, && \forall i, j \in N \\ \text{其余约束} & \text{同第一阶段} \\ \end{aligned}
s.t.yij其余约束≥yij∗,同第一阶段∀i,j∈N
需要保持第一阶段的结果作为约束，确保不会增加冲突MR数。
第三阶段：最小化模3干扰MR数
在前两个阶段的解决方案基础上，增加优化目标
y
i
,
j
≥
y
i
,
j
∗
∗
∀
i
,
j
y_{i,j} \geq y_{i,j}^{**} \quad \forall i, j
yi,j≥yi,j∗∗∀i,j，维持第一二阶段的解，确保之前设定的
y
i
,
j
y_{i,j}
yi,j 不被改变。
：
Minimize
F
mod3
=
∑
i
=
0
N
−
1
∑
j
=
i
+
1
N
−
1
(
c
i
,
j
+
c
j
,
i
)
z
i
,
j
\text{Minimize} \quad F_{\text{mod3}} = \sum_{i=0}^{N-1} \sum_{j=i+1}^{N-1} (c_{i,j}+c_{j,i})\ z_{i,j} \\
MinimizeFmod3=i=0∑N−1j=i+1∑N−1(ci,j+cj,i) zi,j
s.t.
y
i
,
j
≥
y
i
,
j
∗
∗
,
∀
i
,
j
∈
N
θ
i
=
k
i
m
o
d
3
,
∀
i
∈
N
z
i
,
j
=
{
1
,
if
θ
i
=
θ
j
,
0
,
otherwise
,
∀
i
,
j
∈
N
其余约束
同第一阶段
\begin{aligned} \text{s.t.} \quad y_{i,j} &\geq y_{i,j}^{**}, && \forall i, j \in N \\ \theta_i &= k_i \mod 3, && \forall i \in N \\ z_{i,j} &= \begin{cases} 1, & \text{if } \theta_i = \theta_j, \\ 0, & \text{otherwise}, \end{cases} && \forall i, j \in N \\ \text{其余约束} & \text{同第一阶段} \\ \end{aligned}
s.t.yi,jθizi,j其余约束≥yi,j∗∗,=kimod3,={1,0,if θi=θj,otherwise,同第一阶段∀i,j∈N∀i∈N∀i,j∈N
问题求解
第一阶段 使用模拟退火优化冲突MR数
在第一阶段，使用模拟退火算法来最小化由冲突矩阵 A A A 定义的冲突MR数。关键步骤包括：
- 转换能量函数：
为了提高计算效率，我们将能量函数转换为矩阵形式：
E conflict = sum ( A ⊙ Y ) E_{\text{conflict}} = \text{sum}(A \odot Y) Econflict=sum(A⊙Y)
其中， A [ i , j ] A[i, j] A[i,j] 表示小区 i i i 和 j j j 分配相同的 PCI 值时引发的冲突 MR 数。 - 初始化：随机生成初始 PCI 配置，确保每个小区获得有效的 PCI 值。
- 迭代过程：
- 选择：随机选择一个小区并尝试更换其PCI值，探索可能的新配置。
- 评估：计算新配置的能量并与当前配置比较。
- 接受或拒绝：如果新配置的能量低于当前配置，根据 Metropolis 准则判断是否接受更高能量的配置。
- 终止条件：当达到预定的迭代次数或系统温度下降到一定阈值。
第二阶段：使用模拟退火优化冲突和混淆MR数
第二阶段在第一阶段的基础上继续进行，目标是最小化由 A + B A+B A+B 矩阵定义的冲突和混淆 MR 数。此阶段考虑了小区间的冲突（由矩阵 A A A 表示）和混淆（由矩阵 B B B 表示），以找到减少这两种干扰的最优 PCI 配置。
- 转换能量函数：
E confusion = sum ( ( A + B ) ⊙ Y ) E_{\text{confusion}} = \text{sum}((A + B) \odot Y) Econfusion=sum((A+B)⊙Y)
其中， A + B A + B A+B 表示冲突和混淆矩阵的总和，每个元素 A [ i , j ] + B [ i , j ] A[i, j] + B[i, j] A[i,j]+B[i,j] 表示小区 i i i 和 j j j 分配相同的 PCI 值时引发的总冲突和混淆 MR 数。 Y Y Y 是一个对称矩阵，其元素 Y [ i , j ] Y[i, j] Y[i,j] 为1表示小区 i i i 和 j j j 分配了相同的 PCI，否则为 0。 - 扰动策略：
- 我们使用矩阵
A
A
A 来约束扰动解的生成过程，使得选择候选 PCI 值时尽可能避免增加冲突MR数。
伪代码：
Algorithm : Perturb Solution for Conflict and Confusion Minimization Input: pcis, A, B Output: candidate_solution function PERTURB_SOLUTION(pcis, A) N ← length of pcis pci_array ← array from pcis candidate_solution ← copy of pci_array idx ← random choice from 1 to N conflict_pcis ← pci_array[where A[idx] == 1] possible_pcis ← array from 0 to 1007 non_conflicting_pcis ← set difference of possible_pcis and conflict_pcis if size of non_conflicting_pcis > 0 then new_pci ← random choice from non_conflicting_pcis candidate_solution[idx] ← new_pci return candidate_solution end function - 我们使用矩阵
A
A
A 来约束扰动解的生成过程，使得选择候选 PCI 值时尽可能避免增加冲突MR数。
- 迭代过程：
- 完全等价于第一阶段，但增加约束使得解在扰动的过程中不增加冲突MR数。
- 降温策略和终止条件：与第一阶段相同，逐步降低温度，细化解空间搜索。
第三阶段：最小化模3干扰MR数（问题抽象化）
灵感来源
在前两个阶段的优化中，已成功将小区之间的冲突和混淆 MR 数降至零，而这个过程中，我发现：具体的 PCI 值并非影响 MR 数的关键因素。这一发现激发了进一步简化问题的思路：将关注点从具体的 PCI 值转移到它们模 3 的结果上，这种抽象不仅简化了问题，还显著缩小了解空间，因为每个小区的 PCI 值现在只需考虑三种可能的模 3 结果（0，1或2）。
故现在将问题进一步抽象化为只关注模 3 结果（0, 1, 2），而不是具体的 PCI 值。
目标函数和约束的重新定义
- 符号定义
- θ i \theta_i θi：表示小区 i i i 的PCI值模3后的结果，取值为 0, 1, 或 2。
- 转换能量函数：现在目标函数仅关注模 3 的结果。
Minimize F m o d 3 = ∑ i = 0 N − 1 ∑ j = 0 N − 1 ( C [ i , j ] ⋅ 1 ( θ i = θ j ) ) \text{Minimize} \quad F_{mod3} = \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} \left(C[i,j] \cdot \mathbf{1}(\theta_i = \theta_j)\right) MinimizeFmod3=i=0∑N−1j=0∑N−1(C[i,j]⋅1(θi=θj))
其中， 1 ( θ i = θ j ) \mathbf{1}(\theta_i = \theta_j) 1(θi=θj) 是一个指示函数，如果 θ i \theta_i θi 和 θ j \theta_j θj 相等则为1，否则为0。 - 定义
θ i = k i m o d 3 \theta_i = \text{k}_i \mod 3 θi=kimod3
θ i \theta_i θi 的取值为 0, 1, 或 2，这代表了 PCI 值模 3 后的三种可能结果。 - 为了将能量函数转为矩阵形式，我们根据
θ
\theta
θ 再次定义矩阵
Z
Z
Z 用于表示小区之间模 3 结果的相同性
Z i , j = { 1 如果 θ i = θ j 0 如果 θ i ≠ θ j Z_{i, j} = \begin{cases} 1 & \text{如果 } \theta_i = \theta_j \\ 0 & \text{如果 } \theta_i \neq \theta_j \end{cases} Zi,j={10如果 θi=θj如果 θi=θj
矩阵形式：
E mod3 = sum ( C ⊙ Z ) E_{\text{mod3}} = \text{sum}(C \odot Z) Emod3=sum(C⊙Z)
求解策略
在模型中，每个小区的模 3 结果被视作遗传算法中的一个基因，而所有小区的模3结果共同构成一个染色体，代表一个完整的模 3 分配方案。这种类比生物 RNA 的模型设定具有以下特征：
- 基因：每个小区的模3结果，相当于遗传算法中解决问题的基本单元。
- 染色体：整个小区群的模3结果组合，表示一个完整的模3分配方案。
遗传算法过程描述
- 选择
- 锦标赛选择法。
- 交叉
- 多点交叉。
- 变异
- 随机变异：在这一过程中，通过以一定的概率（即变异率）随机改变个体中的一个或多个基因（在此模型中为某小区的模3结果）。
在完整的求解过程中，我们发现仅需要找到使得干扰 MR 数最小的那个序列，就可以依据当前序列和 A + B A+B A+B 矩阵来扰动解使得整体的目标函数收敛于最小干扰 MR 总数，在完整的建模周期中，所有找寻到的序列都可以通过该策略将其余 MR 数降为 0。
也就是说，对于该类型的题，我们可以通过先解决第三阶段，然后再解决一，二阶段。第四题便是使用的反向思路来解决。
最终的三个结果（暂定）：
| 指标 | (冲突，混淆，干扰) |
|---|---|
| MR数 | (0，0，27241020) |
核心代码
目标函数（旧版）
这里是旧版，结果没什么差别，只是新版用的是objective_function(pcis, AB, C)，也就是第一题那个（应该是后面重构代码的时候更新的），过了太久了懒得重新梳理最新的代码，所以放个第二版的矩阵目标函数在这里，差别在于推理速度（应该是慢 1/2），可以自行替换。
代码仅供参考学习。
def objective_function(pcis, A, C, B):
# 将pcis数组转换为NumPy数组
pci_array = np.array(pcis)
# 计算PCI相同的掩码
pci_conflict_mask = pci_array[:, None] == pci_array
# 计算模3相同的掩码
mod3_conflict_mask = pci_array[:, None] % 3 == pci_array % 3
# 计算总冲突、混淆和干扰MR数
total_conflicts_confusions = np.sum((A + B) * pci_conflict_mask)
total_interferences = np.sum(C * mod3_conflict_mask)
return int((total_conflicts_confusions + total_interferences)/2)
def conflicts_objective_function(pcis, A):
pci_array = np.array(pcis)
pci_conflict_mask = pci_array[:, None] == pci_array
total_conflicts = np.sum(A * pci_conflict_mask)
return int(total_conflicts / 2)
def confusion_objective_function(pcis, B):
pci_array = np.array(pcis)
pci_conflict_mask = pci_array[:, None] == pci_array
total_confusions = np.sum(B * pci_conflict_mask)
return int(total_confusions / 2)
def interference_objective_function(pcis, C):
pci_array = np.array(pcis)
mod3_conflict_mask = pci_array[:, None] % 3 == pci_array % 3
total_interferences = np.sum(C * mod3_conflict_mask)
return int(total_interferences / 2)
# 得到分开的分数
def get_score(pcis, A, C, B):
# 将pcis数组转换为NumPy数组
pci_array = np.array(pcis)
# 计算PCI相同的掩码
pci_conflict_mask = pci_array[:, None] == pci_array
# 计算模3相同的掩码
mod3_conflict_mask = pci_array[:, None] % 3 == pci_array % 3
# 计算总冲突、混淆和干扰MR数
total_conflicts = np.sum(A * pci_conflict_mask)
total_confusions = np.sum(B * pci_conflict_mask)
total_interferences = np.sum(C * mod3_conflict_mask)
return int(total_conflicts/2), int(total_confusions/2), int(total_interferences/2)
第一阶段
# 用M来指导可以生成具有更多可能性的解
def perturb_solution(pcis, M):
N = len(pcis)
pci_array = np.array(pcis)
candidate_solution = pci_array.copy()
# 随机选择一个小区进行PCI更改尝试
idx = np.random.choice(N)
# 找出当前不会引起冲突的所有可能的PCI
possible_pcis = np.arange(1008)
conflict_pcis = pci_array[M[idx] == 1] # 从M中找出与idx冲突的PCI
non_conflicting_pcis = np.setdiff1d(possible_pcis, conflict_pcis)
# 如果存在非冲突的PCI，随机选择一个进行更改
if non_conflicting_pcis.size > 0:
candidate_solution[idx] = np.random.choice(non_conflicting_pcis)
return candidate_solution
def simulated_annealing_conflicts(A, C, B, M, num_cells=1613, temp=1e3, initial_temp=1e3, cooling_rate=0.99, max_steps=10000, current_solution=None): # temp=1000, initial_temp=1000, cooling_rate=0.999, max_steps=1000
# 随机初始化PCI分配
if current_solution is None:
current_solution = np.random.randint(0, 1008, size=num_cells)
current_cost = conflicts_objective_function(current_solution, A)
best_solution = current_solution.copy()
best_cost = current_cost
costs = [] # 用于存储每一步的成本
temperatures = [] # 用于存储每一步的温度
acceptance_rate = [] # 用于记录接受率
accepted_changes = 0 # 用于计算接受率
# 如果后续代码有加什么新东西，可以都列表记录下来
last_improvement_step = 0
with tqdm(total=max_steps, desc="Simulated Annealing for Confusion Optimization", unit="step") as pbar:
for step in range(max_steps):
if best_cost == 0:
break
# 产生新的候选解
# candidate_solution = perturb_solution(current_solution.copy(), M)
base_change_rate = 0.01 # 初始变化率
min_change_rate = 0.001 # 最小变化率
change_rate = max(min_change_rate, base_change_rate * (temp / initial_temp))
num_changes = int(np.ceil(num_cells * change_rate))
change_indices = np.random.choice(num_cells, num_changes, replace=False)
# 产生新的候选解
candidate_solution = current_solution.copy()
candidate_solution[change_indices] = np.random.randint(0, 1008, size=num_changes)
# 计算新解的成本
candidate_cost = conflicts_objective_function(candidate_solution, A,)
# 接受概率，以便于跳出局部最小值
if candidate_cost < current_cost or np.random.rand() < np.exp((current_cost - candidate_cost) / temp):
current_solution, current_cost = candidate_solution, candidate_cost
accepted_changes += 1
# 更新最优解和最低成本
if candidate_cost < best_cost:
best_solution = candidate_solution.copy()
best_cost = candidate_cost
last_improvement_step = step
costs.append(current_cost)
acceptance_rate.append(accepted_changes / (step + 1))
temperatures.append(temp)
# 自适应调整冷却率（简易版，可以自行修改）
# max_temp = initial_temp
if step - last_improvement_step > 100:
temp *= 0.99999 # 减缓冷却
else:
temp *= cooling_rate # 正常冷却
pbar.set_description(f"(Cost: {current_cost}, Temp: {temp:.2f})")
pbar.update(1)
return best_solution, best_cost, costs, temperatures, acceptance_rate
def plot_simulation_results(costs, temperatures, acceptance_rate):
fig, axs = plt.subplots(3, 1, figsize=(10, 12))
axs[0].plot(costs, label='Cost')
axs[0].set_title('Cost Over Time')
axs[0].set_xlabel('Iteration')
axs[0].set_ylabel('Cost')
axs[0].legend()
axs[1].plot(temperatures, label='Temperature', color='red')
axs[1].set_title('Temperature Over Time')
axs[1].set_xlabel('Iteration')
axs[1].set_ylabel('Temperature')
axs[1].legend()
axs[2].plot(acceptance_rate, label='Acceptance Rate', color='green')
axs[2].set_title('Acceptance Rate Over Time')
axs[2].set_xlabel('Iteration')
axs[2].set_ylabel('Acceptance Rate')
axs[2].legend()
plt.tight_layout()
plt.show()
temp = initial_temp = 1e3 # 不用写那么多0，e后面是几就等于几个0
cooling_rate = 0.99
max_steps = 10000
results = simulated_annealing_conflicts(A, C, B, M, temp=temp, initial_temp=initial_temp, cooling_rate=cooling_rate, max_steps=max_steps,current_solution=None)
可视化
# 查找具有最低成本的运行
best_solution, best_cost, best_costs, best_temperatures, best_acceptance_rate = results
# 打印当前初始解中的最佳成本和对应的解
print("Best solution:", best_solution)
print("Minimum cost:", best_cost)
# 画图
plot_simulation_results(best_costs, best_temperatures, best_acceptance_rate)
第二阶段
def simulated_annealing_confusion(A, C, B, M, num_cells=1613, temp=1e3, initial_temp=1e3, cooling_rate=0.99, max_steps=10000, current_solution=None):
# 随机初始化PCI分配
if current_solution is None:
current_solution = np.random.randint(0, 1008, size=num_cells)
current_cost = confusion_objective_function(current_solution, B)
best_solution = current_solution.copy()
best_cost = current_cost
costs = [] # 用于存储每一步的成本
temperatures = [] # 用于存储每一步的温度
acceptance_rate = [] # 用于记录接受率
accepted_changes = 0 # 用于计算接受率
last_improvement_step = 0
with tqdm(total=max_steps, desc="Simulated Annealing for Confusion Optimization", unit="step") as pbar:
for step in range(max_steps):
# 产生新的候选解
candidate_solution = perturb_solution(current_solution.copy(), M)
# 计算新解的成本
candidate_cost = confusion_objective_function(candidate_solution, B)
# 接受概率，以便于跳出局部最小值
if candidate_cost < current_cost or np.random.rand() < np.exp((current_cost - candidate_cost) / temp):
current_solution, current_cost = candidate_solution, candidate_cost
accepted_changes += 1
# 更新最优解和最低成本
if candidate_cost < best_cost:
best_solution = candidate_solution.copy()
best_cost = candidate_cost
last_improvement_step = step
costs.append(current_cost)
acceptance_rate.append(accepted_changes / (step + 1))
temperatures.append(temp)
# 自适应调整冷却率（简易版，可以自行修改）
# max_temp = initial_temp
if step - last_improvement_step > 100:
temp *= 0.99999 # 减缓冷却
else:
temp *= cooling_rate # 正常冷却
pbar.set_description(f"(Cost: {current_cost}, Temp: {temp:.2f})")
pbar.update(1)
if best_cost == 0:
break
return best_solution, best_cost, costs, temperatures, acceptance_rate
results2 = simulated_annealing_confusion(A, C, B, M, temp=temp, initial_temp=initial_temp, cooling_rate=cooling_rate, max_steps=max_steps,current_solution=best_solution)
可视化
best_solution2, best_cost2, best_costs2, best_temperatures2, best_acceptance_rate2 = results2
# 打印当前初始解中的最佳成本和对应的解
print("Best solution:", best_solution2)
print("Minimum cost:", best_cost2)
# 画图
plot_simulation_results(best_costs2, best_temperatures2, best_acceptance_rate2)
第三阶段
当前是常规解决方案，依然是模拟退火。
抽象化后的新思路是抛弃前两个阶段的模拟退火算法，利用遗传算法替代： interference_objective_function() 作为目标函数，令其最小化，然后将这个模 3 的序列丢给第三阶段的模拟退火算法，最后将会得到（0, 0, optimal score)。
def create_mod3_interference_matrix(num_pcis=1008):
Mod3M = np.zeros((num_pcis, num_pcis), dtype=int)
for i in range(num_pcis):
for j in range(num_pcis):
if i % 3 == j % 3:
Mod3M[i, j] = 1
return Mod3M
def perturb_solution_for_interference(pcis, Mod3M, M):
N = len(pcis)
pci_array = np.array(pcis)
candidate_solution = pci_array.copy()
idx = np.random.choice(N)
possible_pcis = np.arange(1008)
conflict_pcis = pci_array[M[idx] == 1] # 从M中找出与idx冲突的PCI
non_conflicting_pcis = np.setdiff1d(possible_pcis, conflict_pcis)
# 过滤出不引起模3冲突的PCI
mod3_conflict_pcis = pci_array[np.where(Mod3M[pci_array[idx] % 3] == 1)[0]]
non_conflicting_mod3_pcis = np.setdiff1d(non_conflicting_pcis, mod3_conflict_pcis)
if non_conflicting_mod3_pcis.size > 0:
candidate_solution[idx] = np.random.choice(non_conflicting_mod3_pcis)
return candidate_solution
def simulated_annealing_interference(A, C, B, Mod3M, M, num_cells=1613, temp=1e3, initial_temp=1e3, cooling_rate=0.99, max_steps=10000, current_solution=None): # temp=1000, initial_temp=1000, cooling_rate=0.999, max_steps=1000
# 随机初始化PCI分配
if current_solution is None:
current_solution = np.random.randint(0, 1008, size=num_cells)
current_cost = interference_objective_function(current_solution, C)
best_solution = current_solution.copy()
best_cost = current_cost
costs = [] # 用于存储每一步的成本
temperatures = [] # 用于存储每一步的温度
acceptance_rate = [] # 用于记录接受率
accepted_changes = 0 # 用于计算接受率
last_improvement_step = 0
with tqdm(total=max_steps, desc="Simulated Annealing for Confusion Optimization", unit="step") as pbar:
for step in range(max_steps):
candidate_solution = perturb_solution_for_interference(current_solution.copy(), Mod3M, M)
# 计算新解的成本
candidate_cost = interference_objective_function(candidate_solution, C)
# 接受概率，以便于跳出局部最小值
if candidate_cost < current_cost or np.random.rand() < np.exp((current_cost - candidate_cost) / temp):
current_solution, current_cost = candidate_solution, candidate_cost
accepted_changes += 1
# 更新最优解和最低成本
if candidate_cost < best_cost:
best_solution = candidate_solution.copy()
best_cost = candidate_cost
last_improvement_step = step
costs.append(current_cost)
acceptance_rate.append(accepted_changes / (step + 1))
temperatures.append(temp)
# 自适应调整冷却率（简易版，可以自行修改）
# max_temp = initial_temp
if step - last_improvement_step > 100:
temp *= 0.99999 # 减缓冷却
else:
temp *= cooling_rate # 正常冷却
pbar.set_description(f"(Cost: {current_cost}, Temp: {temp:.2f})")
pbar.update(1)
if best_cost == 0:
break
return best_solution, best_cost, costs, temperatures, acceptance_rate
def parallel_simulated_annealing_interference(args):
A, C, B, Mod3M, M, temp, initial_temp, cooling_rate, max_steps, initial_solution = args
return simulated_annealing_interference(A, C, B, Mod3M, M, temp=temp, initial_temp=initial_temp, cooling_rate=cooling_rate, max_steps=max_steps,current_solution=initial_solution)
def create_conflict_mask(A):
M = (A > 0).astype(int)
return M
Mod3M = create_mod3_interference_matrix()
M = create_conflict_mask(A+B)
results3 = simulated_annealing_interference(A, C, B, Mod3M, M, temp=temp, initial_temp=initial_temp, cooling_rate=cooling_rate, max_steps=100000,current_solution=best_solution2)
可视化
best_solution3, best_cost3, best_costs3, best_temperatures3, best_acceptance_rate3 = results3
# 打印当前初始解中的最佳成本和对应的解
print("Best solution:", best_solution3)
print("Minimum cost:", best_cost3)
# 画图
plot_simulation_results(best_costs3, best_temperatures3, best_acceptance_rate3)
A_score, B_score, C_score = get_score(best_solution3, A, C, B)
print(f"冲突MR数：{A_score}\n混淆分数：{B_score}\n干扰MR数：{C_score}")
写在前面，题目 3，4 为 1，2 题的简单拓展，如果仅出于学习目的，代码写到这里就够了，后面只需理解思路。
题目3
给这 2067 个小区重新分配 PCI，使得所有可能被影响到的小区间的冲突 MR 数、混淆 MR 数和模 3 干扰 MR 数的总和最少。
建模
要求在考虑更多小区的情况下重新分配 PCI，使得涉及到的所有小区之间的冲突 MR 数、混淆 MR 数和模 3 干扰 MR 数的总和最小。
Minimize
F
obj
=
∑
(
i
,
j
)
∈
S
(
A
i
,
j
y
i
,
j
+
B
i
,
j
y
i
,
j
+
C
i
,
j
z
i
,
j
)
\text{Minimize} \quad F_{\text{obj}} = \sum_{(i, j) \in S} \left( A_{i,j} y_{i,j} + B_{i,j} y_{i,j} + C_{i,j} z_{i,j} \right)
MinimizeFobj=(i,j)∈S∑(Ai,jyi,j+Bi,jyi,j+Ci,jzi,j)
s.t.
∑
k
=
1
M
x
i
,
k
=
1
,
∀
i
∈
N
,
y
i
,
j
=
∑
k
=
1
M
x
i
k
x
j
,
k
,
∀
(
i
,
j
)
∈
S
,
θ
i
=
k
i
m
o
d
3
,
∀
i
∈
N
,
z
i
j
=
{
1
if
θ
i
=
θ
j
,
0
if
θ
i
≠
θ
j
,
∀
(
i
,
j
)
∈
S
,
x
i
,
k
∈
{
0
,
1
}
,
∀
i
∈
N
,
∀
k
∈
{
0
,
…
,
1007
}
,
y
i
,
j
∈
{
0
,
1
}
,
∀
(
i
,
j
)
∈
S
,
z
i
,
j
∈
{
0
,
1
}
,
∀
(
i
,
j
)
∈
S
.
\begin{align*} \text{s.t.} \quad & \sum_{k=1}^{M} x_{i,k} = 1, && \forall i \in N, \\ & y_{i,j} = \sum_{k=1}^{M} x_{ik} x_{j,k}, && \forall (i, j) \in S, \\ & \theta_i = k_i \mod 3, && \forall i \in N, \\ & z_{ij} = \begin{cases} 1 & \text{if } \theta_i = \theta_j, \\ 0 & \text{if } \theta_i \neq \theta_j, \end{cases} && \forall (i, j) \in S, \\ & x_{i,k} \in \{0, 1\}, && \forall i \in N, \forall k \in \{0, \ldots, 1007\}, \\ & y_{i,j} \in \{0, 1\}, && \forall (i, j) \in S, \\ & z_{i,j} \in \{0, 1\}, && \forall (i, j) \in S. \end{align*}
s.t.k=1∑Mxi,k=1,yi,j=k=1∑Mxikxj,k,θi=kimod3,zij={10if θi=θj,if θi=θj,xi,k∈{0,1},yi,j∈{0,1},zi,j∈{0,1},∀i∈N,∀(i,j)∈S,∀i∈N,∀(i,j)∈S,∀i∈N,∀k∈{0,…,1007},∀(i,j)∈S,∀(i,j)∈S.
其中 S S S 表示所有可能被影响的小区对的集合， A A A, B B B, 和 C C C 分别代表冲突、混淆和模 3 干扰的 MR 数矩阵。
问题求解
实际问题求解过程中，我们依旧可以先针对模 3 干扰进行优化，因为这个是干扰数大头。
pcis 是我们实际的优化对象，对于目标函数来说， A B C ABC ABC 都是固定的，只有 Y Y Y, Z Z Z 是变的，而 Y Y Y, Z Z Z 由 pcis 直接创建，我们想找到一个最优的解就得找到其对应分配的 pcis，对于 3，4 问来讲，我们仍然只需要对优化小区进行 PCI 分配，其他相关小区的 PCI 是固定死的。
数据预处理
数据观察：64345, 63409 小区也和所有小区无冲突。
为了给 2067 个小区重新分配 PCI，使得所有可能被影响的小区间总 MR 数量最少。 我们首先筛选优化小区的相关数据，并且排除了所有频数为 156510 的小区。
在数据预处理方面与问题 1 和问题 2 有所不同。具体表现在以下几个方面:
- 获取与优化小区(1613，因为排除掉了频数为 156510 的小区)所有相关的的小区编号，同时去掉两个不参与分配 PCI 的小区(编号为 64345 与 63409)，因为 这两个小区的 MR 数为 0。因此最后的矩阵维度为 2890×2890。
- 在 PCI 分配方面，我们的数据矩阵首先按问题 1 当中的优化小区排序，然后再将其余小区扩展到矩阵当中，因此，PCI 分配的长度也需要为 2890。由于我们是按顺序对生成矩阵，并且只需要对优化的小区进行分配 PCI，于是我们在进行寻找最优解的过程当中只截取优化小区的长度进行 PCI 的分配，而并不是对整个 2890 长度的列表进行计算，这样极大减少了计算资源的消耗以及计算时间。
经过上述预处理后，与问题 1 一致，使用模拟退火算法寻找最优解。
题目4
考虑冲突、混淆和干扰的不同优先级，给这 2067 个小区重新 分配 PCI，也是考虑所有可能被影响到的小区间的冲突、混淆和模 3 干扰。首先保证冲突的MR 数降到最低，在此基础上保证混淆的MR 数降到最低，最后尽量降低模 3 干扰的 MR 数。
问题求解
这里我们没有采取问题二的逐层优化顺序，而是根据抽象后的问题，直接对随机初始解进行了模 3 处理，并先通过遗传算法找到序列优解， 再通过类似于问题二的约束调整，确保在不增加模 3 干扰的前提下，将冲突和混淆 MR 数降降到 0。
| 指标 | (冲突，混淆，干扰) |
|---|---|
| MR数 | (0，0，31205749) |
4389
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
