Modeling error for sure

[lucazav]

Hi guys,

I'm trying this great package to solve the following MIP problem:

Here my data:

warehouse_supply
> [1] 50000 30000 40000 55000

country_demands
> [1] 40000 15000 25000 45000 25000 25000

cost_matrix
>      Italy France Germany Japan China USA
> [1,]     8     18      14    40    40  25
> [2,]    12     10       8    18    40  18
> [3,]    34     32      30    10    33  35
> [4,]    25     20      18    35    30  10

Here the model I'm trying to "translate" in ompr language:

model <- MIPModel() %>% 
    add_variable( x[i, j], i = 1:n_warehouses, j = 1:n_countries, type = "integer" ) %>% 
    set_objective( sum_expr(cost_matrix[i, j] * x[i, j], i = 1:n_warehouses, j = 1:n_countries), sense = 'min' ) %>% 
    add_constraint( sum_expr(x[i, j], j = 1:n_countries) <= warehouse_supply[i], i = 1:n_warehouses ) %>% 
    add_constraint( sum_expr(x[i, j], i = 1:n_warehouses) >= country_demands[j], j = 1:n_countries )

If I try to solve it using glpk, the solution is infeasible:

result <- model %>% 
    solve_model(with_ROI(solver = 'glpk'))
result
> Status: infeasible
> Objective value: 0

As I know the problem is feasible, there is something wrong in the way I'm writing constraints.

Could you help me in finding the error, please?

[lucazav]

Wow, I just discovered that adding the lower bound constraints to the x variables fixed it:

model <- MIPModel() %>% 
    add_variable( x[i, j], i = 1:n_warehouses, j = 1:n_countries, type = "integer", lb = 0 ) %>% 
    set_objective( sum_expr(cost_matrix[i, j] * x[i, j], i = 1:n_warehouses, j = 1:n_countries), sense = 'min' ) %>% 
    add_constraint( sum_expr(x[i, j], j = 1:n_countries) <= warehouse_supply[i], i = 1:n_warehouses ) %>% 
    add_constraint( sum_expr(x[i, j], i = 1:n_warehouses) >= country_demands[j], j = 1:n_countries )

[dirkschumacher]

That is interesting, that adding a lower bound makes the problem feasible and looks like a bug. Could you create a reprex using the reprex package?

[lucazav]

Here a working code:

n_countries <- 6
n_warehouses <- 4

country_demands <- c(40000, 15000, 25000, 45000, 25000, 25000)

warehouse_supply <- c(50000, 30000, 40000, 55000)

cost_matrix <- structure(c(8, 12, 34, 25, 18, 10, 32, 20, 14, 8, 30, 18, 40, 
            18, 10, 35, 40, 40, 33, 30, 25, 18, 35, 10), .Dim = c(4L, 6L), .Dimnames = list(
                c("Warehouse ITA", "Warehouse DEU", "Warehouse JPN", "Warehouse USA"
                ), c("Italy", "France", "Germany", "Japan", "China", "USA"
                )))


model <- MIPModel() %>% 
    add_variable( x[i, j], i = 1:n_warehouses, j = 1:n_countries, type = "integer" ) %>% 
    set_objective( sum_expr(cost_matrix[i, j] * x[i, j], i = 1:n_warehouses, j = 1:n_countries), sense = 'min' ) %>% 
    add_constraint( sum_expr(x[i, j], j = 1:n_countries) <= warehouse_supply[i], i = 1:n_warehouses ) %>% 
    add_constraint( sum_expr(x[i, j], i = 1:n_warehouses) >= country_demands[j], j = 1:n_countries )

result <- model %>% 
    solve_model(with_ROI(solver = 'glpk'))

result

[dirkschumacher]

The problem is that your linear program is unbounded without the lower bound on the variables. If the linear program is unbounded, then the solution code is infeasible and thus you cannot get a result. I hope this helps. Below is the output of GLPK signalling that the LP is unbounded.

library(ompr)
library(ompr.roi)
library(ROI.plugin.glpk)
n_countries <- 6
n_warehouses <- 4

country_demands <- c(40000, 15000, 25000, 45000, 25000, 25000)

warehouse_supply <- c(50000, 30000, 40000, 55000)

