Задача планирования Job-Shop: модели смешанно-целочисленного программирования

Математическое моделирование и реализация Python классической задачи секвенирования с использованием Pyomo

Задача планирования рабочих мест (JSSP) — это широко изучаемая задача оптимизации, имеющая несколько промышленных применений. Цель состоит в том, чтобы определить, как свести к минимуму период времени, необходимый для выделения общих ресурсов (машин) для выполнения конкурирующих действий (работ). Что касается других задач оптимизации, смешанно-целочисленное программирование может быть эффективным инструментом для получения хороших решений, хотя для крупномасштабных случаев, вероятно, следует прибегать к эвристикам.

В этой статье можно найти две самые распространенные формулировки смешанно-целочисленного программирования для JSSP с реализацией на Python с использованием библиотеки pyomo (Bynum et al., 2021). Те, кому интересны подробности, могут ознакомиться с полным кодом, доступным в этом git-репозитории.

Если вы незнакомы со смешанно-целочисленным программированием или оптимизацией в целом, вы можете получить больше опыта после прочтения этого введения по этому вопросу.

Введение в линейное программирование смешанных целых чисел: задача о рюкзаке
Узнайте, как решать задачи оптимизации в Python с использованием scipy и pyomotowardsdatascience.com

Теперь давайте погрузимся!

Постановка задачи

Предположим, что набор заданий J необходимо обработать на наборе машин M, каждое в заданном порядке. Например, задание № 1 может потребоваться для обработки на машинах (1, 4, 3, 2), а задание № 2 — на (2, 3, 4, 1). В этом случае, прежде чем попасть на машину 4, задание 1 должно быть отправлено на машину 1. Аналогично, перед тем как перейти на машину 1, задание 2 должно быть обработано на машине 4.

Каждая машина может одновременно обрабатывать только одно задание. Операции определяются парами (машина, задание), и каждая из них имеет определенное время обработки p. Таким образом, общий срок службы зависит от того, как распределяются ресурсы для выполнения задач.

На рисунке ниже показана оптимальная последовательность операций для простого экземпляра с 5 машинами и 4 заданиями. Обратите внимание, что каждый компьютер одновременно обрабатывает только одно задание, и каждое задание обрабатывается одновременно только одним компьютером.

Что касается других задач оптимизации, мы должны преобразовать эти правила в математические уравнения, чтобы получить разумное распределение ресурсов. Поэтому в следующем разделе мы рассмотрим две обычные формулировки JSSP.

Смешанно-целочисленные модели

После исследования Ku & Beck (2016) будут представлены две формулировки JSSP: дизъюнктивная модель (Manne, 1960) и модель временного индекса (Bowman, 1959; Kondili, 1993). Желающие могут обратиться к Вагнеру (1959) за третьей формулировкой (ранговой). Дизъюнктивная модель, несомненно, является наиболее эффективной из трех для исходной задачи. Однако с другими может быть проще справиться при включении новых ограничений, которые могут возникнуть в реальных задачах.

В дизъюнктивной модели рассмотрим набор J заданий и набор M машин. Каждое задание j должно соответствовать порядку обработки (σʲ₁, σʲ ₂, …, σʲₖ), и каждая операция (m, j) имеет время обработки п. Рассматриваемые переменные решения — это время запуска задания j на машине m, xₘⱼ; двоичный файл, который отмечает приоритет задания i перед j на машине m, zₘᵢⱼ; и общий срок эксплуатации, C, что само по себе является целью минимизации.

Нам нужно создать ограничения, чтобы гарантировать, что:

Время начала задания j на машине m должно быть больше времени начала предыдущей операции задания j плюс время его обработки.
Каждая машина обрабатывает только одно задание за раз. Для этого мы утверждаем, что если i предшествует j в машине m, время начала задания j в машина m должна быть больше или равна времени запуска задания i плюс время его обработки.
Из каждой пары заданий i, j один элемент должен предшествовать другому для каждой машины m в M.
Общий срок изготовления больше, чем время начала каждой операции плюс время ее обработки.

И получаем следующую формулировку:

В котором V — произвольно большое значение (большое M) ограничения «или-или».

