OR-Tools: Obtaining every optimal solution

Question

I am using OR-Tools to solve a MIP with SCIP. OR-Tools returns optimal values for continuous and integer (binary) variables of my problem.

When I then fix the binary variables of this MIP to the optimal values returned by OR-Tools (for the MIP) and I solve the corresponding LP with GLOP, OR-Tools returns new values for the optimal values of the continuous variables.

My understanding is that the initial problem does not have a unique solution (in terms of optimal values of variables).

So my question is: How can I make OR-tools return every optimal solution, and not just one?

Please find the code bellow:

from ortools.linear_solver import pywraplp

#Fixed parameters
K = 4
L = 3
J = {}
J[0,0] = 3
J[0,1] = 2
J[0,2] = 1
J[1,0] = 3
J[1,1] = 1
J[1,2] = 2
J[2,0] = 2
J[2,1] = 2
J[2,2] = 2
J[3,0] = 1
J[3,1] = 1
J[3,2] = 1

S_up = {}
S_lw = {}
U = {}

S_up[0,0,0] = 20
S_up[0,1,0] = 40
S_up[0,2,0] = 60
S_up[0,0,1] = 15
S_up[0,1,1] = 40
S_up[0,0,2] = 50

S_lw[0,0,0] = 0
S_lw[0,1,0] = 21
S_lw[0,2,0] = 41
S_lw[0,0,1] = 5
S_lw[0,1,1] = 16
S_lw[0,0,2] = 10

U[0,0,0] = 5
U[0,1,0] = 7
U[0,2,0] = 8
U[0,0,1] = 2
U[0,1,1] = 5
U[0,0,2] = 6

S_up[1,0,0] = 25
S_up[1,1,0] = 35
S_up[1,2,0] = 50
S_up[1,0,1] = 30
S_up[1,0,2] = 30
S_up[1,1,2] = 30

S_lw[1,0,0] = 5
S_lw[1,1,0] = 26
S_lw[1,2,0] = 36
S_lw[1,0,1] = 10
S_lw[1,0,2] = 5
S_lw[1,1,2] = 31

U[1,0,0] = 6
U[1,1,0] = 8
U[1,2,0] = 9
U[1,0,1] = 3
U[1,0,2] = 5
U[1,1,2] = 7

S_up[2,0,0] = 40
S_up[2,1,0] = 60
S_up[2,0,1] = 60
S_up[2,1,1] = 80
S_up[2,0,2] = 40
S_up[2,1,2] = 60

S_lw[2,0,0] = 5
S_lw[2,1,0] = 41
S_lw[2,0,1] = 5
S_lw[2,1,1] = 61
S_lw[2,0,2] = 5
S_lw[2,1,2] = 41

U[2,0,0] = 5
U[2,1,0] = 6
U[2,0,1] = 4
U[2,1,1] = 5
U[2,0,2] = 5
U[2,1,2] = 6

S_up[3,0,0] = 100
S_up[3,0,1] = 90
S_up[3,0,2] = 90

S_lw[3,0,0] = 0
S_lw[3,0,1] = 0
S_lw[3,0,2] = 0

U[3,0,0] = 5
U[3,0,1] = 4
U[3,0,2] = 5

D = [100,100,90]
P = [50,30,20]

Q_up = [0,1,1,1]
Q_lw = [0,0,0.1,0]

#Declare MIP solver
solver_mip = pywraplp.Solver.CreateSolver('SCIP')

#Define variables
infinity = solver_mip.infinity()
y = {}
for k in range(K):
    for l in range(L):
        for j in range(J[k,l]):
            y[k, j, l] = solver_mip.NumVar(0, infinity, '')

x = {}
for k in range(K):
    for l in range(L):
        for j in range(J[k,l]):
            x[k, j, l] = solver_mip.IntVar(0, 1, '')

print('Number of variables =', solver_mip.NumVariables())

#Define constraints
for k in range(K):
    for l in range(L):
        for j in range(J[k,l]):
            solver_mip.Add(y[k, j, l] <= x[k, j, l]*S_up[k, j, l])
            solver_mip.Add(x[k, j, l]*S_lw[k, j, l] <= y[k, j, l])

for k in range(K):
    for l in range(L):
        solver_mip.Add(sum([x[k, j, l] for j in range(J[k,l])]) <= 1)

for l in range(L):
    solver_mip.Add(sum([sum([y[k, j, l] for j in range(J[k,l])]) for k in range(K)]) == D[l]) 

for k in range(K):
    solver_mip.Add(sum([sum([y[k, j, l]*P[l] for j in range(J[k,l])]) for l in range(L)]) <= Q_up[k]*sum([D[l]*P[l] for l in range(L)]))
    solver_mip.Add(Q_lw[k]*sum([D[l]*P[l] for l in range(L)]) <= sum([sum([y[k, j, l]*P[l] for j in range(J[k,l])]) for l in range(L)]))

print('Number of constraints =', solver_mip.NumConstraints())

#Define objective
solver_mip.Minimize(sum([sum([sum([y[k,j,l]*U[k,j,l] for j in range(J[k,l])]) for k in range(K)]) for l in range(L)]))

#Call MIP solver
status = solver_mip.Solve()

