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