Gekko not solving Integer Programming Problem

Question

I'm trying to use the integer programming option to find a magic square. The algorithm finds a solution if I remove the requirement the entries are unique. I can add up to six requirements for entries not being equal but fails to find a solution when I add the seventh. As I know there is a solution to this 3x3 case, is there a way to change the solver settings, or reframe the constraints to get a solution?

from gekko import GEKKO
m=GEKKO(remote=False)
m.options.SOLVER=1

n=3 #size of magic square
S=15 #sum of each row, col, diag

M = m.Array(m.Var,(n,n),lb=1,ub=9,integer=True)

#constraint that no two entries are the same: (hard-coded to find where it fails)
m.Equation(abs(M[0,0]-M[0,1])>=1)
m.Equation(abs(M[0,0]-M[0,2])>=1)
m.Equation(abs(M[0,0]-M[1,0])>=1)
m.Equation(abs(M[0,0]-M[1,1])>=1)
m.Equation(abs(M[0,0]-M[1,2])>=1)
m.Equation(abs(M[0,0]-M[2,0])>=1)
m.Equation(abs(M[0,0]-M[2,1])>=1)
#m.Equation(abs(M[0,0]-M[2,2])>=1)
#m.Equation(abs(M[0,1]-M[0,2])>=1)
#m.Equation(abs(M[0,1]-M[1,0])>=1)
#m.Equation(abs(M[0,1]-M[1,1])>=1)
#m.Equation(abs(M[0,1]-M[1,2])>=1)
#m.Equation(abs(M[0,1]-M[2,0])>=1)
#m.Equation(abs(M[0,1]-M[2,1])>=1)
#m.Equation(abs(M[0,1]-M[2,2])>=1)
#m.Equation(abs(M[0,2]-M[1,0])>=1)
#m.Equation(abs(M[0,2]-M[1,1])>=1)
#m.Equation(abs(M[0,2]-M[1,2])>=1)
#m.Equation(abs(M[0,2]-M[2,0])>=1)
#m.Equation(abs(M[0,2]-M[2,1])>=1)
#m.Equation(abs(M[0,2]-M[2,2])>=1)
#m.Equation(abs(M[1,0]-M[1,1])>=1)
#m.Equation(abs(M[1,0]-M[1,2])>=1)
#m.Equation(abs(M[1,0]-M[2,0])>=1)
#m.Equation(abs(M[1,0]-M[2,1])>=1)
#m.Equation(abs(M[1,0]-M[2,2])>=1)
#m.Equation(abs(M[1,1]-M[1,2])>=1)
#m.Equation(abs(M[1,1]-M[2,0])>=1)
#m.Equation(abs(M[1,1]-M[2,1])>=1)
#m.Equation(abs(M[1,1]-M[2,2])>=1)
#m.Equation(abs(M[1,2]-M[2,0])>=1)
#m.Equation(abs(M[1,2]-M[2,1])>=1)
#m.Equation(abs(M[1,2]-M[2,2])>=1)
#m.Equation(abs(M[2,0]-M[2,1])>=1)
#m.Equation(abs(M[2,0]-M[2,2])>=1)
#m.Equation(abs(M[2,1]-M[2,2])>=1)
   

#Each col sums to S
for i in range(n):
    m.Equation(m.sum([M[i,j] for j in range(n)])==S)
        
#Each row sum to S        
for j in range(n):
    m.Equation(m.sum([M[i,j] for i in range(n)])==S)
               
#Diagonals sum to S
m.Equation(m.sum([M[i,i] for i in range(n)])==S)
m.Equation(m.sum([M[i,n-1-i] for i in range(n)])==S)

m.solve()

print(M)

I've tried rewriting the constraints as abs(M[0,0] - M[0,1])>0. I've also changed to squares. It seems the issue is I've hit an upper bound on the number of constraints--it doesn't mention that in the output, but wondering if that is the issue.

Answer 1

Reinderien is correct that it is often easier to translate to a LP problem, even for the MINLP solver (APOPT) in GEKKO. If you do want to use nonlinear models, use m.abs3() or m.abs2() instead of abs() that doesn't work well with gradient-based solvers. Here is an MILP version of the program in Gekko.

