Реализация, решение и визуализация задачи коммивояжера с использованием Python

Решение задачи коммивояжера с использованием Python реализация и визуализация

Переведите модель оптимизации из математики в Python, оптимизируйте ее и визуализируйте решение, чтобы быстро получить обратную связь по ошибкам моделирования

👁️ Это статья № 3 серии о проекте “Интеллектуальная система поддержки принятия решений для туризма на Python”. Рекомендую вам ознакомиться с ней, чтобы получить общее представление о всем проекте. Если вас интересует только реализация модели TSP на Python, то вы все равно находитесь в правильном месте: эта статья самодостаточна, и я проведу вас через все этапы – установка зависимостей, анализ, код, интерпретация результатов и отладка модели.

Данная статья продолжает наше путешествие там, где мы остановились в sprint 2. Здесь мы берем математическую модель, сформулированную в предыдущей статье, и реализуем ее на Python с использованием Pyomo и следуя хорошим практикам. Затем модель оптимизируется, а решения визуализируются и анализируются. Из педагогических соображений мы обнаруживаем, что исходная формулировка модели неполна, поэтому я покажу, как можно получить, из первых принципов, ограничения, необходимые для исправления формулировки. Эти новые ограничения добавляются в модель Pyomo, и новые решения снова анализируются и проверяются.

Оглавление

1. Реализация модели на Python с использованием Pyomo

• 1.1. Установка зависимостей
• 1.2. Математика становится кодом

2. Решение и проверка модели

• 2.1. Решение модели
• 2.2. Визуализация результатов
• 2.3. Анализ результатов

3. Исправление формулировки

• 3.1. Мотивирующая идея
• 3.2. Выражение мотивирующей идеи в виде логических импликаций
• 3.3. Формулирование логических импликаций как линейных ограничений
• 3.4. Получение правильного значения для “большой M”

4. Реализация и проверка новой формулировки

• 4.1. Добавление новой формулировки в модель Pyomo
• 4.2. Визуализация обновленного решения модели

5. Заключение (для следующего спринта)

📖 Если вы не читали предыдущую статью, не беспокойтесь. Математическая формулировка также приведена (но не производная) здесь, с каждым компонентом модели рядом с ее кодовой реализацией. Если вы не понимаете, откуда что берется или что означает, пожалуйста, прочитайте статью “sprint 2”, и если вам нужно больше контекста, прочитайте статью “sprint 1”.

1. Реализация модели на Python с использованием Pyomo

Python используется, так как это основной язык в области науки о данных, а Pyomo является одной из лучших (открытых) библиотек, которая эффективно работает с моделями большого масштаба.

• Hands-On Deep Q-Learning’ (Практическое освоение глубокого Q-обучения)
• Длинные и короткие стороны пропорциональная актуальность для захвата семантики документа от начала до конца
• Ложный пророк Сравнение регрессионной модели с Пророком Меты

В этом разделе я рассмотрю каждый компонент модели, определенный в формулировке, и объясню, как он переводится в код Pyomo. Я постарался не оставить пробелов, но если у вас возникнут вопросы, оставьте их в комментариях.

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

1.1. Установка зависимостей

Для спешащих

Установите (или убедитесь, что у вас уже установлены) библиотеки pyomo, networkx и pandas, а также пакет glpk.

📝 Пакет glpk содержит решатель GLPK, который является (открытым исходным кодом) внешним решателем, который мы будем использовать для оптимизации создаваемых нами моделей. Pyomo используется для создания моделей проблемы и их передачи GLPK, который запускает алгоритмы, выполняющие процесс оптимизации. GLPK затем возвращает решения в виде атрибутов модели Pyomo, чтобы мы могли удобно использовать их без выхода из Python.

Рекомендуется устанавливать GLPK через conda, чтобы Pyomo мог легко найти решатель GLPK. Для одновременной установки всех зависимостей выполните следующую команду:

conda install -y -c conda-forge pyomo=6.5.0 pandas networkx glpk

Для организованных людей

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

name: ttp  # проблема путешествующего туристаchannels:  - conda-forgedependencies:  - python=3.9  - pyomo=6.5.0  - pandas  - networkx  - glpk  # внешний решатель, используемый для оптимизации моделей  - jupyterlab  # закомментируйте эту строку, если не будете использовать Jupyter Lab в качестве IDE