#Display solution
if status == pywraplp.Solver.OPTIMAL:
    print('Solution of MIP:')
    print('Objective value =', solver_mip.Objective().Value())
    x_opt = {} #store optimal values of binary variable
    for k in range(K):
        for l in range(L):
            for j in range(J[k,l]):
                x_opt[k,j,l] = x[k,j,l].solution_value()
                if x[k,j,l].solution_value() == 1:
                    print('y[',k,',',j,',',l,']=',y[k,j,l].solution_value())
else:
    print('The problem does not have an optimal solution.')

print('\nAdvanced usage:')
print('Problem solved in %f milliseconds' % solver_mip.wall_time())
print('Problem solved in %d iterations' % solver_mip.iterations())
print('Problem solved in %d branch-and-bound nodes' % solver_mip.nodes())

#Primal problem with fixed binary variables to optimal value becomes a LP
#Declare LP solver
solver_lp = pywraplp.Solver.CreateSolver('GLOP')

##Quantity variable
y_fixed_binary = {}
for k in range(K):
    for l in range(L):
        for j in range(J[k,l]):
            y_fixed_binary[k, j, l] = solver_lp.NumVar(0, infinity, '')

#Define constraints
##Quantity should be in bounds defined by reinsurer
for k in range(K):
    for l in range(L):
        for j in range(J[k,l]):
            solver_lp.Add(y_fixed_binary[k, j, l] <= x_opt[k, j, l]*S_up[k, j, l])
            solver_lp.Add(x_opt[k, j, l]*S_lw[k, j, l] <= y_fixed_binary[k, j, l])

for l in range(L):
    solver_lp.Add(sum([sum([y_fixed_binary[k, j, l] for j in range(J[k,l])]) for k in range(K)]) == D[l]) 

for k in range(K):
    solver_lp.Add(sum([sum([y_fixed_binary[k, j, l]*P[l] for j in range(J[k,l])]) for l in range(L)]) <= Q_up[k]*sum([D[l]*P[l] for l in range(L)]))
    solver_lp.Add(Q_lw[k]*sum([D[l]*P[l] for l in range(L)]) <= sum([sum([y_fixed_binary[k ,j ,l]*P[l] for j in range(J[k,l])]) for l in range(L)]))

#Define objective
solver_lp.Minimize(sum([sum([sum([y_fixed_binary[k,j,l]*U[k,j,l] for j in range(J[k,l])]) for k in range(K)]) for l in range(L)]))

status = solver_lp.Solve()

if status == pywraplp.Solver.OPTIMAL:
    print('Solution of LP:')
    print('Objective value =', solver_mip.Objective().Value())
    for k in range(K):
        for l in range(L):
            for j in range(J[k,l]):
                if x_opt[k,j,l] == 1:
                    print('y[',k,',',j,',',l,']=',y_fixed_binary[k,j,l].solution_value())
else:
    print('The LP does not have an optimal solution.')

print('Advanced usage:')
print('Problem solved in %f milliseconds' % solver_lp.wall_time())

And the output is

Solution of MIP:
Objective value = 1280.0
y[ 1 , 0 , 1 ]= 30.0
y[ 2 , 0 , 0 ]= 13.999999999999998
y[ 2 , 0 , 1 ]= 5.0
y[ 2 , 0 , 2 ]= 5.0
y[ 3 , 0 , 0 ]= 86.0
y[ 3 , 0 , 1 ]= 55.0
y[ 3 , 0 , 2 ]= 85.0

Solution of LP:
Objective value = 1280.0
y[ 1 , 0 , 1 ]= 30.0
y[ 2 , 0 , 0 ]= 40.0
y[ 2 , 0 , 1 ]= 60.0
y[ 2 , 0 , 2 ]= 40.0
y[ 3 , 0 , 0 ]= 60.0
y[ 3 , 0 , 1 ]= 0.0
y[ 3 , 0 , 2 ]= 50.0

Please send us the model. Mps is fine, code too. Which solver are you using ? — Laurent Perron, Feb 20 '21 at 19:41
What does MPS mean? I'll clean the model (it is part of a bigger program) and post it. I've updated my original post with the solvers. — DeepNet, Feb 20 '21 at 19:57

score 2 · Answer 1 · answered Feb 22 '21 at 18:26

2

You can get all optimal solutions by calling NextSolution() after the first Solve()

answered Feb 22 '21 at 18:26

Laurent Perron

8,594
1
8
22

Thanks for your answer. I added ```solver_mip.NextSolution()``` after printing the values of the solutions of the MIP in the code above. It turns out that the value is ```False```, which is in contradiction with the solution of the LP. – DeepNet Feb 23 '21 at 13:56
PS: I understand via this post https://stackoverflow.com/questions/63291512/multiple-milp-solutions-in-ortools-python that it is in fact not possible with SCIP but only with Gurobi. Is it still the case? – DeepNet Feb 23 '21 at 15:10
No, it is implemented with SCIP and Gurobi in 8.1. – Laurent Perron Feb 23 '21 at 16:57
I am sorry but I don't understand your answer. When I use 'SCIP' as a solver, NextSolution() returns False but when I use 'GUROBI_MIXED_INTEGER_PROGRAMMING' is returns True (however Gurobi returns a higher objective value). – DeepNet Feb 23 '21 at 17:00

OR-Tools: Obtaining every optimal solution

1 Answers1