from gekko import GEKKO

n = 3  # size of the magic square
S = n * (n*n + 1) // 2

# Create a Gekko model
m = GEKKO(remote=False)

# Create matrix of decision variables
selectors = m.Array(m.Var,(n,n,n*n),lb=0,ub=1,integer=True)

# Sum of each row, column must equal S
for i in range(n):
    m.Equation(m.sum([selectors[i][j][k]*(k+1) for j in range(n) for k in range(n*n)]) == S)  # Rows
    m.Equation(m.sum([selectors[j][i][k]*(k+1) for j in range(n) for k in range(n*n)]) == S)  # Columns

m.Equation(m.sum([selectors[i][i][k]*(k+1) for i in range(n) for k in range(n*n)]) == S)  # Diagonal
m.Equation(m.sum([selectors[i][n-i-1][k]*(k+1) for i in range(n) for k in range(n*n)]) == S)  # Anti-diagonal

# Each number must be used exactly once
for k in range(n*n):
    m.Equation(sum(selectors[i][j][k] for i in range(n) for j in range(n)) == 1)
    
# Each number can only be used once for each entry
# Avoid solutions for 3x3 such as:
# [4,   2, 9]
# [7+3, 5, 0]
# [1,   8, 6]
for i in range(n):
    for j in range(n):
        m.Equation(sum(selectors[i][j][k] for k in range(n*n)) == 1)

# Solve the model
m.options.SOLVER=1
m.solve(disp=True)

# Extract and print the solution
solution = [[sum((k+1) * selectors[i][j][k].value[0] for k in range(n*n)) for j in range(n)] for i in range(n)]
for row in solution:
    print(row)

It produces the following solutions:

3x3 Magic Square

 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.01 NLPi:    2 Dpth:    0 Lvs:    3 Obj:  0.00E+00 Gap:       NaN
Iter:     2 I: -1 Tm:      0.00 NLPi:    1 Dpth:    1 Lvs:    2 Obj:  0.00E+00 Gap:       NaN
Iter:     3 I:  0 Tm:      0.01 NLPi:    2 Dpth:    1 Lvs:    4 Obj:  0.00E+00 Gap:       NaN
Iter:     4 I: -1 Tm:      0.00 NLPi:    1 Dpth:    2 Lvs:    3 Obj:  0.00E+00 Gap:       NaN
--Integer Solution:   0.00E+00 Lowest Leaf:   0.00E+00 Gap:   0.00E+00
Iter:     5 I:  0 Tm:      0.01 NLPi:    2 Dpth:    2 Lvs:    3 Obj:  0.00E+00 Gap:  0.00E+00
 Successful solution

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.05400000000000005 sec
 Objective      :  0.
 Successful solution
 ---------------------------------------------------

[8.0, 1.0, 6.0]
[3.0, 5.0, 7.0]
[4.0, 9.0, 2.0]

4x4 Magic Square

Iter:   179 I: 0 Tm: 0.01 NLPi: 2 Dpth: 11 Lvs: 8 Obj: 0.00E+00 Gap: 0.00E+00
 Successful solution

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  29.3463 sec
 Objective      :  0.
 Successful solution
 ---------------------------------------------------

[11.0, 13.0, 6.0, 4.0]
[8.0, 2.0, 9.0, 15.0]
[1.0, 7.0, 16.0, 10.0]
[14.0, 12.0, 3.0, 5.0]

The 3x3 solution is very fast, but the 4x4 (and likely larger problems) will take much longer to run with the APOPT solver. A dedicated MILP solver may work better for larger problems with number of binary variables:

• 3x3 Magic Square = 3^4 = 81 binary variables,
• 4x4 Magic Square = 4^4 = 256 binary variables
• 5x5 Magic Square = 5^4 = 625 binary variables
• 6x6 Magic Square = 6^4 = 1296 binary variables

For the 4x4 problem, I tried without using m.sum() and used the Python sum() function with faster results (29 sec vs 113 sec).

from gekko import GEKKO

n = 3  # size of the magic square
S = n * (n*n + 1) // 2

