Operation Research

Linear Programming

Posted by Pritish Roy on July 12, 2020

Linear programming (LP, also called linear optimization) is a method to achieve the best outcome (such as maximum profit or lowest cost) in a mathematical model whose requirements are represented by linear relationships. Linear programming is a special case of mathematical programming (also known as mathematical optimization).

Linear programming (LP) is a powerful framework for describing and solving optimization problems. It allows you to specify a set of decision variables, and a linear objective and a set of linear constraints on these variables. To give a simple and widely used example, consider the problem of minimizing the cost of a selection of foods that meets all the recommended daily nutrient guidelines. The LP model would have a set of decision variables that capture the amount of each food to buy, a linear objective that minimizes the total cost of purchasing the chosen foods, and a linear constraint for each nutrient, requiring that the chosen foods together contain a sufficient quantity of that nutrient. Using linear algebra notation, a linear program can be described as follows: 

Objective:     minimize c^(T) x 
Constraints:  A x = b (linear constraints)
l <= x <=  u (bound constraints)

When described in this form, the vector x represents the decision variables, the vector c captures the linear objective function, the matrix equation Ax = b specifies the linear constraints on x, and the vectors l and u give the lower and upper bounds on x. The set of applications of linear programming is literally too long to list. It includes everything from production scheduling to web advertising optimization to clothing manufacturing. LP touches nearly every commercial industry in some way.

Demonstration with an example

We represent the small nation of Pacific Paradise where we have been working to improve our response to the devastating cyclones that are forecast for the coming season. We would like your assistance in planning our disaster relief distribution strategy (DRDS).
We have modeled our distribution network using three types of nodes: the central disaster relief depot (CDRD); the forward supply depots (FSD); and the disaster demand points (DDP). Our plan is to transport relief goods from the CDRD to the five FSDs prior to a disaster. This will mean that the relief goods will be closer to where they are needed, but also means that if one or more of the depots are destroyed in the cyclone then other depots can be used for supply. Relief goods cost us 145 per tonne. We have estimated the cost of transporting one tonne of relief goods from the CDRD to each of five FSDs and FSDs to 10 DDPs.
In the event of a cyclone, we estimate demands (tonnes) at each of the disasters. We’ve also realized that this plan will not work for us since we have a maximum capacity at each FSD of 790 tonnes. It is good that we will now be able to fit the supplies in our FSDs. However, we’re concerned that some of the DDPs are too reliant on specific FSDs. We’re now thinking it would be safer to make sure that no more than half of the demand at each DDP is coming from a single FSD. We have received some advice from our weather bureau outlining possible scenarios for the path of Cyclone Cyril. For example, they have told us that there is a 10% chance that the path it takes will pass over FSDs 0 and 3, destroying any relief goods we had moved there. 


I         set of FSDs depots
J         set of DDPs depots
S         set of scenarios


cost_cdrd_fsd_{i}             cost of transport from cdrd to fsd for i in I
cost_fsd_ddp_{ij}             cost of transport from fsd for i in I to ddp for j in J
relief_good_cost              cost of relief goods
demand_{j}                    demand at ddp node for j in J
maximum_capacity         maximum capacity limit
scenario_probability_{s}      probability of scenarios for s in S
FSDs_destroyed_{s}            FSDs destroyed in scenario for s in S

C_F_{i}      Amount of supplies to move from CDRD to FSD for i in I
F_D_{ijs}     Amount of supplies to move from FSD for i in I to DDP for j in J


Minimize sum (cost_cdrd_fsd_{i} + relief_good_cost) * C_F_{i} + sum cost_fsd_ddp_{ij} * scenario_probability_{s} * F_D_{ijs} 


  • The Sum of demand at the DDP should equal FSD nodes for every FSD node and scenarios.
sum F_D_{ijs} = C_F_{i}  for every i in I, for every s in S
  • Capacity at each FSD node should not cross the maximum capacity
C_F_{i} <= maximum_capacity   for every i in I
  • The Sum of the amount of supplies to move from FSD nodes should be greater than the demand at each DDP node in each scenario.
sum F_D_{ijs} >= demand_{j} for every j in J, for every s in S
  • At least half the demand should be coming from the FSD nodes.
