PMCM-2019-F题-问题1（3）

1. 以后写博文会在同名公众号“ManTou馒头”上同步更新啦~ 主要原因是：
   （1）博文容易被搬运
   （2）博客回复大家的问题不及时(当然大家博客有问题还是可以回复和私信一起讨论！ManTou会尽量一起讨论解决的！)
   （3）CSDN广告多，影响观感，公众号没有粉丝就不用担心有广告植
   入了哦~
2. 代码要结合前几次作业的内容，需要整合一下~
3. 每周的事情比较多，数学模型作业的质量可能不能让大家满意~ 大家多理解支持吖~
4. 公众号回复“PMCM2019F”，获取源程序 拜托不要直接copy，拜托了拜托了

作业要求

1. 进一步剪枝
   在 作业（3）的基础上，进一步减少网络邻接矩阵中的有向边的数量。
2. 建立模型
   以 航迹长度尽可能小 和 经过矫正区域进行校正的次数尽可能少 为目标建立模型。

进一步剪枝

1. 第二次剪枝（ 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))

2. 第三次剪枝

• 剪枝的标准
  转弯飞行的波动不太大：删去距离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))

3. 剪枝的结果
   剪枝可以有效减少计算量，可能对最终结果没有太大的影响。

建立数学模型

1. 模型的基本认识

• 模型应该是最短路0-1整数规划模型，$x_{ij}$是一个0-1变量，A点的入度为0出度为1，B点的出度为0入度为1。其他点入度等于出度（有向图的出度和入度）,K是顶点项的惩罚系数（惩罚函数的作用）
  $\min {\sum }_{(i,j)\in E}{w}_{ij}{x}_{ij}+K{\sum }_{(i,j)\in E}{x}_{ij}$
  s.t. $\{\begin{array}{l}{\sum }_{(i,j)\in E}{x}_{ij}={\sum }_{(j,k)\in E}{x}_{jk},\forall j\ne A,B\\ {\sum }_{(A,j)\in E}{x}_{Aj}=1,{\sum }_{(j,A)\in E}{x}_{jA}=0\\ {\sum }_{(j,B)\in E}{x}_{jB}=1,{\sum }_{(B,j)\in E}{x}_{Bj}=0\\ x_{ij}=0\text{ or }1,\forall i,j\end{array}$
• 水平和垂直校正点误差不超过上限
  $v_i\le {\alpha }_1,h_i\le {\alpha }_2,i\in V$
  $v_i\le {\beta }_1,h_i\le {\beta }_2,i\in H$
  $v_B\le \theta ,h_B\le \theta$
• 相邻校正点误差之间的关系
  $\left(I_i{h}_i+\delta w_{ij}-h_j\right)x_{ij}=0,(i,j)\in E$
  $\left(\left(1-I_i\right)v_i+\delta w_{ij}-v_j\right)x_{ij}=0,(i,j)\in E$
• 采用大M法将限制条件线性化（大M法）
  $I_i{h}_i+\delta w_{ij}-h_j\le \left(1-x_{ij}\right)M,(i,j)\in E$
  $\left(1-I_i\right)v_i+\delta w_{ij}-v_j\le \left(1-x_{ij}\right)M,(i,j)\in E$

2. Gurobi：Pyhon API
   具体函数的使用方法，参考我的这篇博文：简单使用 Gurobi：Python API
3. 基本步骤

• 添加变量 x
• 添加变量h和v
• 添加关键的约束
• 添加目标函数

4. 代码实现

# 建立模型
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()

建模及优化结果

结果如图所示~