基于Google OR-Tools的无人机VRP能耗最优求解技术咨询

一、OR-Tools VRP核心逻辑梳理 #

OR-Tools路由模块默认通过启发式算法（含禁忌搜索、局部搜索、遗传算法等，可配置）求解VRP，核心是用成本函数引导搜索方向——默认目标是路径长度，你需要将成本函数替换为自定义的能耗模型。

二、关键修改：替换成本函数为能耗模型 #

你的需求核心是「直线飞行能耗低，转弯能耗高」，因此需要将两点间的移动成本定义为「直线飞行能耗 + 转弯能耗（若存在转弯）」。

1. 明确能耗计算规则 #

基于你提到的论文思路，先落地能耗公式：

• 直线飞行能耗：与飞行距离成正比，公式为 E_straight = k1 × 距离(i,j)，其中k1是单位距离能耗系数
• 转弯能耗：当无人机从点i→j→k时，转弯角度θ（向量ji与jk的夹角）越大能耗越高，公式为 E_turn = k2 × θ，其中k2是单位角度能耗系数；路径首尾的移动（起点→首个节点、最后节点→起点）无转弯，不计入该能耗

2. 在OR-Tools中实现自定义成本 #

OR-Tools的RoutingModel支持自定义弧成本（对应直线飞行能耗）和过渡成本（对应转弯能耗，依赖前序节点）。

步骤1：设置弧成本（直线飞行能耗）

用SetArcCostEvaluatorOfAllVehicles替换默认的距离成本为直线能耗：

import math
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp

# 假设x、y是各节点的坐标列表，k1是直线能耗系数
x = [0, 1, 2, 3]
y = [0, 1, 0, 1]
k1 = 0.5

def straight_energy_cost(from_index, to_index):
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    # 计算欧氏距离
    distance = math.hypot(x[to_node] - x[from_node], y[to_node] - y[from_node])
    # OR-Tools要求成本为整数，按精度缩放
    return int(k1 * distance * 100)

# 初始化管理器和路由模型
manager = pywrapcp.RoutingIndexManager(len(x), 1, 0)
routing = pywrapcp.RoutingModel(manager)

# 注册直线能耗成本函数
transit_callback_index = routing.RegisterTransitCallback(straight_energy_cost)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

步骤2：添加转弯能耗的过渡成本

转弯能耗依赖前序节点，需通过**维度（Dimension）**跟踪路径状态，再计算过渡成本：

k2 = 2.0  # 转弯能耗系数

# 创建维度用于跟踪路径中的前序节点
turn_dimension = routing.AddDimension(
    transit_callback=lambda i,j: 0,  # 基础成本为0，仅用于状态跟踪
    slack_max=0,
    capacity=0,
    fix_start_cumul_to_zero=True,
    name='TurnTracking'
)

# 定义过渡回调函数，计算转弯能耗
def transition_cost(from_node, to_node, vehicle_id):
    # 起点无前序节点，无转弯能耗
    if from_node == routing.Start(vehicle_id):
        return 0
    # 获取前序节点
    prev_node = routing.Start(vehicle_id) if from_node == routing.Start(vehicle_id) else from_node
    # 计算两个向量的夹角（弧度）
    vec1 = (x[from_node] - x[prev_node], y[from_node] - y[prev_node])
    vec2 = (x[to_node] - x[from_node], y[to_node] - y[from_node])
    
    # 避免除以0的边界情况
    if (vec1[0] == 0 and vec1[1] == 0) or (vec2[0] == 0 and vec2[1] == 0):
        angle = 0.0
    else:
        dot_product = vec1[0] * vec2[0] + vec1[1] * vec2[1]
        mag1 = math.hypot(vec1[0], vec1[1])
        mag2 = math.hypot(vec2[0], vec2[1])
        # 限制余弦值范围，避免数值误差导致acos报错
        cos_theta = max(min(dot_product / (mag1 * mag2), 1.0), -1.0)
        angle = math.acos(cos_theta)
    
    # 缩放为整数后返回
    return int(k2 * angle * 100)

# 绑定过渡成本到维度
turn_dimension.SetTransitionCallback(transition_cost)
# 将转弯能耗加入总优化目标
routing.AddVariableMinimizedByFinalizer(turn_dimension.CumulVar(routing.End(0)))

三、启用遗传算法求解 #

OR-Tools默认用禁忌搜索，你可以指定遗传算法作为求解器：

search_params = pywrapcp.DefaultRoutingSearchParameters()
# 设置为遗传算法
search_params.local_search_metaheuristic = routing_enums_pb2.LocalSearchMetaheuristic.GENETIC_ALGORITHM
# 调整算法参数
search_params.time_limit.seconds = 60  # 超时时间
search_params.population_size = 100    # 种群大小
search_params.mutation_rate = 0.1      # 突变率
search_params.crossover_rate = 0.8     # 交叉率

# 求解
solution = routing.SolveWithParameters(search_params)

四、关键注意事项 #

• OR-Tools要求成本为整数，需将能耗值缩放（如乘以100）后取整，避免浮点数精度问题
• 转弯能耗依赖路径历史节点，必须通过维度跟踪状态，无法仅用弧成本实现
• 若无人机有续航、载重等约束，可添加对应维度（如累计能耗不超过最大续航阈值）

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