THE WAY YOU MOVE746 2024-06-24 22:33 采纳率: 0%

# 运筹优化，gurobi,python

from gurobipy import Model, GRB, quicksum, or_
import numpy as np
import gurobipy as gp

# 输入数据
B = 3 # 电动公交车辆数量

L = 20# 行程任务数量
F = L  # 充电节点数量
M = 4  # 分割点集合

τ_p =6 # 每位乘客收入
E_p=28
T_p=15
τ_e = 2  # 电量成本
c_null = 1  # 空驶电量成本
W = 100000  # 大数值
SOC_full = 100  # 初始电量
t_min, t_max = 10, 40  # 最小最大发车间隔，放宽限制
Z_initiate, Z_end = 0, 450# 首末班车发车时刻
Y = 200  # 固定成本参数
γ = 0.1  # 固定成本参数
t_age, β, r_0, S_dr = 10, 0.9, 50, 30  # 固定成本参数
θ = 0.3  # 加权延误时间成本参数

# 修改后的节点顺序：交替出现任务节点和充电节点
task_nodes = list(range(1, L + 1))
charging_nodes = list(range(L + 1, 2 * L + 1))
nodes = [val for pair in zip(task_nodes, charging_nodes) for val in pair]
#e = np.random.uniform(30,38, (H))  # 电量消耗
#p = np.random.randint(20, 50, (H))  # 乘客数量
p=[13, 50, 20, 86, 39, 31]
e=[13.28, 21.54, 18.21, 21.88, 19.75, 14.71]
ρ = [0.33,0.06, 0.33, 0.01, 0.13,0.14]  # 情景发生概率
t = [26.81659693, 87.24644013,40.70677871, 126.112585, 47.89251773, 30.89393045]

H = len(ρ)  # 情景数量
o, d = 0, 2 * L + 1  # 出发和到达节点
N = [o] + nodes + [d]
print(N)
# 初始化模型
model = Model("ElectricBusScheduling")
c_points = [0, 10, 30,50]  # Example custom points for charging time
a_points = [0, 30, 90, 100]  # Example custom points for SOC
print(type(c_points[0]))

# 定义后继节点字典
successors = {nodes[i]: nodes[i+1] for i in range(0, len(nodes)-1)}
successors.update({nodes[i+1]: nodes[i+2] if i+2 < len(nodes) else d for i in range(0, len(nodes)-1, 2)})
successors[o] = nodes[0]
successors[d] = [d]

# 定义变量
Z3 = model.addVars(charging_nodes, B, vtype=GRB.INTEGER, name="Z3")
Z4 = model.addVars(charging_nodes, B, vtype=GRB.INTEGER, name="Z4")
#E3 = model.addVars(charging_nodes, B,lb=0,ub=100, vtype=GRB.CONTINUOUS, name="SOC1")
#E4 = model.addVars(charging_nodes, B,lb=0,ub=100, vtype=GRB.CONTINUOUS, name="SOC2")
x = model.addVars(nodes, nodes, B, vtype=GRB.BINARY, name="x")
u = model.addVars(L, B, vtype=GRB.BINARY, name="u")

y1 = model.addVars(charging_nodes, B, M, vtype=GRB.CONTINUOUS, name="y1")
y2 = model.addVars(charging_nodes, B, M, vtype=GRB.CONTINUOUS, name="y2")
w1 = model.addVars(charging_nodes, B, M, vtype=GRB.BINARY, name="w1")
w2 = model.addVars(charging_nodes, B, M, vtype=GRB.BINARY, name="w2")
# 弧集合
arcs = set()
arcs1 = set()
# 添加起点到任务节点的弧
# 添加任务节点到充电节点的弧，充电节点到其他任务节点或终点的弧
for j in charging_nodes:
if j==i+20:
for j in charging_nodes:
if r+21>=j:
def has_arc_connection(node1, node2, arc_set):
return (node1, node2) in arc_set
print(arcs1)
#      if has_arc_connection(i, j, arcs1):
#       print(f"There is an arc connection between {i} and {j}")
#     else:
#      print(f"There is no arc connection between {i} and {j}")

#sorted_arcs = sorted(list(arcs1))
#print(sorted_arcs)
def calculate_charging_time(SOC_start, SOC_end, charging_rate):
# Placeholder logic for calculating charging time
if charging_rate == 0:
return 0
return (SOC_end - SOC_start) / charging_rate

# Define variables
charge_time = model.addVars(charging_nodes, B, vtype=GRB.CONTINUOUS, name="charge_time")
charging_rate = 10  # Example constant charging rate
# Define charging time calculation within the constraints

