В своей первой публикации мне хочется рассказать о том, как можно быстро и просто решить задачу линейного программирования с помощью замечательной библиотеки scipy. Для подобных задач в python есть так же pulp, но для новичков в scipy более понятный синтаксис.

Зачем может понадобиться линейное программирование на практике? Как правило, с его помощью решают задачу минимизации функции f(x) (или обратную задачу максимизации для — f(x) ).

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

Итак, задача.

У нас есть 8 фабрик, которые каждую неделю производят некоторое количество продукции. Нам нужно распределить продукцию по 13 магазинам так, чтобы максимизировать суммарную прибыль, при этом разрешается закрывать нерентабельные магазины.

Нам известны:

  1. Производительность фабрик (поле supply в кг продукции за неделю), их координаты (X, Y)

              fact	      x	          y	   supply
    0	F_01	59.250276	59.871183	389
    1	F_02	84.739320	14.179336	409
    2	F_03	42.397937	42.474530	124
    3	F_04	19.539202	13.714643	70
    4	F_05	41.280669	37.860993	386
    5	F_06	37.159066	41.353602	196
    6	F_07	96.890453	64.420010	394
    7	F_08	86.267499	81.662811	365

  2. Производительность магазинов (поле demand в кг продукции за неделю), их координаты (X, Y)

         shop           	x	    y	         demand
    0	S_01	13.490869	73.269974	200
    1	S_02	85.435435	66.637250	20
    2	S_03	28.578297	8.997380	320
    3	S_04	31.324145	91.839907	360
    4	S_05	40.338575	15.487028	360
    5	S_06	41.642451	42.121572	120
    6	S_07	53.983692	20.950457	360
    7	S_08	75.761895	87.067552	60
    8	S_09	81.836739	36.799647	80
    9	S_10	54.260517	25.920108	100
    10	S_11	67.918105	68.108601	340
    11	S_12	92.200710	10.898110	360
    12	S_13	19.966539	39.046271	60
    

  3. Так же мы знаем, что затраты на производство конфет: 200 руб/кг; выручка от продажи конфет: 800 руб/кг; стоимость доставки конфет: 20 руб/кг/км; недельные расходы на магазин: 10 000 руб.

Теперь наложим первое логические ограничение:

  1. общий supply <= общего demand, то есть фабрики не могут отгрузить в магазин больше продукции, чем те могут сбыть .

Введем обозначения:

  1. costkg = затраты на производство;
  2. revkg = выручка от производства (цена сбыта за кг);
  3. devkg = затраты на логистику;
  4. costweek = затраты на обслуживание магазина.

Для работы с данными сделаем cross join фабрик с магазинами, чтобы получить все комбинации поставок (фабрика — магазин).

Получим таблицу такого вида:

	fact	supply	shop	demand	distance	cost_kg	rev_kg	del_kg	cost_week
99	F_08	365	S_09	80	44.863164	200	800	20	10000
100	F_08	365	S_10	100	55.742703	200	800	20	10000
101	F_08	365	S_11	340	13.554211	200	800	20	10000
102	F_08	365	S_12	360	70.764701	200	800	20	10000
103	F_08	365	S_13	60	42.616540	200	800	20	10000

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

$profit = volume* ( revkg - costkg - devkg *distance)$

Здесь profit — прибыль в рублях от отгрузок в конкретный магазин, а volume — это объем отгрузок фабрики в данный магазин.

Всего у нас будет 104 слагаемых profit (т.к. у нас 104 комбинации фабрик и магазинов).

Конечный вид функции будет выглядеть как:

$Profit = sum(profit_1...profit_n) - 13*10000$

13 * 10000 — это затраты на содержание магазинов.

Теперь переходим к самому линейному программированию. Описание модуля scipy.optimize вы можете найти тут. В общем виде:

res = linprog(c, A_ub=A, b_ub=b, bounds=(x0_bounds, x1_bounds), options={"disp": True})

Здесь:

  1. с — это коэффициенты при наших volume (то есть издержки). У нас они будут отрицательные, т.к. мы решаем задачу максимизации (минимизация- f(x));
  2. A_ub — это значения при volume для накладываемого нами условия, b_ub — значения ограничений;
  3. x0_bounds and x1_bounds нижние и верхние границы на volume

Ограничения следующие:

  1. volume для каждой строки из нашей таблицы не превышает supply;
  2. сумма volume для фабрики по строкам не превышает supply;
  3. сумма volume для магазина по строкам не превышает demand;
  4. volume принимает значения от нуля до max(supply) для конкретной фабрики.

Код ниже:

## ff - table with factories and shops
coefs = []
for f in ff.iterrows():
    coefs.append((f[1]['rev_kg'] - f[1]['cost_kg'] - f[1]['del_kg']*f[1]['distance'])*(-1))
