В своей первой публикации мне хочется рассказать о том, как можно быстро и просто решить задачу линейного программирования с помощью замечательной библиотеки 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

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

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

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

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

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

$\$inline\$Profit\; =\; sum(profit\_1...profit\_n)\; -\; 13*10000\$inline\$$

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 принимает значения от нуля до одного.

Код ниже:

## 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% спроса.

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

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

https://habrahabr.ru/post/332866/