search

ERGO-Code / HiGHS

Linear optimization software
MIT License
931 stars 173 forks source link

Feasible solution does not satisfy constraints. #1931

Closed lucidfrontier45 closed 1 day ago

lucidfrontier45 commented 1 day ago

Hi.

I implemented TSP with MTZ subtour elimination approach. I found that Highs reported feasible solutions but they violate the constraint and were actually infeasible. I implemented the model by two frameworks.

  1. Python + ORTools 9.11+ Highs 1.7.2
  2. Rust + good_lp 1.81 + Highs 1.6.1

Both of them show the same behaviour so I think it is not the problem of those frameworks but of Highs itself. In Python code if I use CP_SAT solver the result is correct so the modeling should also be valid.

Here are my one-file source codes. They are implemented exactly in the same way.

import datetime
from typing import Self

from ortools.math_opt.python import mathopt

class DistanceCalculator:
    def __init__(self, mat: list[list[int]]) -> None:
        self.mat = mat

    def __call__(self, i: int, j: int) -> int:
        return self.mat[i][j]

    def size(self) -> int:
        return len(self.mat)

    @classmethod
    def from_coordinates(cls, coords: list[tuple[float, float]]) -> Self:
        mat = [[0] * len(coords) for _ in range(len(coords))]
        for i, (x1, y1) in enumerate(coords):
            for j, (x2, y2) in enumerate(coords):
                dist = int(abs(x1 - x2) + abs(y1 - y2))
                mat[i][j] = dist
        return cls(mat)

def solve_tsp(dc: DistanceCalculator):
    n = dc.size()
    model = mathopt.Model()
    x = [[model.add_binary_variable() for _ in range(n)] for _ in range(n)]
    u = [model.add_integer_variable(lb=0, ub=n - 1) for _ in range(n)]

    obj = sum(dc(i, j) * x[i][j] for i in range(n) for j in range(n))
    model.minimize_linear_objective(obj)

    constraints = []

    # edge constraints
    for i in range(n):
        constraints.append(x[i][i] == 0)
        constraints.append(sum(x[i][j] for j in range(n)) == 1)
        constraints.append(sum(x[j][i] for j in range(n)) == 1)

    # start and end at 0
    constraints.append(u[0] == 0)
    for i in range(1, n):
        constraints.append(u[i] >= 1)

    # subtour elimination by MTZ formulation
    for i in range(1, n):
        for j in range(1, n):
            if i != j:
                constraints.append(u[i] - u[j] <= n * (1 - x[i][j]) - 1)

    # add all constraints to the model
    for c in constraints:
        model.add_linear_constraint(c)

    params = mathopt.SolveParameters(
        time_limit=datetime.timedelta(seconds=10), enable_output=True
    )
    result = mathopt.solve(model, mathopt.SolverType.HIGHS, params=params)

    assert result.termination.reason in (
        mathopt.TerminationReason.OPTIMAL,
        mathopt.TerminationReason.FEASIBLE,
    )

    x_vals = [
        [int(result.variable_values(x[i][j])) for j in range(n)] for i in range(n)
    ]

    for i in range(n):
        # error !!
        assert sum(x_vals[i][j] for j in range(n)) == 1
        assert sum(x_vals[j][i] for j in range(n)) == 1

    tour = [0]
    current = 0
    while True:
        j = x_vals[current].index(1)
        if j == 0:
            break
        tour.append(j)
        current = j

    total_distance = int(result.objective_value())

    return total_distance, tour

if __name__ == "__main__":
    coords = [
        (38.24, 20.42),
        (39.57, 26.15),
        (40.56, 25.32),
        (36.26, 23.12),
        (33.48, 10.54),
        (37.56, 12.19),
        (38.42, 13.11),
        (37.52, 20.44),
        (41.23, 9.10),
        (41.17, 13.05),
        (36.08, -5.21),
        (38.47, 15.13),
        (38.15, 15.35),
        (37.51, 15.17),
        (35.49, 14.32),
        (39.36, 19.56),
    ]

    dc = DistanceCalculator.from_coordinates(coords)

    total_distance, tour = solve_tsp(dc)
    print(f"Total distance: {total_distance}")
    print(f"Tour: {tour}")

use good_lp::{constraint, variable, Expression, ProblemVariables, Solution, SolverModel};

struct DistanceCalculator {
    mat: Vec<Vec<i32>>,
}

impl DistanceCalculator {
    pub fn call(&self, i: usize, j: usize) -> i32 {
        self.mat[i][j]
    }

    pub fn size(&self) -> usize {
        self.mat.len()
    }