A = []
b = []
for f in ff['fact'].values:
    A.append((ff['fact'] == f)*1)
    b.append(ff[ff['fact'] ==f]['supply'].max())
    

for f in ff['shop'].values:
    A.append((ff['shop'] == f)*1)
    b.append(ff[ff['shop'] ==f]['demand'].max())

x0_bounds = []
x1_bounds = []

for f in ff.iterrows():
    x0_bounds.append(0)
    x1_bounds.append(f[1]['demand'])

x0_bounds = tuple(x0_bounds)
x1_bounds = tuple(x1_bounds)
A.append(coefs)
b.append(-10000*13)

res = linprog(coefs, A_ub=A, b_ub=b,  bounds=list(zip(x0_bounds, x1_bounds)), options={"disp": True,   'maxiter'  : 50000}
             )

Output:

Optimization terminated successfully.
Current function value: -948302.914122
Iterations: 20


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

ff['supply_best'] = res.x
ff['stay_opened'] = (ff['supply_best'] > 0)*1
ff['profit'] = (ff['supply_best']*(ff['rev_kg']- ff['cost_kg'] - ff['distance'] * ff['del_kg']))*ff['stay_opened']
net_profit = ff['profit'].sum() - ff[ff['stay_opened']==1]['shop'].nunique()*10000

grouped = ff.groupby(['fact', 'shop'])['supply_best'].sum().reset_index()

f = {'supply_best': 'sum', 'supply': 'max'}
ff.groupby('fact')['supply_best', 'supply'].agg(f)

Читая прибыль — 828302.9141220199 рублей

Поставки каждой фабрики:

	fact	shop	supply_best
0	F_01	S_01	166.0
17	F_02	S_05	360.0
24	F_02	S_12	49.0
31	F_03	S_06	120.0
38	F_03	S_13	4.0
50	F_04	S_12	70.0
58	F_05	S_07	346.0
61	F_05	S_10	40.0
73	F_06	S_09	80.0
74	F_06	S_10	60.0
77	F_06	S_13	56.0
78	F_07	S_01	34.0
79	F_07	S_02	20.0
88	F_07	S_11	340.0
94	F_08	S_04	305.0
98	F_08	S_08	60.0

Поставки в магазин:

	
	supply_best	demand
shop		
S_01	200.0	200
S_02	20.0	20
S_03	0.0	320
S_04	305.0	360
S_05	360.0	360
S_06	120.0	120
S_07	346.0	360
S_08	60.0	60
S_09	80.0	80
S_10	100.0	100
S_11	340.0	340
S_12	119.0	360
S_13	60.0	60

Мы видим, что S_03 стоит закрыть, а в S_12 будет поставляться около 30% спроса.

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

Если у вас будут вопросы и замечания, буду рад ответить.
Поделиться с друзьями
-->

Комментарии (9)


  1. lash05
    11.07.2017 22:44

    Вопрос по библиотеке в целом — не пробовали её использовать для целочисленных задач?


    1. kxx
      12.07.2017 00:11

      Если вы имеете ввиду mixed integer programming, то, насколько мне известно, в scipy нет солвера для такого типа задач. Посмотрите на cvxpy.


    1. paveltro
      12.07.2017 00:14

      Согласно документации, на данный момент в scily.linprog поддерживается только simplex. Для ILP (integer linear programming) подходят pulp, cvxopt


  1. oleg-x2001
    12.07.2017 00:15

    А почему бы не использовать cvxopt? Это, пожалуй, самый мощный пакет для выпуклой оптимизации.


    1. paveltro
      12.07.2017 00:17

      Вы правы, cvxopt более мощный пакет. Scipy был просто под рукой)


  1. AC130
    12.07.2017 14:41

    Статья интересная. Я даже и не знал что scipy можно использовать для решения задач линейного программирования.

    Пара непонятностей. Видится мне что ограничение 1 можно убрать (ибо оно выполняется автоматически при наличии ограничения 2). Ограничение 4 вовсе неясно, почему именно от нуля до одного? Может быть, просто больше нуля?

    И почему ваша модель не учитывает возможность закрытия магазинов, коль такое упомянуто в задаче? Оно понятно что ILP всунуть не получится, но ведь можно решить ослабленную задачу и показать какой получается ответ. Ведь не факт что решение ILP даст те же выводы, что и ваше решение.


    1. paveltro
      12.07.2017 18:40

      Ограничение 4 поправил.

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

      Спасибо за замечания


      1. AC130
        12.07.2017 19:53

        А сколько времени требует один обсчёт решения? Если порядка секунды, то можно просто перебрать 2^13 вариантов. Займёт сутки примерно.


        1. paveltro
          13.07.2017 12:40

          Порядка секунды как раз)