可怕的！数学建模——2019F题详细过程解_数模亚太赛2019真题解析 python-CSDN博客

本文链接：https://blog.csdn.net/weixin_45288557/article/details/109564710

数学建模2019F题

1 数学模型
2 Python代码—利用Python进行数据筛选、模型建立及求解全过程
3 输出结果
★佐佑思维公众号提供更多学习学术服务★ 二维码如下：
如有问题欢迎公众号交流

1 数学模型

决策变量：

$x_{ij}$ 代表在或不在最短轨迹上；非负变量 $v_{i}$ 、 $h_{i}$ 代表误差。

$x_{ij}=\left\{ \begin{aligned} 1 && (i,j)\in P_{AB} \\ 0 && (i,j)\notin P_{AB} \end{aligned} \right. ; I_{i}=\left\{ \begin{aligned} 1 && i\in V \\ 0 && i\in H \end{aligned} \right. ; \left\{ \begin{aligned} v_i &\geq& 0 \\ h_i &\geq& 0 \end{aligned} \right.$

约束条件：

1） $v_{A}=0, h_{A}=0, v_{B} \leq \theta, h_{B} \leq \theta$ ：

$\begin{array}{ll} v_{i} \leq & \left(1-I_{i}\right) \times 20+I_{i} \times 25 \\ h_{i} \leq & \left(1-I_{i}\right) \times 25+I_{i} \times 15 \end{array}$

$I_i$ =1， $v_i \le 25$ ， $h_i \le 15$ ，垂直误差校正前需要满足的条件； $I_i$ =0， $v_i \le 20$ ， $h_i \le 25$ ，水平误差校正前需要满足的条件。

2）限制最短路径上前一点和后一点的关系式满足 $I_{i} h_{i}+\delta d_{i j}=h_{j}、(1-I_{i}) v_{i}+\delta d_{i j}=v_{j}$ ：

$\begin{array}{l} I_{i} h_{i}+\delta d_{i j}-h_{j} \leq\left(1-x_{i j}\right) M \\ \left(1-I_{i}\right) v_{i}+\delta d_{i j}-v_{j} \leq\left(1-x_{i j}\right) M \end{array}$

$dist=(d_{ij})$ ， $d_{ij}$ 对应 $i$ 点到 $j$ 点的飞行距离；M表示一个非常大的数，用来增加或释放约束， $x_{ij}=0$ ，表示释放上式， $x_{ij}=1$ ，表示约束上式，右边为0。

3）限制同类顶点之间的有向边距离小于等于10 km

$\begin{array}{l} d_{ij}\leq 10000，i\in V,j \in V \\ d_{ij}\le 10000，i\in H,j \in H \end{array}$

4）在最短路径上，限制A点为起点，那么只一出不进；限制B点为终点，那么只一进不出；其他点为一进一出。

$\left\{ \begin{aligned} \sum_{j=0}^{613}x_{ij} =1 ,i=0（A点），\sum_{i=0}^{613}x_{ij} =0 ,j=0（A点）\\ \sum_{i=0}^{613}x_{ij} =1 ,j=613（B点），\sum_{j=0}^{613}x_{ij} =0 ,i=613（B点）\\ \sum_{i=0}^{613}\sum_{j=1}^{612}x_{ij} =\sum_{j=1}^{612}\sum_{k=0}^{613}x_{jk} \end{aligned} \right.$

目标函数：

$\quad z_1=\sum_{i=0}^{613}\sum_{j=0}^{613}d_{ij}x_{ij}$

$\quad z_2=\sum_{i=0}^{613}\sum_{j=0}^{613}x_{ij}$

可以发现 $z_1$ 和 $z_2$ 相差一个 $d_{ij}$ ，为综合考虑两个目标，将两个目标放在几乎同一个数量级上，并对两个目标函数进行加权组合，权值各为0.5
$\quad 0.5z_1+0.5(max(d_{ij}))z_2$

2 Python代码—利用Python进行数据筛选、模型建立及求解全过程

1）之前的工作—生成dist矩阵

#======================生成dist矩阵======================
# 导入库
import pandas as pd
import numpy as np

# 数据读取
io = r'数据.xlsx' #io，Excel的存储路径
sheet_name=0 #sheet_name，要读取的工作表名称：可以是整型数字、列表名或SheetN，也可以是上述三种组成的列表
usecols=range(0,5) #需要读取哪些列
skiprows=1 #跳过前n行；skiprows = [a, b, c]，跳过第a+1,b+1,c+1行（索引从0开始）

data = pd.read_excel(io=io, sheet_name=sheet_name, usecols=usecols, skiprows=skiprows)
print(data.head())#展示部分结果

#校正列
I=data.iloc[:,4].values

# 编号列
N=list(data.iloc[:,0].values)

# 确定校正点坐标
V=[];#垂直校正点
H=[];#水平校正点

for i in N:
    if data.iloc[i,4]==1:
        V.append(N[i]);#垂直校正点        
    if data.iloc[i,4]==0:       
        H.append(N[i]);#水平校正点
        
# data删掉编号列
data=data.iloc[:,1:5]
print('data:\n',data)
print('----------------------------------------------------')
print('\ndata的数据类型',type(data))

# 定义dist初始矩阵
dist=np.zeros((613,613))

# 定义欧式距离函数
def Ed_Dis(m, n):
    return np.linalg.norm(m - n)

# 计算距离并存入dist
for h in range(613):
    m=data.iloc[h,:3].values    
    for l in range(613):
        n=data.iloc[l,:3].values
        
        dist[h,l]=Ed_Dis(m, n)

2）进行边权值的筛选

#======================将dist中不同点，即不为0的权值以编号的形式存储======================
#这里注意，只对编号列表进行删除，保留dist，可以避免例如元素数据为0等造成不必要的麻烦（毕竟普通的矩阵计算对于计算机来说非常简单）
Boundary=[]
for i in range(613):
    for j in range(613):
        if dist[i,j]!=0:
            Boundary.append([i,j])
print('\n初始权值数：',len(Boundary))

#======================（1）约束：删掉以AB为轴，半径10000外的点======================
#定义A→B向量x
A=data.iloc[0,:3].values
B=data.iloc[-1,:3].values
x_AB=B-A

cut_dian1=[]
for i in range(1,612):
    y_i=data.iloc[i,:3].values
    iB=B-y_i
    if dist[i, 612]!=0:
        
        if dist[i, 612] * np.sqrt(1 - ((np.dot(iB,x_AB)) / (dist[0, 612] * dist[i, 612])) ** 2) >= 10000:
            cut_dian1.append(i)
#————————————————————————————————————————————————————————————————————————————————————
#实现删除
import copy
Boundary1 = copy.copy(Boundary)
#print(len(Boundary1))

for [i, j] in Boundary1:    
    if i in cut_dian1 or j in cut_dian1:
        Boundary.remove([i, j])
print("删掉以AB为轴，半径10000外的点之后：", len(Boundary))

#======================（2）约束：、利用校正条件筛选======================

cut_dian2=[]

for i in range(1,612):
    # 垂直误差不大于25个单位，水平误差不大于15个单位时才能进行垂直误差校正
    for j in V:
        if (0.001*dist[i,j])<=15:
            continue # 跳出本次循环，进入下一次循环
        else:
            cut_dian2.append([i,j]) # 不满足条件不可进行垂直误差调整
    
    # 垂直误差不大于20个单位，水平误差不大于25个单位时才能进行水平误差校正
    for j in H:
        if (0.001*dist[i,j])<=20:
            continue # 跳出本次循环，进入下一次循环
        else:
            cut_dian2.append([i,j]) # 不满足条件不可进行水平误差调整
    
    # 到达终点时垂直误差和水平误差均应小于30个单位
    for j in [612]:
        if (0.001*dist[i,j])<30:
            continue # 跳出本次循环，进入下一次循环
        else:
            cut_dian2.append([i,j]) # 不满足条件不能按照规划路径到达B点
#————————————————————————————————————————————————————————————————————————————————————
#实现删除
import copy
Boundary2 = copy.copy(Boundary) 
#print(len(Boundary2))

for [i, j] in Boundary2:
    if [i,j] in cut_dian2:
        
        Boundary.remove([i, j])
print("校正条件筛选之后：", len(Boundary))

#======================（3）约束：删去距离大于 10 km 的同类顶点之间的有向边======================

cut_dian3=[]

# 删去距离大于 10 km 的同为垂直校正顶点之间的有向边；
for i in V:    
    for j in V:        
        if dist[i,j]>10000:
            cut_dian3.append([i,j])
        else:            
            continue # 跳出本次循环，进入下一次循环

# 删去距离大于 10 km 的同为水平校正顶点之间的有向边；
for i in H:
    for j in H:        
        if dist[i,j]>10000:
            cut_dian3.append([i,j])
        else:            
            continue # 跳出本次循环，进入下一次循环

#————————————————————————————————————————————————————————————————————————————————————
#实现删除
import copy
Boundary3 = copy.copy(Boundary)  # 复制一个
#print(len(Boundary3))

for [i, j] in Boundary3:
    if [i,j] in cut_dian3:
        
        Boundary.remove([i, j])
print("删去距离大于 10 km 的同类顶点之间的有向边之后：", len(Boundary))

#======（4）约束：删除方向与前进方向相反的有向边（标准：在这里定义与A→B向量与任意向量夹角大于90°认为与前进方向相反）======

# 定义计算夹角函数 cos值在[0,1]之间，代表角度在0～90°；cos值在[-1,0]之间，代表角度在90°～180°
# 所以这里只需要计算cos值即可判断角度的范围
def Angle_cos(x, y):
    
    # 分别计算两个向量的模：
    l_x=np.sqrt(x.dot(x))
    l_y=np.sqrt(y.dot(y))
    #print('向量的模=',l_x,l_y)
    
    # 计算两个向量的点积
    dian=x.dot(y)
    #print('向量的点积=',dian)
    
    # 计算夹角的cos值：
    cos_=dian/(l_x*l_y)
    #print('夹角的cos值=',cos_)
    
    return cos_
#————————————————————————————————————————————————————————————————————————————————————
#定义A→B向量x
A=data.iloc[0,:3].values
B=data.iloc[-1,:3].values
x_AB=B-A

cut_dian4=[]
for i in range(613):
    y_i=data.iloc[i,:3].values
    
    for j in range(613):
        y_j=data.iloc[j,:3].values
        
        if i!=j:            
            # i→j向量y
            y_ij=y_j-y_i
            
            Ang_cos=Angle_cos(x_AB, y_ij)  # cos值在[0,1]之间，代表角度在0～90°；cos值在[-1,0]之间，代表角度在90°～180°
            if Ang_cos<=0:
                cut_dian4.append([i,j])
            else:
                continue # 跳出本次循环，进入下一次循环
#————————————————————————————————————————————————————————————————————————————————————
#实现删除
import copy
Boundary4 = copy.copy(Boundary)  # 复制一个
#print(len(Boundary4))

for [i, j] in Boundary4:
    if [i,j] in cut_dian4:
        
        Boundary.remove([i, j])
print("删除方向与前进方向相反的有向边之后：", len(Boundary))

3）模型计算前需要进行的一些准备工作

Boundary_ab = copy.copy(Boundary)  # 复制一下，以防后面对Boundary还有什么改变方便操作

I[0]=0 # 由于I中起点和终点信息为字符无法进行运算，故在这里将起点的字符改变为0，从后面的约束条件可以发现对结果没有影响，而终点不会用到，故不用更改
#————————————————————————————————————————————————————————————————————————————————————
#根据Gurobi数据类型的要求对dist进行格式更换
# 将dist转换为tupledict类型
from gurobipy import *
dist_tup = {}
for i in range(613):
    for j in range(613):
        dist_tup[i,j] = dist[i,j]
print(type(dist_tup))
dist_tup = tupledict(dist_tup)
print(type(dist_tup))
#————————————————————————————————————————————————————————————————————————————————————
#寻找dij的最大值
max_dij=np.max(dist)
print('max_dij=',max_dij)

max_d = {}
for i in range(613):
    for j in range(613):
        max_d[i,j] = max_dij
print(type(max_d))
max_d = tupledict(max_d)
print(type(max_d))

4）建立模型并求解

from gurobipy import *

# Create a new model
m = Model("feiji")

# Create variables vtype(可选):新变量的数据类型（GRB.CONTINUOUS,GRB.BINARY,GRB.INTEGER,GRB.SEMICONT,GRB.SEMIINT）
x = m.addVars(613, 613, vtype=gurobipy.GRB.BINARY, name="x")
v = m.addVars(613, vtype=gurobipy.GRB.CONTINUOUS, name="v")
h = m.addVars(613, vtype=gurobipy.GRB.CONTINUOUS, name="h")

# 更新变量环境
m.update()

# Set objective
m.setObjective(0.5*dist_tup.prod(x) + 0.5*max_d.prod(x), GRB.MINIMIZE)

# Add constraint1:
for i in range(613):
    if i == 0:
        # 起点的垂直和水平误差为0
        m.addConstr(h[i] == 0)
        m.addConstr(v[i] == 0)
    elif 0 < i < 612:
        # 中间点的误差约束条件
        m.addConstr(v[i]<=20*(1-I[i])+25*I[i])
        m.addConstr(h[i]<=25*(1-I[i])+15*I[i])
    else:
        # 终点前的垂直和水平误差约束条件
        m.addConstr(h[i] <= 30)
        m.addConstr(v[i] <= 30)

## Add constraint2:
M = 10000
for (i, j) in Boundary_ab:

    m.addConstr(I[i] * h[i] + 0.001 * dist[i, j] - h[j] <= M * (1 - x[i, j]))
    m.addConstr((1 - I[i]) * v[i] + 0.001 * dist[i, j] - v[j] <= M * (1 - x[i, j]))
            
#Add constraint3: 在最短路径上，限制A点为起点，那么只出不进；限制B点为终点，那么只进不出；其他点为一进一出
out_degree = [0] * 613
in_degree = [0] * 613
# out_degree[i]表示第i个点的出度

for (i, j) in Boundary_ab:
    out_degree[i] = out_degree[i] + x[i, j]
# in_degree[i]表示第i个点的入度
for (j, i) in Boundary_ab:
    in_degree[i] = in_degree[i] + x[j, i]

for i in range(613):
    if i == 0:
        # 起点A的出度为1，入度为0
        m.addConstr(out_degree[i] == 1)
        m.addConstr(in_degree[i] == 0)
    elif 0 < i < 612:
        # 中间点的出度等于入度
        m.addConstr(out_degree[i] == in_degree[i])
    else:
        # 终点B的入度为1，出度为0
        m.addConstr(out_degree[i] == 0)
        m.addConstr(in_degree[i] == 1)


# 执行模型
m.optimize()

#computeIIS()用于检查哪些约束是互相矛盾的，即去掉这些矛盾约束剩下的约束构成的问题是可行的
#m.computeIIS()
#m.write("model3.ilp")

#输出数据
for v in m.getVars():
    if v.x == 1:
        print('%s %g' % (v.varName, v.x))

print('Obj: %g' % m.objVal) # 查看单目标规划模型的目标函数值

3 输出结果

#========================print(data.head())========================
编号    X坐标（单位: m）     Y坐标（单位:m）   Z坐标（单位: m） 校正点类型
0   0      0.000000  50000.000000  5000.000000   A 点
1   1  33070.825831   2789.480761  5163.524680     0
2   2  54832.887019  49179.219108  1448.303791     1
3   3  77991.545911  63982.175286  5945.823038     0
4   4  16937.179535  84714.341903  5360.293002     0

#=======================print('data:\n',data)=======================
data:
         X坐标（单位: m）     Y坐标（单位:m）   Z坐标（单位: m） 校正点类型
0         0.000000  50000.000000  5000.000000   A 点
1     33070.825831   2789.480761  5163.524680     0
2     54832.887019  49179.219108  1448.303791     1
3     77991.545911  63982.175286  5945.823038     0
4     16937.179535  84714.341903  5360.293002     0
..             ...           ...          ...   ...
608   45789.201490  21191.207906   440.001872     1
609   94917.734859  82958.725321  6169.657707     0
610   14870.601682  95939.173362  8248.836185     0
611   93009.567446   4549.333260  7882.608297     1
612  100000.000000  59652.343380  5022.001164    B点

[613 rows x 4 columns]
----------------------------------------------------

#=======================print('\ndata的数据类型',type(data))=======================
data的数据类型 <class 'pandas.core.frame.DataFrame'>

#=======================#将dist中不同点，即不为0的权值以边的形式存储为Boundary列表后===========

初始权值数： 375156

删掉以AB为轴，半径10000外的点之后： 17292

校正条件筛选之后： 5097

删去距离大于 10 km 的同类顶点之间的有向边之后： 3737

删除方向与前进方向相反的有向边之后： 1809

#===========================将dist转换为tupledict类型、以及最大值max_dij的tupledict类型===================
<class 'dict'>
<class 'gurobipy.tupledict'>
max_dij= 139286.85817998415
<class 'dict'>
<class 'gurobipy.tupledict'>
#===========================模型求解结果===========================
Gurobi Optimizer version 9.0.3 build v9.0.3rc0 (win64)
Optimize a model with 5459 rows, 376995 columns and 13889 nonzeros
Model fingerprint: 0xf473aac1
Variable types: 1226 continuous, 375769 integer (375769 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+04]
  Objective range  [7e+04, 1e+05]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+04]
Presolve removed 3077 rows and 375344 columns
Presolve time: 0.28s
Presolved: 2382 rows, 1651 columns, 10685 nonzeros
Variable types: 121 continuous, 1530 integer (1530 binary)
Found heuristic solution: objective 2520005.5056

Root relaxation: objective 5.386343e+05, 195 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 538634.323    0   14 2520005.51 538634.323  78.6%     -    0s
H    0     0                    1246717.7584 538634.323  56.8%     -    0s
H    0     0                    963471.90541 538634.323  44.1%     -    0s
H    0     0                    891937.45079 538634.323  39.6%     -    0s
H    0     0                    820316.25261 538944.757  34.3%     -    0s
     0     0 539292.457    0   16 820316.253 539292.457  34.3%     -    0s
     0     0 539344.334    0   18 820316.253 539344.334  34.3%     -    0s
H    0     0                    749581.67384 539344.334  28.0%     -    0s
H    0     0                    679412.93773 539344.334  20.6%     -    0s
     0     0 539776.400    0   14 679412.938 539776.400  20.6%     -    0s
H    0     0                    679367.18791 539776.400  20.5%     -    0s
     0     0 539809.111    0   14 679367.188 539809.111  20.5%     -    0s
     0     0 540309.127    0   16 679367.188 540309.127  20.5%     -    0s
     0     0 540361.414    0   20 679367.188 540361.414  20.5%     -    0s
     0     0 540768.325    0   14 679367.188 540768.325  20.4%     -    0s
     0     0 540974.035    0   14 679367.188 540974.035  20.4%     -    0s
     0     0 608968.758    0    8 679367.188 608968.758  10.4%     -    0s
     0     0 609166.422    0    8 679367.188 609166.422  10.3%     -    0s
     0     0 610045.789    0   14 679367.188 610045.789  10.2%     -    0s
     0     0 610045.789    0   10 679367.188 610045.789  10.2%     -    0s
     0     0 610045.789    0   17 679367.188 610045.789  10.2%     -    0s
     0     0 613720.884    0   20 679367.188 613720.884  9.66%     -    0s
     0     0 633510.455    0   18 679367.188 633510.455  6.75%     -    0s
     0     0 678556.660    0    6 679367.188 678556.660  0.12%     -    0s
*    0     0               0    679221.39430 679221.394  0.00%     -    0s

Explored 1 nodes (792 simplex iterations) in 0.69 seconds
Thread count was 12 (of 12 available processors)

Solution count 9: 679221 679367 679413 ... 2.52001e+06

Optimal solution found (tolerance 1.00e-04)
Best objective 6.792213943009e+05, best bound 6.792213943009e+05, gap 0.0000%

#===========================模型输出数据========================
x[0,503] 1
x[91,607] 1
x[250,340] 1
x[277,612] 1
x[294,91] 1
x[340,277] 1
x[503,294] 1
x[540,250] 1
x[607,540] 1
Obj: 679221