    pub fn from_coordinates(coordinates: Vec<(f64, f64)>) -> Self {
        let n = coordinates.len();
        let mut mat = vec![vec![0; n]; n];
        for (i, (x1, y1)) in coordinates.iter().enumerate() {
            for (j, (x2, y2)) in coordinates.iter().enumerate() {
                let dist = (x1 - x2).abs() + (y1 - y2).abs();
                mat[i][j] = dist as i32;
            }
        }
        Self { mat }
    }
}

#[allow(clippy::needless_range_loop)]
fn solve_tsp(dc: &DistanceCalculator) -> (i32, Vec<usize>) {
    let n = dc.size();
    let mut variables = ProblemVariables::new();

    let mut x = vec![vec![]; n];
    for i in 0..n {
        for _ in 0..n {
            x[i].push(variables.add(variable().binary()));
        }
    }
    let mut u = vec![];
    for _ in 0..n {
        u.push(variables.add(variable().integer().clamp(0, n as i32 - 1)));
    }
    let mut obj: Expression = 0.into();
    for i in 0..n {
        for j in 0..n {
            obj += dc.call(i, j) * x[i][j];
        }
    }

    let mut constraints = vec![];

    // edge constraints
    for i in 0..n {
        constraints.push(constraint! {x[i][i] == 0});
        let s_row = (0..n).map(|j| x[i][j]).sum::<Expression>();
        constraints.push(constraint! {s_row == 1 });
        let s_col = (0..n).map(|j| x[j][i]).sum::<Expression>();
        constraints.push(constraint! {s_col == 1 });
    }

    // start and end at 0
    constraints.push(constraint! {u[0] == 0});
    for i in 1..n {
        constraints.push(constraint! {u[i] >= 1});
    }

    let mut problem = variables
        .minimise(obj.clone())
        .using(good_lp::solvers::highs::highs)
        .set_time_limit(10.0);
    problem.set_verbose(true);

    // subtour elimination by MTZ formulation
    for i in 1..n {
        for j in 1..n {
            if i != j {
                constraints.push(constraint! {
                    u[i] - u[j] <= n as i32 * (1 - x[i][j]) - 1
                });
            }
        }
    }

    // add all constraints to the problem
    for c in constraints {
        problem.add_constraint(c);
    }

    let solution = problem.solve().unwrap();

    let mut x_vals = vec![vec![0; n]; n];
    for i in 0..n {
        for j in 0..n {
            x_vals[i][j] = solution.value(x[i][j]) as i32;
        }
    }

    for i in 0..n {
        let mut s_row = 0;
        for j in 0..n {
            s_row += x_vals[i][j];
        }
        // error !!
        assert_eq!(s_row, 1);

        let mut s_col = 0;
        for j in 0..n {
            s_col += x_vals[j][i];
        }
        assert_eq!(s_col, 1);
    }

    // tour by following x_vals
    let mut tour = vec![0];
    let mut current = 0;
    loop {
        let next = x_vals[current]
            .iter()
            .enumerate()
            .find(|(_, &x)| x == 1)
            .unwrap()
            .0;
        if next == 0 {
            break;
        }

        tour.push(next);
        current = next;
    }

    assert_eq!(tour.len(), n);

    let obj_val = solution.eval(obj) as i32;

    (obj_val, tour)
}

fn main() {
    let coordinates = vec![
        (38.24, 20.42),
        (39.57, 26.15),
        (40.56, 25.32),
        (36.26, 23.12),
        (33.48, 10.54),
        (37.56, 12.19),
        (38.42, 13.11),
        (37.52, 20.44),
        (41.23, 9.10),
        (41.17, 13.05),
        (36.08, -5.21),
        (38.47, 15.13),
        (38.15, 15.35),
        (37.51, 15.17),
        (35.49, 14.32),
        (39.36, 19.56),
    ];
    let dc = DistanceCalculator::from_coordinates(coordinates);
    let (total_distance, tour) = solve_tsp(&dc);
    dbg!(total_distance, tour);
}
odow commented 1 day ago

I cannot reproduce from JuMP.

julia> using JuMP, HiGHS