Следующей исследуемой формулировкой будет модель с временной индексацией. Она ограничена в том смысле, что можно учитывать только целочисленное время обработки, и можно заметить, что она создает матрицу ограничений с несколькими ненулевыми элементами, что делает ее вычислительно более затратной, чем дизъюнктивную модель. Кроме того, по мере увеличения времени обработки увеличивается и количество переменных решения.

В модели с временной индексацией мы будем рассматривать те же наборы заданий J и машин M, кроме набора дискретных интервалов T. Выбор размера T может быть ориентирован так же, как и определение V: сумма всех времен обработки. Будут использоваться те же параметры порядка заданий и времени обработки. Однако при таком подходе мы рассматриваем только двоичные переменные, которые отмечают, что задание j начинается на машине m в момент времени t, xₘⱼₜ, кроме действительного (или целого) значения, которое составляет C.

Сформулируем ограничения:

Каждое задание j на машине m запускается только один раз.
Убедитесь, что каждый компьютер обрабатывает только одно задание за раз. И это сложно в подходе с временной индексацией. Для этого мы утверждаем, что не более одного задания j может начаться на машине m в течение промежутка времени между текущим временем t и pₘⱼ предыдущий раз. Для каждой машины и момента времени.
Время начала задания j на машине m должно быть больше времени начала предыдущей операции задания j плюс время его обработки.
Общий срок изготовления больше, чем время начала каждой операции плюс время ее обработки.

Выполнение

Прежде чем погрузиться в модели, давайте создадим несколько служебных классов для обработки параметров этой задачи. Первым будет JobSequence дочерний класс Python list с методами для поиска предыдущих и следующих элементов в последовательности. Это будет полезно при обращении к последовательности машин для каждого задания.

class JobSequence(list):
    
    def prev(self, x):
        if self.is_first(x):
            return None
        else:
            i = self.index(x)
            return self[i - 1]
    
    def next(self, x):
        if self.is_last(x):
            return None
        else:
            i = self.index(x)
            return self[i + 1]
    
    def is_first(self, x):
        return x == self[0]
    
    def is_last(self, x):
        return x == self[-1]
    
    def swap(self, x, y):
        i = self.index(x)
        j = self.index(y)
        self[i] = y
        self[j] = x
    
    def append(self, __object) -> None:
        if __object not in self:
            super().append(__object)
        else:
            pass

Теперь давайте создадим класс белой метки для параметров. Он должен хранить набор заданий J, набор машин M, последовательность операций каждого задания j в словаре JobSequences и время обработки каждой пары m, j в индексном словаре кортежей p_times.

class JobShopParams:
    
    def __init__(self, machines, jobs, p_times, seq):
        self.machines = machines
        self.jobs = jobs
        self.p_times = p_times
        self.seq = seq

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

import numpy as np


class JobShopRandomParams(JobShopParams):
    
    def __init__(self, n_machines, n_jobst, t_span=(1, 20), seed=None):

        self.t_span = t_span
        self.seed = seed
        
        machines = np.arange(n_machines, dtype=int)
        jobs = np.arange(n_jobs, dtype=int)
        p_times = self._random_times(machines, jobs, t_span)
        seq = self._random_sequences(machines, jobs)
        super().__init__(machines, jobs, p_times, seq)
    
    def _random_times(self, machines, jobs, t_span):
        np.random.seed(self.seed)
        t = np.arange(t_span[0], t_span[1])
        return {
            (m, j): np.random.choice(t)
            for m in machines
            for j in jobs
        }
    
    def _random_sequences(self, machines, jobs):
        np.random.seed(self.seed)
        return {
            j: JobSequence(np.random.permutation(machines))
            for j in jobs
        }

Теперь мы можем легко создавать случайные проблемы для проверки наших моделей.

На следующих шагах мы создадим три класса, которые наследуются от ConcreteModel от pyomo. Первым будет класс white label для моделей MIP. Второй и третий будут классами дизъюнктивной и индексированной по времени модели соответственно.

import pyomo.environ as pyo