# 目标函数
income = quicksum((p[h] * T_p * ρ[h] * o[i, b]) for i in task_nodes for h in range(H) for b in range(B))
cost = quicksum((SOC2[a,b]-SOC1[a,b])* E_p * ρ[h] * o[a, b] + c_null * 5 * x[i, a, b] for i in task_nodes for a in charging_nodes for b in range(B) for h in range(H))
fixed_cost = (Y * (1 - γ) / (t_age * β * 365) + r_0 + S_dr) / M
continue  # 跳过最后一个点
delay_cost = quicksum((θ * p[h] * (Z[i+1] - Z[i]) * ρ[h])  for h in range(H) for b in range(B))
model.setObjective(income - cost - fixed_cost - delay_cost, GRB.MAXIMIZE)

for b in range(B):
a=successors[i]
model.addConstr(SOC2[i, b] >= 40 * (1 - x[i, a, b]))
model.addConstr(SOC1[a, b] >= 40 * (1 - o[a, b]))

#model.addConstrs((quicksum(Oa[i, b, h] for b in range(B)) == 1 for h in range(H)), name=f"TaskAssigned1_{i}")
#model.addConstrs((Oa[i, (i % B), h] == 1 for h in range(H)), name=f"TaskAssignmentOrder_{i}")
model.addConstr((Z[i + 1] - Z[i] >= t_min), name=f"DepartureInterval_LowerBound2_{i}")
model.addConstr((Z[i + 1] - Z[i] <= t_max), name=f"DepartureInterval_UpperBound2_{i}")

for a in charging_nodes:
#model.addConstrs((Oa[i, (i % B), h] == 1 for h in range(H)), name=f"TaskAssignmentOrder_{i}")

if i != j:
if i != r:
model.addConstr(((quicksum(x[j, i, b] for b in range(B))  == quicksum(x[i, r, b] for b in range(B)))), name=f"TaskAssigned1_{i}_{j}_{r}")
for a in charging_nodes:

model.addConstr(SOC[1, b] == SOC_full, name=f"SOC_Full1_{b}") #满电约束
for b in range(B):
#任务时刻更新
model.addConstr(Z[i] <= Z1[i, b]+W * (1 - o[i, b]) )
model.addConstr(Z[i] >= Z1[i, b]-W * (1 - o[i, b]) )

#开始时间电量更新
model.addConstr((SOC1[i, b] - E_p +W * (1 - o[i, b])>=SOC2[i,b]),name=f"SOC_UpperBound6_{i}_{b}")
model.addConstr((SOC1[i, b] -E_p -W * (1 - o[i, b] )<= SOC2[i,b]),name=f"SOC_LowerBound7_{i}_{b}")
model.addConstr((Z1[i,b] +15 +W * (1 - o[i, b]) >=Z2[i,b] ),name=f"Z2_LowerBound8_{i}_{b}")
model.addConstr((Z1[i,b] +15- W * (1 - o[i, b]) <= Z2[i,b]),name=f"Z2_UpperBound9_{i}_{b}")

if has_arc_connection(i, j, arcs1):
#print(b,i,j)
model.addConstr((Z2[i,b] +1-W * (1 - x[i, j, b]) <=Z[j] ),name=f"Z2_LowerBound8_{i}")
model.addConstr((Z2[i,b]+1+ W * (1 - x[i, j, b])>=Z[j]  ), name=f"Z2_LowerBound9_{i}")
model.addConstr((SOC2[i, b]-5+W * (1 - x[i, j, b])>=SOC1[j,b]),name=f"SOC2_LowerBound8_{i}")
model.addConstr((SOC2[i, b]-5- W * (1 - x[i, j, b]) <=SOC1[j,b]), name=f"SOC2_LowerBound9_{i}")

model.addConstr((Z2[i, b] + 1 - W * (1 - x[i, j, b]) <= Z1[j,b]), name=f"Z2_LowerBound8_{i}")
model.addConstr((Z2[i, b] + 1 + W * (1 - x[i, j, b]) >= Z1[j,b]), name=f"Z2_LowerBound9_{i}")
model.addConstr((SOC2[i, b] - 5 + W * (1 - x[i, j, b]) >= SOC1[j, b]),name=f"SOC2_LowerBound8_{i}")
model.addConstr((SOC2[i, b] - 5 - W * (1 - x[i, j, b]) <= SOC1[j, b]), name=f"SOC2_LowerBound9_{i}")