# Create a Gekko model
m = GEKKO(remote=False)

# Create matrix of decision variables
selectors = m.Array(m.Var,(n,n,n*n),lb=0,ub=1,integer=True)

# Sum of each row, column must equal S
for i in range(n):
    m.Equation(sum(selectors[i][j][k]*(k+1) for j in range(n) for k in range(n*n)) == S)  # Rows
    m.Equation(sum(selectors[j][i][k]*(k+1) for j in range(n) for k in range(n*n)) == S)  # Columns

m.Equation(sum(selectors[i][i][k]*(k+1) for i in range(n) for k in range(n*n)) == S)  # Diagonal
m.Equation(sum(selectors[i][n-i-1][k]*(k+1) for i in range(n) for k in range(n*n)) == S)  # Anti-diagonal

# Each number must be used exactly once
for k in range(n*n):
    m.Equation(sum(selectors[i][j][k] for i in range(n) for j in range(n)) == 1)

for i in range(n):
    for j in range(n):
        m.Equation(sum(selectors[i][j][k] for k in range(n*n)) == 1)

# Solve the model
m.options.SOLVER=1
m.solve(disp=True)

# Extract and print the solution
solution = [[sum((k+1) * selectors[i][j][k].value[0] for k in range(n*n)) for j in range(n)] for i in range(n)]
for row in solution:
    print(row)

The m.sum() is used for larger problems if a model equation exceeds 15,000 characters.

Answer 2

If I understand how Gekko works correctly, it's only somewhat of an accident that abs() works at all. In a regular linear program, abs() is never supported directly and needs to be constructed on inference by auxiliary variables. However, since the auxiliary absolute variable would be the subject of a minimum bound, I cannot think of a way to construct an LP where this would work, since the auxiliary variable would just take some arbitrary value above its constraint but not reflect the actual absolute value. Long story short, under the hood your problem as currently constructed is not linear.

The most straightforward way to implement this is through a traditional integer linear program. I demonstrate this with a less-clever PuLP which does not have Gekko's non-linear support. This works via binary assignment variables, one for each cell-value pair.

import pulp


n = 3   # size of magic square
target_total = n*(n*n + 1)//2

values = range(1, 1 + n*n)

selectors = pulp.LpVariable.matrix(
    name='sel',
    indices=(range(n), range(n), values),
    cat=pulp.LpBinary,
)

entries = [
    [
        pulp.lpDot(elm, values)
        for elm in row
    ]
    for row in selectors
]

# No optimization sense - there is no objective
prob = pulp.LpProblem(name='magic_square')

# Each element must exclusively select one value
for i, row in enumerate(selectors):
    for j, elm in enumerate(row):
        prob.addConstraint(
            name=f'sel_{i}_{j}_k',
            constraint=pulp.lpSum(elm) == 1,
        )

# Each value must be unique
for k, value in enumerate(values):
    prob.addConstraint(
        name=f'sel_ij{value}',
        constraint=pulp.lpSum(
            elm[k]
            for row in selectors
            for elm in row
        ) == 1,
    )

# Segment sum constraints
for i, row in enumerate(entries):
    prob.addConstraint(
        name=f'row_{i}j',
        constraint=pulp.lpSum(row) == target_total,
    )
for j in range(n):
    prob.addConstraint(
        name=f'col_i{j}',
        constraint=pulp.lpSum(
            row[j] for row in entries
        ) == target_total,
    )
prob.addConstraint(
    name='diagonal_ij',
    constraint=pulp.lpSum(
        row[ij]
        for ij, row in enumerate(entries)
    ) == target_total,
)
prob.addConstraint(
    name='antidiagonal_ij',
    constraint=pulp.lpSum(
        row[-i]
        for i, row in enumerate(entries, start=1)
    ) == target_total,
)

print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal

for row in entries:
    print([
        int(elm.value())
        for elm in row
    ])

magic_square:
MINIMIZE
None
SUBJECT TO
sel_0_0_k: sel_0_0_1 + sel_0_0_2 + sel_0_0_3 + sel_0_0_4 + sel_0_0_5
 + sel_0_0_6 + sel_0_0_7 + sel_0_0_8 + sel_0_0_9 = 1