F_D_{ijs} <= demand_{j} * 1/2 for i in I, for j in J, for s in S
  • The DDP nodes should be equal to 0 if the scenario is in FSDs_destroyed for every FSD, DDP, scenarios.
F_D_{ijs} = 0  for j in J, for i in I, for s in S | if i in FSDs_destroyed[s]
  • And lastly the non-negativity constraints
C_F_{i}, F_D_{ijs} >= 0, for i in I, for j in J, for s in S

The Gurobi Optimizer is a commercial optimization solver for linear programming (LP), quadratic programming (QP), quadratically constrained programming (QCP), mixed integer linear programming (MILP), mixed-integer quadratic programming (MIQP), and mixed-integer quadratically constrained programming (MIQCP). 

Gurobi Optimizer is one of the fastest solvers until now. The license fee of Gurobi optimizer is quite expensive though. If you have an educational email, you can get a license till the end of your educational email validity.

Full Code:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
Created on Wed Mar 11 12:21:12 2020

@author: Pritish Roy

from gurobipy import *

# Sets
I = range(0, 5)

J = range(0, 10)

S = range(0, 6)

# Data
demand = [112, 254, 117, 272, 135, 125, 133, 223, 257, 201]

Cost_CDRD_FSD = [91, 71, 94, 52, 30]

Cost_FSD_DDP = [
    [50, 46, 35, 85, 151, 98, 140, 54, 31, 121],
    [67, 57, 42, 79, 144, 90, 132, 54, 50, 115],
    [70, 67, 57, 103, 168, 115, 157, 74, 48, 139],
    [74, 56, 41, 37, 98, 45, 86, 31, 70, 71],
    [100, 81, 67, 47, 94, 47, 80, 55, 95, 72]

relief_good_cost = 145

maximum_capacity = 790

scenario_probability = [0.10, 0.09, 0.13, 0.18, 0.24, 0.26]

FSDs_destroyed = [
    [0, 3],
    [0, 1],
    [1, 2],
    [1, 4],

# Model
m = Model("Pacific Paradise")

# Variables
FSD_DDP = {}
for i in I:
    for j in J:
        for s in S:
            CDRD_FSD[i] = m.addVar()
            FSD_DDP[i, j, s] = m.addVar()

# Objective
        (Cost_CDRD_FSD[i] + relief_good_cost) * CDRD_FSD[i] for i in I
    ) +
        Cost_FSD_DDP[i][j] * FSD_DDP[i, j, s] * scenario_probability[s]
        for i in I for j in J for s in S

# Set constraints
for i in I:
    m.addConstr(CDRD_FSD[i] <= maximum_capacity)
    for s in S:
            quicksum(FSD_DDP[i, j, s] for j in J) <= CDRD_FSD[i]

for j in J:
    for s in S:
            FSD_DDP[i, j, s] for i in I) == demand[j]
        for i in I:
            m.addConstr(FSD_DDP[i, j, s] <= 1 / 2 * demand[j])

for j in J:
    for i in I:
        for s in S:
            if i in FSDs_destroyed[s]:
                    FSD_DDP[i, j, s] == 0

# Optimize

# Answer
print(f'Optimal cost of Transport + Purchase: {m.objVal}')

print('CDRD ----> FSD')
for i in I:
    print(f"{i} ---- {[round(CDRD_FSD[i].x, 1)]}")


print('\t\t\tFSD ----> DDP')
for i in I:
    print(f"{i} ---- {[round(FSD_DDP[i, j, s].x, 1) for j in J]}")


Optimal cost of Transport + Purchase: 716436.7299999997
CDRD ----> FSD
0 ---- [249.0]
1 ---- [249.0]
2 ---- [790.0]
3 ---- [790.0]
4 ---- [790.0]
			FSD ----> DDP
0 ---- [56.0, 64.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 128.5, 0.0]
1 ---- [0.0, 62.5, 58.5, 0.0, 0.0, 0.0, 0.0, 0.0, 128.0, 0.0]
2 ---- [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
3 ---- [56.0, 127.0, 58.5, 136.0, 67.5, 62.5, 66.5, 111.5, 0.5, 100.5]
4 ---- [0.0, 0.0, 0.0, 136.0, 67.5, 62.5, 66.5, 111.5, 0.0, 100.5]