и сохраните его в файл YAML с именем environment.yml. Откройте Anaconda Prompt в том же месте и выполните команду

conda env create --file environment.yml

Через несколько минут будет создано окружение с установленными всеми зависимостями. Запустите conda activate ttp, чтобы “попасть внутрь” окружения, запустите Jupyter Lab (выполнив команду jupyter lab в терминале) и приступите к написанию кода!

1.2. Математика становится кодом

Сначала убедимся, что решатель GLPK доступен для Pyomo

### =====  Кодовый блок 3.1  ===== ###import pandas as pdimport pyomo.environ as pyoimport pyomo.versionimport networkx as nxsolver = pyo.SolverFactory("glpk")solver_available = solver.available(exception_flag=False)print(f"Solver '{solver.name}' доступен: {solver_available}")if solver_available:    print(f"Версия решателя: {solver.version()}")print("Версия pyomo:", pyomo.version.__version__)print("Версия networkx:", nx.__version__)

Solver 'glpk' доступен: TrueВерсия решателя: (5, 0, 0, 0)Версия pyomo: 6.5Версия networkx: 2.8.4

⛔ Если вы получили сообщение 'glpk' доступен: False, решатель не установлен должным образом. Попробуйте исправить проблемы, выполнив одно из следующих действий:

– повторно выполните установку осторожно

– выполните команду conda install -y -c conda-forge glpk в базовом (по умолчанию) окружении

– попробуйте установить другой решатель, который работает для вас

Затем прочитайте CSV-файл с данными о расстояниях

### =====  Блок кода 3.2  ===== ###path_data = (    "https://raw.githubusercontent.com/carlosjuribe/"    "traveling-tourist-problem/main/"    "Paris_sites_spherical_distance_matrix.csv")df_distances = pd.read_csv(path_data, index_col="site")df_distances

Теперь мы переходим к “этапу 4” рабочего процесса [Гибкая оптимизация исследований], отмеченному зеленым блоком ниже:

Задача состоит в том, чтобы взять математическую модель, созданную ранее, и реализовать ее в коде точно так, как она была определена математически.

👁️ Мы можем создавать столько объектов Python, сколько хотим, если это облегчит реализацию модели, но нам запрещено изменять основную модель во время написания кода. Это приведет к расхождению математической модели и компьютерной модели, что затруднит отладку модели позже.

Мы создаем пустую модель Pyomo, в которой будем хранить компоненты модели в виде атрибутов:

model = pyo.ConcreteModel("TSP")

1.2.1. Наборы

Для создания набора сайтов 𝕊 = {Лувр, Эйфелева башня, …, отель} мы извлекаем их имена из индекса таблицы данных и используем их для создания Pyomo Set с именем sites:

### =====  Блок кода 3.3  ===== ###list_of_sites = df_distances.index.tolist()model.sites = pyo.Set(initialize=list_of_sites,                       domain=pyo.Any,                       doc="набор всех сайтов для посещения (𝕊)")

Для создания производного набора

Мы сохраняем фильтр 𝑖 ≠ 𝑗 внутри правила конструирования (функции Python _rule_domain_arcs), и передаем это правило в ключевой аргумент filter при инициализации Set. Обратите внимание, что этот фильтр будет применяться к декартову произведению сайтов (𝕊 × 𝕊) и будет отсеивать те элементы декартового произведения, которые не удовлетворяют правилу.

### =====  Блок кода 3.4  ===== ###def _rule_domain_arcs(model, i, j):    """ Все возможные дуги, соединяющие сайты (𝔸) """    # создаем пару (i, j) только если сайты i и j разные    return (i, j) if i != j else None rule = _rule_domain_arcsmodel.valid_arcs = pyo.Set(    initialize=model.sites * model.sites,  # 𝕊 × 𝕊    filter=rule, doc=rule.__doc__)

1.2.2. Параметры

Параметр

𝐷ᵢⱼ ≔ Расстояние между сайтом 𝑖 и сайтом 𝑗

создается с помощью конструктора pyo.Param, который принимает в качестве первого (позиционного) аргумента домен 𝔸 (model.valid_arcs), который индексирует его, а в качестве ключевого аргумента initialize другое правило конструирования, _rule_distance_between_sites, которое вычисляется для каждой пары (𝑖, 𝑗) ∈ 𝔸. При каждом вычислении числовое значение расстояния извлекается из dataframe df_distances и связывается внутренне с парой (𝑖, 𝑗):