for a in charging_nodes:
if has_arc_connection(i, a, arcs):
print(b, i, a)
model.addConstr((SOC2[i, b] - 5 + W * (1 - x[i, a, b]) >= SOC1[a, b]), name=f"SOC_UpperBound6_{i}")
model.addConstr((SOC2[i, b] - 5 - W * (1 - x[i, a, b]) <= SOC1[a, b]), name=f"SOC_LowerBound7_{i}")
model.addConstr((Z2[i, b] + 1 - W * (1 - x[i, a, b]) <= Z3[a, b]), name=f"Z2_LowerBound8_{i}")
model.addConstr((Z2[i, b] + 1 + W * (1 - x[i, a, b]) >= Z3[a, b]), name=f"Z2_LowerBound9_{i}")

model.addConstr(Z2[i, b] + 1 + W * (1 - x[i, a, b]) >= Z[a], name=f"Z_a_LowerBound7_{a}")
model.addConstr(Z2[i, b] + 1 - W * (1 - x[i, a, b]) <= Z[a], name=f"Z_a_UpperBound8_{a}")
model.addConstr(SOC2[i, b] - 5 + W * (1 - x[i, a, b]) >= SOC[a, b], name=f"SOC_a_LowerBound7_{a}_{b}")
model.addConstr(SOC2[i, b] - 5 - W * (1 - x[i, a, b]) <= SOC[a, b], name=f"SOC_a_UpperBound8_{a}_{b}")

model.addConstr(SOC1[a, b] + 30 + W * (1 - o[a, b]) >= SOC2[a, b], name=f"SOC_a_LowerBound7_{a}_{b}")
model.addConstr(SOC1[a, b] + 30 - W * (1 - o[a, b]) <= SOC2[a, b], name=f"SOC_a_LowerBound7_{a}_{b}")
model.addConstr(Z3[a, b] + 10 + W * (1 - o[a, b]) >= Z4[a, b], name=f"SOC_a_LowerBound7_{a}_{b}")
model.addConstr(Z3[a, b] + 10 - W * (1 - o[a, b]) <= Z4[a, b], name=f"SOC_a_LowerBound7_{a}_{b}")

model.addConstr(SOC[a,b] + 30 + W * (1 - o[a, b]) >= SOC2[a, b], name=f"SOC_a_LowerBound7_{a}_{b}")
model.addConstr(SOC[a,b] + 30 - W * (1 - o[a, b]) <= SOC2[a, b], name=f"SOC_a_LowerBound7_{a}_{b}")
model.addConstr(Z[a] + 10 + W * (1 - o[a, b]) >= Z4[a, b], name=f"SOC_a_LowerBound7_{a}_{b}")
model.addConstr(Z[a] + 10 - W * (1 - o[a, b]) <= Z4[a, b], name=f"SOC_a_LowerBound7_{a}_{b}")

if has_arc_connection(a, k, arcs):
print(b,i,a,k)
model.addConstr((SOC2[a, b] - 5 + W * (1 - x[a, k, b]) >= SOC1[k, b]),name=f"SOC_k_Update1_{a}_{b}")
model.addConstr((SOC2[a, b]- 5 - W * (1 - x[a,k ,b]) <= SOC1[k, b]),name=f"SOC_k_Update2_{a}_{b}")
model.addConstr((Z4[a, b] + 1 + W * (1 - x[a, k,b]) >= Z1[k, b]),name=f"Z_k_Update1_{a}_{b}")
model.addConstr((Z4[a, b] + 1 - W * (1 - x[a,k, b]) <= Z1[k, b]),name=f"Z_k_Update2_{a}_{b}")

#.addConstrs((SOC[a, a % B] - 5 + charge_time[a, a % B] * 2 + W * (1 - x[a, k, a % k % B]) >= SOC[k, k % B]),name=f"SOC1_k_UpperBound_{a}")
#model.addConstrs((SOC[a, a % B] - 5 + charge_time[a, a % B] * 2 - W * (1 - x[a, k, a % k % B]) <= SOC[k, k % B]),name=f"SOC1_k_LowerBound_{a}")
#model.addConstrs((Z[a] + 1 + charge_time[a, a % B] + W * (1 - x[a, k, a % k % B]) >= Z[k]),name=f"Z1_k_LowerBound_{a}")
#model.addConstrs((Z[a] + 1 + charge_time[a, a % B] - W * (1 - x[a, k, a % k % B]) <= Z[k]),name=f"Z1_k_UpperBound_{a}")

# Assuming you need to calculate and constrain charge_time based on SOC1, SOC2, and charging_rate
#model.addConstr(charge_time[a, a % B] == calculate_charging_time(SOC1[a, a % B], SOC2[a, a % B], charging_rate),name=f"ChargingTime_{a}_{b}")

# Optimize the model
model.optimize()

