Pyomo: Простое ограничение неравенства требует неоправданного времени

Question

Я сталкиваюсь с проблемой при кодировании целочисленной линейной программы в pyomo.Модель, которую я хочу реализовать, должна найти оптимальный график для ряда задач.В частности, мы хотим запланировать обучение сотрудников на несколько периодов вперед.

Одним из ключевых ограничений, которое следует учитывать, является то, что не существует неограниченного количества инструкторов / тренеров, доступных каждый период для проведения соответствующих тренингов.Т.е. мы должны наложить ограничение мощности в каждом периоде.Чтобы гарантировать, что решатель всегда может найти решение, я также включаю переменную слабины для ограничения емкости и наказываю за это в задаче.Небольшой, воспроизводимый пример включен ниже.

Проблема : профилирование приведенного ниже кода показало, что при построении модели pyomo тратит большую часть времени, необходимого для построения модели, на наращивание потенциалаограничение (последнее ограничение в фрагменте кода ниже).При масштабировании до размера реального приложения время для создания именно этого ограничения становится непомерно большим.То есть, исключая время, затрачиваемое на решение проблемы, простое построение модели занимает несколько (читай:> 8) часов.Если исключить ограничение по емкости, это сокращается до десяти минут.

Вопросы: 1. Существует ли очевидная причина, по которой время сборки для этого ограничения намного выше, чем для любого издругие?Меня смущает тот факт, что для моей заявки число сотрудников (E ниже) составит ~ 2000, а количество отделов ~ 5 и количество периодов ~ 30.Следовательно, число строк, добавленных первым и вторым ограничением к общей матрице ограничений, намного выше, чем из-за ограничения емкости.По общему признанию, ограничение емкости должно суммироваться по большому количеству элементов, но я ожидаю, что это будет дешевая операция.Я ошибаюсь?2. Существует ли более эффективный синтаксис для кодирования этого ограничения?2а.Это поможет перейти к AbstractModel ()?2b.Поможет ли мне переместить K (словарь Python) в какой-нибудь набор Pyomo?

Вот небольшой пример, который, я надеюсь, поможет разобраться в проблеме:

Воспроизводимый пример:

from pyomo.environ import *

#set of employees
E = ['emp0', 'emp1']
#set of trainings per employee
J = ['train0', 'train1']
#periods in our planning horizon
T = [0, 1, 2]
#departments that have to provide trainers
DP = ['dept1', 'dept2']


#value for the company from each employee taking a training in particular period (want to maximize the sum over this)
v = {('emp0', 'train0', 0): 3,
     ('emp0', 'train0', 1): 0,
     ('emp0', 'train0', 2): 3,
     ('emp0', 'train1', 0): 9,
     ('emp0', 'train1', 1): 4,
     ('emp0', 'train1', 2): 8,
     ('emp1', 'train0', 0): 3,
     ('emp1', 'train0', 1): 0,
     ('emp1', 'train0', 2): 8,
     ('emp1', 'train1', 0): 7,
     ('emp1', 'train1', 1): 0,
     ('emp1', 'train1', 2): 7}

#M is a penalty for using the slack variable that leads the algorithm to violate the capacity constraint
M = 50000

#departments capacity in each period
K = {('dept1', 0): 10000,
     ('dept1', 1): 10000,
     ('dept1', 2): 10000,
     ('dept2', 0): 10000,
     ('dept2', 1): 10000,
     ('dept2', 2): 10000}

#number of coaches neeeded per employee-training-department-period capacity needed
a = {('emp0', 'train0', 'dept1', 0): 2,
     ('emp0', 'train0', 'dept2', 0): 2,
     ('emp0', 'train0', 'dept1', 1): 2,
     ('emp0', 'train0', 'dept2', 1): 1,
     ('emp0', 'train0', 'dept1', 2): 0,
     ('emp0', 'train0', 'dept2', 2): 2,
     ('emp0', 'train1', 'dept1', 0): 0,
     ('emp0', 'train1', 'dept2', 0): 1,
     ('emp0', 'train1', 'dept1', 1): 1,
     ('emp0', 'train1', 'dept2', 1): 0,
     ('emp0', 'train1', 'dept1', 2): 0,
     ('emp0', 'train1', 'dept2', 2): 1,
     ('emp1', 'train0', 'dept1', 0): 2,
     ('emp1', 'train0', 'dept2', 0): 1,
     ('emp1', 'train0', 'dept1', 1): 2,
     ('emp1', 'train0', 'dept2', 1): 0,
     ('emp1', 'train0', 'dept1', 2): 1,
     ('emp1', 'train0', 'dept2', 2): 0,
     ('emp1', 'train1', 'dept1', 0): 0,
     ('emp1', 'train1', 'dept2', 0): 1,
     ('emp1', 'train1', 'dept1', 1): 2,
     ('emp1', 'train1', 'dept2', 1): 1,
     ('emp1', 'train1', 'dept1', 2): 2,
     ('emp1', 'train1', 'dept2', 2): 2}

