PAO.Pyomo model sets two variables x and y to be unique

Question

The background is I am trying to implement a bilevel optimization program using Python-based package for Adversarial Optimization (PAO), which is based on Pyomo. I am new to Pyomo, though I have used many solvers such as Gurobi and Cplex previously. To solve the bilevel optimization program, I need to use the pao.pyomo.MIBS solver. However, I found my Pyomo program can't even solve a naive unequal constraint z[0] != z[1] correctly. My code is as follows, which is trying to assign values to z[0], z[1] \in {0, 1} that satisfy constraint z[0] != z[1] :

import pyomo.environ as pyo
from pyomo.environ import *
from pao.pyomo import *

N = 2
M=1024

model = ConcreteModel()
model.z = Var(range(N), domain=Integers)
for i in range(N):
    model.z[i].setlb(0)
    model.z[i].setub(N-1)

model.delta = Var(range(N), range(N), domain=Binary)

model.unique_z = ConstraintList()
for i in range(N):
    for j in range(i+1, N):
        model.unique_z.add(model.z[i] - model.z[j] + M*(1-model.delta[i, j]) >= 1)
        model.unique_z.add(model.z[i] - model.z[j] - M*model.delta[i, j] <= -1)

model.obj = Objective(expr=sum(model.z[i] for i in range(N)), sense=maximize)
model.L = SubModel(fixed=model.z)

with Solver('pao.pyomo.MIBS') as solver:
    results = solver.solve(model)

print("obj =", model.obj.expr(), "z =", [model.z[i].value for i in range(N)])

So the basic constraint is like here, which expresses the not equal constraint z[0] != z[1] as a linear constraint with the help of binary variables delta. An obvious solution is either z[0] or z[1] is 1 and the other variable is 0.

However, I have printed the model's output and found the model sets both z[0] and z[1] to be 0.

>>> model.z.pprint()
z : Size=2, Index=z_index
    Key : Lower : Value : Upper : Fixed : Stale : Domain
      0 :     0 :     0 :     1 : False : False : Integers
      1 :     0 :     0 :     1 : False : False : Integers

>>> model.delta.pprint()
delta : Size=4, Index=delta_index
    Key    : Lower : Value : Upper : Fixed : Stale : Domain
    (0, 0) :     0 :  None :     1 : False :  True : Binary
    (0, 1) :     0 :     0 :     1 : False : False : Binary
    (1, 0) :     0 :  None :     1 : False :  True : Binary
    (1, 1) :     0 :  None :     1 : False :  True : Binary

>>> model.unique_z.pprint()
unique_z : Size=2, Index=unique_z_index, Active=True
    Key : Lower : Body                                : Upper : Active
      1 :   1.0 : z[0] - z[1] + 1024*(1 - delta[0,1]) :  +Inf :   True
      2 :  -Inf :       z[0] - z[1] - 1024*delta[0,1] :  -1.0 :   True

So when z[0], z[1], and delta(0,1) are 0, doesn't this violate the constraint z[0] - z[1] - 1024*delta[0,1] <= -1? I don't understand why this can be feasible.

I'm a beginner with Pyomo. If there are any misconceptions in the model I've built or in the Pyomo outputs, I would greatly appreciate any clarification or suggestions. Thank you.

---

Edit: Thanks to @AirSquid, it seems the problem is with my 'pao.pyomo.MIBS' solver. Could anyone give some suggestions on what's wrong with my current implementation of 'pao.pyomo.MIBS' solver? I have to use this solver for my bilevel optimization program.

Answer 1

I suspect there is something wrong with your solving code and you are looking at an un-solved or erroneous solve of some kind. I'm not familiar with the code sequence you are using there. I usually:

• get a solver instance from the SolverFactory
• catch the solve result
• inspect the solve result. <-- must do to ensure optimal or result is junk.

Your math is correct. I tweaked your model a tiny bit and used a different solver, but it's working. You are double importing. You should be able either use your first or second import line. I prefer the pyo prefix, but either way. I commented out the pao as it isn't used in my example.

Note that in pyomo, pprint shows the construction w/ variable names, and display can be used to fill in (or show) the solved values

Code:

# import pyomo.environ as pyo
from pyomo.environ import *
# from pao.pyomo import *


N = 6
M=1024

model = ConcreteModel()
model.z = Var(range(N), domain=Integers)
for i in range(N):
    model.z[i].setlb(0)
    model.z[i].setub(N-1)

model.delta = Var(range(N), range(N), domain=Binary)

model.unique_z = ConstraintList()
for i in range(N):
    for j in range(i+1, N):
        model.unique_z.add(model.z[i] - model.z[j] + M*(1-model.delta[i, j]) >= 1)
        model.unique_z.add(model.z[i] - model.z[j] - M*model.delta[i, j] <= -1)

model.obj = Objective(expr=sum(model.z[i] for i in range(N)), sense=maximize)
# model.L = SubModel(fixed=model.z)

solver = SolverFactory('cbc')
results = solver.solve(model)

# DO THIS.  Check for "OPTIMAL"
print(results)

print("obj =", model.obj.expr(), "z =", [model.z[i].value for i in range(N)])

# easy way to inspect the output of a variable (or the whole model)
model.z.display()
# model.display()

Output:

Problem: 
- Name: unknown
  Lower bound: 15.0
  Upper bound: 15.0
  Number of objectives: 1
  Number of constraints: 30
  Number of variables: 21
  Number of binary variables: 15
  Number of integer variables: 21
  Number of nonzeros: 6
  Sense: maximize
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.27
  Wallclock time: 0.46
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 270
      Number of created subproblems: 270
    Black box: 
      Number of iterations: 4106
  Error rc: 0
  Time: 0.47368288040161133
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

obj = 15.0 z = [2.0, 4.0, 0.0, 1.0, 3.0, 5.0]
z : Size=6, Index=z_index
    Key : Lower : Value : Upper : Fixed : Stale : Domain
      0 :     0 :   2.0 :     5 : False : False : Integers
      1 :     0 :   4.0 :     5 : False : False : Integers
      2 :     0 :   0.0 :     5 : False : False : Integers
      3 :     0 :   1.0 :     5 : False : False : Integers
      4 :     0 :   3.0 :     5 : False : False : Integers
      5 :     0 :   5.0 :     5 : False : False : Integers