https://habrahabr.ru/post/332866/В своей первой публикации мне хочется рассказать о том, как можно быстро и просто решить задачу линейного программирования с помощью замечательной библиотеки
scipy. Для подобных задач в
python есть так же
pulp, но для новичков в
scipy более понятный синтаксис.
Зачем может понадобиться линейное программирование на практике? Как правило, с его помощью решают задачу минимизации функции f(x) (или обратную задачу максимизации для — f(x) ).
Здесь я не буду приводить теоретические выкладки (можно посмотреть
тут), а рассмотрю конкретный пример.
Итак, задача.
У нас есть 8 фабрик, которые каждую неделю производят некоторое количество продукции. Нам нужно распределить продукцию по 13 магазинам так, чтобы максимизировать суммарную прибыль, при этом разрешается закрывать нерентабельные магазины.
Нам известны:
- Производительность фабрик (поле 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
- Производительность магазинов (поле 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
- Так же мы знаем, что затраты на производство конфет: 200 руб/кг; выручка от продажи конфет: 800 руб/кг; стоимость доставки конфет: 20 руб/кг/км; недельные расходы на магазин: 10 000 руб.
Теперь наложим первое логические ограничение:
- общий supply <= общего demand, то есть фабрики не могут отгрузить в магазин больше продукции, чем те могут сбыть .
Введем обозначения:
- costkg = затраты на производство;
- revkg = выручка от производства (цена сбыта за кг);
- devkg = затраты на логистику;
- 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 — это объем отгрузок фабрики в данный магазин.
Всего у нас будет 104 слагаемых profit (т.к. у нас 104 комбинации фабрик и магазинов).
Конечный вид функции будет выглядеть как:
13 * 10000 — это затраты на содержание магазинов.
Теперь переходим к самому линейному программированию. Описание модуля
scipy.optimize вы можете найти
тут. В общем виде:
res = linprog(c, A_ub=A, b_ub=b, bounds=(x0_bounds, x1_bounds), options={"disp": True})
Здесь:
- с — это коэффициенты при наших volume (то есть издержки). У нас они будут отрицательные, т.к. мы решаем задачу максимизации (минимизация- f(x));
- A_ub — это значения при volume для накладываемого нами условия, b_ub — значения ограничений;
- x0_bounds and x1_bounds нижние и верхние границы на volume
Ограничения следующие:
- volume для каждой строки из нашей таблицы не превышает supply;
- сумма volume для фабрики по строкам не превышает supply;
- сумма volume для магазина по строкам не превышает demand;
- 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% спроса.
Таким образом, с помощью нехитрых манипуляций мы решили хоть и базовую и сильно упрощенную, но часто встречающуюся на практике задачу оптимизации цепочки поставок.
Если у вас будут вопросы и замечания, буду рад ответить.