CP-SAT求解器进度与最优性分析

本文详细阐述了如何准确衡量 CP-SAT 求解器的优化进度，特别是通过 `ObjectiveValue` 和 `BestObjectiveBound` 计算最优性间隙。文章分析了简单比率法的局限性，并引入了适用于正负目标值的通用间隙计算公式，同时提供了代码示例和关键注意事项，帮助用户更专业地评估求解器性能。

理解求解器进度与最优性间隙

在约束规划（CP）和混合整数规划（MIP）领域，评估求解器的优化进度是至关重要的。通常，我们希望了解当前找到的最佳解距离理论上的最优解还有多远。CP-SAT 求解器提供了 ObjectiveValue()（当前找到的最佳可行解的目标值）和 BestObjectiveBound()（目标值的最佳理论界限）这两个关键指标来辅助这一评估。

初始方法的局限性

一种直观的进度衡量方法是计算目标值与最佳界限的比率，例如 100 * ObjectiveValue() / BestObjectiveBound()。然而，这种方法存在显著局限性：

1. 仅适用于正值目标函数： 当目标函数值始终为正时（例如，最大化布尔变量乘以正系数的和），此比率可能提供一个粗略的进度指示。
2. 负系数或负目标值： 一旦引入负系数或允许目标值变为负数，此比率将失效。例如，当目标值为负时，一个接近 0 的负目标值可能代表一个糟糕的解，但其比率可能看起来“更好”。当最佳界限为负时，比率的解释也会变得复杂且不直观。
3. 零值问题： 如果 BestObjectiveBound() 或 ObjectiveValue() 为零，比率计算将导致除以零错误或无意义的结果。

这些局限性表明，需要一种更通用、更鲁棒的方法来计算求解器进度，即最优性间隙（Optimality Gap）。

最优性间隙的定义与计算

最优性间隙是衡量当前最佳可行解与最佳理论界限之间差距的标准指标。它通常表示为百分比，反映了当前解距离最优解的相对距离。为了处理各种目标函数值（正、负、零），MIP 求解器通常采用以下通用公式来定义间隙：

通用最优性间隙公式： 对于最小化问题，目标是找到最小的 ObjectiveValue，而 BestObjectiveBound 是目标值的下限。 对于最大化问题，目标是找到最大的 ObjectiveValue，而 BestObjectiveBound 是目标值的上限。

一个鲁棒的间隙计算方法需要考虑目标值的符号和零值情况。参考商业求解器（如 CPLEX）的做法，一个广泛接受的、对符号和零值健壮的间隙公式如下：

$$ \text{Gap} = \frac{|\text{BestObjective} - \text{BestBound}|}{10^{-10} + |\text{BestObjective}|} $$

其中：

• BestObjective 指的是 solver.ObjectiveValue()，即当前找到的最佳可行解的目标值。
• BestBound 指的是 solver.BestObjectiveBound()，即目标值的最佳理论界限。
• 10^{-10} (或一个小的正数 epsilon) 是为了防止当 BestObjective 为零时出现除以零的错误，同时避免在 BestObjective 接近零时导致间隙值过大。

解释：

• 分子 |BestObjective - BestBound| 表示当前最佳解与最佳界限之间的绝对差值。
• 分母 10^{-10} + |BestObjective| 将这个绝对差值归一化，使其成为一个相对百分比。使用 |BestObjective| 作为基准是因为它代表了当前找到的最佳解的“大小”，能够更好地反映相对误差。

CP-SAT 中的应用： 在 CP-SAT 中，solver.ObjectiveValue() 对应于 BestObjective，solver.BestObjectiveBound() 对应于 BestBound。无论目标是最小化还是最大化，这个公式都能提供一个一致的间隙度量。

示例代码

以下 Python 代码片段展示了如何在 CP-SAT 求解过程中计算并监控最优性间隙：

from ortools.sat.python import cp_model

