Linear programming in python using scipy library

In my first publication, I want to talk about how you can quickly and easily solve the linear programming problem using the wonderful scipy library . For similar tasks in python there is also pulp , but for beginners in scipy a more understandable syntax.

Why might linear programming be needed in practice? As a rule, with its help they solve the problem of minimizing the function f (x) (or the inverse problem of maximizing for - f (x)).

Here I will not give theoretical calculations (you can see here ), but I will consider a specific example.

So the challenge.

We have 8 factories that produce a certain amount of products every week. We need to distribute products in 13 stores in such a way as to maximize total profit, while it is allowed to close unprofitable stores.

We know:

1. Productivity of factories (supply field in kg of products per week), their coordinates (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. Productivity of stores (demand field in kg of products per week), their coordinates (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. We also know that the cost of producing sweets: 200 rubles / kg; proceeds from the sale of chocolates: 800 rubles / kg; cost of delivery of sweets: 20 rubles / kg / km; Weekly store expenses: 10,000 rubles.
   

Now we impose the first logical restriction:

1. total supply <= total demand, that is, factories cannot ship more products to the store than they can sell.

We introduce the following notation:

1. costkg = production costs;
2. revkg = revenue from production (sales price per kg);
3. devkg = logistics costs;
4. costweek = store maintenance costs.

To work with the data, we will make cross join factories with stores in order to get all combinations of supplies (factory - store).

We get a table of this kind:

	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

Now let's look at the function that we will optimize (so far without taking into account the content of stores).

$$

Here profit is the profit in rubles from shipments to a particular store, and volume is the volume of factory shipments to this store.

In total, we will have 104 terms of profit (since we have 104 combinations of factories and stores).

The final form of the function will look like:

$$

13 * 10000 is the cost of maintaining stores.

Now we pass to the most linear programming. The description of the scipy.optimize module can be found here . In general:

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

Here:

1. c are the coefficients of our volume (i.e., the costs). In our country they will be negative, because we solve the maximization problem (minimization - f (x));
2. A_ub - these are values ​​at volume for the condition imposed by us, b_ub - values ​​of restrictions;
3. x0_bounds and x1_bounds lower and upper bounds on volume

The limitations are as follows:

1. volume for each row from our table does not exceed supply;
2. the amount of volume for the factory line by line does not exceed supply;
3. the amount of volume for the store line by line does not exceed demand;
4. volume takes values ​​from zero to max (supply) for a particular factory.

Code below:

## 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

Calculate the net profit and see how much each factory will supply and which stores are worth closing:

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)

Reading profit - 828302.9141220199 rubles

Deliveries of each factory:

	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

Deliveries to the store:

	
	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

We see that S_03 is worth closing, and about 30% of demand will be supplied to S_12.

Thus, with the help of simple manipulations, we solved the basic and very simplified, but often encountered in practice problem of optimizing the supply chain.

If you have questions or comments, I will be happy to answer.