class JobShopModel(pyo.ConcreteModel):
    
    def __init__(self, params, **kwargs):
        super().__init__()
        self.params = params
        self._construct_sets()
        self._construct_params()
    
    @property
    def seq(self):
        return self.params.seq
    
    def _construct_sets(self):
        self.M = pyo.Set(initialize=self.params.machines)
        self.J = pyo.Set(initialize=self.params.jobs)
    
    def _construct_params(self):
        self.p = pyo.Param(self.M, self.J, initialize=self.params.p_times)
        self.V = sum(self.p[key] for key in self.p)

Можно заметить, что наборы заданий J и машин M хранятся в атрибутах экземпляра с одинаковым именем. Атрибут p содержит время обработки, а V — разумный верхний предел времени выполнения.

Давайте теперь создадим дизъюнктивную модель, класс DisjModel.

class DisjModel(JobShopModel):

    def __init__(self, params, **kwargs):
        super().__init__(params, **kwargs)
        self._create_vars()
        self._create_cstr()
        self.obj = pyo.Objective(rule=self.C, sense=pyo.minimize)

    def _create_vars(self):
        self.x = pyo.Var(self.M, self.J, within=pyo.NonNegativeReals)
        self.z = pyo.Var(self.M, self.J, self.J, within=pyo.Binary)
        self.C = pyo.Var(within=pyo.NonNegativeReals)

    def _create_cstr(self):
        self.cstr_seq = pyo.Constraint(self.M, self.J, rule=cstr_seq)
        self.cstr_precede = pyo.Constraint(self.M, self.J, self.J, rule=cstr_precede)
        self.cstr_comp_precede = pyo.Constraint(self.M, self.J, self.J, rule=cstr_comp_precede)
        self.cstr_total_time = pyo.Constraint(self.J, rule=cstr_total_time)

Экземпляры DisjModel содержат атрибуты x, z и C — ранее описанных переменных. Цель довольно проста: минимизировать одну из переменных решения: C. Обратите внимание, что нам все еще нужно определить правила для ограничений. Они определяются в том же порядке, в котором они были перечислены ранее при введении модели. Давайте теперь напишем их в стиле pyomo.

def cstr_seq(model, m, j):
    o = model.seq[j].prev(m)
    if o is not None:
        return model.x[o, j] + model.p[o, j] <= model.x[m, j]
    else:
        return pyo.Constraint.Skip


def cstr_precede(model, m, j, k):
    if j != k:
        return model.x[m, j] + model.p[m, j] <= model.x[m, k] + model.V * (1 - model.z[m, j, k])
    else:
        return pyo.Constraint.Skip


def cstr_comp_precede(model, m, j, k):
    if j != k:
        return model.z[m, j, k] + model.z[m, k, j] == 1
    else:
        return model.z[m, j, k] == 0


def cstr_total_time(model, j):
    m = model.seq[j][-1]
    return model.x[m, j] + model.p[m, j] <= model.C

И мы готовы решать проблемы JSSP с подходом дизъюнктивной модели. Определим также индексированную по времени модель.

class TimeModel(JobShopModel):

    def __init__(self, params, **kwargs):
        super().__init__(params, **kwargs)
        self.T = pyo.RangeSet(self.V)
        self._create_vars()
        self._create_cstr()
        self.obj = pyo.Objective(rule=self.C, sense=pyo.minimize)

    def _create_vars(self):
        self.x = pyo.Var(self.M, self.J, self.T, within=pyo.Binary)
        self.C = pyo.Var(within=pyo.NonNegativeReals)

    def _create_cstr(self):
        self.cstr_unique_start = pyo.Constraint(self.M, self.J, rule=cstr_unique_start)
        self.cstr_unique_machine = pyo.Constraint(self.M, self.T, rule=cstr_unique_machine)
        self.cstr_precede = pyo.Constraint(self.M, self.J, rule=cstr_precede)
        self.cstr_total_time = pyo.Constraint(self.J, rule=cstr_total_time)

    def _get_elements(self, j):
        machines = [x.index()[0] for x in self.x[:, j, :] if x.value == 1]
        starts = [x.index()[2] for x in self.x[:, j, :] if x.value == 1]
        spans = [self.p[m, j] for m in machines]
        return machines, starts, spans

