07_Column_Generation_Algorithm.ipynb

{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "provenance": [],
   "toc_visible": true,
   "authorship_tag": "ABX9TyNpNgVw8+/QQqgLY5+jR04G",
   "include_colab_link": true
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "view-in-github",
    "colab_type": "text"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/Lozkins/mos/blob/master/07_Column_Generation_Algorithm.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Алгоритм генерации столбцов (Column Generation)\n",
    "\n",
    "*Aleksejs Lozkins,*\n",
    "[*la-optimization*](https://la-optimization.ru/)\n",
    "\n",
    "![Схема алгоритма генерации столбцов](https://drive.google.com/uc?export=view&id=1Ice5Bn3mQLHRivBFeHgZUaDGeZabUR3P)\n",
    "\n",
    "Генерация столбцов - подход к решению задач смешанного линейного программирования (MIP) с большим кол-вом переменных или столбцов.\n",
    "\n",
    "В статье представил теоретическую предпосылку, схему алгоритма и python реализацию подхода. Практическая часть состоит из решения двух задач: задача планирования расписания смен экипажей и задача раскроя.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ],
   "metadata": {
    "id": "YYZBQtZOMEbg"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Математическая база\n",
    "\n",
    "Погружение в алгоритм генерации столбцов начнем с симплекс метода. Описание всего [симплекс метода](https://habr.com/ru/articles/474286/) оставлю за рамками этой статьи, здесь рассмотрим только определенный этап симплекс метода - определение ведущего столбца для добавления в базис (pricing problem).\n",
    "\n",
    "\n",
    "\n"
   ],
   "metadata": {
    "id": "wo27tfgV9252"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Вектор приведенной стоимости (reduced cost)\n",
    "\n",
    "Пусть имеем задачу линейного программирования в канонической форме. Для определенности, рассмотрим задачу минимизации:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\min cx  \\\\\n",
    "Ax = b & & (1)\\\\\n",
    "x \\in \\mathbb{R}^{m+} &  & \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "где $c$ вектор коэффициентов целевой функции размерности $n$, $x$ вектор вещественных переменных размерности $m$, $A$ и $b$ матрица коэффициентов при ограничениях и вектор свободных членов, соответственно.\n",
    "\n",
    "**Определения:** Пусть $(B \\quad N)$ разбиение столбцов матрицы $A$ таким образом, что $B$ имеет размерность $m \\times m$ и определитель матрицы отличен от нуля. Тогда $x = (x_B \\quad x_N)$ и $c=(c_B \\quad c_N)$ соответствующее разбиение для $x$ и $c$. Будем называть допустимое решение $x = (x_B, x_N)$ - **базисным (основным)**, если $x_N = 0$. Переменные $x_B$ называются **базисными (основными) переменными**, а $x_N$ - **свободными**.\n",
    "\n",
    "В процессе решения задачи линейного программирования, на каждой итерации симплекс метода происходит оценка (pricing) свободных переменных и выбор нового кандидата для добавления в базис (ведущий столбец). Этот шаг играет ключевую роль в алгоритме генерации столбцов, так и в самом симплекс методе.\n",
    "\n",
    "$$Ax = (B \\quad N)\n",
    "\\begin{pmatrix}\n",
    "  x_B \\\\\n",
    "  0\n",
    "\\end{pmatrix} = b, \\text{ тогда } x_B = B^{-1}b.$$\n",
    "\n",
    "Если $x$ является допустимым базисным решением, тогда значение целевой функции (стоимость) будет:\n",
    "$$cx = (c_B \\quad c_N)\n",
    "\\begin{pmatrix}\n",
    "  x_B \\\\\n",
    "  0\n",
    "\\end{pmatrix} = c_BB^{-1}b.$$\n",
    "\n",
    "Куда уйти с текущей угловой точки (текущего базисного решения) чтобы улучшить целевую функцию? Предположим, что свободные переменные теперь могут принимать не нулевые значения, обозначим такой вектор $x' = (x'_B \\quad x'_N)$, тогда\n",
    "\n",
    "$$Ax' = (B \\quad N)\n",
    "\\begin{pmatrix}\n",
    "  x'_B \\\\\n",
    "  x'_N\n",
    "\\end{pmatrix} = b.$$\n",
    "\n",
    "До множим правую и левую часть на $B^{-1}$\n",
    "$$B^{-1}Ax' = (I \\quad B^{-1}N)\n",
    "\\begin{pmatrix}\n",
    "  x'_B \\\\\n",
    "  x'_N\n",
    "\\end{pmatrix} = B^{-1}b,$$\n",
    "\n",
    "$$x'_B + B^{-1}Nx'_N = B^{-1}b,$$\n",
    "\n",
    "$$x'_B = B^{-1}b - B^{-1}Nx'_N.$$\n",
    "\n",
    "Подставим значения вектора $x'$ в целевую функцию:\n",
    "\n",
    "$$cx'= c_Bx'_B+c_Nx'_N = c_B(B^{-1}b - B^{-1}Nx'_N) + c_Nx'_N = c_BB^{-1}b + (c_N - c_BB^{-1}N)x'_N.$$\n",
    "\n",
    "Как изменилась целевая функция при допущении, что свободные переменные не равны нулю? Запишем изменение целевой функции как $cx' - cx$ что будет соответствовать\n",
    "\n",
    "$$c_BB^{-1}b + (c_N - c_BB^{-1}N)x'_N - c_BB^{-1}b= (c_N-c_BB^{-1}N)x'_N.$$\n",
    "\n",
    "Мы можем видеть что если $x_N$ увеличивается ($x_N \\ge 0$), функция  затрат растет или убывает в зависимости от знака $$\\bar{c} = c_N-c_BB^{-1}N.$$\n",
    "\n",
    "Вектор $\\bar{c}$ называют приведенной стоимостью (reduced cost). Если существует отрицательная компонента в векторе $\\bar{c}$, то соответствующая ей свободная переменная может уменьшить значение целевой функции - добавляем её в базис.\n",
    "\n",
    "В симплекс методе в базис добавляется одна переменная с наиболее отрицательной компонентой $\\bar{c}$, т.к. она доставляет наискорейшее уменьшение целевой функции. Если в векторе $\\bar{c}$ нет отрицательных компонент, то мы не можем улучшить целевую функцию, следовательно нашли оптимальное решение.\n",
    "\n",
    "На этом закончим погружение в симплекс метод и перейдем к алгоритму генерации столбцов."
   ],
   "metadata": {
    "id": "_095xOjEtNGs"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Усеченная задача и генерация столбцов\n",
    "\n",
    "Рассмотрим задачу смешанного линейного программирования (MIP) в стандартной форме следующего вида:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\min \\sum_{j \\in J} c_{j} x_{j} \\\\\n",
    "\\sum_{j \\in J} \\textbf{a}_j x_j \\ge \\textbf{b} & & (2)\\\\\n",
    "x_{j} \\in \\mathbb{R}^+|\\mathbb{N}^+ & \\quad \\forall j \\in J & \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Будем называть задачу (2) основной (координирующей) сокращенно MP (Master Problem). Размерность задачи $|J|=n$ столбцов (переменных) и $|\\textbf{b}|=m$ строк (ограничений). Предполагается, что кол-во переменных в задаче $n$ экспоненциально больше кол-ва ограничений $m$ ($n \\gg m$). Здесь $c_j$ коэффициенты при переменных в целевой функции, $x_j$ решающие переменные, матрица коэффициентов при ограничениях $A=\\{\\textbf{a}_1, \\dots, \\textbf{a}_n\\}$ и вектор свободных членов $\\textbf{b}$.\n",
    "\n",
    "**Расслабленной задачей** к задаче (2) будем называть задачу с ослабленным ограничением на целочисленность переменных, т.е. в такой задаче все переменные являются целочисленными. **Двойственная задача** к прямой расслабленной задаче (2) запишется следующим образом:\n",
    "\n",
    "\\begin{align}\n",
    "\\max \\textbf{u}^\\intercal \\textbf{b} \\\\\n",
    "\\textbf{u}^\\intercal A \\le \\textbf{c} & & (3)\\\\\n",
    "\\textbf{u} \\in \\mathbb{R}^m & \\\\\n",
    "\\end{align}\n",
    "\n",
    "здесь $\\textbf{u}$ вектор двойственных переменных размерности $m$.\n",
    "\n"
   ],
   "metadata": {
    "id": "60cIJqZoLoH_"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Приведенная стоимость через двойственные переменные\n",
    "\n",
    "Когда обсуждали вектор приведенной стоимости уже получали выражение $x_B = B^{-1}b$ для базисных переменных. Для этого базисного решения можем сформулировать соответствующую двойственную задачу: $B^\\intercal u = c_B^\\intercal$. Применим операцию транспонирования к правой и левой части равенства, до множим на $B^{-1}$ справа: $u^\\intercal = c_B B^{-1}$. Подставим полученный результат в формулу расчета приведенной стоимости (reduced cost), получим:\n",
    "$$\\bar{c} = c_N-c_BB^{-1}N=c_N-u^\\intercal N.$$\n",
    "\n",
    "Пусть $J_N = \\{j_1, \\dots, j_{n-m}\\}$ множество индексов столбцов матрицы $A$, которые соответствуют индексам свободных переменных и составляют матрицу $N$. По симплекс методу, добавляем в базис переменную, которая имеет самую отрицательную компоненту $\\bar{c}_j$ из небазисных переменных, т.е.\n",
    "\n",
    "$$j_B = \\text{arg} \\min \\{\\bar{c_j} = c_j-u^\\intercal a_j|j \\in J_N\\}$$\n",
    "\n",
    "Вполне правомерным можно назвать использование вместо $J_N$ множество всех столбцов $J$. Это оправдывается тем, что reduced cost для базисных переменных равен нулю: $c_B-c_BB^{-1}B = 0$. Отвязываемся от деления матрицы $A$ на базисную и свободную часть:\n",
    "\n",
    "$$j_B = \\text{arg} \\min \\{\\bar{c_j} = c_j-u^\\intercal a_j|j \\in J\\}.$$\n",
    "\n",
    "Ранее делали предположение что кол-во строк в матрице $A$ значительно меньше кол-ва столбцов. Громоздкость матрицы усложняет вычисления. На практике достаточно решить значительно меньшую задачу для получения того же оптимального решения."
   ],
   "metadata": {
    "id": "5qLIdovetgB1"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "\n",
    "### Идея алгоритма генерации столбцов\n",
    "\n",
    "Рассмотрим некоторую подзадачу (Restricted Master Problem, RMP) исходной задачи $J' \\subset J$ при ($|J'| \\ge m$), которые отличаются между собой только количеством столбцов в матрице $A$. Воспользуемся идеей симплекс метода добавления ведущего столбца в базис (pricing problem), но будем добавлять столбец в RMP. Таким образом, если для всех оставшихся столбцов приведенная стоимость (reduced cost) $\\bar{c_j} \\ge 0, \\forall j \\in J \\text{\\\\} J'$, то значение целевой функции RMP улучшить невозможно за счет добавления новых столбцов. Это означает, что в этом случае оптимальное решение подзадачи RMP и полной задачи MP совпадают.\n",
    "\n",
    "**Комментарий:** Для расчета $\\bar{c_j}$ в RMP используем $\\bar{u}$ вектор оптимальных значений двойственных переменных расслабленной задачи RMP.\n",
    "\n"
   ],
   "metadata": {
    "id": "zkY2I4WGt7Su"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Явное и неявное задание столбцов\n",
    "\n",
    "Обратим внимание на множество столбцов $J$, которое ранее предполагали получать в явном виде - определенный конечный набор столбцов. В некоторых задачах множество столбцов может быть безумно большим или бесконечным (например, столбцы могут описывать все возможные маршруты, паттерны, подмножества или перестановки). В этом случае можно попробовать прибегнуть к неявному способу задания множества столбцов через набор ограничений $a_j \\in \\mathcal{A}$. Тогда поиск столбца будет сводится к оптимизационной задаче:\n",
    "\n",
    "$$\\min \\{c(a) - \\bar{u}^\\intercal a| a \\in \\mathcal{A} \\}.$$\n",
    "\n",
    "Базовым примером применения такого множества ограничений $\\mathcal{A}$ является задача оптимального раскроя, где $a_j$ представляет собой паттерн разреза, а ограничения $\\mathcal{A}$ с целевой функцией приведенной стоимости представляют собой задачу о рюкзаке. На каждой итерации алгоритма генерации столбцов решаем задачу о рюкзаке для построения паттерна. Затем добавляем этот паттерн в RMP."
   ],
   "metadata": {
    "id": "RkUuuI84vKcJ"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Схема алгоритма\n",
    "\n",
    "Зафиксируем основные этапы алгоритма генерации столбцов. Приведу для случая явного задания столбцов:\n",
    "\n",
    "1. Фиксируем начальный набор столбцов $J'$ для инициализации алгоритма;\n",
    "2. Формулируем расслабленную задачу RMP. В исходной матрице $A$ и векторе $c$ оставляем столбцы соответствующие зафиксированному набору столбцов $J'$;\n",
    "3. Строим двойственную задачу к RMP и решаем ее;\n",
    "4. Pricing problem. Рассчитываем приведенную стоимость (reduced cost) и выбираем кандидата для добавления в RMP. Критерий: $j_B = \\text{arg} \\min \\{\\bar{c_j} = c_j-u^\\intercal \\textbf{a}_j|j \\in J\\}.$\n",
    "5. Проверяем условия:\n",
    "\n",
    "  5.1. Если приведенная стоимость (reduced cost) для кандидата из п.4 - отрицательная, то добавляем его в $J'$ и переходим в п.2.;\n",
    "\n",
    "  5.2. Если приведенная стоимость для $\\forall j \\in J$ неотрицательная, то останавливаем алгоритм (stop). Оптимальное решение RMP будет совпадать с оптимальным решением полной задачи MP;\n",
    "\n",
    "7. Решаем прямую задачу RMP (на множестве столбцов $J'$) с учетом целочисленности переменных. Полученное решение является оптимальным.\n",
    "\n",
    "![Схема алгоритма генерации столбцов](https://drive.google.com/uc?export=view&id=1t1PnVlQhmkjgkpsNREkQQJuDcIk49PNu)\n",
    "\n",
    "Дополнительным критерием остановки алгоритма рекомендую добавлять условие $J = J'$. В случае, когда все столбцы добавлены в RMP мы получаем полную исходную задачу. Этот сценарий должен наводить на сомнение в корректности реализации алгоритма или его применимости к решению задачи в общем."
   ],
   "metadata": {
    "id": "dVyayF5Av9Oc"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Практическая реализация\n",
    "\n",
    "В технической части статьи рассмотрим 2 случая использования алгоритма генерации столбцов. В первом случае, множество столбцов задается в явном виде, для реалистичности примера будем решать задачу планирования сменных графиков экипажей поездов в пассажира-перевозках. Второй случай, столбцы заданы в неявном виде через набор ограничений, речь пойдет о задаче раскроя.\n",
    "\n",
    "По стеку технологий: python, ortools, солверы SCIP и GLOP, для работы с данными - numpy."
   ],
   "metadata": {
    "id": "w6K1Lo_URUOF"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Задача планирования смен экипажей в пассажира-перевозках\n",
    "\n",
    "Метод генерации столбцов широко используется для решения задачи планирования смен экипажей ([Crew Scheduling Problem](https://en.wikipedia.org/wiki/Crew_scheduling), CSP), особенно актуально в сфере ЖД пассажира-перевозок. Конечно, это далеко не единственная задача, где уместно применять метод генерации столбцов. Здесь с большим успехом можно также рассмотреть задачи о назначениях, задачи покрытия, задачи маршрутизации или др.\n",
    "\n",
    "Задача заключается в планировании графиков работы экипажей на заранее известном наборе рейсов. На входе имеем план работ в виде участков движения поездов откуда-куда, время начала движения по участку, требования к квалификации экипажа или его составу и т.д. Цель задачи найти минимальный набор экипажей, который обеспечит выполнение всех запланированных рейсов.\n",
    "\n",
    "Одним из вариантов моделирования задачи является ее представление в виде [задачи покрытия](ru.wikipedia.org/wiki/Задача_о_покрытии_множества). На вход модели подаются готовые варианты графиков работы экипажей (допустимые), а результатом является выбор минимального кол-ва таких графиков (соответствует кол-ву экипажей), которые покрывают все рейсы.\n",
    "\n",
    "![Пример сменных графиков экипажей](https://drive.google.com/uc?export=view&id=1lZAaOp7Qg9qaYOEeE4YIBmxahqinMb7m)\n",
    "\n",
    "Кол-во возможных сменных графиков $n$ значительно больше кол-ва рейсов $m$, потому что графики состоят из различных допустимых наборов рейсов."
   ],
   "metadata": {
    "id": "kDtvO0mG5Xyk"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Математическая постановка\n",
    "$i \\in I$ - множество запланированных рейсов;\n",
    "\n",
    "$j \\in J$ - множество допустимых сменных графиков, каждый из которых может выполнить один экипаж;\n",
    "\n",
    "$\\textbf{a}_j$ - булевый $|I|$-мерный вектор констант, где $i$-ая компонента равна 1, если соответствующий рейс находится в сменном графике $j$. В противном случае, 0. Тот самый сменный график в векторной форме;\n",
    "\n",
    "$x_j$ - бинарная переменная, принимает значение 1, если сменный график $j$ выбран к исполнению, 0 в противном случае;\n",
    "\n",
    "$\\sum_{j \\in J} a_{ij}x_j \\ge 1, \\quad \\forall i \\in I$ - каждый рейс должен быть покрыт как минимум одним экипажем;\n",
    "\n",
    "$\\min \\sum_{j \\in J} x_j$ - минимизируем кол-во сменных графиков для обеспечения выполнения всех запланированных рейсов (кол-во сменных графиков = кол-во экипажей)."
   ],
   "metadata": {
    "id": "OFXwXuoO5UX6"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Python реализация задачи CSP"
   ],
   "metadata": {
    "id": "6cdVyUa8KqBk"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Устанавливаем ortools\n",
    "%pip install ortools"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "0yvkLV35nFeX",
    "outputId": "fda1afe8-55ed-4fbc-c555-cfe5a1fd69cc"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "Данные для задачи CSP сгенерируем случайным образом. По умолчанию стоит маленький кейс. Чтобы получить расширенную задачу, то необходимо указать `is_small = False`. Поиск оптимального решения полной задачи на расширенных данных займет приличное время."
   ],
   "metadata": {
    "id": "PWlEPmw-mmTh"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Импортируем необходимые библиотеки\n",
    "from ortools.linear_solver import pywraplp\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import random\n",
    "import time\n",
    "\n",
    "# Подготовка искусственного сценария данных\n",
    "is_small = True  # Переключатель размерности данных\n",
    "trips_cnt = 30 if is_small else 200 # Кол-во рейсов\n",
    "duties_cnt = 200 if is_small else 4000 # Кол-во сменных графиков\n",
    "seed = 3\n",
    "\n",
    "rng = np.random.default_rng(seed=seed)\n",
    "df = rng.choice(100, size=(trips_cnt, duties_cnt))\n",
    "\n",
    "columns = [f\"duty_{i + 1}\" for i in range(duties_cnt)]\n",
    "indexes = [f\"trip_{i}\" for i in range(trips_cnt)]\n",
    "\n",
    "# Вектор весов целевой функции.\n",
    "c = np.array([1 for i in range(duties_cnt)])\n",
    "\n",
    "# Инициализируем матрицу А\n",
    "A = pd.DataFrame(df, columns=columns, index=indexes)\n",
    "A = A / 100 - 0.3 if is_small else A / 100 - 0.48  # Прореживаем сменные графики\n",
    "A = A.round().astype(int)\n",
    "A = A.to_numpy()\n",
    "\n",
    "# Проверим, что для каждого рейса есть хотя бы один сменный график, который его покрывает\n",
    "sat = A.sum(axis=1)\n",
    "sat[sat[0] == 0]"
   ],
   "metadata": {
    "id": "5rp6G3Z6uEyZ",
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "outputId": "56d26b2d-7ebe-489d-f576-4e1dad6e36ca"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Двойственная задача к расслабленной задаче RMP\n",
    "def optimize_dual(A, c, indexes, selected_x_keys=None):\n",
    "    \"\"\"\n",
    "      A - исходная матрица задачи\n",
    "      с - вектор коэффициентов целевой функции\n",
    "      indexes - список рейсов к выполнению\n",
    "      selected_x_keys - список индексов сменных графиков J'\n",
    "    \"\"\"\n",
    "    model = pywraplp.Solver.CreateSolver('GLOP')  # Солвер LP\n",
    "    A_T = A.T\n",
    "\n",
    "    # Инициализация вещественных переменных\n",
    "    y_vars = []\n",
    "    for i, name in enumerate(indexes):\n",
    "        y_vars.append(model.NumVar(0, 2, name=f\"{name}\"))  # Remark: ограничение сверху будет пояснено далее\n",
    "\n",
    "    # Добавляем ограничения в модель\n",
    "    selected_x_keys = selected_x_keys or list(range(len(A_T)))\n",
    "    for i in selected_x_keys:  # Строим ограничения только для сменных графиков из J'\n",
    "        model.Add(A_T[i] @ y_vars <= c[i])\n",
    "\n",
    "    # Целевая функция двойственной задачи\n",
    "    model.Maximize(sum(y_vars))\n",
    "    model_str = model.ExportModelAsLpFormat(False)  # Задача в формате .lp\n",
    "\n",
    "    # Поиск решения задачи\n",
    "    status = model.Solve()\n",
    "\n",
    "    print(\"Найдено оптимальное решение: \", status == pywraplp.Solver.OPTIMAL)\n",
    "    print('Целевая функция = %f' % model.Objective().Value())\n",
    "\n",
    "    # Извлекаем и выводим решение задачи\n",
    "    solution = {indexes[i]: var.solution_value()\n",
    "                for i, var in enumerate(y_vars)}\n",
    "    print(solution)\n",
    "\n",
    "    # Вектор значений двойственных переменных (для расчета приведенной стоимости)\n",
    "    solution_vect = np.array([var.solution_value() for var in y_vars])\n",
    "    return solution_vect, model_str"
   ],
   "metadata": {
    "id": "iVY58cfmOCxA"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "%%time\n",
    "result, dual_model_str = optimize_dual(A, c, indexes)"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "4Rb7rdluOvjM",
    "outputId": "0fa25dac-8c83-4388-bbec-7c781951404a"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Прямая MIP задача\n",
    "def optimize_primal_MIP(A, c, columns, selected_x_keys=None):\n",
    "    \"\"\"\n",
    "      A - исходная матрица задачи\n",
    "      с - вектор коэффициентов целевой функции\n",
    "      columns - список сменных графиков\n",
    "      selected_x_keys - список индексов сменных графиков J'\n",
    "    \"\"\"\n",
    "    selected_x_keys = selected_x_keys or list(range(A.shape[1]))\n",
    "\n",
    "    model = pywraplp.Solver.CreateSolver('SCIP')\n",
    "\n",
    "    # Инициализация бинарных переменных выбора сменного графика:\n",
    "    # принимает значение 1 если сменный график выбран к исполнению, 0 в противном случае.\n",
    "    x_vars = []\n",
    "    for i in selected_x_keys:\n",
    "        name = columns[i]\n",
    "        x_vars.append(model.BoolVar(name=f\"{name}\"))\n",
    "\n",
    "    # Добавляем ограничения в модель\n",
    "    # Ограничение покрытия формируем исходя из доступного набора сменных графиков J'\n",
    "    A_tmp = A[:, selected_x_keys]\n",
    "    for i in range(len(A_tmp)):\n",
    "        constr = model.Add(A_tmp[i] @ x_vars >= 1)\n",
    "\n",
    "    # Целевая функция\n",
    "    model.Minimize(c[selected_x_keys] @ x_vars)\n",
    "    model_str = model.ExportModelAsLpFormat(False) # Задача в формате .lp\n",
    "\n",
    "    # Поиск решения задачи\n",
    "    model.EnableOutput()\n",
    "    status = model.Solve()\n",
    "\n",
    "    print(\"Найдено оптимальное решение: \", status == pywraplp.Solver.OPTIMAL)\n",
    "    print('Целевая функция = %f' % model.Objective().Value())\n",
    "\n",
    "    # Извлекаем и выводим решение задачи\n",
    "    solution = {columns[selected_x_keys[i]]: var.solution_value()\n",
    "                for i, var in enumerate(x_vars)\n",
    "               if var.solution_value() }\n",
    "    print(solution)\n",
    "\n",
    "    # Вектор значений решающих переменных\n",
    "    solution_vect = np.array([var.solution_value() for var in x_vars])\n",
    "    return solution_vect, model_str"
   ],
   "metadata": {
    "id": "o15LU9PvOROQ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "%%time\n",
    "res, model_str = optimize_primal_MIP(A, c, columns)"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "i3R_y0kEOoWc",
    "outputId": "61c86edb-efca-43cf-96a4-99b949761f5a"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "Найдено оптимальое решение:  True\n",
    "Целевая функция = 43.000000\n",
    "{'duty_3175': 1.0, 'duty_3795': 1.0, 'duty_2441': 1.0, 'duty_3032': 1.0, 'duty_86': 1.0, 'duty_3929': 1.0, 'duty_928': 1.0, 'duty_1069': 1.0, 'duty_3820': 1.0, 'duty_373': 1.0, 'duty_2810': 1.0, 'duty_54': 1.0, 'duty_813': 1.0, 'duty_1837': 1.0, 'duty_2044': 1.0, 'duty_45': 1.0, 'duty_440': 1.0, 'duty_1337': 1.0, 'duty_1338': 1.0, 'duty_1746': 1.0, 'duty_2463': 1.0, 'duty_2801': 1.0, 'duty_2869': 1.0, 'duty_58': 1.0, 'duty_550': 1.0, 'duty_3973': 1.0, 'duty_1922': 1.0, 'duty_2134': 1.0, 'duty_1148': 1.0, 'duty_233': 1.0, 'duty_2642': 1.0, 'duty_2673': 1.0, 'duty_224': 1.0, 'duty_2499': 1.0, 'duty_2623': 1.0, 'duty_2914': 1.0, 'duty_692': 1.0, 'duty_2100': 1.0, 'duty_1960': 1.0, 'duty_1020': 1.0, 'duty_1521': 1.0, 'duty_3174': 1.0, 'duty_3723': 1.0}\n",
    "CPU times: user 3min 21s, sys: 146 ms, total: 3min 21s\n",
    "Wall time: 3min 22s"
   ],
   "metadata": {
    "id": "bOW0DTUcrLo5"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "stop_bound = -1e-5  # Для остановки будем использовать число близкое к нулю\n",
    "def column_generation_by_dual(A, c, indexes, columns):\n",
    "    # indexes - dual variables\n",
    "    # columns - primal variables\n",
    "\n",
    "    # Инициализируем начальное множество J'\n",
    "    selected_keys = [0]\n",
    "    print(\"Инициализируем J': \", selected_keys, '\\n')\n",
    "\n",
    "    min_val = -1  # инициализация минимального значения приведенной стоимости\n",
    "    iter_count = 0  # счетчик итераций\n",
    "    # Критерии остановки алгоритма:\n",
    "    # 1. В J/J' не осталось столбцов, которые уменьшают целевую функцию\n",
    "    # 2. |J'| == |J|\n",
    "    while min_val < stop_bound and len(selected_keys) < len(columns):\n",
    "        iter_count += 1\n",
    "        print('Итерация: ', iter_count)\n",
    "        dual_solution, model_str = optimize_dual(A, c, indexes, selected_keys)\n",
    "\n",
    "        reduced_cost = 1 - A.T @ dual_solution\n",
    "        print(\"Приведенная стоимость: \", reduced_cost)\n",
    "\n",
    "        best_key = reduced_cost.argmin()\n",
    "        min_val = reduced_cost[best_key]\n",
    "        if min_val < stop_bound:\n",
    "            selected_keys.append(best_key)\n",
    "            print('Индекс столбца с минимальной приведенной стоимостью: ', selected_keys, '\\n')\n",
    "\n",
    "    m, n = len(selected_keys), len(columns)\n",
    "    print(f\"Кол-во столбцов в RMP {m}/{n}. Степень сжатия задачи {round(100 * m / n, 3)}%\")\n",
    "    return selected_keys"
   ],
   "metadata": {
    "id": "KFzB4HyMObpF"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "%%time\n",
    "selected_keys = column_generation_by_dual(A, c, indexes, columns)"
   ],
   "metadata": {
    "id": "KqZVOju-aeSA"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "%%time\n",
    "# Общий запуск\n",
    "selected_keys = column_generation_by_dual(A, c, indexes, columns)\n",
    "res, model_str = optimize_primal_MIP(A, c, columns)"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "VuErtApuO7ui",
    "outputId": "919a1c21-025e-4931-db74-4291344cfa13"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "res"
   ],
   "metadata": {
    "id": "ea2_JxvIa4lB",
    "outputId": "da77d247-84bd-437f-ecd7-25c57a8e47d9",
    "colab": {
     "base_uri": "https://localhost:8080/"
    }
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "В численном эксперименте рассмотрел два набора данных, обозначим их как усеченный сценарий (`is_small=True`) и расширенный сценарий (`is_small=False`). Оценку эффекта по производительности для алгоритма генерации столбцов предлагаю сделать по отношению к решению полной задачи для этого оценим следующие показатели: степень сжатия модели - кол-во столбцов достаточный для получения оптимального решения исходной задачи; выражается в виде процента \"полезных\" столбцов к общему кол-ву столбцов в задаче. Скорость выявления этого множества столбцов - время работы алгоритма генерации столбцов (CG). Время поиска оптимального решения задачи RMP MIP на множестве \"полезных\" столбцов. С целью оценки преимуществ решения задачи с помощью CG, найдем оптимальное решение для полной задачи, где весь набор столбцов подается на вход модели. Резюмируем отношением скорости нахождения оптимального решения полной задачи к времени расчета CG + время решения RMP MIP (Performance).\n",
    "\n",
    "<!-- - Степень сжатия модели - отношение кол-ва \"полезных\" столбцов в финальной RMP к кол-ву столбцов в MP;\n",
    "уменьшенное кол-во столбцов достаточное для получения оптимального решения исходной задачи. Выражается в виде процента \"полезных\" столбцов к общему кол-ву столбцов в задаче;\n",
    "- Время работы алгоритма генерации столбцов (CG) - длительность всех итераций алгоритма до оставноки;\n",
    "- Время решения RMP MIP - длительность решения задачи RMP на финальном множестве столбцов\n",
    "- С целью оценки преимуществ решения задачи с помощью CG, найдем оптимальное решение для полной задачи, где весь набор столбцов подается на вход модели.\n",
    "-Резюмируем отношением скорости нахождения оптимального решения полной задачи к времени расчета CG + время решения RMP MIP (Performance).   -->\n",
    "\n",
    "\n",
    "\n",
    "<!-- Метрики по расчетам на усеченной и расширенном наборах данных представлены ниже. Оценивал такие метрики как степень сжатия модели - кол-во столбцов достаточный для получения оптимального решения исходной задачи; выражается в виде процента \"полезных\" столбцов к общему кол-ву столбцов в задаче. Скорость выявления этого множества столбцов - время работы алгоритма генерации столбцов (CG). Время поиска оптимального решения задачи RMP MIP на множестве \"полезных\" столбцов. С целью оценки преимуществ решения задачи с помощью CG, найдем оптимальное решение для полной задачи, где весь набор столбцов подается на вход модели. Резюмируем отношением скорости нахождения оптимального решения полной задачи к времени расчета CG + время решения RMP MIP (Performance).   -->\n",
    "\n",
    "**Усеченный кейс**\n",
    "- Степень сжатия задачи: 19.5% или 39 столбцов из 200;\n",
    "- Время работы алгоритма CG: ~0.565сек;\n",
    "- Время решения RMP MIP: ~0.373сек;\n",
    "- Время решения полной задачи MIP: ~2.9сек.\n",
    "- Performance: 2.9 / (0.565 + 0.373) = x3.\n",
    "\n",
    "\n",
    "**Расширенный кейс**\n",
    "- Степень сжатия задачи: 6.125% или 245 столбцов из 4000;\n",
    "- Время работы алгоритма CG: ~90сек;\n",
    "- Время решения RMP MIP: ~200сек;\n",
    "- Время решения полной задачи MIP: ~1160сек.\n",
    "- Performance: 1160 / (200 + 90) = x4."
   ],
   "metadata": {
    "id": "BdpZXK_TOdKn"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Комментарии\n",
    "\n",
    "Обозначу некоторые нюансы в работе с CG. Конечно, это не все нюансы. Если есть опыт в решении других вопросов - комментарии открыты.\n",
    "\n",
    "**Проблема инициализации.** В нашем случае, в качестве начального множества столбцов взяли нулевой (`selected_keys = [0]`), но настолько бездумно поступать недопустимо при инициализации алгоритма. Необходимо сформировать такой набор столбцов, который в совокупности покрывает каждый из рейсов (строку) хотя бы один раз - чтобы RMP была решабельная (feasible). В своей реализации добавил в двойственную задачу ограничение сверху на двойственные переменные (`y_vars.append(model.NumVar(0, 2, name=f\"{name}\"))`). Это позволило избежать неограниченности (unbounded) двойственной постановки (но RMP по-прежнему infeasible, что при высоких моральных ценностях - недопустимо). В случае, когда изначальный набор сменных графиков покрывает весь набор рейсов (RMP - feasible), такой трюк не потребуется.\n",
    "\n",
    "**Добавление группы столбцов.** С количеством итераций решение двойственной задачи становится затратнее по времени несмотря на то, что решаем LP задачу. Чтобы сократить кол-во решений двойственных задач иногда добавляют в RMP сразу группу сменных графиков. Можно сортировать столбцы по приведенной стоимости и добавлять topN столбцов или все столбцы с отрицательной приведенной стоимостью. Но это сказывается на качестве компрессии задачи. Кроме того, можно использовать эвристики или концепцию доминирования, чтобы формировать группу из \"различных\" по некоторой оценке набор столбцов для добавления в RMP на каждой итерации CG.\n",
    "\n",
    "**Прямая расслабленная RMP.** В большинстве своем солверы позволяют извлечь оптимальное решение из прямой задачи. Поэтому реализовывать и решать двойственную задачу явно не требуется. Ниже участок кода, который позволяет извлекать значения двойственных переменных из решения прямой задачи."
   ],
   "metadata": {
    "id": "-TrZA_oZOdHs"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Прямая расслабленная задача RMP\n",
    "def optimize_primal_LP(A, c, columns, selected_x_keys=None):\n",
    "    \"\"\"\n",
    "      A - исходная матрица задачи\n",
    "      с - вектор коэффициентов целевой функции\n",
    "      columns - список сменных графиков\n",
    "      selected_x_keys - список индексов сменных графиков J'\n",
    "    \"\"\"\n",
    "    selected_x_keys = selected_x_keys or list(range(len(A)))\n",
    "\n",
    "    model = pywraplp.Solver.CreateSolver('GLOP')\n",
    "\n",
    "    # Инициализация бинарных переменных выбора сменного графика\n",
    "    x_vars = []\n",
    "    for i in selected_x_keys:\n",
    "        name = columns[i]\n",
    "        x_vars.append(model.NumVar(0, model.infinity(), name=f\"{name}\"))\n",
    "\n",
    "    # Добавляем ограничения в модель\n",
    "    # Ограничение покрытия формируем исходя из доступного набора сменных графиков J'\n",
    "    A_tmp = A[:, selected_x_keys]\n",
    "    constraints = []\n",
    "    for i in range(len(A_tmp)):\n",
    "        constr = model.Add(A_tmp[i] @ x_vars >= 1)\n",
    "        constraints.append(constr)\n",
    "\n",
    "    # Целевая функция\n",
    "    model.Minimize(c[selected_x_keys] @ x_vars)\n",
    "    model_str = model.ExportModelAsLpFormat(False)\n",
    "\n",
    "    model.EnableOutput()\n",
    "    status = model.Solve()\n",
    "\n",
    "    print(\"Найдено оптимальое решение: \", status == pywraplp.Solver.OPTIMAL)\n",
    "    print('Целевая функция = %f' % model.Objective().Value())\n",
    "\n",
    "    solution = {columns[selected_x_keys[i]]: var.solution_value()\n",
    "                for i, var in enumerate(x_vars)\n",
    "               if var.solution_value() }\n",
    "    print(solution)\n",
    "\n",
    "    # Вектор значений решающих переменных\n",
    "    solution_vect = np.array([var.solution_value() for var in x_vars])\n",
    "    # Вектор значений двойственных переменных\n",
    "    dual_solution_vect = np.array([constr.dual_value() for constr in constraints])\n",
    "    return solution_vect, dual_solution_vect, model_str"
   ],
   "metadata": {
    "id": "YwYtv08R2CnQ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "%%time\n",
    "res, res_dual, model_str = optimize_primal_LP(A, c, columns)"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "dTN54VL23seX",
    "outputId": "9679ebfc-e7e1-48e2-fb27-ee234999e454"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Задача раскроя\n",
    "\n",
    "В [википедии](https://ru.wikipedia.org/wiki/%D0%97%D0%B0%D0%B4%D0%B0%D1%87%D0%B0_%D1%80%D0%B0%D1%81%D0%BA%D1%80%D0%BE%D1%8F) есть достаточно хорошее описание задачи, поэтому воспользуюсь уже готовым текстом.\n",
    "\n",
    "\"Представим себе, что вы работаете на целлюлозно-бумажном предприятии, и у вас имеется некоторое количество рулонов бумаги фиксированной ширины, но различным заказчикам нужны различные количества рулонов различной ширины. Как разрезать бумагу, чтобы минимизировать отходы?\n",
    "\n",
    "Согласно данным Конфедерации европейских производителей бумаги[en][1], в 2012 году 1331 бумагоделательная машина в регионе производит в среднем отходов на 56 млн евро (примерно 73 млн долларов США) каждая. Экономии даже 1 % будет очень существенной.\"\n",
    "\n",
    "Отмечу, что на хабре есть статья посвященная задаче раскроя, которую опубликовал мой коллега Михаил Красильников. С ней можете ознакомится [здесь](https://habr.com/ru/articles/733804/)."
   ],
   "metadata": {
    "id": "izlHK58m5wHk"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Математическая постановка\n",
    "$i \\in I$ - множество заказов;\n",
    "\n",
    "$j \\in J$ - множество допустимых схем раскроя, которые можно применить к фиксированной ширине производимого рулона;\n",
    "\n",
    "$\\textbf{a}_j$ - целочисленный $|I|$-мерный вектор констант, где $i$-ая компонента соответствует кол-ву рулонов с шириной заказа $i$, которое получаем при схеме раскроя $j$;\n",
    "\n",
    "$q_i$ - целочисленная константа, объем заказа $i$. Кол-во рулонов определенной ширины, требуемых по заказу $i$;\n",
    "\n",
    "$x_j$ - целочисленная переменная, кол-во раз применения схемы раскроя $j$;\n",
    "\n",
    "$\\sum_{j \\in J} a_{ij}x_j \\ge q_i, \\quad \\forall i \\in I$ - необходимо выполнить все заказы;\n",
    "\n",
    "$\\min \\sum_{j \\in J} x_j$ - минимизируем кол-во рулонов фиксированной ширины к производству.\n",
    "\n",
    "-----------------------------\n",
    "\n",
    "Постановка представленная выше предполагает наличие всех возможных схем раскроя в явном виде. На входе известна ширина производимого рулона (фиксированная ширина производства) и ширина рулона требуемого по заказу $i$. Тогда допустимые схемы раскроя можно задать следующий неравенством:\n",
    "$\\sum_{i \\in I} a_{ij} w_i \\le W$, где $w_i$ ширина рулона по заказу $i$, а $W$ ширина фиксированного рулона. Неявный вид схем раскроя можно использовать на этапе расчета приведенной стоимости (reduced cost).\n",
    "\n",
    "Выбор схемы раскроя для добавления в RMP в явном виде:\n",
    "\n",
    "$j_B =  \\text{arg} \\min \\{\\bar{c_j} = 1-u^\\intercal \\textbf{a}_j|j \\in J\\}$.\n",
    "\n",
    "Теперь представим, что вектор (схему раскроя) $\\textbf{a}_j$ мы можем построить по некоторому критерию, а именно, минимизируя приведенную стоимость. Обозначим $z_i$ - целочисленная переменная, которая соответствует кол-ву рулонов по заказу $i$, тогда схема раскроя для добавления в RMP является оптимальным решением следующей задачи:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\min 1-\\sum_{i\\in I} u_i z_i \\\\\n",
    "\\sum_i w_i z_i \\leqslant W & & \\\\\n",
    "z_i \\in \\mathbb{Z}  \\quad \\forall i \\in I & \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "![Пример схем раскроя](https://drive.google.com/uc?export=view&id=1etM2aP5owBP8sXXr7EKys3d1Iv65K6aU)\n",
    "\n",
    "Оптимальное решение задачи целочисленного программирования будет тем самым столбцом $\\textbf{a}_j$ для добавления в модель RMP. Отмечу, что полученная задача является задачей о рюкзаке.\n",
    "\n",
    "В нашей постановке минимизируем кол-во рулонов к производству. В случае минимизации остатков (отходов) потребуется ввести еще одно слагаемое $c_j = (W - \\sum_i w_i z_i) p$, где $p$ величина убытка за потерю единицы бумаги."
   ],
   "metadata": {
    "id": "jbYUhdX46Prx"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Python реализация задачи раскроя\n",
    "\n",
    "Набор входных данных для эксперимента взят из [википедии](https://ru.wikipedia.org/wiki/%D0%97%D0%B0%D0%B4%D0%B0%D1%87%D0%B0_%D1%80%D0%B0%D1%81%D0%BA%D1%80%D0%BE%D1%8F). Набор небольшой, но в качестве примера - подойдет.\n",
    "\n",
    "В комментария к предыдущей задаче указывал, что можно решать прямую задачу, а из солвера извлечь значения двойственных переменных. В этот раз так и поступим."
   ],
   "metadata": {
    "id": "QahCFgiFK9Zm"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Импортируем необходимые библиотеки\n",
    "from ortools.linear_solver import pywraplp\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import random\n",
    "import time\n",
    "\n",
    "# Входные данные\n",
    "W = 5600  # Фиксированная ширина рулона, результат единицы производства\n",
    "# Ширина рулона по заказам и их кол-во\n",
    "data = [(1380, 22),\n",
    "        (1520, 25),\n",
    "        (1560, 12),\n",
    "        (1710, 14),\n",
    "        (1820, 18),\n",
    "        (1880, 18),\n",
    "        (1930, 20),\n",
    "        (2000, 10),\n",
    "        (2050, 12),\n",
    "        (2100, 14),\n",
    "        (2140, 16),\n",
    "        (2150, 18),\n",
    "        (2200, 20)]\n",
    "\n",
    "data_df = pd.DataFrame(data, columns=['w', 'q'])\n",
    "w = np.array(data_df['w'])  # Ширина рулона по каждому из заказов\n",
    "q = np.array(data_df['q'])  # Кол-во рулонов  по заказам"
   ],
   "metadata": {
    "id": "goIkJ-DUKT8M"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Расслабленная задача RMP\n",
    "def cutting_stock_relaxed(A, q, is_relaxed=True):\n",
    "    \"\"\"\n",
    "      A - текущая матрица схем раскроя\n",
    "      q - кол-во рулонов  по заказам\n",
    "      is_relaxed - параметр переключения между LP и MIP задачами\n",
    "    \"\"\"\n",
    "    m, n = A.shape\n",
    "\n",
    "    if is_relaxed:\n",
    "      model = pywraplp.Solver.CreateSolver('GLOP')\n",
    "    else:\n",
    "      model = pywraplp.Solver.CreateSolver('SCIP')  # На случай RMP MIP\n",
    "\n",
    "    # Инициализация переменных\n",
    "    x_vars = []\n",
    "    for i in range(n):\n",
    "        if is_relaxed:\n",
    "          x_vars.append(model.NumVar(0, model.infinity(), name=f\"x_{i}\"))\n",
    "        else:\n",
    "          x_vars.append(model.IntVar(0, 100, name=f\"x_{i}\"))\n",
    "\n",
    "    # Добавление ограничений в модель\n",
    "    constraints = []\n",
    "    for i in range(m):\n",
    "        constr = model.Add(A[i] @ x_vars >= q[i])\n",
    "        constraints.append(constr)\n",
    "\n",
    "    # Целевая функция\n",
    "    model.Minimize(sum(x_vars))\n",
    "    model_str = model.ExportModelAsLpFormat(False)\n",
    "\n",
    "    status = model.Solve()\n",
    "\n",
    "    print('Значение целевой функции = %f' % model.Objective().Value())\n",
    "\n",
    "    solution = {f\"x_{i}\": var.solution_value()\n",
    "                for i, var in enumerate(x_vars)\n",
    "               if var.solution_value() }\n",
    "\n",
    "    print(solution)\n",
    "    solution_vect = np.array([var.solution_value() for var in x_vars])\n",
    "    if is_relaxed:\n",
    "      dual_solution_vect = np.array([constr.dual_value() for constr in constraints])\n",
    "    else:\n",
    "      dual_solution_vect = []\n",
    "    return solution_vect, dual_solution_vect, model_str"
   ],
   "metadata": {
    "id": "aapOQs0iRwPW"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Оптимизация схемы раскроя для добавления в модель\n",
    "def pricing_problem(u, w, W):\n",
    "    \"\"\"\n",
    "      u - вектор значений двойственных переменных\n",
    "      w - ширина рулонов в зависимости от заказа\n",
    "      W - ширина производимого рулона для раскроя\n",
    "    \"\"\"\n",
    "    m = len(w)\n",
    "\n",
    "    model = pywraplp.Solver.CreateSolver('SCIP')\n",
    "\n",
    "    # Инициализация переменных\n",
    "    z_vars = []\n",
    "    for i in range(m):\n",
    "        z_vars.append(model.IntVar(0, model.infinity(), name=f\"z_{i}\"))\n",
    "\n",
    "    # Добавление ограничений в модель\n",
    "    model.Add(w @ z_vars <= W)\n",
    "\n",
    "    # Целевая функция\n",
    "    model.Minimize(1 - u @ z_vars)\n",
    "\n",
    "    status = model.Solve()\n",
    "\n",
    "    print('Значение целевой функции = %f' % model.Objective().Value())\n",
    "\n",
    "    solution = {f\"z_{i}\": var.solution_value()\n",
    "                for i, var in enumerate(z_vars)\n",
    "               if var.solution_value() }\n",
    "    print(solution)\n",
    "    solution_vect = np.array([var.solution_value() for var in z_vars])\n",
    "    return solution_vect, model.Objective().Value()"
   ],
   "metadata": {
    "id": "Fr5V2J-pQD_B"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Алгоритм генерации столбцов для задачи раскроя в случае неявного задания схем раскроя\n",
    "def cutting_stock_col_gen(q, w, W):\n",
    "    # Инициализация решения\n",
    "    m = q.shape[0]  # кол-во ограничений\n",
    "    # Инициализируем начальный набор схем: одна схема - один заказ\n",
    "    A = np.eye(m, dtype=int)\n",
    "\n",
    "    iter_count = 0\n",
    "    while True:\n",
    "        iter_count += 1\n",
    "        print('Итерация', iter_count)\n",
    "\n",
    "        sol, dual_sol, m_str = cutting_stock_relaxed(A, q)\n",
    "        new_pattern, obj = pricing_problem(dual_sol, w, W)\n",
    "        if obj >= 0:\n",
    "            break\n",
    "        new_pattern = new_pattern.astype(int)\n",
    "        print('Новая схема: ', new_pattern)\n",
    "\n",
    "        new_pattern_matr = new_pattern.reshape((-1, 1))\n",
    "        A = np.hstack((A, new_pattern_matr))\n",
    "\n",
    "    return A"
   ],
   "metadata": {
    "id": "dK2JsvPfQtqE"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "%%time\n",
    "A_schemes = cutting_stock_col_gen(q, w, W)\n",
    "res_vec, _, _ = cutting_stock_relaxed(A_schemes, q, is_relaxed=False)"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "GZQX7NS1Psq9",
    "outputId": "824f49dd-5155-4e2c-a25e-0857f41ab723"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "В результате алгоритма генерации столбцов получим матрицу `A_schemes`, которая содержит достаточное кол-во схем раскроя чтобы получить оптимальное решение полной задачи.\n",
    "\n",
    "Провели самую очевидную инициализацию алгоритма: в начальный набор схем раскроя добавили схему одного разреза, когда весь произведенный рулон отдаем под один заказ. Конечно, более умная эвристика позволит сократить время генерации столбцов и ускорит процесс решения.\n",
    "\n",
    "**Показатели алгоритма:**\n",
    "- Кол-во итераций CG: 37шт;\n",
    "- Время работы CG: ~1.42сек;\n",
    "- Время оптимизации RMP: ~0.05сек;\n",
    "- Оптимальное кол-во рулонов: 73шт.\n",
    "\n",
    "Интегральная проверка решения. По заказам требуется $\\sum_{i \\in I} w_iq_i=407160$м бумаги, запланировано к производству 5600*73=408800м, остаток 1640м или в среднем 1640 / 73 ~ 22.47м с рулона. Кроме этого, результат сходится с [википедией](https://ru.wikipedia.org/wiki/%D0%97%D0%B0%D0%B4%D0%B0%D1%87%D0%B0_%D1%80%D0%B0%D1%81%D0%BA%D1%80%D0%BE%D1%8F).\n",
    "\n",
    "Вы можете убедится в достоверности полученного решения путем реализации своей версии решения задачи раскроя или эвристики инициализации алгоритма CG. Пишите в комментариях, обсудим производительность!\n",
    "\n"
   ],
   "metadata": {
    "id": "djDKx7i0VtOa"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Финальное слово\n",
    "\n",
    "Алгоритм генерации столбцов (Column Generation) инструмент уменьшения размерности задачи. Сильным сигналом к применению данного метода является вытянутая матрица $A$ по кол-ву столбцов.\n",
    "\n",
    "В статье кратко изложил математическое обоснование подхода. Схема работы алгоритма позаимствована у симплекс метода, поэтому в статье присутствует часть описания симплекс метода.\n",
    "\n",
    "На практических примерах продемонстрировал способы применения алгоритма для явного и неявного случаев задания столбцов. Рассмотрел такие задачи как планирование графиков смен экипажей в жд пассажира-перевозках и задачу раскроя.\n",
    "\n",
    "Метод генерации столбцов предоставляет широкие возможности применения эвристик: инициализация алгоритма, размер батчей столбцов для добавления в RMP, отбор столбцов не только по приведенной стоимости, но и с учетом доминирования, стабилизация алгоритма генерации столбцов и др.\n",
    "\n",
    "Существуют более специфические модификации алгоритма, например: строко-столбцовая генерация; генерация столбцов зависимых от строк; использование на этапе branch-and-bound.\n",
    "\n",
    "На практике мне доводилось несколько раз применять этот подход (в BIA Technologies, Huawei и ПГК). Всегда всплывали свои нюансы. Не смотря на это, кратный рост производительности алгоритмов оптимизации оправдывал использование генерации столбцов.\n"
   ],
   "metadata": {
    "id": "Swt-8mtqBPoK"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Ссылки\n",
    "\n",
    "- Ссылка на Jupyter Notebook;\n",
    "- Обзор задачи [CSP](https://www.sciencedirect.com/science/article/abs/pii/S0377221719304916);\n",
    "- Задача [раскроя](https://habr.com/ru/articles/733804/);\n",
    "- Обзор по [алгоритму генерации столбцов](https://pubsonline.informs.org/doi/10.1287/opre.1050.0234?cookieSet=1);\n",
    "- Формализованное описание [алгоритма](https://optimization.cbe.cornell.edu/index.php?title=Column_generation_algorithms#:~:text=Column%20generation%20is%20a%20way,the%20entire%20problem%20is%20achieved).\n",
    "\n",
    "# См. также\n",
    "\n",
    "- Моделирование нелинейных функций - [обзор](https://habr.com/ru/articles/765116/).\n"
   ],
   "metadata": {
    "id": "H2Rx_IQmOdE-"
   }
  }
 ]
}