julia> function solve_tsp(dc)
           n = size(dc, 1)
           model = Model(HiGHS.Optimizer)
           @variable(model, x[1:n, 1:n], Bin)
           @variable(model, 1 <= u[1:n] <= n - 1)
           @objective(model, Min, sum(dc .* x))
           for i in 1:n
               @constraint(model, x[i, i] == 0)
               @constraint(model, sum(x[i, :]) == 1)
               @constraint(model, sum(x[:, i]) == 1)
           end
           fix(u[1], 0; force = true)
           for i in 2:n, j in 2:n
               if i != j
                   @constraint(model, u[i] - u[j] <= n * (1 - x[i, j]) - 1)
               end
           end
           optimize!(model)
           x_vals = round.(Int, value.(x))
           for i in 1:n
               @assert sum(x_vals[i, :]) == 1
               @assert sum(x_vals[:, i]) == 1
           end
           tour = [1]
           current = 1
           while true
               j = findfirst(==(1), x_vals[current, :])
               if j === 1
                   break
               end
               push!(tour, j)
               current = j
           end
           total_distance = round(Int, objective_value(model))
           return total_distance, tour
       end
solve_tsp (generic function with 1 method)

julia> coords = [
           (38.24, 20.42),
           (39.57, 26.15),
           (40.56, 25.32),
           (36.26, 23.12),
           (33.48, 10.54),
           (37.56, 12.19),
           (38.42, 13.11),
           (37.52, 20.44),
           (41.23, 9.10),
           (41.17, 13.05),
           (36.08, -5.21),
           (38.47, 15.13),
           (38.15, 15.35),
           (37.51, 15.17),
           (35.49, 14.32),
           (39.36, 19.56),
       ]
16-element Vector{Tuple{Float64, Float64}}:
 (38.24, 20.42)
 (39.57, 26.15)
 (40.56, 25.32)
 (36.26, 23.12)
 (33.48, 10.54)
 (37.56, 12.19)
 (38.42, 13.11)
 (37.52, 20.44)
 (41.23, 9.1)
 (41.17, 13.05)
 (36.08, -5.21)
 (38.47, 15.13)
 (38.15, 15.35)
 (37.51, 15.17)
 (35.49, 14.32)
 (39.36, 19.56)

julia> function from_coordinates(coords)
           mat = zeros(length(coords), length(coords))
           for (i, (x1, y1)) in enumerate(coords)
               for (j, (x2, y2)) in enumerate(coords)
                   mat[i, j] = round(Int, abs(x1 - x2) + abs(y1 - y2))
               end
           end
           return mat
       end
from_coordinates (generic function with 1 method)

julia> dc = from_coordinates(coords)
16×16 Matrix{Float64}:
  0.0   7.0   7.0   5.0  15.0   9.0   7.0   1.0  14.0  10.0  28.0   6.0   5.0   6.0   9.0   2.0
  7.0   0.0   2.0   6.0  22.0  16.0  14.0   8.0  19.0  15.0  35.0  12.0  12.0  13.0  16.0   7.0
  7.0   2.0   0.0   7.0  22.0  16.0  14.0   8.0  17.0  13.0  35.0  12.0  12.0  13.0  16.0   7.0
  5.0   6.0   7.0   0.0  15.0  12.0  12.0   4.0  19.0  15.0  29.0  10.0  10.0   9.0  10.0   7.0
 15.0  22.0  22.0  15.0   0.0   6.0   8.0  14.0   9.0  10.0  18.0  10.0   9.0   9.0   6.0  15.0
  9.0  16.0  16.0  12.0   6.0   0.0   2.0   8.0   7.0   4.0  19.0   4.0   4.0   3.0   4.0   9.0
  7.0  14.0  14.0  12.0   8.0   2.0   0.0   8.0   7.0   3.0  21.0   2.0   3.0   3.0   4.0   7.0
  1.0   8.0   8.0   4.0  14.0   8.0   8.0   0.0  15.0  11.0  27.0   6.0   6.0   5.0   8.0   3.0
 14.0  19.0  17.0  19.0   9.0   7.0   7.0  15.0   0.0   4.0  19.0   9.0   9.0  10.0  11.0  12.0
 10.0  15.0  13.0  15.0  10.0   4.0   3.0  11.0   4.0   0.0  23.0   5.0   5.0   6.0   7.0   8.0
 28.0  35.0  35.0  29.0  18.0  19.0  21.0  27.0  19.0  23.0   0.0  23.0  23.0  22.0  20.0  28.0
  6.0  12.0  12.0  10.0  10.0   4.0   2.0   6.0   9.0   5.0  23.0   0.0   1.0   1.0   4.0   5.0
  5.0  12.0  12.0  10.0   9.0   4.0   3.0   6.0   9.0   5.0  23.0   1.0   0.0   1.0   4.0   5.0
  6.0  13.0  13.0   9.0   9.0   3.0   3.0   5.0  10.0   6.0  22.0   1.0   1.0   0.0   3.0   6.0
  9.0  16.0  16.0  10.0   6.0   4.0   4.0   8.0  11.0   7.0  20.0   4.0   4.0   3.0   0.0   9.0
  2.0   7.0   7.0   7.0  15.0   9.0   7.0   3.0  12.0   8.0  28.0   5.0   5.0   6.0   9.0   0.0