print(x[1,20,2].X)
print(o[7,1].X)
print(o[10,1].X)
# Check for constraint violations
def check_constraints():
violations = []
if model.status == GRB.OPTIMAL:
for constr in model.getConstrs():
slack = constr.getAttr(GRB.Attr.Slack)
if slack < 0:
violations.append(constr.ConstrName)
else:
print("Model did not reach an optimal solution. Status:", model.status)
return violations
if model.status == GRB.OPTIMAL:
for b in range(B):
for a in charging_nodes:
#print(a,b)
if x[i,a, b].X > 0.5:
charge_start = Z3[a,b].X
charge_end = Z4[a, b].X
soc_start = SOC1[a, b].X
soc_end = SOC2[a, b].X
print(f"  Charging at node {a}: Start time = {charge_start}, End time = {charge_end}, SOC start = {soc_start}, SOC end = {soc_end}")
#
# Output results if optimal solution is found
if model.status == GRB.OPTIMAL:
for b in range(B):
print(f"Bus {b }:")
if x[i,j, b].X > 0.5:
print(f"  Task {i}: Departure time = {Z1[i,b].X},{Z2[i,b].X},SOC:={SOC1[i,b].X},{SOC2[i,b].X}")

for i in charging_nodes:
print(Z[i])

# Save the model

model.computeIIS()
model.write(r"D:/model.ilp")
model.write(r"D:/electric_bus_scheduling_with_charging.lp")

#Check violated constraints()
#violated_constraints = check_constraints()
#if violated_constraints:
#    print(f"Violated constraints: {violated_constraints}")
# Check for infeasibility or unboundedness
if model.status == GRB.UNBOUNDED:
print("Model is unbounded. Checking constraints and variables...")
model.computeIIS()
model.write(r"D:/model.ilp")
print("Irreducible inconsistent subsystem written to 'model.ilp'")
elif model.status == GRB.INFEASIBLE:
print("Model is infeasible. Checking for conflicts...")
model.computeIIS()
model.write(r"D:/model.ilp")
print("Irreducible inconsistent subsystem written to 'model.ilp'")
elif model.status == GRB.OPTIMAL:
print("Optimal solution found.")
else:
print(f"Optimization ended with status {model.status}")


• 写回答

#### 26条回答默认 最新

• 阿里嘎多学长 2024-06-24 23:16
关注
获得0.30元问题酬金

以下内容由CHATGPT及阿里嘎多学长共同生成、有用望采纳：

您好！我很高兴能帮助您解决使用Python和Gurobi求解电动公交车最优行车链问题。首先，我会对您的代码进行一些基本的检查，并提供一些可能的改进方向。

1. 代码结构：您的代码结构看起来比较清晰，变量定义和模型构建逻辑也比较直观。这是一个很好的开始。

2. 数据输入：您已经定义了所有必要的输入数据，这是构建优化模型的基础。

3. 变量定义：您定义了很多决策变量，包括Z, SOC, x等。请确保每个变量的含义和作用都是清晰的，并且在模型中正确使用。

4. 约束条件：您的代码中包含了大量的约束条件。请仔细检查每个约束条件是否正确反映了问题的实际需求。例如，发车间隔的约束、电量更新的逻辑等。

5. 目标函数：您已经定义了目标函数，包括收入、成本和固定成本等。请确保目标函数的计算是正确的，并且符合您的问题设定。

6. 模型求解：您使用了model.optimize()来求解模型。请检查求解过程是否有错误，比如是否有未定义的变量或约束条件导致求解失败。

7. 结果检查：您的代码中包含了一些打印语句来输出结果，这是检查模型是否正确求解的好方法。请确保输出的结果与您的预期一致。

8. 调试和验证：作为初学者，调试和验证模型的正确性是非常重要的。您可以通过设置一些简单的测试案例来验证模型的输出是否合理。

9. 代码优化：您的代码中可能存在一些可以优化的地方，比如使用循环来减少重复代码，或者使用更高效的数据结构来存储和处理数据。

10. 参考链接：以下是一些可能对您有帮助的参考链接：

请注意，这些链接可能需要您根据具体问题进行筛选和深入研究。

最后，如果您在代码中遇到具体的错误或问题，欢迎您提供更详细的信息，我会尽力为您提供更具体的帮助。祝您在运筹优化的学习道路上一切顺利！

评论 编辑记录

• 系统已结题 7月2日
• 创建了问题 6月24日

#### 悬赏问题

• ¥15 咨询一下有关于王者荣耀赢藏战绩
• ¥100 求购一套带接口实现实习自动签到打卡
• ¥50 MacOS 使用虚拟机安装k8s