### =====  Блок кода 3.5  ===== ###def _rule_distance_between_sites(model, i, j):    """ Расстояние между сайтами i и j (𝐷𝑖𝑗) """    return df_distances.at[i, j]  # получить расстояние из dataframe-аргумент = _rule_distance_between_sitesmodel.distance_ij = pyo.Param(model.valid_arcs,                               initialize=rule,                               doc=rule.__doc__)

1.2.3. Решающие переменные

Поскольку 𝛿ᵢⱼ имеет ту же “индексную область” что и 𝐷ᵢⱼ, способ построения этой компоненты очень похож, за исключением того, что здесь не требуется правило создания.

model.delta_ij = pyo.Var(model.valid_arcs, within=pyo.Binary,                          doc="Перейти ли с сайта i на сайт j (𝛿𝑖𝑗)")

Первым позиционным аргументом pyo.Var зарезервирован для его индексного множества 𝔸, а “тип” переменной указывается с помощью именованного аргумента within. Под “типом переменной” я имею в виду диапазон значений, которые переменная может принимать. Здесь диапазон 𝛿ᵢⱼ равен только 0 и 1, поэтому он имеет бинарный тип. Математически, мы бы написали 𝛿ᵢⱼ ∈ {0, 1}, но вместо создания отдельных ограничений, чтобы указать это, мы можем указать это непосредственно в Pyomo, установив within=pyo.Binary в момент создания переменной.

1.2.4. Целевая функция

Для построения целевой функции мы можем «сохранить» выражение в функции, которая будет использоваться в качестве правила создания для цели. Эта функция принимает только один аргумент – модель, который используется для получения любого компонента модели, необходимого для построения выражения.

### =====  Блок кода 3.6  ===== ###def _rule_total_distance_traveled(model):    """ общая пройденная дистанция """    return pyo.summation(model.distance_ij, model.delta_ij)rule = _rule_total_distance_traveledmodel.obj_total_distance = pyo.Objective(rule=rule,                                          sense=pyo.minimize,                                          doc=rule.__doc__)

Обратите внимание на параллелизм между математическим выражением и оператором return функции. Мы указываем, что это проблема минимизации с помощью ключевого слова sense.

1.2.5. Ограничения

Если вы помните из предыдущей статьи, удобный способ обеспечить, чтобы каждый сайт был посещен только один раз, заключается в том, чтобы обеспечить, чтобы каждый сайт был “введен” только один раз и “вышел” только один раз одновременно.

Каждый сайт входит только один раз

def _rule_site_is_entered_once(model, j):    """ каждый сайт j должен быть посещен только с одного другого сайта """    return sum(model.delta_ij[i, j] for i in model.sites if i != j) == 1rule = _rule_site_is_entered_oncemodel.constr_each_site_is_entered_once = pyo.Constraint(                                          model.sites,                                           rule=rule,                                           doc=rule.__doc__)

“`html

Каждый сайт покидается только один раз

def _rule_site_is_exited_once(model, i):    """ каждый сайт i должен уходить только в один другой сайт """    return sum(model.delta_ij[i, j] for j in model.sites if j != i) == 1rule = _rule_site_is_exited_oncemodel.constr_each_site_is_exited_once = pyo.Constraint(                                          model.sites,                                           rule=rule,                                           doc=rule.__doc__)

1.2.6. Завершающая проверка модели

Реализация модели завершена. Чтобы увидеть, как она выглядит в целом, мы должны выполнить model.pprint() и немного пронавигировать, чтобы увидеть, пропустили ли мы какие-то объявления или допустили ошибки.

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

def print_model_info(model):    print(f"Name: {model.name}",           f"Num variables: {model.nvariables()}",          f"Num constraints: {model.nconstraints()}", sep="\n- ")print_model_info(model)#[Out]:# Name: TSP#  - Num variables: 72#  - Num constraints: 18

У модели меньше 100 ограничений или переменных, это проблема небольшого размера, и ее можно оптимизировать относительно быстро.

2. Решение и проверка модели

2.1. Решение модели

Следующий шаг в [AOR flowchart] – оптимизация модели и проверка решений:

res = solver.solve(model)  # оптимизация моделиprint(f"Найдено оптимальное решение: {pyo.check_optimal_termination(res)}")# [Out]: Найдено оптимальное решение: True

Хорошие новости! Солвер нашел оптимальное решение для этой проблемы! Давайте проверим его, чтобы знать, по какому маршруту идти, когда мы приедем в Париж!

Для очень быстрой проверки мы можем выполнить model.delta_ij.pprint(), который напечатает все (оптимальные) значения переменных 𝛿ᵢⱼ, которые могут быть либо 0, либо 1:

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

Для действительного понимания решения нам нужно его визуализировать.

2.2. Визуализация результатов

Изображение стоит тысячу записей

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

Функция extract_solution_as_arcs принимает решенную модель Pyomo и извлекает “выбранные дуги” из решения. Затем функция plot_arcs_as_graph хранит список активных дуг в объекте Graph для более простого анализа и строит этот график так, чтобы отель был единственным красным узлом в качестве ориентира. Наконец, функция plot_solution_as_graph вызывает вышеупомянутые две функции, чтобы отобразить решение заданной модели в виде графика.

“`

