Pyomo Objective trying to maximize the sum of the product of consecutive indexed variables

I made some code for a scheduling task using Pyomo in SolverStudio. The overall goal of the solve is to maximise the amount of consecutive weeks (model.Dates) that are worked while being constrained by the amount of people working per week (doctorsequalteamsRule) and the amount of weeks in the time period that a person can work (docsworktocontractRule). There may be a varying number of teams and a person can only work on one team at a time (onedocatatimeRule). The model.NamesLarge is a stupid way to index the worker's names and the various (up to 10) teams they might work on (1-10 is worker 1, 11-20 is worker 2, 21-30...).

from pyomo.environ import * # For Pyomo 4.0 & later
# from coopr.pyomo import * # For earlier versions
model = AbstractModel()
## Define sets
model.NamesSmall = Set()
model.NamesLarge = Set()
model.Dates = Set()
model.Teams = Set()
## Define parameters
#how many weeks a doc works in the year
model.doctorsperweek = Param(model.NamesSmall)
#how many teams/docs are on each week
model.numberworked = Param(model.Dates)
model.abletowork = Param(model.NamesSmall, model.Dates)
model.mustwork = Param(model.NamesSmall, model.Dates)
model.maxweeksinarow = Param(model.NamesSmall, model.Dates)
model.teamsperweek = Param(model.Teams, model.Dates)
## Define variables
model.variables = Var(model.NamesLarge, model.Dates, within = Binary, initialize=0)
## Define Objective Function
def costRule(model):
return sum(model.variables[nm,qm]*model.variables[nm,qm+1] for nm in model.NamesLarge for qm in range(1,len(model.Dates)))
model.SolverResult = Objective(rule=costRule, sense=maximize)
## Constraints
#Maybe need a constraint that if a var = 1 then the next date other var in the block of 10 cant be 1
#Need a constraint that the different teams are specifically reflected in the variables
#Ensures the docs are working the amount of weeks in their contract
blahg = []
def docsworktocontractRule(model,tr):
for mr in model.NamesSmall:
blah = sum(model.variables[(10*mr)-p,bo] for bo in model.Dates for p in range(10))
blahg.append(blah)
return blahg[tr-1] == model.doctorsperweek[tr]
model.docsworktocontractConstraint = Constraint(model.NamesSmall, rule=docsworktocontractRule)
#Ensures same amount of doctors as teams per week
def doctorsequalteamsRule(model,i):
return  sum(model.variables[m,i] for m in model.NamesLarge) == value(model.numberworked[i])
model.doctorsequalteamsConstraint = Constraint(model.Dates, rule=doctorsequalteamsRule)
def onedocatatimeRule(model,fe,da):
grapesof = []
grape = sum(model.variables[(10*fe)-p,da] for p in range(0,10))
grapesof.append(grape)
return sum(grapesof[yh-1] for yh in range(len(grapesof))) <= 1
model.onedocatatimeConstraint = Constraint(model.NamesSmall, model.Dates, rule=onedocatatimeRule)

The constraints work, but the objective takes an incredible amount of time to solve even only for 6 people, 2 teams and 10 weeks (I'll need to solve for ~40 people, 52 weeks and 7-8 teams).

Is there any way to make the objective solve more efficiently?

I was thinking maybe instead to using * (where 0x0 = 0, 0x1 = 0, 1x0 = 0 and 1x1 = 1) I could use + and then filter out to only sum the sum of variables that == 2, but this doesn't work as the variables are unknown at the time of instance creation?

I'm at a loss, maybe I'm missing something obvious or have I completely set the whole thing up inefficiently?

You have (inadvertently?) made this a non-linear model by multiplying variables together. It can be made linear, which is much simpler to solve. If the linear/non-linear/integer issue are confusing, you should consult a linear programming text or tutorial on how those types of models are solved and why non-linear is quite different.

That aside, here is an example that shows a linear construct for counting sequential weeks by doing a little integer math on the assignment variables.

A side note: If you are going to post your code and ask for help, you might want to take a few more minutes and clean it up. A variable called variables ? Inconsistent (and non-standard) capitalization, variables called blah and grape? C'mon.... show some pride in the work.

Code:

import pyomo.environ as pyo
### DATA
weeks = 8
doc_limits = { 	'bob':		3,	# doc : max assignments
'sam':  	4,
'sally': 	4,
'xavier':	2}
requirements = {	1: 2,     # week : number of docs
2: 2,
3: 1,
4: 0,
5: 2,
6: 2,
7: 2,
8: 2}
### MODEL
m = pyo.ConcreteModel('doc_sked')
m.W = pyo.Set(initialize=range(1, weeks+1), ordered=True, doc='week')
m.W2 = pyo.Set(initialize=range(2, weeks+1), ordered=True, doc='week 2 plus')
m.D = pyo.Set(initialize=doc_limits.keys(), doc='doc')
m.assign = pyo.Var(m.D, m.W, domain=pyo.Binary)
m.consecutive_assign = pyo.Var(m.D, m.W2, domain=pyo.Binary)
### OBJ
m.obj = pyo.Objective(expr=pyo.sum_product(m.consecutive_assign), sense=pyo.maximize)
### CONSTRAINTS
# don't bust max assignments
def max_assignments(m, d):
return sum(m.assign[d, w] for w in m.W) <= doc_limits[d]
m.C1 = pyo.Constraint(m.D, rule=max_assignments)
# don't bust requirement
def requirement_limit(m, w):
return sum(m.assign[d, w] for d in m.D) <= requirements[w]
m.C2 = pyo.Constraint(m.W, rule=requirement_limit)
# link assignment to consecutive assignment indicator
def link(m, d, w):
return m.consecutive_assign[d, w] <= (m.assign[d, w] + m.assign[d, m.W.prev(w)])/2
m.C3 = pyo.Constraint(m.D, m.W2, rule=link)
m.pprint()
solver = pyo.SolverFactory('cbc')
res = solver.solve(m)
print(res)
m.assign.display()