Strictly positive limits of a variable are not respected

[JosepFanals]

What version of OR-Tools and what language are you using?
Version: 9.8.3296
Language: Python 3.10

Which solver are you using (e.g. CP-SAT, Routing Solver, GLOP, BOP, Gurobi)
Three solvers have been tested: CBC, SCIP, GLOP

What operating system (Linux, Windows, ...) and version?
Windows 10

What did you do?
Together with @SanPen we are optimizing the operation of a hydro system coupled to an electrical grid. Each hydro plant has a variable called "spillage" that symbolizes the lost water exiting the system. This variable has to be positive as it makes no sense to have a magically positive inflow.

What did you expect to see
We were expecting to get non-negative spillages.

What did you see instead?
Solving the problem with OR-Tools yields a value of -90.0 for the spillage at node 3 at the first time step (NodeSpillage_3_0). We tried reproducing the same problem with PuLP and we got non-negative results which we deem as correct. For comparison purposes, we are attaching both LP and MPS models from OR-Tools and PuLP.
Spillage_LP_MPS.zip

[SanPen]

@lperron how is this not an issue?

[lperron]

we check solutions against the model. So the model is wrong.

[lperron]

so, apart from the fact that they are badly conditioned, the ortools and pulp models do not produce the same optimal objective value.

[lperron]

File        : '/Users/lperron/Downloads/Spillage_LP_MPS/spill_ortools.mps'
Solver      : GLOP_LINEAR_PROGRAMMING
Parameters  : 
Dimension   : 400 x 530

Initial problem: 400 rows, 530 columns, 1365 entries with magnitude in [1.000000e+00, 1.086960e+02]
Objective stats: 280 non-zeros, range [-1.020000e+02, 1.200000e+03]
Bounds stats: 930 non-zeros, range [-1.000000e+20, 1.000000e+20]
Parameters: use_preprocessing: true log_search_progress: true

Starting presolve...
FixedVariablePreprocessor                    : 400(0) rows, 510(-20) columns, 1276(-89) entries. (0.000009s)
SingletonPreprocessor                        : 400(0) rows, 450(-80) columns, 1216(-149) entries. (0.000035s)
ImpliedFreePreprocessor                      : 400(0) rows, 450(-80) columns, 1216(-149) entries. (0.000039s)
DoubletonFreeColumnPreprocessor              : 399(-1) rows, 450(-80) columns, 1215(-150) entries. (0.000019s)
DoubletonEqualityRowPreprocessor             : 359(-41) rows, 410(-120) columns, 1135(-230) entries. (0.000021s)
SingletonPreprocessor                        : 359(-41) rows, 409(-121) columns, 1134(-231) entries. (0.000025s)
ImpliedFreePreprocessor                      : 359(-41) rows, 409(-121) columns, 1134(-231) entries. (0.000030s)
DoubletonFreeColumnPreprocessor              : 358(-42) rows, 409(-121) columns, 1133(-232) entries. (0.000017s)
DoubletonEqualityRowPreprocessor             : 348(-52) rows, 399(-131) columns, 1113(-252) entries. (0.000012s)
SingletonPreprocessor                        : 348(-52) rows, 398(-132) columns, 1112(-253) entries. (0.000019s)
ImpliedFreePreprocessor                      : 348(-52) rows, 398(-132) columns, 1112(-253) entries. (0.000028s)
Reached fixed point after presolve pass #3
EmptyConstraintPreprocessor                  : 338(-62) rows, 398(-132) columns, 1112(-253) entries. (0.000005s)
ProportionalRowPreprocessor                  : 258(-142) rows, 398(-132) columns, 872(-493) entries. (0.000034s)
SingletonColumnSignPreprocessor              : 258(-142) rows, 398(-132) columns, 872(-493) entries. (0.000002s)
ScalingPreprocessor                          : 258(-142) rows, 398(-132) columns, 872(-493) entries. (0.000028s)

Presolved problem: 258 rows, 398 columns, 872 entries with magnitude in [3.207161e-02, 1.000000e+00]
Objective stats: 239 non-zeros, range [-3.100272e+01, 1.629456e+02]
Bounds stats: 775 non-zeros, range [-1.000000e+21, 1.000000e+21]

Starting basis: create from scratch.
Trying to remove 120 fixed variables from the initial basis.
The primal residual of the initial basis is above the tolerance, 2970.02 vs. 1e-08
The matrix with slacks has 258 rows, 656 columns, 1130 entries.
Number of basic infeasible variables: 30
Number of basic slack variables: 58
Number of basic variables at bound: 127
Number of basic fixed variables: 0
Number of basic free variables: 71
Number of super-basic variables: 0

Primal feasibility phase, iteration # 0, sum_primal_infeasibilities = 1.746502864883356E+05 [norms]
Primal feasibility phase, iteration # 65, sum_primal_infeasibilities = 1.081028648833557E+04
Primal feasibility phase, iteration # 130, sum_primal_infeasibilities = 1.081095939528713E+04
Primal feasibility phase, iteration # 170, sum_primal_infeasibilities = 1.081125816864771E+04 [check]
Current status: PRIMAL_INFEASIBLE
Primal infeasibility (bounds) = 1201.02
Primal residual |A.x - b| = 2970.02
Dual infeasibility (reduced costs) = 2.83632e-15
Dual residual |c_B - y.B| = 4.44089e-16


Final unscaled solution:
Primal objective (before moving primal/dual values) = 2.046684482774433E+09
Dual objective (before moving primal/dual values) = 1.999987647824382E+09
Primal objective (after moving primal/dual values) = 2.046684482774433E+09
Max. rhs perturbation = 16384.3
Max. cost perturbation = 2095.22
Max. primal infeasibility = 16384.3
Max. dual infeasibility = 0
Objective error <= 2046.72
status: IMPRECISE
objective: 2.04668e+09
iterations: 170
time: 0.001488
deterministic_time: 0.000815637

Status      : MPSOLVER_ABNORMAL
Objective   : 0.000000000000000e+00
BestBound   : 0.000000000000000e+00
StatusString: 
Time        : 0.001877 s

[SanPen]

There is something to investigate indeed.

I have just run spill_ortools.lp and spill_pulp.lp with highs using the commands:

highs spill_pulp.lp --solution_file spill_pulp.highs_solution
highs spill_ortools.lp --solution_file spill_ortools.highs_solution

Both solutions respect the bounds and both objective function values are 191.

To clarify our experiment:

We do have an interface that abstracts the use of ortools / pulp so that we can choose the framework and our formulation doesn't change. The steps have been:

1. Run the optimization with ortools -> We observe that the variable NodeSpillage_0_3 =-90 when this variable is defined in [0, 1e20]. The objective function value is -1488.26

2. We save the .lp and .mps files generated with ortools.

3. We change the OR framework to pulp and run the optimization -> we observe everything is ok.

4. We save the .lp and .mps files generated with pulp.

Now, all the .lp and .mps files sent by josep solve with highs to an objective function value of 191.

spill_models_and_solutions.zip

[SanPen]

@lperron any feedback on this?

[lperron]

So GLOP assumes that you use infinity when you mean it and not any numerical proxy.
You use 1e20 (which is generally infinity for other solvers). And GLOP treats it as 1e20, and performs presolve on it.

It happens with at least one constraint between -1e20..1e20. Meaning the constraint can be removed.