PMCM-2019-F题-问题1（3）
- 作业要求
- 进一步剪枝
- 建立数学模型
- 绘制优化结果图
- 建模及优化结果
- 以后写博文会在同名公众号“ManTou馒头”上同步更新啦~ 主要原因是：
（1）博文容易被搬运
（2）博客回复大家的问题不及时(当然大家博客有问题还是可以回复和私信一起讨论！ManTou会尽量一起讨论解决的！)
（3）CSDN广告多，影响观感，公众号没有粉丝就不用担心有广告植
入了哦~- 代码要结合前几次作业的内容，需要整合一下~
- 每周的事情比较多，数学模型作业的质量可能不能让大家满意~ 大家多理解支持吖~
- 公众号回复“PMCM2019F”，获取源程序 拜托不要直接copy，拜托了拜托了
作业要求
- 进一步剪枝
在 作业（3）的基础上，进一步减少网络邻接矩阵中的有向边的数量。 - 建立模型
以 航迹长度尽可能小 和 经过矫正区域进行校正的次数尽可能少 为目标建立模型。
进一步剪枝
- 第二次剪枝（ PMCM-2019-F题-问题1（2）中的剪枝称为第一次剪枝）
- 剪枝的标准
不走回头路：删除方向和前进方向相反的有向边。 - 采用的方法
向量乘积小于零则说明是钝角，也就说明是回头路，这样的边就剪掉 - 代码实现
# 第二次剪枝
AB = [x_ab[1] - x_ab[0], y_ab[1] - y_ab[0], z_ab[1] - z_ab[0]]
edge1 = copy.copy(edge) # 建立一个新的edge
for (i, j) in edge1:
if A_B[0] * (x_coord[j] - x_coord[i]) \
+ A_B[1] * (y_coord[j] - y_coord[i]) \
+ A_B[1] * (z_coord[j] - z_coord[i]) <= 0:
edge.remove((i, j))
print("第二次剪枝之后的边数：", len(edge))
- 第三次剪枝
- 剪枝的标准
转弯飞行的波动不太大：删去距离AB路线太远的校正点，以AB为轴，R为半径的圆柱之外的点 - 采用的方法
计算校正点到直线的距离（参考链接），大于R的剪去。 - 代码实现
# 第三次剪枝
# 剪去以AB为轴，R为半径的圆柱之外的点
# https://zhuanlan.zhihu.com/p/24760577
R = 10000
third_cut = []
for i in range(1, shape - 1):
if dist[i, shape - 1] * np.sqrt(1 - (((x_coord[shape - 1] - x_coord[i]) * (x_coord[shape - 1] - x_coord[0])
+ (y_coord[shape - 1] - y_coord[i]) * (y_coord[shape - 1] - y_coord[0])
+ (z_coord[shape - 1] - z_coord[i]) * (z_coord[shape - 1] - z_coord[0]))
/ (dist[0, shape - 1] * dist[i, shape - 1])) ** 2) >= R:
third_cut.append(i)
edge2 = copy.copy(edge) # 建立一个新的edge
for (i, j) in edge2:
if i in third_cut or j in third_cut:
edge.remove((i, j))
print("第三次剪枝之后的边数：", len(edge))
- 剪枝的结果
剪枝可以有效减少计算量，可能对最终结果没有太大的影响。
建立数学模型
- 模型的基本认识
- 模型应该是最短路0-1整数规划模型， x i j {x_{ij}} xij是一个0-1变量，A点的入度为0出度为1，B点的出度为0入度为1。其他点入度等于出度（有向图的出度和入度）,K是顶点项的惩罚系数（惩罚函数的作用）
min ∑ ( i , j ) ∈ E w i j x i j + K ∑ ( i , j ) ∈ E x i j \min \sum_{(i, j) \in E} w_{i j} x_{i j}+K \sum_{(i, j) \in E} x_{i j} min∑(i,j)∈Ewijxij+K∑(i,j)∈Exij
s.t. { ∑ ( i , j ) ∈ E x i j = ∑ ( j , k ) ∈ E x j k , ∀ j ≠ A , B ∑ ( A , j ) ∈ E x A j = 1 , ∑ ( j , A ) ∈ E x j A = 0 ∑ ( j , B ) ∈ E x j B = 1 , ∑ ( B , j ) ∈ E x B j = 0 x i j = 0 or 1 , ∀ i , j \left\{\begin{array}{l}\sum_{(i, j) \in E} x_{i j}=\sum_{(j, k) \in E} x_{j k}, \forall j \neq A, B \\ \sum_{(A, j) \in E} x_{A j}=1, \quad \sum_{(j, A) \in E} x_{j A}=0 \\ \sum_{(j, B) \in E} x_{j B}=1, \sum_{(B, j) \in E} x_{B j}=0 \\ x_{i j}=0 \text { or } 1, \quad \forall i, j\end{array}\right. ⎩⎪⎪⎨⎪⎪⎧∑(i,j)∈Exij=∑(j,k)∈Exjk,∀j=A,B∑(A,j)∈ExAj=1,∑(j,A)∈ExjA=0∑(j,B)∈ExjB=1,∑(B,j)∈ExBj=0xij=0 or 1,∀i,j - 水平和垂直校正点误差不超过上限
v i ≤ α 1 , h i ≤ α 2 , i ∈ V v_{i} \leq \alpha_{1}, h_{i} \leq \alpha_{2}, i \in V vi≤α1,hi≤α2,i∈V
v i ≤ β 1 , h i ≤ β 2 , i ∈ H v_{i} \leq \beta_{1}, h_{i} \leq \beta_{2}, i \in H vi≤β1,hi≤β2,i∈H
v B ≤ θ , h B ≤ θ v_{B} \leq \theta, h_{B} \leq \theta vB≤θ,hB≤θ - 相邻校正点误差之间的关系
( I i h i + δ w i j − h j ) x i j = 0 , ( i , j ) ∈ E \left(I_{i} h_{i}+\delta w_{i j}-h_{j}\right) \quad x_{i j}=0,(i, j) \in E (Iihi+δwij−hj)xij=0,(i,j)∈E
( ( 1 − I i ) v i + δ w i j − v j ) x i j = 0 , ( i , j ) ∈ E \left(\left(1-I_{i}\right) v_{i}+\delta w_{i j}-v_{j}\right) x_{i j}=0,(i, j) \in E ((1−Ii)vi+δwij−vj)xij=0,(i,j)∈E - 采用大M法将限制条件线性化（大M法）
I i h i + δ w i j − h j ≤ ( 1 − x i j ) M , ( i , j ) ∈ E I_{i} h_{i}+\delta w_{i j}-h_{j} \leq\left(1-x_{i j}\right) M,(i, j) \in E Iihi+δwij−hj≤(1−xij)M,(i,j)∈E
( 1 − I i ) v i + δ w i j − v j ≤ ( 1 − x i j ) M , ( i , j ) ∈ E \left(1-I_{i}\right) v_{i}+\delta w_{i j}-v_{j} \leq\left(1-x_{i j}\right) M,(i, j) \in E (1−Ii)vi+δwij−vj≤(1−xij)M,(i,j)∈E
- Gurobi：Pyhon API
具体函数的使用方法，参考我的这篇博文：简单使用 Gurobi：Python API - 基本步骤
- 添加变量 x
- 添加变量h和v
- 添加关键的约束
- 添加目标函数
- 代码实现
# 建立模型
model1 = Model()
# 添加新变量x[i,j], i=0 to shape-1 ,j = 0 to shape-1
x = model1.addVars(shape, shape, vtype=GRB.BINARY, name='x')
# 添加新变量h[i],v[i],i = 0 to shape-1
h = model1.addVars(shape, vtype=GRB.CONTINUOUS, name='h')
v = model1.addVars(shape, vtype=GRB.CONTINUOUS, name='v')
# 添加限制条件 x[i,j]==0 根据矩阵信息进行剪枝
for i in range(shape):
for j in range(shape):
if C[i, j] == 0:
model1.addConstr(x[i, j] == 0)
# 添加限制条件的起始点，中间点，终点的出度和入度条件
# 对于有向图来说，顶点的出边条数称为该顶点的出度，顶点的入边条数称为该顶点的入度
out_degree = [0] * shape
in_degree = [0] * shape
# out_degree[i]表示第i个点的出度
for (i, j) in edge:
out_degree[i] = out_degree[i] + x[i, j]
# in_degree[i]表示第i个点的入度
for (j, i) in edge:
in_degree[i] = in_degree[i] + x[j, i]
for i in range(shape):
if i == 0:
# 起点A的出度为1，入度为0
model1.addConstr(out_degree[i] == 1)
model1.addConstr(in_degree[i] == 0)
elif 0 < i < shape - 1:
# 中间点的出度等于入度
model1.addConstr(out_degree[i] == in_degree[i])
else:
# 终点B的入度为1，出度为0
model1.addConstr(out_degree[i] == 0)
model1.addConstr(in_degree[i] == 1)
# 添加限制条件h，v进行约束，M = 100000
M = 10000
for (i, j) in edge:
model1.addConstr(prop[i] * h[i] + delta * dist[i, j] - h[j] <= M * (1 - x[i, j]))
model1.addConstr((1 - prop[i]) * v[i] + delta * dist[i, j] - v[j] <= M * (1 - x[i, j]))
for i in range(shape):
if i == 0:
# 起点的垂直和水平误差为0
model1.addConstr(h[i] == 0)
model1.addConstr(v[i] == 0)
elif 0 < i < shape - 1:
# 中间点的误差约束条件
if prop[i] == 1:
# 垂直校正点前的误差约束条件
model1.addConstr(h[i] <= alpha2)
model1.addConstr(v[i] <= alpha1)
else:
# 水平校正点前的误差约束条件
model1.addConstr(h[i] <= beta2)
model1.addConstr(v[i] <= beta1)
else:
# 终点前的垂直和水平误差约束条件
model1.addConstr(h[i] <= theta)
model1.addConstr(v[i] <= theta)
# 添加目标函数
K = 10000
dicct_modify = { }
for i in range(shape):
for j in range(shape):
dicct_modify[i, j] = K
dicct_modify = tupledict(dicct_modify)
model1.setObjective(dict_dist.prod(x) + dicct_modify.prod(x), GRB.MINIMIZE)
model1.update()
model1.optimize()
for s in model1.getVars():
if s.x == 1:
print('%s' '%g' % (s.varName, s.x))
print('Obj:%g' % model1.objVal)
绘制优化结果图
# 将优化的路径点取出绘图
node = [0, 503, 294, 91, 607, 540, 250, 340, 277, 612]
# 定义空数组用来存放校正点
x_coord_node = []
y_coord_node = []
z_coord_node = []
for j in node:
x_coord_node.append(x_coord[j])
y_coord_node.append(y_coord[j])
z_coord_node.append(z_coord[j])
ax.plot(x_coord_node, y_coord_node, z_coord_node, linewidth=2, c='r')
# 保存绘制的图形
plt.savefig("FigData1.jpg", bbox_inches='tight')
plt.show()
建模及优化结果
结果如图所示~