### =====  Блок кода 3.7  ===== ###
def extract_solution_as_arcs(model):
    """ Извлечение списка активных (выбранных) дуг из решенной модели,
    путем оставления только тех дуг, у которых бинарные переменные равны 1 """
    active_arcs = [(i, j)
                   for i, j in model.valid_arcs
                   if model.delta_ij[i, j].value == 1]
    return active_arcs

def plot_arcs_as_graph(tour_as_arcs):
    """ Взять список кортежей, представляющий дуги, преобразовать его
    в граф networkx и нарисовать его """
    G = nx.DiGraph()
    G.add_edges_from(tour_as_arcs)  # сохранение решения в виде графа
    node_colors = ['red' if node == 'hotel' else 'skyblue'
                   for node in G.nodes()]
    nx.draw(G, node_color=node_colors, with_labels=True, font_size=6,
            node_shape='o', arrowsize=5, style='solid')

def plot_solution_as_graph(model):
    """ Построение решения данной модели в виде графа """
    print(f"Общее расстояние: {model.obj_total_distance()}")
    active_arcs = extract_solution_as_arcs(model)
    plot_arcs_as_graph(active_arcs)

Теперь мы можем увидеть, как на самом деле выглядит решение:

plot_solution_as_graph(model)

Ну, это явно не то, что мы ожидали! Нет одного тура, проходящего через все достопримечательности и возвращающегося в отель. Верно, все достопримечательности посещаются, но только как часть небольших отдельных групп достопримечательностей. Технически, да, ограничения, которые мы указали, тщательно соблюдаются: каждая достопримечательность посещается только один раз и покидается только один раз, но общий результат – не один тур, как мы хотели, а группа подтурмов. Это означает, что предположение, сделанное в предыдущей статье, неверно, поэтому в модели что-то еще отсутствует, что бы закодировать требование “подтурма в решении не допускаются”.

2.3. Анализ результатов

Что пошло не так?

Когда решение модели не имеет смысла, есть только одно возможное объяснение: модель неверна¹.

Солвер дал нам действительно оптимальное решение модели, но мы дали ему модель, которая не соответствует проблеме, которую мы хотим решить. Задача сейчас заключается в том, чтобы выяснить, почему и где мы сделали ошибку. По размышлению, очевидным кандидатом является сомнительное предположение, сделанное в последних двух абзацах раздела “4.4. Создание ограничений” в предшествующей статье, где мы разработали математическую модель. Мы (неправильно, как мы теперь знаем) полагали, что образование одного тура предполагается естественным образом из двух ограничений, которые обеспечивают, что каждая достопримечательность посещается только один раз. Но, как мы только что визуализировали, это не так. Почему?

Коренная причина ошибки заключается в том, что я называю “невысказанным здравым смыслом”: знания, которые мы имеем о мире, которые настолько очевидны, что мы забываем указать их в модели

Мы знали, неявно, что телепортация невозможна при посещении достопримечательностей, но мы неявно не сказали модели об этом. Вот почему мы наблюдаем эти маленькие подтурмы, соединяющие некоторые из достопримечательностей, но не все. Модель “считает”, что перемещение с одной достопримечательности на другую допустимо, достаточно лишь находясь на достопримечательности, ее покидать один раз и входить один раз (см. рисунок 3.2 снова). Если мы видим подтурмы, это всего лишь потому, что мы сказали модели минимизировать расстояние маршрута, и случайно оказалось, что телепортация полезна для сокращения расстояния.

