pyomo example

pyomo_test.py

import pyomo.environ as pyo
from pyomo.opt import SolverFactory

S = [
  (2,2), (5,2), (8,2),
  (2,5), (5,5), (8,5),
  (2,8), (5,8), (8,8)
]

V = [
  (-1,-1), (-1,0), (-1,1),
  (0,-1), (0,0), (0,1),
  (1,-1), (1,0), (1,1)
]

W = [
  ((1,1), 4), ((3,1), 5), ((5,1), 2), ((8,1), 7), ((9,1), 1),
  ((5,2), 5), ((9,2), 9),
  ((4,3), 4), ((5,3), 7), ((6,3), 3), ((7,3), 2),
  ((2,4), 4), ((3,4), 9), ((4,4), 6),
  ((3,5), 2),
  ((1,6), 3), ((5,6), 1), ((7,6), 7), ((8,6), 6),
  ((1,7), 1), ((8,7), 2),
  ((1,8), 5), ((7,8), 6), ((9,8), 3),
  ((3,9), 7), ((4,9), 3), ((6,9), 2), ((7,9), 9), ((8,9), 1), ((9,9), 5)
]

model = pyo.ConcreteModel()
N = 9
model.I = pyo.RangeSet(1,N)
model.x = pyo.Var(model.I, model.I, model.I, within=pyo.Binary)
model.constraints = pyo.ConstraintList()

def rule1(m, i, n):
  return sum(m.x[i,j,n] for j in m.I) == 1
model.Rule1 = pyo.Constraint(model.I, model.I, rule=rule1)

def rule2(m, j, n):
  return sum(m.x[i,j,n] for i in m.I) == 1
model.Rule2 = pyo.Constraint(model.I, model.I, rule=rule2)

for s in S:
  for n in model.I:
    model.constraints.add(sum(model.x[s[0]+v[0], s[1]+v[1], n] for v in V) == 1)

for w in W:
  coord = w[0]
  num = w[1]
  model.constraints.add(model.x[coord[0], coord[1], num] == 1)

for i in model.I:
  for j in model.I:
    model.constraints.add(sum(model.x[i,j,n] for n in model.I) == 1)

model.objective = pyo.Objective(expr=0, sense=pyo.minimize)
solver = SolverFactory("cbc")
result = solver.solve(model)

if result.solver.status == pyo.SolverStatus.ok and result.solver.termination_condition == pyo.TerminationCondition.optimal:
  for i in model.I:
    for j in model.I:
      item = [n for n in model.I if model.x[i,j,n] == 1]
      if len(item) == 1:
        sep = ''
        if j == 9:
          sep = '\n'
        elif j%3 == 0:
          sep = " | "
        else:
          sep = " "
      
        print("{}{}".format(item[0], sep), end="")
    if i%3 == 0 and i != 9:
      print("----------------------")