CP-SAT求解带容量VRP时车辆冗余问题及优化方案问询

我正在用CP-SAT求解器解决一个含6辆可用车辆、17个节点的带容量约束车辆路径问题（VRP）。目前求解结果用了全部6辆车，但从节点总需求和车辆容量来看，明明用更少车辆就能满足需求。请问该如何调整模型，选择合适的车辆数量来最小化目标函数？

现有求解代码 #

from ortools.sat.python import cp_model as cp

locations = range(17)

dist_matrix = [
        [
            0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354,
            468, 776, 662
        ],
        [
            548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674,
            1016, 868, 1210
        ],
        [
            776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164,
            1130, 788, 1552, 754
        ],
        [
            696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822,
            1164, 560, 1358
        ],
        [
            582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708,
            1050, 674, 1244
        ],
        [
            274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628,
            514, 1050, 708
        ],
        [
            502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856,
            514, 1278, 480
        ],
        [
            194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320,
            662, 742, 856
        ],
        [
            308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662,
            320, 1084, 514
        ],
        [
            194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388,
            274, 810, 468
        ],
        [
            536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764,
            730, 388, 1152, 354
        ],
        [
            502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114,
            308, 650, 274, 844
        ],
        [
            388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194,
            536, 388, 730
        ],
        [
            354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0,
            342, 422, 536
        ],
        [
            468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536,
            342, 0, 764, 194
        ],
        [
            776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274,
            388, 422, 764, 0, 798
        ],
        [
            662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730,
            536, 194, 798, 0
        ],
]

distance_matrix = {}
for i in locations:
    for j in locations:
        distance_matrix[(i, j)] = dist_matrix[i][j]

num_vehicles = 6
demands = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
vehicle_capacity = 30
depot = 0

model = cp.CpModel()

# Build decision variables :
dv = {}
for i in locations:
    for j in locations:
        for k in range(num_vehicles):
            dv[(i, j, k)] = model.NewBoolVar("from_%i_to_%i_vehc_%i" % (i, j, k))

# Objective function
total_distance_route = sum(dv[i, j, k] * distance_matrix[i, j] for i in locations for j in locations for k in range(num_vehicles))

model.Minimize(total_distance_route)

# Capacity constraint
for k in range(num_vehicles):
    model.Add(sum(dv[i, j, k] * demands[j] for i in locations for j in locations[1:]) <= vehicle_capacity)

# no travel from a node to itself
for k in range(num_vehicles):
    model.Add(dv[0, 0, k] == 0)

# which nodes were covered in a tour by a vehicle 
dv_vehc_node_bin = {}
for k in range(num_vehicles):
    lst = {}
    for i in locations[1:]:
        lst[i] = model.NewBoolVar("vehc_%i_node_%i" % (k, i))
    dv_vehc_node_bin[k] = lst

# constraint: each vertex is in exactly one tour
for v in locations[1:]:
    model.add(sum(dv_vehc_node_bin[i][v] for i in range(num_vehicles)) == 1)

# sub-tour elimination
arcs_dict = {}
for k in range(num_vehicles):
    arcs_list = []
    for i in locations:
        for j in locations:
            if (i == j) & (i == 0):
                arcs_list.append([0, 0, dv[0, 0, k]])
            elif (i == j) & (i > 0):
                arcs_list.append([i, i, ~dv_vehc_node_bin[k][i]])
                
            else:
                arcs_list.append([i, j, dv[(i, j, k)]])
    arcs_dict[k] = arcs_list

for a in arcs_dict:
    model.AddCircuit(arcs_dict[a])
    

# solve the model
solver = cp.CpSolver() 
solver.parameters.num_search_workers = 8
solver.parameters.log_search_progress = True
solver.log_callback = print
status = solver.Solve(model)

当前求解结果 #

得到的无环路径解决方案如下，但使用了全部6辆车辆：

0
[(0, 12), (1, 0), (3, 4), (4, 1), (11, 15), (12, 11), (15, 3)]

1
[(0, 14), (2, 6), (6, 8), (8, 0), (10, 2), (14, 16), (16, 10)]

2
[(0, 13), (13, 0)]

3
[(0, 5), (5, 0)]

4
[(0, 7), (7, 0)]

5
[(0, 9), (9, 0)]

调整方案：最小化车辆使用数量 #

1. 给车辆启用添加惩罚成本（推荐） #

在目标函数中加入车辆启用的固定惩罚，让求解器优先选择用更少的车辆，即使总距离略有增加。具体操作：

• 为每辆车定义布尔变量标记是否启用
• 给启用的车辆设置足够高的惩罚系数（需大于单辆车空载往返仓库的最大距离，确保节省车辆的收益远大于路径合并增加的成本）

修改后的代码片段：

# 新增车辆启用标记变量
used_vehicle = [model.NewBoolVar(f"used_vehicle_{k}") for k in range(num_vehicles)]

# 关联车辆启用状态与行驶路径：只要车辆有行驶动作，就标记为启用
for k in range(num_vehicles):
    model.Add(sum(dv[i, j, k] for i in locations for j in locations) >= 1).OnlyEnforceIf(used_vehicle[k])
    model.Add(sum(dv[i, j, k] for i in locations for j in locations) == 0).OnlyEnforceIf(used_vehicle[k].Not())

# 调整目标函数：总行驶距离 + 车辆启用惩罚（系数可根据实际距离规模调整）
vehicle_penalty = 10000
total_cost = total_distance_route + vehicle_penalty * sum(used_vehicle)
model.Minimize(total_cost)

2. 先计算理论最小车辆数，再逐步验证 #

先通过总需求计算理论最小车辆数，再从该数值开始逐步尝试求解，直到找到可行解：

• 总需求：sum(demands) = 60
• 车辆容量30，理论最小车辆数为 ceil(60/30) = 2
• 从2辆开始尝试，逐步增加车辆数，直到模型可行，再在可行的车辆数量中选择总距离最小的方案

3. 优化子回路约束（可选） #

当前子回路消除约束可以简化，同时可以添加约束避免车辆空载往返（除非没有节点需要服务），但核心优化还是通过惩罚成本推动求解器选择更少车辆。

内容的提问来源于stack exchange，提问作者Bhartendu Awasthi