<!–Таким образом, нам необходимо предотвратить образование таких подтур, чтобы получить реалистическое решение. Мы должны разработать некоторые новые ограничения, которые “говорят модели”, что подтуры запрещены, или, что равносильно, что решение должно быть единственным туром. Будем придерживаться последнего и выведем, из первых принципов, один набор ограничений, который интуитивен и отлично справляется с задачей.

3. Исправление формулировки

Ссылаясь на [рабочий процесс Agile Operations Research], мы сейчас находимся в фазе переформулировки модели. Переформулировка модели может иметь целью ее улучшение или исправление. Нашей целью будет исправление.

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

3.1. Подстегивающая идея

У нас есть 𝑁 объектов для “прохождения”, включая отель. Поскольку мы начинаем с отеля, это означает, что осталось посетить 𝑁 − 1 объектов. Если мы отслеживаем “хронологический порядок” посещения этих объектов, такой, что первому пункту назначения (после отеля) присваивается номер 1, второму пункту назначения присваивается номер 2 и так далее, то последний пункт назначения перед возвращением в отель будет иметь номер 𝑁 − 1. Если мы назовем эти номера, используемые для отслеживания порядка посещений, “рангами”, то образующийся в туре шаблон будет таков, что ранг любого объекта (кроме отеля) всегда на 1 больше, чем ранг объекта, предшествующего ему в туре. Если мы могли бы сформулировать набор ограничений, которые навязывают такой шаблон любому допустимому решению, мы, говоря образно, вводим требование “хронологического порядка” в модель. И оказывается, что мы действительно можем это сделать.

3.2. Выражение подстегивающей идеи в виде логических импликаций

💡 Вот “шаблон”, который мы хотим, чтобы любое допустимое решение удовлетворяло:

Ранг любого объекта (кроме отеля) всегда должен быть на 1 больше ранга объекта, предшествующего ему в туре

Мы можем переформулировать этот шаблон как логическую импликацию, так: “ранг объекта 𝑗 должен быть на 1 больше ранга объекта 𝑖, если и только если 𝑗 посещается сразу после 𝑖, для всех дуг (𝑖, 𝑗), не включающих отель 𝐻”. Эту формулировку можно математически записать так:

где 𝑟ᵢ – новые переменные, которые нам нужно отслеживать (еще неизвестный) порядок посещения. Для отличия их от решающих переменных давайте назовем их “ранговыми” переменными. Правая сторона читается как “для всех 𝑖 и 𝑗, принадлежащих множеству всех объектов, исключая отель”. Для удобства введем новое множество 𝕊*, в котором хранятся все объекты за исключением отеля (обозначенного 𝐻):

что позволяет нам кратко определить переменные ранга следующим образом:

👁️ Очень важно, чтобы в отеле не было связанной переменной ранга, потому что он одновременно будет являться начальной и конечной точкой для любого тура, что нарушит шаблон «постепенного роста переменных ранга в туре». Таким образом, ранг каждого места всегда вынужден увеличиваться с каждым новым соединением, что гарантирует, что замкнутые петли запрещены, за исключением случая, когда петли замыкаются по единственному месту, у которого нет переменной ранга: отелю

Границы 𝑟ᵢ выводятся из его описания: ранг начинается с 1 и монотонно увеличивается до тех пор, пока не будут посещены все места в 𝕊*, достигая значения | 𝕊* | (размера набора мест, не являющихся отелем). Кроме того, мы позволяем им принимать любое положительное вещественное значение:

3.3. Формулировка логических следствий в виде линейных ограничений

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

Одним из способов сделать это является так называемый метод Большой-М, который заключается в объявлении ограничения таким образом, что, когда выполняется условие, которое вас интересует, ограничение применяется (становится активным), а когда условие, которое вас интересует, не выполняется, ограничение становится избыточным (становится неактивным). Техника называется «большой-М», потому что она использует постоянное значение 𝑀, которое достаточно велико, чтобы, если оно появится в ограничении, сделать его избыточным для каждого случая. Когда 𝑀 не появляется в ограничении, ограничение является «активным» в том смысле, что оно обеспечивает желаемое следствие.