#duration of training per employee
D = {('emp0', 'train0'): 4,
     ('emp0', 'train1'): 4,
     ('emp1', 'train0'): 4,
     ('emp1', 'train1'): 3}

# Model setup:
model = ConcreteModel() 

model.J = Set(initialize = J)
model.E = Set(initialize = E)
model.T = Set(initialize = T)
model.EJ = Set(initialize = model.E*model.J)
model.EJT = Set(initialize = model.EJ*model.T)

model.x = Var(model.EJT, within=Binary)
model.y = Var(model.EJT, within=Binary)
model.R = Var(T, DP, within = NonNegativeIntegers, initialize = 0)

#objective: maximize the total value of trainings minus the penalty for violating the capacity constraint
def obj_rule(m):
    return sum(v[e,j,t]*m.x[e,j,t] for e in E for j in J for t in T) \
            - sum(m.R[t,dp] for dp in DP for t in T)*M 
model.obj = Objective(rule=obj_rule, sense = maximize)

#start each employee-training once
def one_per_emp_rule(m, e, j, t):
    return sum(m.x[e,j,t] for t in m.T)==1
model.one_per_emp_rule = Constraint(E, J, T, rule=one_per_emp_rule)

#for the duration of each training keep y=1
def y1_for_duration(m, e, j, t):
    return m.y[e,j,t] == sum(m.x[e,j,t] for t in T[max((t-D[e,j]+1), 0):t+1])         
model.y1_for_duration = Constraint(E, J, T, rule=y1_for_duration)    

# sum up all the consumed capacity across all employees per period, to check if departments can handle it
def period_capacity_dept(m, e, j, t, dp):
    return sum(a[e, j, dp, t]*m.y[e,j,t] for (e,j) in model.EJ)<= K[dp,t] + m.R[t,dp] 
model.period_capacity_dept = Constraint(E, J, T, DP, rule=period_capacity_dept)

# hand over to solver:

opt = SolverFactory('glpk')
results = opt.solve(model, tee=True)
results.write()
model.solutions.load_from(results)

for v in model.component_objects(Var, active=True):
    print ("Variable",v)
    varobject = getattr(model, str(v))
    for index in varobject:
        print ("   ",index, varobject[index].value)

Answer 1

В уравнении

def period_capacity_dept(m, e, j, t, dp):
    return sum(a[e, j, dp, t]*m.y[e,j,t] for (e,j) in model.EJ)<= K[dp,t] + m.R[t,dp] 
model.period_capacity_dept = Constraint(E, J, T, DP, rule=period_capacity_dept)

вы должны отбросить аргументы функции e, j.Теперь вы повторяете одно и то же ограничение много раз.Вы можете увидеть это, выполнив model.period_capacity_dept.pprint().Это показывает:

