OR-Tools：获取每一个最优解

Posted 2023-02-16

【中文标题】OR-Tools：获取每一个最优解

【英文标题】：OR-Tools: Obtaining every optimal solution

【发布时间】：2021-02-20 18:33:34

【问题描述】：

我正在使用 OR-Tools 通过 SCIP 解决 MIP。 OR-Tools 为我的问题的连续和整数（二进制）变量返回最佳值。

然后，当我将此 MIP 的二进制变量固定为 OR-Tools（对于 MIP）返回的最佳值并使用 GLOP 求解相应的 LP 时，OR-Tools 为连续变量的最佳值返回新值.

我的理解是初始问题没有唯一的解决方案（就变量的最优值而言）。

所以我的问题是：如何让 OR-tools 返回每个最佳解决方案，而不仅仅是一个？

请在下面找到代码：

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())

输出是

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

【问题讨论】：

请将模型发给我们。 Mps很好，代码也很好。您使用的是哪个求解器？

MPS 是什么意思？我将清理模型（它是更大程序的一部分）并发布它。我已经用求解器更新了我的原始帖子。

@LaurentPerron 我已按要求添加了代码。谢谢！

【参考方案1】：

在第一个Solve()之后调用NextSolution()可以获得所有最优解

【讨论】：

感谢您的回答。我在上面的代码中打印了 MIP 的解决方案的值后添加了solver_mip.NextSolution()。原来该值为False，与LP的解相矛盾。

PS：我通过这篇帖子 ***.com/questions/63291512/… 了解到，实际上 SCIP 不可能，而只有 Gurobi 才能实现。还是这样吗？

不，它是在 8.1 中使用 SCIP 和 Gurobi 实现的。

对不起，我不明白你的回答。当我使用“SCIP”作为求解器时，NextSolution() 返回 False，但是当我使用“GUROBI_MIXED_INTEGER_PROGRAMMING”时返回 True（但是 Gurobi 返回更高的目标值）。