Python与C#中Google OrTools求解结果不一致问题排查

在使用Google OR-Tools LinearSolver时，Python和C#实现的主问题求解逻辑，输入的fsi、zsl、ol数据已逐行验证完全一致，但求解得到的objectiveValue、alphas和duals结果却不同。

数据详情 #

• fsi：300×300的对角矩阵（主对角线为1，其余为0），类型为numpy.ndarray
• zsl：300×322的numpy.ndarray矩阵
• ol：长度为322的numpy.ndarray数组

Python实现代码 #

def master_problem(fsi, zsl, ol, tlim=None, relaxation=True, enable_output=False):
    """
    Solve the master problem, either in its full version (MP)
    or in its relaxed version (RMP). Returns the following:
    - Objective value: minimization of sum(alpha[l] * h[l]), with h heights and l layer
    - Alpha values: alpha[l] represents layer selection
    - [RMP] Duals: one dual for each item
    """
    #logger.info("RMP defining variables and constraints")
    
    # Solver
    if relaxation:
        slv = pywraplp.Solver("RMP", pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
    else:
        slv = pywraplp.Solver("MP", pywraplp.Solver.BOP_INTEGER_PROGRAMMING)

    # Enable verbose output from solver
    if enable_output:
        slv.EnableOutput()

    # Utility
    infinity = slv.infinity()
    n_layers = 300
    n_items = fsi.shape[-1]

    # Variables
    if relaxation:
        al = [slv.NumVar(0, infinity, f"alpha_{l}") for l in range(n_layers)]
    else:
        al = [slv.BoolVar(f"alpha_{l}") for l in range(n_layers)]

    # Constraints
    constraints = []
    coefficients = np.matmul(fsi.T, zsl)

    # Select each item at least once
    # sum(al[l] * zsl[s, l] * fsi[s, i])
    for i in range(n_items):
        c = slv.Constraint(1, infinity, f"c_{i}")
        for l in range(n_layers):
            if coefficients[i, l] > 0:
                c.SetCoefficient(al[l], float(coefficients[i, l]))
        constraints += [c]

    # Objective
    obj = slv.Objective()
    for l, h in enumerate(ol):
        #logger.warning(str(ol) + " | " + str(h))
        obj.SetCoefficient(al[l], float(h))
    obj.SetMinimization()

    # Set a time limit in milliseconds
    if tlim is not None:
        slv.SetTimeLimit(1000 * tlim)

    # Solve
    #logger.debug(f"RMP variables: {slv.NumVariables()}")
    #logger.debug(f"RMP constraints: {slv.NumConstraints()}")
    status = slv.Solve()
    #logger.debug(f"RMP iterations: {slv.iterations()}")

    # Extract results
    duals, alphas = None, None
    objective = float("inf")
    logger.debug(f"RMP alphas: {alphas}")
    if status in (slv.OPTIMAL, slv.FEASIBLE):
        logger.info(f"RMP solved")

        # Extract alpha values
        alphas = [al[l].solution_value() for l in range(n_layers)]
        logger.debug(f"RMP alphas: {alphas}")
        if not all(alphas[l] in (0, 1) for l in range(n_layers)):
            logger.debug("RMP solution not feasible (at least one alpha value is not binary)")

        # Extract objective value
        objective = slv.Objective().Value()
        #logger.debug(f"RMP objective: {objective}")

        # Extract duals
        if relaxation:
            duals = np.array([c.DualValue() for c in constraints])
            #logger.debug(f"RMP duals: {duals}")
    else:
        logger.warning("RMP unfeasible")
    
    return objective, alphas, duals

C#实现代码 #

public static (double? objective, List<double>? alphas, List<double>? duals) MasterProblem(double[,]fsi, double[,]zsl, double[]ol, double? tlim = null, bool relaxation = true, bool enableOutput = false)
{
            // Prepare data
            var nLayers = layerPool.Count;
            var nItems = fsi.GetLength(0);

            /*Solve the master problem, either in its full version (MP)
            or in its relaxed version (RMP). Returns the following:
            - Objective value: minimization of sum(alpha[l] * h[l]), with h heights and l layer
            - Alpha values: alpha[l] represents layer selection
            - [RMP] Duals: one dual for each item*/

            // Create solver
            Solver solver;
            if (relaxation)
            {
                solver = Solver.CreateSolver("GLOP");
            }
            else
            {
                solver = Solver.CreateSolver("BOP");
            }

            // Enable verbose output from solver
            if (enableOutput)
            {
                solver.EnableOutput();
            }

            // Utility
            // infinity = slv.infinity() <- In CSharp anders gelöst. Wird direkt in Variable verwendet (z.B. double.PositiveInfinity)

            // Variables
            var al = new Variable[nLayers];
            if (relaxation)
            {
                for (int l = 0; l < nLayers; l++)
                {
                    al[l] = solver.MakeNumVar(0.0, double.PositiveInfinity, $"alpha_{l}");
                }
            }
            else
            {
                for (int l = 0; l < nLayers; l++)
                {
                    al[l] = solver.MakeBoolVar($"alpha_{l}");
                }

            }

            // Constraints
            var constraints = new List<Google.OrTools.LinearSolver.Constraint>();
            var fsiMatrix = Matrix<double>.Build.DenseOfArray(fsi);
            var zslMatrix = Matrix<double>.Build.DenseOfArray(zsl);
            var coefficients = fsiMatrix.Transpose() * zslMatrix;
          
            for (int i = 0; i < nItems; i++)
            {
                var c = solver.MakeConstraint(1, double.PositiveInfinity, $"c_{i}");
                for (int l = 0; l < nLayers; l++)
                {
                    if (coefficients[i, l] > 0)
                    {
                        c.SetCoefficient(al[l], coefficients[i, l]);
                    }
                }
                constraints.Add(c);
            }

            // Objective
            var objective = solver.Objective();
            for (int l = 0; l < ol.Count(); l++)
            {
                objective.SetCoefficient(al[l], ol[l]);
            }
            objective.SetMinimization();
            //var objective = solver.Objective();
            //for (int l = 0; l < nLayers; l++)
            //{
            //    objective.SetCoefficient(al[l], ol[l]);
            //}
            //objective.SetMinimization();

            // Set time limit
            if (tlim.HasValue)
            {
                solver.SetTimeLimit((long)(tlim.Value * 1000));
            }

            // Solve
            //Console.WriteLine($"RMP variables: {solver.NumVariables()}");
            //Console.WriteLine($"RMP constraints: {solver.NumConstraints()}");
            var status = solver.Solve();
            //Console.WriteLine($"RMP iterations: {solver.Iterations()}");

            // Extract results
            List<double> alphas = null;
            List<double> duals = null;
            double objectiveValue = double.PositiveInfinity;

            if (status == Solver.ResultStatus.OPTIMAL || status == Solver.ResultStatus.FEASIBLE)
            {
                //Console.WriteLine("RMP solved");

                // Extract alpha values
                alphas = al.Select(a => a.SolutionValue()).ToList();

                // Log alphas
                //Console.WriteLine($"RMP alphas: {string.Join(", ", alphas)}");

                // Check if all alphas are 0 or 1
                if (!alphas.All(alpha => alpha == 0 || alpha == 1))
                {
                    Console.WriteLine("RMP solution not feasible (at least one alpha value is not binary)");
                }

                // Extract objective value
                objectiveValue = objective.Value();

                // Extract duals
                if (relaxation)
                {
                    duals = constraints.Select(c => c.DualValue()).ToList();
                }
            }
            else
            {
                Console.WriteLine("RMP unfeasible");
            }

            return (objectiveValue, alphas, duals);
        }

核心排查方向 #

1. 维度匹配与索引错误 #

• Python中n_layers固定为300，但目标函数遍历的是长度322的ol，会出现索引越界风险（当l>=300时），需确认实际运行时ol长度是否为300；C#中用ol.Count()循环，若nLayers是300而ol长度322，同样会触发索引越界，这是代码中的明显问题，优先验证。
• C#中nLayers = layerPool.Count，必须确认layerPool.Count是否严格等于300，否则变量数量与Python不一致，直接导致结果不同。

2. 矩阵计算差异 #

• Python用np.matmul(fsi.T, zsl)计算系数矩阵，C#用MathNet的fsiMatrix.Transpose() * zslMatrix，需确认：
  • 矩阵乘法的维度是否正确：fsi是300×300，zsl是300×322，转置后fsi.T为300×300，matmul结果应为300×322；C#中矩阵乘法是否得到相同维度的结果。
  • 浮点数精度差异：两种语言的浮点数计算精度可能导致系数矩阵中微小数值差异，比如接近0的数在Python中被过滤，C#中被保留，可将判断条件改为abs(coefficients[i,l]) > 1e-9规避。

3. 求解器配置与版本 #

• 确认Python和C#使用的OR-Tools版本是否一致，不同版本的GLOP求解器可能存在行为差异。
• 时间限制设置：Python和C#的时间限制单位是否均为毫秒，若tlim不为null，需验证转换逻辑一致。

4. 约束与变量的构建顺序 #

• 检查约束的添加顺序是否完全一致（i从0到n_items的循环顺序），确保dual值的提取顺序匹配。
• 变量创建顺序是否与Python一致，保证alpha值的提取顺序无偏差。

内容的提问来源于stack exchange，提问作者I.Weix