sel_0_1_k: sel_0_1_1 + sel_0_1_2 + sel_0_1_3 + sel_0_1_4 + sel_0_1_5
 + sel_0_1_6 + sel_0_1_7 + sel_0_1_8 + sel_0_1_9 = 1

sel_0_2_k: sel_0_2_1 + sel_0_2_2 + sel_0_2_3 + sel_0_2_4 + sel_0_2_5
 + sel_0_2_6 + sel_0_2_7 + sel_0_2_8 + sel_0_2_9 = 1

sel_1_0_k: sel_1_0_1 + sel_1_0_2 + sel_1_0_3 + sel_1_0_4 + sel_1_0_5
 + sel_1_0_6 + sel_1_0_7 + sel_1_0_8 + sel_1_0_9 = 1

sel_1_1_k: sel_1_1_1 + sel_1_1_2 + sel_1_1_3 + sel_1_1_4 + sel_1_1_5
 + sel_1_1_6 + sel_1_1_7 + sel_1_1_8 + sel_1_1_9 = 1

sel_1_2_k: sel_1_2_1 + sel_1_2_2 + sel_1_2_3 + sel_1_2_4 + sel_1_2_5
 + sel_1_2_6 + sel_1_2_7 + sel_1_2_8 + sel_1_2_9 = 1

sel_2_0_k: sel_2_0_1 + sel_2_0_2 + sel_2_0_3 + sel_2_0_4 + sel_2_0_5
 + sel_2_0_6 + sel_2_0_7 + sel_2_0_8 + sel_2_0_9 = 1

sel_2_1_k: sel_2_1_1 + sel_2_1_2 + sel_2_1_3 + sel_2_1_4 + sel_2_1_5
 + sel_2_1_6 + sel_2_1_7 + sel_2_1_8 + sel_2_1_9 = 1

sel_2_2_k: sel_2_2_1 + sel_2_2_2 + sel_2_2_3 + sel_2_2_4 + sel_2_2_5
 + sel_2_2_6 + sel_2_2_7 + sel_2_2_8 + sel_2_2_9 = 1

sel_ij1: sel_0_0_1 + sel_0_1_1 + sel_0_2_1 + sel_1_0_1 + sel_1_1_1 + sel_1_2_1
 + sel_2_0_1 + sel_2_1_1 + sel_2_2_1 = 1

sel_ij2: sel_0_0_2 + sel_0_1_2 + sel_0_2_2 + sel_1_0_2 + sel_1_1_2 + sel_1_2_2
 + sel_2_0_2 + sel_2_1_2 + sel_2_2_2 = 1

sel_ij3: sel_0_0_3 + sel_0_1_3 + sel_0_2_3 + sel_1_0_3 + sel_1_1_3 + sel_1_2_3
 + sel_2_0_3 + sel_2_1_3 + sel_2_2_3 = 1

sel_ij4: sel_0_0_4 + sel_0_1_4 + sel_0_2_4 + sel_1_0_4 + sel_1_1_4 + sel_1_2_4
 + sel_2_0_4 + sel_2_1_4 + sel_2_2_4 = 1

sel_ij5: sel_0_0_5 + sel_0_1_5 + sel_0_2_5 + sel_1_0_5 + sel_1_1_5 + sel_1_2_5
 + sel_2_0_5 + sel_2_1_5 + sel_2_2_5 = 1

sel_ij6: sel_0_0_6 + sel_0_1_6 + sel_0_2_6 + sel_1_0_6 + sel_1_1_6 + sel_1_2_6
 + sel_2_0_6 + sel_2_1_6 + sel_2_2_6 = 1

sel_ij7: sel_0_0_7 + sel_0_1_7 + sel_0_2_7 + sel_1_0_7 + sel_1_1_7 + sel_1_2_7
 + sel_2_0_7 + sel_2_1_7 + sel_2_2_7 = 1