И снова ограничения были определены в том же порядке, в котором они были представлены ранее. Напишем их тоже в стиле pyomo.

def cstr_unique_start(model, m, j):
    return sum(model.x[m, j, t] for t in model.T) == 1


def cstr_unique_machine(model, m, t):
    total = 0
    start = model.T.first()
    for j in model.J:
        duration = model.p[m, j]
        t0 = max(start, t - duration + 1)
        for t1 in range(t0, t + 1):
            total = total + model.x[m, j, t1]
    return total <= 1


def cstr_precede(model, m, j):
    o = model.seq[j].prev(m)
    if o is None:
        prev_term = 0
    else:
        prev_term = sum(
            (t + model.p[o, j]) * model.x[o, j, t]
            for t in model.T
        )
    current_term = sum(
        t * model.x[m, j, t]
        for t in model.T
    )
    return prev_term <= current_term


def cstr_total_time(model, j):
    m = model.seq[j][-1]
    return sum((t + model.p[m, j]) * model.x[m, j, t] for t in model.T) <= model.C

И мы готовы проверить, как эти модели работают в некоторых случайно сгенерированных задачах!

Полученные результаты

Давайте создадим случайную задачу 4x3 (JxM) и посмотрим, как работают наши модели.

params = JobShopRandomParams(3, 4, seed=12)
disj_model = DisjModel(params)
time_model = TimeModel(params)

Чтобы решить эти примеры, я буду использовать решатель CBC с открытым исходным кодом. Вы можете скачать бинарники CBC с AMPL или по этой ссылке. Вы также можете найти учебник по установке здесь. Поскольку исполняемый файл CBC включен в переменную PATH моей системы, я могу создать экземпляр решателя без указания пути к исполняемому файлу. Если у вас нет, проанализируйте аргумент ключевого слова исполняемый файл с путем к вашему исполняемому файлу.

В качестве альтернативы для решения этой проблемы можно было использовать GLPK (или любой другой решатель, совместимый с pyomo). Последнюю доступную версию GLPK можно найти здесь, а исполняемые файлы Windows можно найти здесь.

# Time limited in 20 seconds
solver = pyo.SolverFactory("cbc", options=dict(cuts="on", sec=20))

# solve
res_disj = solver.solve(disj_model, tee=False)
res_time = solver.solve(time_model, tee=False)

Решатель без труда нашел оптимальное решение для дизъюнктивной модели и доказал оптимальность менее чем за одну секунду.

Однако мы видим, что даже для этой простой задачи решатель не смог найти оптимальное решение для временной модели в пределах 20 секунд.

Удивительно видеть разницу в производительности двух моделей с одной и той же идеей, просто переставляя математические уравнения.

Кстати, желающие могут найти полный код (включая сюжеты) в этом репозитории.

дальнейшее чтение

Для более крупных экземпляров из-за комбинаторных аспектов этой проблемы даже высокопроизводительные коммерческие решатели, такие как Gurobi или Cplex, могут столкнуться с трудностями при предоставлении качественных решений и доказательстве их оптимальности. В этом контексте метаэвристика может быть интересной альтернативой. Я бы посоветовал заинтересованному читателю поискать статьи Параллельный GRASP с перекомпоновкой путей для планирования рабочих мест (Aiex et al., 2003) и Расширенный графический метод Акерса со смещенным генетическим алгоритмом со случайным ключом для планирование рабочих мест» (Goncalves & Resende, 2014). Недавно я пытался реализовать упрощенные версии этих алгоритмов и получил некоторые интересные результаты, хотя реализация на чистом Python по-прежнему требует больших затрат времени. Вы можете найти их в этом хранилище.

Выводы

В этой статье были реализованы и решены два разных подхода смешанно-целочисленного программирования для задачи планирования рабочих мест (JSSP) с использованием библиотеки Python pyomo и решателя CBC с открытым исходным кодом. Дизъюнктивная модель оказалась лучшей альтернативой исходному JSSP, хотя более сложные модели реального мира могут иметь сходство с индексированной по времени формулировкой для включения дополнительных правил. Полный код, использованный в этих примерах, доступен для дальнейшего использования.