cost_matrix <- structure(c(
  8, 12, 34, 25, 18, 10, 32, 20, 14, 8, 30, 18, 40,
  18, 10, 35, 40, 40, 33, 30, 25, 18, 35, 10
), .Dim = c(4L, 6L), .Dimnames = list(
  c("Warehouse ITA", "Warehouse DEU", "Warehouse JPN", "Warehouse USA"), c("Italy", "France", "Germany", "Japan", "China", "USA")
))


model <- MIPModel() |>
  add_variable(x[i, j], i = 1:n_warehouses, j = 1:n_countries, type = "integer") |>
  set_objective(sum_expr(cost_matrix[i, j] * x[i, j], i = 1:n_warehouses, j = 1:n_countries), sense = "min") |>
  add_constraint(sum_expr(x[i, j], j = 1:n_countries) <= warehouse_supply[i], i = 1:n_warehouses) |>
  add_constraint(sum_expr(x[i, j], i = 1:n_warehouses) >= country_demands[j], j = 1:n_countries)

result <- model |>
  solve_model(with_ROI(solver = "glpk", verbose = TRUE))
#> <SOLVER MSG>  ----
#> GLPK Simplex Optimizer, v4.65
#> 10 rows, 24 columns, 48 non-zeros
#>       0: obj =   0.000000000e+00 inf =   1.750e+05 (6)
#>       8: obj =   2.325000000e+06 inf =   0.000e+00 (0)
#> LP HAS UNBOUNDED PRIMAL SOLUTION
#> GLPK Integer Optimizer, v4.65
#> 10 rows, 24 columns, 48 non-zeros
#> 24 integer variables, none of which are binary
#> glp_intopt: optimal basis to initial LP relaxation not provided
#> <!SOLVER MSG> ----

result
#> Status: infeasible
#> Objective value: 0

Adding a lb = 0 to add_variable makes the problem feasible, as you noticed.

library(ompr)
library(ompr.roi)
library(ROI.plugin.glpk)
n_countries <- 6
n_warehouses <- 4

country_demands <- c(40000, 15000, 25000, 45000, 25000, 25000)

warehouse_supply <- c(50000, 30000, 40000, 55000)

cost_matrix <- structure(c(
  8, 12, 34, 25, 18, 10, 32, 20, 14, 8, 30, 18, 40,
  18, 10, 35, 40, 40, 33, 30, 25, 18, 35, 10
), .Dim = c(4L, 6L), .Dimnames = list(
  c("Warehouse ITA", "Warehouse DEU", "Warehouse JPN", "Warehouse USA"), c("Italy", "France", "Germany", "Japan", "China", "USA")
))


model <- MIPModel() |>
  add_variable(x[i, j], i = 1:n_warehouses, j = 1:n_countries, type = "integer", lb = 0) |>
  set_objective(sum_expr(cost_matrix[i, j] * x[i, j], i = 1:n_warehouses, j = 1:n_countries), sense = "min") |>
  add_constraint(sum_expr(x[i, j], j = 1:n_countries) <= warehouse_supply[i], i = 1:n_warehouses) |>
  add_constraint(sum_expr(x[i, j], i = 1:n_warehouses) >= country_demands[j], j = 1:n_countries)

result <- model |>
  solve_model(with_ROI(solver = "glpk", verbose = TRUE))
#> <SOLVER MSG>  ----
#> GLPK Simplex Optimizer, v4.65
#> 10 rows, 24 columns, 48 non-zeros
#>       0: obj =   0.000000000e+00 inf =   1.750e+05 (6)
#>       8: obj =   2.325000000e+06 inf =   0.000e+00 (0)
#> *    13: obj =   2.270000000e+06 inf =   0.000e+00 (0)
#> OPTIMAL LP SOLUTION FOUND
#> GLPK Integer Optimizer, v4.65
#> 10 rows, 24 columns, 48 non-zeros
#> 24 integer variables, none of which are binary
#> Integer optimization begins...
#> Long-step dual simplex will be used
#> +    13: mip =     not found yet >=              -inf        (1; 0)
#> +    13: >>>>>   2.270000000e+06 >=   2.270000000e+06   0.0% (1; 0)
#> +    13: mip =   2.270000000e+06 >=     tree is empty   0.0% (0; 1)
#> INTEGER OPTIMAL SOLUTION FOUND
#> <!SOLVER MSG> ----

result
#> Status: optimal
#> Objective value: 2270000