sel_ij8: sel_0_0_8 + sel_0_1_8 + sel_0_2_8 + sel_1_0_8 + sel_1_1_8 + sel_1_2_8
 + sel_2_0_8 + sel_2_1_8 + sel_2_2_8 = 1

sel_ij9: sel_0_0_9 + sel_0_1_9 + sel_0_2_9 + sel_1_0_9 + sel_1_1_9 + sel_1_2_9
 + sel_2_0_9 + sel_2_1_9 + sel_2_2_9 = 1

row_0j: sel_0_0_1 + 2 sel_0_0_2 + 3 sel_0_0_3 + 4 sel_0_0_4 + 5 sel_0_0_5
 + 6 sel_0_0_6 + 7 sel_0_0_7 + 8 sel_0_0_8 + 9 sel_0_0_9 + sel_0_1_1
 + 2 sel_0_1_2 + 3 sel_0_1_3 + 4 sel_0_1_4 + 5 sel_0_1_5 + 6 sel_0_1_6
 + 7 sel_0_1_7 + 8 sel_0_1_8 + 9 sel_0_1_9 + sel_0_2_1 + 2 sel_0_2_2
 + 3 sel_0_2_3 + 4 sel_0_2_4 + 5 sel_0_2_5 + 6 sel_0_2_6 + 7 sel_0_2_7
 + 8 sel_0_2_8 + 9 sel_0_2_9 = 15

row_1j: sel_1_0_1 + 2 sel_1_0_2 + 3 sel_1_0_3 + 4 sel_1_0_4 + 5 sel_1_0_5
 + 6 sel_1_0_6 + 7 sel_1_0_7 + 8 sel_1_0_8 + 9 sel_1_0_9 + sel_1_1_1
 + 2 sel_1_1_2 + 3 sel_1_1_3 + 4 sel_1_1_4 + 5 sel_1_1_5 + 6 sel_1_1_6
 + 7 sel_1_1_7 + 8 sel_1_1_8 + 9 sel_1_1_9 + sel_1_2_1 + 2 sel_1_2_2
 + 3 sel_1_2_3 + 4 sel_1_2_4 + 5 sel_1_2_5 + 6 sel_1_2_6 + 7 sel_1_2_7
 + 8 sel_1_2_8 + 9 sel_1_2_9 = 15

row_2j: sel_2_0_1 + 2 sel_2_0_2 + 3 sel_2_0_3 + 4 sel_2_0_4 + 5 sel_2_0_5
 + 6 sel_2_0_6 + 7 sel_2_0_7 + 8 sel_2_0_8 + 9 sel_2_0_9 + sel_2_1_1
 + 2 sel_2_1_2 + 3 sel_2_1_3 + 4 sel_2_1_4 + 5 sel_2_1_5 + 6 sel_2_1_6
 + 7 sel_2_1_7 + 8 sel_2_1_8 + 9 sel_2_1_9 + sel_2_2_1 + 2 sel_2_2_2
 + 3 sel_2_2_3 + 4 sel_2_2_4 + 5 sel_2_2_5 + 6 sel_2_2_6 + 7 sel_2_2_7
 + 8 sel_2_2_8 + 9 sel_2_2_9 = 15

col_i0: sel_0_0_1 + 2 sel_0_0_2 + 3 sel_0_0_3 + 4 sel_0_0_4 + 5 sel_0_0_5
 + 6 sel_0_0_6 + 7 sel_0_0_7 + 8 sel_0_0_8 + 9 sel_0_0_9 + sel_1_0_1
 + 2 sel_1_0_2 + 3 sel_1_0_3 + 4 sel_1_0_4 + 5 sel_1_0_5 + 6 sel_1_0_6
 + 7 sel_1_0_7 + 8 sel_1_0_8 + 9 sel_1_0_9 + sel_2_0_1 + 2 sel_2_0_2
 + 3 sel_2_0_3 + 4 sel_2_0_4 + 5 sel_2_0_5 + 6 sel_2_0_6 + 7 sel_2_0_7
 + 8 sel_2_0_8 + 9 sel_2_0_9 = 15