period_capacity_dept : Size=24, Index=period_capacity_dept_index, Active=True
    Key                            : Lower : Body                                                                                                   : Upper : Active
    ('emp0', 'train0', 0, 'dept1') :  -Inf :                                         2*y[emp1,train0,0] + 2*y[emp0,train0,0] - (10000 + R[0,dept1]) :   0.0 :   True
    ('emp0', 'train0', 0, 'dept2') :  -Inf :     y[emp0,train1,0] + y[emp1,train0,0] + y[emp1,train1,0] + 2*y[emp0,train0,0] - (10000 + R[0,dept2]) :   0.0 :   True
    ('emp0', 'train0', 1, 'dept1') :  -Inf : y[emp0,train1,1] + 2*y[emp1,train0,1] + 2*y[emp1,train1,1] + 2*y[emp0,train0,1] - (10000 + R[1,dept1]) :   0.0 :   True
    ('emp0', 'train0', 1, 'dept2') :  -Inf :                                             y[emp1,train1,1] + y[emp0,train0,1] - (10000 + R[1,dept2]) :   0.0 :   True
    ('emp0', 'train0', 2, 'dept1') :  -Inf :                                           y[emp1,train0,2] + 2*y[emp1,train1,2] - (10000 + R[2,dept1]) :   0.0 :   True
    ('emp0', 'train0', 2, 'dept2') :  -Inf :                      y[emp0,train1,2] + 2*y[emp1,train1,2] + 2*y[emp0,train0,2] - (10000 + R[2,dept2]) :   0.0 :   True
    ('emp0', 'train1', 0, 'dept1') :  -Inf :                                         2*y[emp1,train0,0] + 2*y[emp0,train0,0] - (10000 + R[0,dept1]) :   0.0 :   True
    ('emp0', 'train1', 0, 'dept2') :  -Inf :     y[emp0,train1,0] + y[emp1,train0,0] + y[emp1,train1,0] + 2*y[emp0,train0,0] - (10000 + R[0,dept2]) :   0.0 :   True
    ('emp0', 'train1', 1, 'dept1') :  -Inf : y[emp0,train1,1] + 2*y[emp1,train0,1] + 2*y[emp1,train1,1] + 2*y[emp0,train0,1] - (10000 + R[1,dept1]) :   0.0 :   True
    ('emp0', 'train1', 1, 'dept2') :  -Inf :                                             y[emp1,train1,1] + y[emp0,train0,1] - (10000 + R[1,dept2]) :   0.0 :   True
    ('emp0', 'train1', 2, 'dept1') :  -Inf :                                           y[emp1,train0,2] + 2*y[emp1,train1,2] - (10000 + R[2,dept1]) :   0.0 :   True
    ('emp0', 'train1', 2, 'dept2') :  -Inf :                      y[emp0,train1,2] + 2*y[emp1,train1,2] + 2*y[emp0,train0,2] - (10000 + R[2,dept2]) :   0.0 :   True
    ('emp1', 'train0', 0, 'dept1') :  -Inf :                                         2*y[emp1,train0,0] + 2*y[emp0,train0,0] - (10000 + R[0,dept1]) :   0.0 :   True
    ('emp1', 'train0', 0, 'dept2') :  -Inf :     y[emp0,train1,0] + y[emp1,train0,0] + y[emp1,train1,0] + 2*y[emp0,train0,0] - (10000 + R[0,dept2]) :   0.0 :   True
    ('emp1', 'train0', 1, 'dept1') :  -Inf : y[emp0,train1,1] + 2*y[emp1,train0,1] + 2*y[emp1,train1,1] + 2*y[emp0,train0,1] - (10000 + R[1,dept1]) :   0.0 :   True
    ('emp1', 'train0', 1, 'dept2') :  -Inf :                                             y[emp1,train1,1] + y[emp0,train0,1] - (10000 + R[1,dept2]) :   0.0 :   True
    ('emp1', 'train0', 2, 'dept1') :  -Inf :                                           y[emp1,train0,2] + 2*y[emp1,train1,2] - (10000 + R[2,dept1]) :   0.0 :   True
    ('emp1', 'train0', 2, 'dept2') :  -Inf :                      y[emp0,train1,2] + 2*y[emp1,train1,2] + 2*y[emp0,train0,2] - (10000 + R[2,dept2]) :   0.0 :   True
    ('emp1', 'train1', 0, 'dept1') :  -Inf :                                         2*y[emp1,train0,0] + 2*y[emp0,train0,0] - (10000 + R[0,dept1]) :   0.0 :   True
    ('emp1', 'train1', 0, 'dept2') :  -Inf :     y[emp0,train1,0] + y[emp1,train0,0] + y[emp1,train1,0] + 2*y[emp0,train0,0] - (10000 + R[0,dept2]) :   0.0 :   True
    ('emp1', 'train1', 1, 'dept1') :  -Inf : y[emp0,train1,1] + 2*y[emp1,train0,1] + 2*y[emp1,train1,1] + 2*y[emp0,train0,1] - (10000 + R[1,dept1]) :   0.0 :   True
    ('emp1', 'train1', 1, 'dept2') :  -Inf :                                             y[emp1,train1,1] + y[emp0,train0,1] - (10000 + R[1,dept2]) :   0.0 :   True
    ('emp1', 'train1', 2, 'dept1') :  -Inf :                                           y[emp1,train0,2] + 2*y[emp1,train1,2] - (10000 + R[2,dept1]) :   0.0 :   True
    ('emp1', 'train1', 2, 'dept2') :  -Inf :                      y[emp0,train1,2] + 2*y[emp1,train1,2] + 2*y[emp0,train0,2] - (10000 + R[2,dept2]) :   0.0 :   True

Если мы используем:

def period_capacity_dept(m, t, dp):
    return sum(a[e, j, dp, t]*m.y[e,j,t] for (e,j) in model.EJ)<= K[dp,t] + m.R[t,dp] 
model.period_capacity_dept = Constraint(T, DP, rule=period_capacity_dept)

, вы увидите:

period_capacity_dept : Size=6, Index=period_capacity_dept_index, Active=True
    Key          : Lower : Body                                                                                                   : Upper : Active
    (0, 'dept1') :  -Inf :                                         2*y[emp1,train0,0] + 2*y[emp0,train0,0] - (10000 + R[0,dept1]) :   0.0 :   True
    (0, 'dept2') :  -Inf :     y[emp1,train0,0] + y[emp0,train1,0] + 2*y[emp0,train0,0] + y[emp1,train1,0] - (10000 + R[0,dept2]) :   0.0 :   True
    (1, 'dept1') :  -Inf : 2*y[emp1,train0,1] + y[emp0,train1,1] + 2*y[emp0,train0,1] + 2*y[emp1,train1,1] - (10000 + R[1,dept1]) :   0.0 :   True
    (1, 'dept2') :  -Inf :                                             y[emp0,train0,1] + y[emp1,train1,1] - (10000 + R[1,dept2]) :   0.0 :   True
    (2, 'dept1') :  -Inf :                                           y[emp1,train0,2] + 2*y[emp1,train1,2] - (10000 + R[2,dept1]) :   0.0 :   True
    (2, 'dept2') :  -Inf :                      y[emp0,train1,2] + 2*y[emp0,train0,2] + 2*y[emp1,train1,2] - (10000 + R[2,dept2]) :   0.0 :   True

Теперь все дублирующие ограничения опущены.