＃＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝生成dist矩阵＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＃ 导入库import pandas as pdimport 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：＇，data）print（＇－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－＇）print（＇data的数据类型＇，type（data））
＃ 定义dist初始矩阵dist＝np．zeros（（613，613））
＃ 定义欧式距离函数def Ed＿Dis（m， n）： return np．linalg．norm（m － n）
＃ 计算距离并存入distfor 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）＃＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝将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（＇初始权值数：＇，len（Boundary））＃＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝（1）删掉以AB为轴，半径10000外的点－这个条件有参考其他博客＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＃定义A→B向量xA＝data．iloc［0，：3］．valuesB＝data．iloc［－1，：3］．valuesx＿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 copyBoundary1 ＝ 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 copyBoundary2 ＝ 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 copyBoundary3 ＝ 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向量xA＝data．iloc［0，：3］．valuesB＝data．iloc［－1，：3］．valuesx＿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 copyBoundary4 ＝ copy．copy（Boundary） ＃ 复制一个＃print（len（Boundary4））
for ［i， j］ in Boundary4： if ［i，j］ in cut＿dian4： Boundary．remove（［i， j］）print（＂删除方向与前进方向相反的有向边之后：＂， len（Boundary））
模型计算前需要进行的一些准备工作：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＿dijprint（type（max＿d））max＿d ＝ tupledict（max＿d）print（type（max＿d））
建立模型并求解：from gurobipy import ＊
＃ Create a new modelm ＝ 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 objectivem．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 ＝ 10000for （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］ ＊ 613in＿degree ＝ ［0］ ＊ 613
for （i， j） in Boundary＿ab： out＿degree［i］ ＝ out＿degree［i］ ＋ x［i， j］
for （j， i） in Boundary＿ab： in＿degree［i］ ＝ in＿degree［i］ ＋ x［j， i］
for i in range（613）： if i ＝＝ 0： ＃ 起点A 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 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） ＃ 查看单目标规划模型的目标函数值附录1：输出结果自然是：node ＝ ［0； 503； 294； 91； 607； 540； 250； 340； 277； 612］
