Pyomo | Couenne solver | restrict indexed variable domain to two integer values

Question

I'm a newbie to Pyomo and struggling with understanding the intuition behind Pyomo's syntax and how it builds its models. This may be why I cannot figure out how to define and solve a 'binary' problem where N variables must take only ±1 values using Pyomo and Couenne solver.

First I tried the Integers domain with bounds=(-1, 1) and tried to add a strict inequality:

import pyomo.environ as pyo
import numpy as np
N = 5
w = np.ones(range(N))
pyoModel = pyo.ConcreteModel('binary model')
pyoModel.o = pyo.Var(range(N), bounds=(-1, 1), within=pyo.Integers, initialize=1)
pyoModel.binaryConstraintP = pyo.Constraint(range(N), rule=strictlyPositive)
pyoModel.binaryConstraintN = pyo.Constraint(range(N), rule=strictlyNegative)
pyoModel.objective = pyo.Objective(expr=pyo.sum_product(pyoModel.o, w, index=range(N)), sense=pyo.maximize)
def strictlyPositive(model, i):
    return model.o[i] > 0
def strictlyNegative(model, i):
    return model.o[i] < 0

This ends up with:

ValueError: Constraint 'binaryConstraintP[0]' encountered a strict inequality expression ('>' or '<'). All constraints must be formulated using using '<=', '>=', or '=='.

Okay, no strict inequalities allowed (don't know why!), I tried to switch to the Binary domain and do a workaround by manipulating the variable in the objective so it lies with in {-1, 1} - i.e., if o ∈ {0, 1} then 2 x o - 1 ∈ {-1, 1}:

import pyomo.environ as pyo
import numpy as np
N = 5
w = np.ones(range(N))
pyoModel = pyo.ConcreteModel('binary model')
pyoModel.o = pyo.Var(range(N), within=pyo.Binary, initialize=1)
pyoModel.objective = pyo.Objective(expr=pyo.sum_product(2 * pyoModel.o - 1, w, index=range(N)), sense=pyo.maximize)

Got:

TypeError: unsupported operand type(s) for *: 'int' and 'IndexedVar'

So I used an array of twos and ones instead of 2 and 1 but got another error about shape broadcasting. I'm sure I'm missing something here because it should be easy to construct a linear equation in the objective right?

I also tried to change the domain into a user-defined one:

...
pyoModel.domain = pyo.Set(initialize=[-1, 1])
...
pyoModel.o = pyo.Var(range(N), domain=pyoModel.domain, initialize=1)
...
with SolverFactory('couenne') as opt:
    results = opt.solve(pyoModel, load_solutions=False)
    ...

and ended up with Couenne error:

TypeError: Invalid domain type for variable with name '%s'. Variable is not continuous, integer, or binary.

I also thought of using SOSs but it was even harder to understand how they work!

Again, I must be missing something in each approach. Any help would be appriciated.

Side note: I simplified the original code as possible to make it easier to read.

Answer 1

The reason why strict inequalities are not allowed is that it might not be supported by the underlying solution theory. I believe this applies to MIP. Did you try setting the limits to very small values within the constraints you defined? E.g.

def strictlyPositive(model, i):
    return model.o[i] >= 0.000000000000001
def strictlyNegative(model, i):
    return model.o[i] <= -0.000000000000001

emulating strict inequalities this way is not ideal as it might cause numerical issues, but it's worth a try.

Answer 2

Your first attempt fails due to use of strict inequalities, which is a no-no. There is theory behind this as solvers of these types of problems work on the "convex hull" of the problem space. For more info, pick up a text on linear programming--it is beyond the scope of a stack overflow answer.

Your second attempt is on the right track. It is completely appropriate to let x be a binary variable and use the linear conversion 2x-1 if you are looking for math equivalent of ±1. Your attempt there is failing because you are declaring your x variable to be an indexed variable by providing an index in the Var() construct, yet you are not using the index in the Objective.

Here is an example that employs an indexed variable. If you only have a singleton, then just remove the indexing set references.

from pyomo.environ import *

some_constants = {  0: 100,
                    1: 200,
                    2: -50,
                    3: 300,
                    4: 50}

m = ConcreteModel('plus & minus ones project')

m.S = Set(initialize=range(5))
m.c = Param(m.S, initialize=some_constants)
m.x = Var(m.S, within=Binary)    # creates {m.x[0], m.x[1], ... , m.x[4]}

# try to maximize the sum of x*c
m.obj = Objective(expr=sum((2*m.x[s] - 1)*m.c[s] for s in m.S), sense=maximize)

# some constraint to limit the number of "+1 picks" to 2...easy with binary vars.
m.C1 = Constraint(expr=sum(m.x[s] for s in m.S) <= 2)

m.pprint()

Yields:

1 Set Declarations
    S : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    5 : {0, 1, 2, 3, 4}

1 Param Declarations
    c : Size=5, Index=S, Domain=Any, Default=None, Mutable=False
        Key : Value
          0 :   100
          1 :   200
          2 :   -50
          3 :   300
          4 :    50

1 Var Declarations
    x : Size=5, Index=S
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :     1 : False :  True : Binary
          1 :     0 :  None :     1 : False :  True : Binary
          2 :     0 :  None :     1 : False :  True : Binary
          3 :     0 :  None :     1 : False :  True : Binary
          4 :     0 :  None :     1 : False :  True : Binary

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : (2*x[0] - 1)*100 + (2*x[1] - 1)*200 + (2*x[2] - 1)*-50 + (2*x[3] - 1)*300 + (2*x[4] - 1)*50

1 Constraint Declarations
    C1 : Size=1, Index=None, Active=True
        Key  : Lower : Body                             : Upper : Active
        None :  -Inf : x[0] + x[1] + x[2] + x[3] + x[4] :   2.0 :   True

5 Declarations: S c x obj C1