col_i1: sel_0_1_1 + 2 sel_0_1_2 + 3 sel_0_1_3 + 4 sel_0_1_4 + 5 sel_0_1_5
 + 6 sel_0_1_6 + 7 sel_0_1_7 + 8 sel_0_1_8 + 9 sel_0_1_9 + sel_1_1_1
 + 2 sel_1_1_2 + 3 sel_1_1_3 + 4 sel_1_1_4 + 5 sel_1_1_5 + 6 sel_1_1_6
 + 7 sel_1_1_7 + 8 sel_1_1_8 + 9 sel_1_1_9 + sel_2_1_1 + 2 sel_2_1_2
 + 3 sel_2_1_3 + 4 sel_2_1_4 + 5 sel_2_1_5 + 6 sel_2_1_6 + 7 sel_2_1_7
 + 8 sel_2_1_8 + 9 sel_2_1_9 = 15

col_i2: sel_0_2_1 + 2 sel_0_2_2 + 3 sel_0_2_3 + 4 sel_0_2_4 + 5 sel_0_2_5
 + 6 sel_0_2_6 + 7 sel_0_2_7 + 8 sel_0_2_8 + 9 sel_0_2_9 + sel_1_2_1
 + 2 sel_1_2_2 + 3 sel_1_2_3 + 4 sel_1_2_4 + 5 sel_1_2_5 + 6 sel_1_2_6
 + 7 sel_1_2_7 + 8 sel_1_2_8 + 9 sel_1_2_9 + sel_2_2_1 + 2 sel_2_2_2
 + 3 sel_2_2_3 + 4 sel_2_2_4 + 5 sel_2_2_5 + 6 sel_2_2_6 + 7 sel_2_2_7
 + 8 sel_2_2_8 + 9 sel_2_2_9 = 15

diagonal_ij: sel_0_0_1 + 2 sel_0_0_2 + 3 sel_0_0_3 + 4 sel_0_0_4 + 5 sel_0_0_5
 + 6 sel_0_0_6 + 7 sel_0_0_7 + 8 sel_0_0_8 + 9 sel_0_0_9 + sel_1_1_1
 + 2 sel_1_1_2 + 3 sel_1_1_3 + 4 sel_1_1_4 + 5 sel_1_1_5 + 6 sel_1_1_6
 + 7 sel_1_1_7 + 8 sel_1_1_8 + 9 sel_1_1_9 + sel_2_2_1 + 2 sel_2_2_2
 + 3 sel_2_2_3 + 4 sel_2_2_4 + 5 sel_2_2_5 + 6 sel_2_2_6 + 7 sel_2_2_7
 + 8 sel_2_2_8 + 9 sel_2_2_9 = 15

antidiagonal_ij: sel_0_2_1 + 2 sel_0_2_2 + 3 sel_0_2_3 + 4 sel_0_2_4
 + 5 sel_0_2_5 + 6 sel_0_2_6 + 7 sel_0_2_7 + 8 sel_0_2_8 + 9 sel_0_2_9
 + sel_1_1_1 + 2 sel_1_1_2 + 3 sel_1_1_3 + 4 sel_1_1_4 + 5 sel_1_1_5
 + 6 sel_1_1_6 + 7 sel_1_1_7 + 8 sel_1_1_8 + 9 sel_1_1_9 + sel_2_0_1
 + 2 sel_2_0_2 + 3 sel_2_0_3 + 4 sel_2_0_4 + 5 sel_2_0_5 + 6 sel_2_0_6
 + 7 sel_2_0_7 + 8 sel_2_0_8 + 9 sel_2_0_9 = 15