julia> total_distance, tour = solve_tsp(dc)
Running HiGHS 1.7.2 (git hash: 5ce7a2753): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 2e+01]
  Cost   [1e+00, 4e+01]
  Bound  [1e+00, 2e+01]
  RHS    [1e+00, 2e+01]
Presolving model
242 rows, 255 cols, 1110 nonzeros  0s
242 rows, 255 cols, 1110 nonzeros  0s
Objective function is integral with scale 1

Solving MIP model with:
   242 rows
   255 cols (240 binary, 0 integer, 0 implied int., 15 continuous)
   1110 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0               inf                  inf        0      0      0         0     0.0s
         0       0         0   0.00%   73.06666667     inf                  inf        0      0      0        44     0.0s
 R       0       0         0   0.00%   80              106               24.53%     1100     61     20       223     0.1s
 L       0       0         0   0.00%   81              91                10.99%     1585     47     20       306     0.4s
 L       0       0         0   0.00%   81              90                10.00%     1642     52     20       925     0.7s

20.0% inactive integer columns, restarting
Model after restart has 200 rows, 207 cols (192 bin., 0 int., 0 impl., 15 cont.), and 888 nonzeros

         0       0         0   0.00%   81              90                10.00%       19      0      0      1828     0.7s
         0       0         0   0.00%   81              90                10.00%       19     17      4      1851     0.7s

4.7% inactive integer columns, restarting
Model after restart has 191 rows, 198 cols (183 bin., 0 int., 0 impl., 15 cont.), and 843 nonzeros

         0       0         0   0.00%   82              90                 8.89%       20      0      0      2841     1.1s
         0       0         0   0.00%   82              90                 8.89%       20     15      2      2873     1.1s
 L       0       0         0   0.00%   82              89                 7.87%      832     21      2      4385     1.9s

Solving report
  Status            Optimal
  Primal bound      89
  Dual bound        89
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    89 (objective)
                    0 (bound viol.)
                    1.13242748512e-14 (int. viol.)
                    0 (row viol.)
  Timing            5.33 (total)
                    0.01 (presolve)
                    0.00 (postsolve)
  Nodes             2085
  LP iterations     27727 (total)
                    7264 (strong br.)
                    2021 (separation)
                    4351 (heuristics)

Nothing seems wrong with the solution:

julia> solution_summary(model)
* Solver : HiGHS

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "kHighsModelStatusOptimal"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : 8.90000e+01
  Objective bound    : 8.90000e+01
  Relative gap       : 0.00000e+00

* Work counters
  Solve time (sec)   : 5.33034e+00
  Simplex iterations : 27727
  Barrier iterations : -1
  Node count         : 2085

and the solution satisfies all constraints to 1e-13 or so:

