Блог

Главная - Вопросы по программированию - Как добавить фрагментарное линейное ограничение в pyomo

Как добавить фрагментарное линейное ограничение в pyomo

#pyomo

Вопрос:

В задаче распределения задач давайте предположим:

  1. нам нужно выделить 4 задачи на 4 дня.
  2. Есть только один работник, и он должен работать не более 8 часов в день.
  3. Мы хотим минимизировать общее рабочее время.
  4. Мы можем объединить задачи и сэкономить время за один день. Необходимое количество часов равно

0 задач —> 0 часов

1 задача -> 6 часов

2 задачи -> 8 часов

Иллюстрация для # задач V.S. # необходимых часов:

введите описание изображения здесь

Синяя пунктирная линия: y = 6 * n

Черная строка: фактическое время = min(6 * n, 6 2 * (n — 1))

Красная строка: дневной лимит

Эта вычислительная функция (для черной линии) может быть определена в Python как

 def compute_hours(n):
    return min(6*n, 6   2 * (n - 1))

Используя функцию, определенную выше в настройке ограничения, pyomo выдавайте ошибки, когда я помещаю это в свой список ограничений.

 > File "<ipython-input-69-57dbda13a71f>", line 23, in compute_hours
>     return min(6*n, 6   2 * (n - 1))   File "pyomocoreexprlogical_expr.pyx", line 183, in
> pyomo.core.expr.logical_expr.InequalityExpression.__bool__
> pyomo.common.errors.PyomoException: Cannot convert non-constant
> expression to bool. This error is usually caused by using an
> expression in a boolean context such as an if statement. For example, 
>     m.x = Var()
>     if m.x <= 0:
>         ... would cause this exception.

Ошибку также легко воспроизвести, используя следующий пример кода:

 # Allocating tasks to days

from pyomo.environ import *
import pandas as pd

m = ConcreteModel()

# We have 4 tasks in total. Each cost 6 hours
m.TASKS = ('task1', 'task2', 'task3', 'task4')

# There are 4 days to allocate to
m.DAYS = (1, 2, 3, 4)

# Hours per task
HOURS_PER_TASK = 6


def compute_hours(n):
    return min(6*n, 6   2 * (n - 1))

# for x in range(0, 4):
#     print(x, compute_hours(x))

# Create a decision array for the allocation
m.ALLOCATION = Var(m.TASKS, m.DAYS, domain=Binary)

# Create an array for total hours computation
m.TOTAL_HOURS = Var(m.DAYS, domain=Integers)

# minimize the number of days that are allocated tasks
m.OBJ = Objective(expr=sum([m.TOTAL_HOURS[day] for day in m.DAYS]), sense=minimize)

m.c = ConstraintList()

# One task is done once only in all days
for task in m.TASKS:
    m.c.add(sum([m.ALLOCATION[task, day] for day in m.DAYS]) == 1)

# Compute Total hours a day
for day in m.DAYS:
    m.c.add(sum([m.ALLOCATION[task, day] for task in m.TASKS])*HOURS_PER_TASK == m.TOTAL_HOURS[day])
    # The following computation does not work
    #m.c.add(compute_hours(sum([m.ALLOCATION[task, day] for task in m.TASKS])*HOURS_PER_TASK) == m.TOTAL_HOURS[day])


# Add max working hours per day
for day in m.DAYS:
    m.c.add(m.TOTAL_HOURS[day] <= 8)

SolverFactory('glpk').solve(m).write()

# Show the results
SCHEDULE_df = pd.DataFrame(index = m.TASKS, columns= m.DAYS)
SCHEDULE_HOURS_df = pd.DataFrame(index = m.TASKS, columns= m.DAYS)

for i in m.TASKS:
    for j in m.DAYS:
        SCHEDULE_df.loc[i,j] = m.ALLOCATION[i,j]()

print('------------------------------------------')
print('ALLOCATION:')
print(SCHEDULE_df)

print('------------------------------------------')
print('Total hours per day:')
print(m.TOTAL_HOURS.get_values())
print('Total hours:')
print(value(m.OBJ))