VARIABLES
0 <= sel_0_0_1 <= 1 Integer
0 <= sel_0_0_2 <= 1 Integer
0 <= sel_0_0_3 <= 1 Integer
0 <= sel_0_0_4 <= 1 Integer
0 <= sel_0_0_5 <= 1 Integer
0 <= sel_0_0_6 <= 1 Integer
0 <= sel_0_0_7 <= 1 Integer
0 <= sel_0_0_8 <= 1 Integer
0 <= sel_0_0_9 <= 1 Integer
0 <= sel_0_1_1 <= 1 Integer
0 <= sel_0_1_2 <= 1 Integer
0 <= sel_0_1_3 <= 1 Integer
0 <= sel_0_1_4 <= 1 Integer
0 <= sel_0_1_5 <= 1 Integer
0 <= sel_0_1_6 <= 1 Integer
0 <= sel_0_1_7 <= 1 Integer
0 <= sel_0_1_8 <= 1 Integer
0 <= sel_0_1_9 <= 1 Integer
0 <= sel_0_2_1 <= 1 Integer
0 <= sel_0_2_2 <= 1 Integer
0 <= sel_0_2_3 <= 1 Integer
0 <= sel_0_2_4 <= 1 Integer
0 <= sel_0_2_5 <= 1 Integer
0 <= sel_0_2_6 <= 1 Integer
0 <= sel_0_2_7 <= 1 Integer
0 <= sel_0_2_8 <= 1 Integer
0 <= sel_0_2_9 <= 1 Integer
0 <= sel_1_0_1 <= 1 Integer
0 <= sel_1_0_2 <= 1 Integer
0 <= sel_1_0_3 <= 1 Integer
0 <= sel_1_0_4 <= 1 Integer
0 <= sel_1_0_5 <= 1 Integer
0 <= sel_1_0_6 <= 1 Integer
0 <= sel_1_0_7 <= 1 Integer
0 <= sel_1_0_8 <= 1 Integer
0 <= sel_1_0_9 <= 1 Integer
0 <= sel_1_1_1 <= 1 Integer
0 <= sel_1_1_2 <= 1 Integer
0 <= sel_1_1_3 <= 1 Integer
0 <= sel_1_1_4 <= 1 Integer
0 <= sel_1_1_5 <= 1 Integer
0 <= sel_1_1_6 <= 1 Integer
0 <= sel_1_1_7 <= 1 Integer
0 <= sel_1_1_8 <= 1 Integer
0 <= sel_1_1_9 <= 1 Integer
0 <= sel_1_2_1 <= 1 Integer
0 <= sel_1_2_2 <= 1 Integer
0 <= sel_1_2_3 <= 1 Integer
0 <= sel_1_2_4 <= 1 Integer
0 <= sel_1_2_5 <= 1 Integer
0 <= sel_1_2_6 <= 1 Integer
0 <= sel_1_2_7 <= 1 Integer
0 <= sel_1_2_8 <= 1 Integer
0 <= sel_1_2_9 <= 1 Integer
0 <= sel_2_0_1 <= 1 Integer
0 <= sel_2_0_2 <= 1 Integer
0 <= sel_2_0_3 <= 1 Integer
0 <= sel_2_0_4 <= 1 Integer
0 <= sel_2_0_5 <= 1 Integer
0 <= sel_2_0_6 <= 1 Integer
0 <= sel_2_0_7 <= 1 Integer
0 <= sel_2_0_8 <= 1 Integer
0 <= sel_2_0_9 <= 1 Integer
0 <= sel_2_1_1 <= 1 Integer
0 <= sel_2_1_2 <= 1 Integer
0 <= sel_2_1_3 <= 1 Integer
0 <= sel_2_1_4 <= 1 Integer
0 <= sel_2_1_5 <= 1 Integer
0 <= sel_2_1_6 <= 1 Integer
0 <= sel_2_1_7 <= 1 Integer
0 <= sel_2_1_8 <= 1 Integer
0 <= sel_2_1_9 <= 1 Integer
0 <= sel_2_2_1 <= 1 Integer
0 <= sel_2_2_2 <= 1 Integer
0 <= sel_2_2_3 <= 1 Integer
0 <= sel_2_2_4 <= 1 Integer
0 <= sel_2_2_5 <= 1 Integer
0 <= sel_2_2_6 <= 1 Integer
0 <= sel_2_2_7 <= 1 Integer
0 <= sel_2_2_8 <= 1 Integer
0 <= sel_2_2_9 <= 1 Integer

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 
...

Result - Optimal solution found

Objective value:                0.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.03
Time (Wallclock seconds):       0.03

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.03   (Wallclock seconds):       0.03

[4, 3, 8]
[9, 5, 1]
[2, 7, 6]