julia> primal_feasibility_report(model)
Dict{Any, Float64} with 29 entries:
  x[1,7] + x[2,7] + x[3,7] + x[4,7] + x[5,7] + x[6,7] + x[7,7] + x[8,7] + x[9,7] + x[10,7] + x[11,7] + x[12,7] + x[13,7] + x[14,7] + x[15,7] + x[16,7] = 1                 => 1.44329e-15
  u[8] ≤ 15                                                                                                                                                                => 8.33111e-13
  x[16,3] binary                                                                                                                                                           => 1.13243e-14
  16 x[14,15] + u[14] - u[15] ≤ 15                                                                                                                                         => 1.06581e-14
  x[1,13] binary                                                                                                                                                           => 1.11022e-14
  x[4,2] binary                                                                                                                                                            => 1.13243e-14
  x[4,8] binary                                                                                                                                                            => 1.13243e-14
  x[1,16] + x[2,16] + x[3,16] + x[4,16] + x[5,16] + x[6,16] + x[7,16] + x[8,16] + x[9,16] + x[10,16] + x[11,16] + x[12,16] + x[13,16] + x[14,16] + x[15,16] + x[16,16] = 1 => 5.55112e-15
  x[6,1] + x[6,2] + x[6,3] + x[6,4] + x[6,5] + x[6,6] + x[6,7] + x[6,8] + x[6,9] + x[6,10] + x[6,11] + x[6,12] + x[6,13] + x[6,14] + x[6,15] + x[6,16] = 1                 => 1.44329e-15
  x[3,2] binary                                                                                                                                                            => 1.13243e-14
  x[1,1] + x[2,1] + x[3,1] + x[4,1] + x[5,1] + x[6,1] + x[7,1] + x[8,1] + x[9,1] + x[10,1] + x[11,1] + x[12,1] + x[13,1] + x[14,1] + x[15,1] + x[16,1] = 1                 => 3.77476e-15
  x[12,1] + x[12,2] + x[12,3] + x[12,4] + x[12,5] + x[12,6] + x[12,7] + x[12,8] + x[12,9] + x[12,10] + x[12,11] + x[12,12] + x[12,13] + x[12,14] + x[12,15] + x[12,16] = 1 => 5.55112e-15
  16 x[15,5] - u[5] + u[15] ≤ 15                                                                                                                                           => 9.23706e-14
  x[2,4] binary                                                                                                                                                            => 1.13243e-14
  x[1,4] + x[2,4] + x[3,4] + x[4,4] + x[5,4] + x[6,4] + x[7,4] + x[8,4] + x[9,4] + x[10,4] + x[11,4] + x[12,4] + x[13,4] + x[14,4] + x[15,4] + x[16,4] = 1                 => 1.13243e-14
  x[1,8] binary                                                                                                                                                            => 1.13243e-14
  x[6,7] binary                                                                                                                                                            => 1.44329e-15
  x[2,3] binary                                                                                                                                                            => 1.13243e-14
  16 x[16,3] - u[3] + u[16] ≤ 15                                                                                                                                           => 7.10543e-15
  16 x[4,8] + u[4] - u[8] ≤ 15                                                                                                                                             => 8.88178e-15
  x[1,13] + x[2,13] + x[3,13] + x[4,13] + x[5,13] + x[6,13] + x[7,13] + x[8,13] + x[9,13] + x[10,13] + x[11,13] + x[12,13] + x[13,13] + x[14,13] + x[15,13] + x[16,13] = 1 => 1.11022e-14
  16 x[3,2] - u[2] + u[3] ≤ 15                                                                                                                                             => 6.03961e-14
  x[16,1] binary                                                                                                                                                           => 1.13243e-14
  x[8,1] binary                                                                                                                                                            => 7.54952e-15
  x[1,1] + x[1,2] + x[1,3] + x[1,4] + x[1,5] + x[1,6] + x[1,7] + x[1,8] + x[1,9] + x[1,10] + x[1,11] + x[1,12] + x[1,13] + x[1,14] + x[1,15] + x[1,16] = 1                 => 2.22045e-16
  16 x[7,12] + u[7] - u[12] ≤ 15                                                                                                                                           => 1.7053e-13
  x[3,1] + x[3,2] + x[3,3] + x[3,4] + x[3,5] + x[3,6] + x[3,7] + x[3,8] + x[3,9] + x[3,10] + x[3,11] + x[3,12] + x[3,13] + x[3,14] + x[3,15] + x[3,16] = 1                 => 1.13243e-14
  x[8,1] + x[8,2] + x[8,3] + x[8,4] + x[8,5] + x[8,6] + x[8,7] + x[8,8] + x[8,9] + x[8,10] + x[8,11] + x[8,12] + x[8,13] + x[8,14] + x[8,15] + x[8,16] = 1                 => 7.54952e-15
  x[12,16] binary                                                                                                                                                          => 5.55112e-15
(89, [1, 13, 14, 15, 5, 11, 9, 10, 6, 7, 12, 16, 3, 2, 4, 8])
odow commented 1 day ago

Ah. @lucidfrontier45, you have a bug in your code:

[int(result.variable_values(x[i][j])) for j in range(n)] for i in range(n)

Note that int in Python rounds down, not to the nearest.

>>> int(0.9999)
0

See the docstring:

>>> help(int)

Help on class int in module builtins:

class int(object)
 |  int([x]) -> integer
 |  int(x, base=10) -> integer
 |  
 |  Convert a number or string to an integer, or return 0 if no arguments
 |  are given.  If x is a number, return x.__int__().  For floating point
 |  numbers, this truncates towards zero.
 |  
 |  If x is not a number or if base is given, then x must be a string,
 |  bytes, or bytearray instance representing an integer literal in the
 |  given base.  The literal can be preceded by '+' or '-' and be surrounded
 |  by whitespace.  The base defaults to 10.  Valid bases are 0 and 2-36.
 |  Base 0 means to interpret the base from the string as an integer literal.
 |  >>> int('0b100', base=0)

The important part is

For floating point numbers, this truncates towards zero.

This is not a bug in HiGHS.

lucidfrontier45 commented 1 day ago

@odow Thank you so much! Now both programs work find.