from gurobipy import Model, GRB, quicksum
import math
# ==========================
# ① 原始数据
# ==========================
V = list(range(1, 11)) # 10 个任务点
N = [0] + V # 含基站 0
U = list(range(1, 11)) # 10 架无人机
R = list(range(1, 11)) # 最多 10 次行程
coords = {0: (0.0, 0.0), 1: (2.4, 7.0), 2: (-4.0, 8.0), 3: (7.6, 4.4),
4: (-8.0, 3.0), 5: (-6.0, -5.6), 6: (5.0, -6.4), 7: (9.6, -2.0),
8: (-2.4, -9.0), 9: (1.0, -7.6), 10: (-9.0, -2.0)}
d_ij = {(i, j): round(math.dist(coords[i], coords[j]) * 1000, 3)
for i in N for j in N if i != j}
# 速度 v_ij
v_ij = {(i,j): 12.0 for i in N for j in N if i!=j}
# 飞行时间 t_ij = d_ij / v_ij
t_ij = {(i,j): d_ij[(i,j)]/v_ij[(i,j)] for i in N for j in N if i!=j}
# 基础能耗 E_ij = 距离（线性能耗）
E_ij = {(i,j): t_ij[(i,j)]/45 for i in N for j in N if i!=j}
# E_max: 单块电池最大能量
E_max = 600
#服务时间
s_i = {i: (0.0 if i==0 else 2.0) for i in N}
# 备用大M
BIGM = 30000
# ==========================
# ② 建立模型
# ==========================
m = Model("weighted_VRP")
# 变量
x = m.addVars(N, N, U, R, vtype=GRB.BINARY, name="x")
y_ur = m.addVars(U, R, vtype=GRB.BINARY, name="y_ur")
y_u = m.addVars(U, vtype=GRB.BINARY, name="y_u")
# t_iur: 到达/开始服务时间（连续）
t_iur = m.addVars(N, U, R, vtype=GRB.CONTINUOUS, lb=0.0, name="t")
# e_iur: 剩余电量（连续）
e_iur = m.addVars(N, U, R, vtype=GRB.CONTINUOUS, lb=0.0, name="e")
# ==========================
# ③ 约束（沿用之前裸骨架 + 联动约束）
# ==========================
# 每个任务点恰好被访问一次
for i in V:
m.addConstr(quicksum(x[j, i, u, r] for j in N if j != i for u in U for r in R) == 1)
m.addConstr(quicksum(x[i, j, u, r] for j in N if j != i for u in U for r in R) == 1)
# 流守恒
for u in U:
for r in R:
for j in V:
m.addConstr(quicksum(x[i, j, u, r] for i in N if i != j) ==
quicksum(x[j, k, u, r] for k in N if k != j))
# 基站出/回与 y_ur 绑定
for u in U:
for r in R:
m.addConstr(quicksum(x[0, j, u, r] for j in V) == y_ur[u, r])
m.addConstr(quicksum(x[i, 0, u, r] for i in V) == y_ur[u, r])
# 弧与行程激活联动
for i in N:
for j in N:
if i == j: continue
for u in U:
for r in R:
m.addConstr(x[i, j, u, r] <= y_ur[u, r])
# 无人机启用与行程关系
R_max = len(R)
for u in U:
m.addConstr(quicksum(y_ur[u, r] for r in R) <= R_max * y_u[u])
# 对称破缺（可选）
for idx in range(len(U)-1):
m.addConstr(y_u[U[idx]] >= y_u[U[idx+1]])
for u in U:
for idx in range(len(R)-1):
m.addConstr(y_ur[u, R[idx]] >= y_ur[u, R[idx+1]])
# 时间顺序约束（含绕飞额外时间）：
# t_jur >= t_iur + s_i + t_ij + sum_m t_ijm * w_ijurm - M*(1 - x_ijur)
for i in N:
for j in N:
if i==j: continue
for u in U:
for r in R:
m.addConstr(
t_iur[j,u,r] >= t_iur[i,u,r] + s_i[i] + t_ij[(i,j)] +
- BIGM * (1 - x[i,j,u,r]),
name=f"time_order_{i}_{j}_u{u}_r{r}"
)
# 返回基站后的再出发时间连贯：
# t_ju(r+1) >= t_0ur + t_0j + sum_m t_0jm * w_0ju(r+1)m - M*(1 - x_0ju(r+1))
for u in U:
for idx in range(len(R)-1):
r = R[idx]
rnext = R[idx+1]
for j in N:
if j == 0: # 再出发的目标不能还是基站
continue
m.addConstr(
t_iur[j,u,rnext] >= t_iur[0,u,r] + t_ij[(0,j)] +
- BIGM * (1 - x[0,j,u,rnext]),
name=f"relaunch_u{u}_r{r}_to_r{rnext}_j{j}"
)
# t_0ur >= t_iur + s_i + t_i0 + sum_m t_i0m w_i0urm - M(1-x_i0ur)
for u in U:
for r in R:
for i in V:
m.addConstr(
t_iur[0,u,r] >= t_iur[i,u,r] + s_i[i] + t_ij[(i,0)] +
- BIGM * (1 - x[i,0,u,r]),
name=f"return_base_time_u{u}_r{r}_from_{i}"
)
# t_(0u(r+1)) ≥ t_0ur - M(1-y_u(r+1))
for u in U:
for idx in range(len(R)-1):
r = R[idx]
rnext = R[idx+1]
m.addConstr(
t_iur[0,u,rnext] >= t_iur[0,u,r] - BIGM * (1 - y_ur[u,rnext]),
name=f"time_base_chain_u{u}_r{r}"
)
# 电量消耗与换电：
# e_jur ≤ e_iur - E_ij - sum_m c_ijm * w_ijurm - SE_i + M(1 - x_ijur)
# 文档中含 SE_i（换电或停靠消耗/补充量），若无可设为0
SE_i = {i: 0.0 for i in N} # 若需要修改请填入
for i in V:
for j in V:
if i==j: continue
for u in U:
for r in R:
m.addConstr(
e_iur[j,u,r] <= e_iur[i,u,r] - E_ij[(i,j)]
- SE_i[i]
+ BIGM * (1 - x[i,j,u,r]),
name=f"energy_consume_{i}_{j}_u{u}_r{r}"
)
# 启动出发时电量满电
for u in U:
for r in R:
m.addConstr(
e_iur[0,u,r] == E_max * y_ur[u,r],
name=f"full_charge_at_base_u{u}_r{r}"
)
# 变量上下界约束：0 ≤ e_iur ≤ E_max
for i in N:
for u in U:
for r in R:
m.addConstr(e_iur[i,u,r] <= E_max, name=f"e_max_{i}_{u}_{r}")
m.addConstr(e_iur[i,u,r] >= 0, name=f"e_min_{i}_{u}_{r}")
# ==========================
# ④ 加权和目标函数
# ==========================
alpha = 10000 # 无人机数量权重（大）
beta = 1 # 距离权重（小）
total_distance = quicksum(E_ij[(i,j)] * x[i, j, u, r]
for i in N for j in N if i != j
for u in U for r in R)
m.setObjective(alpha * quicksum(y_u[u] for u in U) + beta * total_distance, GRB.MINIMIZE)
# ==========================
# ⑤ 求解
# ==========================
m.Params.OutputFlag = 1
m.Params.TimeLimit = 120
m.optimize()
# ==========================
# ⑥ 结果打印
# ==========================
if m.status in (GRB.OPTIMAL, GRB.TIME_LIMIT):
print("\n===== 最终结果 =====")
print("启用无人机数 :", int(sum(y_u[u].X > 0.5 for u in U)))
print("总行程数 :", int(sum(y_ur[u, r].X > 0.5 for u in U for r in R)))
print("总路径距离（米）:", sum(d_ij[(i, j)] * x[i, j, u, r].X
for i in N for j in N if i != j
for u in U for r in R))
for u in U:
if y_u[u].X < 0.5:
continue
print(f"\n=== UAV {u} ===")
for r in R:
if y_ur[u, r].X < 0.5:
continue
print(f"\n --- 行程 r = {r} ---")
print(" 使用的弧：")
for i in N:
for j in N:
if i != j and x[i, j, u, r].X > 0.5:
print(f" {i} -> {j}")
else:
print("未求得可行解，写冲突文件")
m.computeIIS()
m.write("conflict.ilp")
这模型有什么问题吗解不出来
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