Этот пример все еще может дать результаты, если мы предположим, что каждая задача всегда занимает 6 часов.

 ------------------------------------------
ALLOCATION:
       1  2  3  4
task1  1  0  0  0
task2  0  1  0  0
task3  0  0  1  0
task4  0  0  0  1
------------------------------------------
Total hours per day:
{1: 6.0, 2: 6.0, 3: 6.0, 4: 6.0}
Total hours:
24.0

Если мы сможем реализовать объединение, то этому работнику нужно будет работать только два дня с рабочей нагрузкой 8 часов в день (вместо работы 6 часов в день в течение 4 дней).

Ответ №1:

Джон.

Pyomo получил кусочно-линейное выражение для этих случаев. ( kernel слой также получил библиотеку кусочных функций)

В вашем случае, на самом деле, легко использовать Piecewise in Pyomo .

Вам нужно создать вспомогательную переменную с тем же индексом, что m.TOTAL_HOURS и (запрошенный pyomo ), назовем его m.TASK_PER_DAY , который будет учитывать количество задач, назначенных на каждый день. Затем этот TASK_PER_DAY параметр будет вычисляться (ограничен) с m.ALLOCATION использованием суммирования задачи каждый день (аналогично тому, как вы пытаетесь ограничить). Затем вы можете использовать pyomo Piecewise для вычисления в TOTAL_HOURS зависимости от TASK_PER_DAY

Редактировать: Дополнительные пояснения к pw_pts и f_rule , которые являются источником ошибок в моделировании Piecewise . pw_pts являются ли точки останова (независимые значения, домен) в Piecewise . f_rule является известным значением в заданных точках pw_pts . Это может быть функция или пара точек. Если вы хотите линейно интерполировать точки (как ваше изображение), вам просто нужно определить такие точки pw_pts=[0,1,2,3,...] и f_rule=[0,6,8,10,...] , но для постоянных кусочных функций вам нужно указать одну и ту же точку дважды, чтобы гарантировать, что для точек домена кусочно возвращает то же значение, что и в пошаговой функции. В вашем моделировании, поскольку переменные есть integers , проблем нет, но лучше устранить проблему.

 m.TASK_PER_DAY = Var(m.DAYS, domain=Integers)
m.PIECEWISE_DAILY_HOURS = Piecewise(
    m.DAYS,
    m.TOTAL_HOURS,
    m.TASK_PER_DAY,
    pw_pts=[0,1,2,3],
    f_rule=[0,6,8,10],
    pw_repn='SOS2',
    pw_constr_type='EQ',
    unbounded_domain_var=True
)
#Constraining TASK_PER_DAY
def Assign_Task(model, day):
    return model.TASK_PER_DAY[day] == sum(model.ALLOCATION[task,day] for task in model.TASKS)
m.assign_task = Constraint(m.DAYS, rule=Assign_Task)

При таком подходе вы не ограничиваете день для выполнения задачи, поэтому назначенный день может быть любой парой дней, но каждый день, безусловно, получает по 2 задачи в день. Если вы хотите ограничить модель для выполнения задачи как можно скорее (в первые дни), вам нужно добавить это новое ограничение.

Следующий код может воспроизвести желаемый результат:

 import pyomo.environ as pyo

m = pyo.ConcreteModel()

#Sets
# We have 4 tasks in total. Each cost 6 hours
m.TASKS = ('task1', 'task2', 'task3', 'task4')
# There are 4 days to allocate to
m.DAYS = (1, 2, 3, 4)
# Hours per task
HOURS_PER_TASK = 6

#Variables
# Create a decision array for the allocation
m.ALLOCATION = pyo.Var(m.TASKS, m.DAYS, domain=Binary)
m.TOTAL_HOURS = pyo.Var(m.DAYS, domain=Integers)
m.TASK_PER_DAY = pyo.Var(m.DAYS, domain=Integers)