Но что определяет, является ли ограничение «активным» или нет? Краткий ответ – значения самих решающих переменных, к которым применяется ограничение. Давайте посмотрим, как это работает.

Желаемое следствие состоит в том, чтобы иметь 𝑟ⱼ = 𝑟ᵢ + 1 только тогда, когда 𝛿ᵢⱼ = 1. Мы можем заменить 1 в выражении на 𝛿ᵢⱼ, что дает

Это отношение, которое мы хотим, чтобы выполнялось, когда 𝛿ᵢⱼ = 1, но не при 𝛿ᵢⱼ = 0. Чтобы «скорректировать» случай, когда 𝛿ᵢⱼ = 0, мы добавляем избыточный член, 𝑅𝑇, который служит для «деактивации ограничения» только тогда, когда 𝛿ᵢⱼ = 0. Поэтому этот избыточный член должен включать переменную 𝛿ᵢⱼ, так как он зависит от нее.

В этом контексте «деактивация ограничения» означает «сделать его избыточным», так как избыточное ограничение имеет такой же эффект, как несуществующее ограничение в модели.

Давайте посмотрим, как мы можем вывести выражение для RT. Выражение для 𝑅𝑇(𝛿ᵢⱼ) должно удовлетворять следующим условиям:

Для удовлетворения пункта (1) нам нужно, чтобы 𝑅𝑇(𝛿ᵢⱼ = 1) = 0, и, таким образом, выражение для 𝑅𝑇 должно содержать множитель (1 − 𝛿ᵢⱼ), так как оно становится равным 0, когда 𝛿ᵢⱼ = 1. Такая форма делает 𝑅𝑇 “исчезающим” при 𝛿ᵢⱼ = 1 или “сводит его к постоянной (давайте назовем ее М), когда 𝛿ᵢⱼ = 0. Таким образом, кандидатом на роль избыточного члена является

где 𝑀 должно быть определено на основе данных задачи (подробнее об этом позже).

Для удовлетворения пункта (2) для всех возможных 𝑖 и 𝑗 нам нужно превратить равенство в выражении (3.11) в неравенство (знак “=” становится “≥”) и найти константу 𝑀, достаточно большую (по абсолютной величине), чтобы, независимо от значений 𝑟ⱼ и 𝑟ᵢ, ограничение всегда выполнялось. Вот откуда берется “большое” в “большом М”.

Когда мы найдем такую достаточно большую константу 𝑀, наше ограничение “логической импликации” примет следующую форму

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

3.4. Вывод правильного значения «большого М»

Так как наша цель – иметь 𝑅𝑇(𝛿ᵢⱼ = 0) = 𝑀, мы можем вывести правильное значение для 𝑀, установив 𝛿ᵢⱼ = 0 в общем ограничении, указанном в Выражении (3.12):

Чтобы 𝑟ⱼ − 𝑟ᵢ ≥ 𝑀 выполнялось для всех сайтов, не являющихся отелем 𝑖, 𝑗, нам нужно, чтобы нижняя граница 𝑟ⱼ − 𝑟ᵢ тоже была больше 𝑀. Нижняя граница (LB) 𝑟ⱼ − 𝑟ᵢ – это минимальное значение, которое может принимать 𝑟ⱼ − 𝑟ᵢ и может быть получено посредством

где 𝑟ᵐⁱⁿ является минимально возможным рангом, а 𝑟ᵐᵃˣ – максимально возможным рангом. Поэтому, для неравенства (1) в Выражении (3.13), чтобы оно было истинным для всех рангов всех сайтов, должно выполняться также следующее неравенство:

Благодаря этому неравенству мы знаем минимальное значение, которое должно принимать 𝑀 для работы метода большого M, и оно равно

А каковы значения 𝑟ᵐⁱⁿ и 𝑟ᵐᵃˣ? В нашем соглашении, мы присвоили первому посещенному сайту ранг 1, который, конечно же, является минимальным рангом (т.е., 𝑟ᵐⁱⁿ = 1). Так как ранг увеличивается на 1 единицу при посещении каждого сайта, последний неотельный сайт в туре будет иметь максимальный ранг, равный количеству всех неотельных сайтов. Так как количество неотельных сайтов может варьироваться, нам нужно общее выражение. Если вы помните, мы определили множество 𝕊*, содержащее все неотельные сайты, поэтому число, которое нас интересует, это размер (т.е., количество элементов) множества 𝕊*, которое в математической нотации записывается как | 𝕊* |. Таким образом, мы выяснили, что 𝑟ᵐⁱⁿ = 1, а 𝑟ᵐᵃˣ = | 𝕊* |. Подставляя в Выражение (3.14), мы наконец получаем правильное значение для M:

Поскольку 𝕊* всегда будет содержать более 2 сайтов для посещения (иначе вообще не возникло бы проблемы выбора в первую очередь), значение «большого M» в этой модели всегда является отрицательным “большим” значением.

📝 Теоретические значения против вычисленных значений

Теоретически мы можем выбрать произвольно «более отрицательные» значения для 𝑀, чем полученное здесь значение, даже придумать огромные числа, чтобы быть уверенными, что они достаточно большие, чтобы избежать этого вычисления — но это нехорошая практика. Если 𝑀 становится слишком большим (по модулю), это может вызвать проблемы с производительностью в алгоритмах решателя или, в худшем случае, даже заставить решатель рассматривать недопустимые решения как допустимые. Поэтому рекомендуется извлекать из данных задачи жесткое, но достаточно большое значение для 𝑀.

Теперь, когда мы получили правильное значение для «большого M», мы сохраним его в новом параметре модели для более удобного его использования. Таким образом, набор ограничений на устранение подтурений готов, и его «полная форма»:

Чтобы сохранить все вещи в перспективе, обратите внимание, что это фактически “ограничительный эквивалент” первого логического следствия, которое мы сформулировали ранее в Выражении (3.6):

Поздравляю! У нас наконец-то есть набор ограничений, которые можно добавить к нашей модели. Их разработка была самой сложной частью. Теперь давайте проверим, что их добавление к модели действительно приводит к исчезновению подциклов.

4. Реализация и проверка новой формулировки

4.1. Расширение модели Pyomo новой формулировкой

Изменение и исправление модели потребовали дополнительного добавления нескольких множеств, параметров, переменных и ограничений. Давайте добавим эти новые компоненты модели Pyomo в той же последовательности, в которой они были в начальной фазе формулирования.

4.1.1. Множества и параметры

• Интересующие места, 𝕊*: множество всех мест за исключением отеля:

Объекты Pyomo Set имеют операции, совместимые с множествами Python, поэтому мы можем сразу сделать разницу между множеством Pyomo и множеством Python:

model.sites_except_hotel = pyo.Set(    initialize=model.sites - {'hotel'},     domain=model.sites,    doc="Интересующие места, т.е. все места, кроме отеля (𝕊*)")

• Большая M, новый параметр для ограничения устранения подциклов:

model.M = pyo.Param(initialize=1 - len(model.sites_except_hotel),                    doc="Большая M, чтобы некоторые ограничения стали избыточными")

4.1.2. Вспомогательные переменные

• Переменные ранга, rᵢ: для отслеживания порядка посещения мест:

model.rank_i = pyo.Var(    model.sites_except_hotel,  # i ∈ 𝕊* (индекс)    within=pyo.NonNegativeReals,  # rᵢ ∈ ℝ₊ (область)    bounds=(1, len(model.sites_except_hotel)),  # 1 ≤ rᵢ ≤ |𝕊*|    doc="Ранг каждого места для отслеживания порядка посещения")

В комментариях вы можете видеть, как красиво элементы полного математического определения переменных отображаются на аргументы функции объявления переменной Pyomo pyo.Var. Надеюсь, это поможет вам оценить ценность иметь ясно определенную математическую модель, прежде чем начинать создавать модель Pyomo. Реализация будет проходить естественным образом, и ошибки будут менее вероятными.

4.1.3. Ограничения

• Решение должно быть единственным туром, начинающимся и заканчивающимся в отеле:

def _rule_path_is_single_tour(model, i, j):    """ Для каждой пары непроживающих отелей (i, j), если сайт j посещается из сайта i, ранг j должен быть строго больше ранга i. """    if i == j:  # если сайты совпадают, пропустить создание ограничения        return pyo.Constraint.Skip        r_i = model.rank_i[i]    r_j = model.rank_i[j]    delta_ij = model.delta_ij[i, j]    return r_j >= r_i + delta_ij + (1 - delta_ij) * model.M# декартово произведение непроживающих отелей для индексации ограниченияnon_hotel_site_pairs = model.sites_except_hotel * model.sites_except_hotelrule = _rule_path_is_single_tourmodel.constr_path_is_single_tour = pyo.Constraint(    non_hotel_site_pairs,    rule=rule,     doc=rule.__doc__)

Модель Pyomo была обновлена. Насколько она увеличилась после добавления ограничений на устранение подциклов?

print_model_info(model)# [Out]:# Название: TSP#  - Количество переменных: 80#  - Количество ограничений: 74

Мы перешли от 72 до 80 переменных и от 18 до 74 ограничений. Очевидно, что эта формулировка требует большего количества ограничений, чем переменных, так как она увеличила количество ограничений в четыре раза. Это “цена”, которую мы платим за то, чтобы модели были более “реалистичными”, поскольку реализм обычно ограничивает количество допустимых решений при неизменных данных.

Как всегда, мы можем изучить структуру модели с помощьюmodel.pprint(). Хотя это “печать” быстро теряет свою ценность по мере увеличения количества компонентов модели, она все же дает нам быстрый обзор того, из чего состоит модель и насколько она большая.

4.2. Графическое представление решения обновленной модели

Давайте решим обновленную модель и построим новое решение. Следим за успехом.

res = solver.solve(model)  # оптимизация моделиsolution_found = pyo.check_optimal_termination(res)print(f"Оптимальное решение найдено: {solution_found}")# [Out]: Оптимальное решение найдено: Trueif solution_found:    plot_solution_as_graph(model)

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

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

👁️ Рисунок графа – это не тур на карте

Обратите внимание, что этот рисунок – это только одно возможное изображение графа, а не географическая траектория. Круги и связи, которые вы видите, не соответствуют реальному пути в географическом пространстве (как они могут, если мы не использовали никакую географическую информацию при их создании?). Вы можете проверить это сами, запустив plot_solution_as_graph(model) несколько раз: каждый раз узлы будут занимать разные позиции. Графы – это абстрактные математические структуры, которые соединяют “точки” через “связи”, представляя отношения любого вида между сущностями любого вида. Мы использовали граф здесь для изучения правильности решения, а не для визуализации реального тура в Париже. Мы делаем это в [этой статье].

5. Заключение (или планирование следующего спринта)

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

💡 Эволюция решений, по одному спринту за раз

Этот POC пока не решает нашу изначальную сложную проблему туризма, но он решает минимальную ценную проблему, которую мы предложили в качестве первого опорного камня на пути к более сложному решению. Таким образом, он доставляет нам доказуемое (то есть основанное на доказательствах) приближение к ценному решению более сложной проблемы. Пользуясь минимальным рабочим примером, мы можем лучше оценить, что нужно сделать, чтобы двигаться в нужном направлении, и что можно временно упростить до более зрелой версии решения. Постепенно развиваясь и имея всегда полезное решение, мы создадим все более ценную систему, пока не будем удовлетворены. В этом суть “гибкой дороги” к эффективным решениям.

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

5.1. Что дальше

В следующей статье (спринте 4) мы будем работать над классом, который автоматически генерирует матрицу расстояний из любого списка сайтов. Эта функциональность, в сочетании с только что построенной моделью, позволит нам, среди прочего, быстро решать множество моделей с разными входными данными и сравнивать их. Кроме того, обобщение решения таким образом упростит нашу жизнь позже, когда мы проведем небольшой анализ чувствительности и сценариев в будущих спринтах. Кроме того, поскольку мы повысим этот прототип до статуса «MVP», мы начнем использовать код, основанный на объектах, для хорошей организации и готовности к расширению. Не теряйте время и переходите прямо сейчас!

Примечания

1. На самом деле, есть еще одна причина неправильного результата: данные, на которых питается модель, могут быть неправильными, а не только формулировка модели. Но, в некотором смысле, если вы подумаете о “модели” как о “экземпляре модели”, то есть о конкретной модели с конкретными данными, то модель, конечно же, будет неправильной, если данные неправильные, что я и имел в виду.

Спасибо за чтение, и увидимся в следующей статье! 📈😊

Не стесняйтесь следовать за мной, задавать вопросы, давать мне обратную связь в комментариях или связаться со мной на LinkedIn.