def solve_and_monitor_gap():
    model = cp_model.CpModel()

    # 声明变量
    x = model.NewIntVar(-10, 10, 'x')
    y = model.NewIntVar(-10, 10, 'y')

    # 添加约束
    model.Add(x + y >= 5)
    model.Add(x - y <= 3)

    # 定义目标函数 (这里以最小化为例，但公式对最大化也适用)
    # 尝试一个可能产生负值或零的目标函数
    model.Minimize(3 * x - 2 * y)

    solver = cp_model.CpSolver()
    solver.parameters.log_search_progress = True # 打印求解进度
    solver.parameters.num_workers = 8 # 可以根据需要调整并行工作数

    # 定义一个自定义的解决方案回调函数来监控进度
    class GapMonitor(cp_model.CpSolverSolutionCallback):
        def __init__(self):
            cp_model.CpSolverSolutionCallback.__init__(self)
            self.__solution_count = 0
            self.__initial_bound = None # 可以选择记录初始界限

        def on_solution_callback(self):
            self.__solution_count += 1
            current_obj = self.ObjectiveValue()
            best_bound = self.BestObjectiveBound()

            # 记录初始界限（可选）
            if self.__initial_bound is None:
                self.__initial_bound = best_bound

            # 计算最优性间隙
            # 使用一个小的epsilon值来避免除以零，并处理目标值接近零的情况
            epsilon = 1e-10
            if abs(current_obj) < epsilon and abs(best_bound) < epsilon:
                gap = 0.0 # 如果两者都接近零，则认为间隙为零
            else:
                gap = abs(current_obj - best_bound) / (epsilon + abs(current_obj))

            print(f"Solution #{self.__solution_count}:")
            print(f"  ObjectiveValue: {current_obj}")
            print(f"  BestObjectiveBound: {best_bound}")
            print(f"  Optimality Gap: {gap:.4f} ({(gap * 100):.2f}%)")
            print("-" * 30)

    # 创建并注册回调
    monitor = GapMonitor()
    status = solver.Solve(model, monitor)

    print("\nSolver finished.")
    print(f"Solver status: {solver.StatusName(status)}")
    if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
        print(f"Final Objective Value: {solver.ObjectiveValue()}")
        print(f"Final Best Objective Bound: {solver.BestObjectiveBound()}")
        final_obj = solver.ObjectiveValue()
        final_bound = solver.BestObjectiveBound()
        epsilon = 1e-10
        if abs(final_obj) < epsilon and abs(final_bound) < epsilon:
            final_gap = 0.0
        else:
            final_gap = abs(final_obj - final_bound) / (epsilon + abs(final_obj))
        print(f"Final Optimality Gap: {final_gap:.4f} ({(final_gap * 100):.2f}%)")
    else:
        print("No solution found or problem is infeasible/unbounded.")

if __name__ == '__main__':
    solve_and_monitor_gap()

注意事项与总结

1. 间隙并非时间预测： 最优性间隙可以很好地衡量当前解的质量，但它并不能准确预测求解器还需要多少时间才能找到最优解或达到更小的间隙。求解器的搜索路径是非线性的，间隙的收敛速度可能在不同阶段差异很大。
2. 最佳界限的动态性： BestObjectiveBound() 也会在求解过程中不断改进。有时，界限的改进可能比目标值的改进更显著。为了更全面地反映进度，可以考虑记录初始的 BestObjectiveBound()，并与当前的界限进行比较，以了解界限本身的收敛情况。
3. 数值稳定性： 在计算间隙时，使用一个小的正数 epsilon (例如 1e-10) 作为分母的一部分至关重要，以避免在 ObjectiveValue() 为零时出现除以零的错误。
4. 最大化与最小化： 上述通用间隙公式对最大化和最小化问题都适用，因为它关注的是目标值和界限之间的绝对差距，并将其相对化。对于最小化问题，ObjectiveValue 应该接近 BestObjectiveBound（从上方逼近）；对于最大化问题，ObjectiveValue 应该接近 BestObjectiveBound（从下方逼近）。

通过理解和正确应用最优性间隙的概念，并使用鲁棒的计算公式，开发者可以更准确、专业地评估 CP-SAT 求解器的性能和优化进度，无论目标函数是正、负还是零。