#Piecewise
m.PIECEWISE_DAILY_HOURS = pyo.Piecewise(
    m.DAYS,
    m.TOTAL_HOURS,
    m.TASK_PER_DAY,
    pw_pts=[0,0,1,1,2,2,3,3],
    f_rule=[0,0,6,6,8,8,10,10],
    pw_repn='SOS2',
    pw_constr_type='EQ',
    unbounded_domain_var=True
)

#Objective
# minimize the number of days that are allocated tasks
m.OBJ = pyo.Objective(expr=sum([m.TOTAL_HOURS[day] for day in m.DAYS]), sense=minimize)

m.c = pyo.ConstraintList()

# One task is done once only in all days
for task in m.TASKS:
    m.c.add(sum([m.ALLOCATION[task, day] for day in m.DAYS]) == 1)

#Compute the number of tasks per day
def Assign_Task(model, day):
    return model.TASK_PER_DAY[day] == sum(model.ALLOCATION[task,day] for task in model.TASKS)

m.assign_task = pyo.Constraint(m.DAYS, rule=Assign_Task)

# Add max working hours per day
for day in m.DAYS:
    m.c.add(m.TOTAL_HOURS[day] <= 8)

pyo.SolverFactory('gurobi').solve(m, tee=True)

Затем время настаивания раствора сокращают до 16 часов.

 >>>m.OBJ()
16.0
>>>m.TOTAL_HOURS.display()
TOTAL_HOURS : Size=4, Index=TOTAL_HOURS_index
    Key : Lower : Value : Upper : Fixed : Stale : Domain
      1 :  None :  -0.0 :  None : False : False : Integers
      2 :  None :   8.0 :  None : False : False : Integers
      3 :  None :  -0.0 :  None : False : False : Integers
      4 :  None :   8.0 :  None : False : False : Integers
>>>m.TASK_PER_DAY.display()
TASK_PER_DAY : Size=4, Index=TASK_PER_DAY_index
    Key : Lower : Value : Upper : Fixed : Stale : Domain
      1 :  None :   0.0 :  None : False : False : Integers
      2 :  None :   2.0 :  None : False : False : Integers
      3 :  None :   0.0 :  None : False : False : Integers
      4 :  None :   2.0 :  None : False : False : Integers
>>>m.ALLOCATION.display()
ALLOCATION : Size=16, Index=ALLOCATION_index
    Key          : Lower : Value : Upper : Fixed : Stale : Domain
    ('task1', 1) :     0 :  -0.0 :     1 : False : False : Binary
    ('task1', 2) :     0 :   1.0 :     1 : False : False : Binary
    ('task1', 3) :     0 :  -0.0 :     1 : False : False : Binary
    ('task1', 4) :     0 :  -0.0 :     1 : False : False : Binary
    ('task2', 1) :     0 :  -0.0 :     1 : False : False : Binary
    ('task2', 2) :     0 :  -0.0 :     1 : False : False : Binary
    ('task2', 3) :     0 :  -0.0 :     1 : False : False : Binary
    ('task2', 4) :     0 :   1.0 :     1 : False : False : Binary
    ('task3', 1) :     0 :  -0.0 :     1 : False : False : Binary
    ('task3', 2) :     0 :   1.0 :     1 : False : False : Binary
    ('task3', 3) :     0 :  -0.0 :     1 : False : False : Binary
    ('task3', 4) :     0 :  -0.0 :     1 : False : False : Binary
    ('task4', 1) :     0 :  -0.0 :     1 : False : False : Binary
    ('task4', 2) :     0 :  -0.0 :     1 : False : False : Binary
    ('task4', 3) :     0 :  -0.0 :     1 : False : False : Binary
    ('task4', 4) :     0 :   1.0 :     1 : False : False : Binary

Комментарии:

1. Большое вам спасибо 🙂 Кстати, знаете ли вы, что какой-либо решатель с открытым исходным кодом также может поддерживать ограничения уровня SOS 2? Я тестировал glpk , и это не работает.

2. Нет, извините. Но если вы измените Piecewise аргумент на pw_repn='DCC' (проверьте другие параметры в документации Pyomo), вы можете использовать glpk

3. Да! Это работает! Спасибо.

Метки: