How to Develop Optimization Models in Python


Determining how to design and operate a system in the best way, under the given circumstances such as allocation of scarce resources, usually requires leveraging on quantitative methods in decision making. Mathematical optimization is one of the main approaches for deciding the best action for a given situation. It consists of maximizing or minimizing the real function by systematically choosing input values from an allowed set and computing the value of the objective function.

Some use cases for optimization are:

  • Product planning and inventory control: planning the issue of orders while avoiding stock-out and exceeding capacity limitations
  • Routing decisions: process of determining the most cost-effective route
  • Packing problems: deciding the packing method without exceeding the capacity while minimizing the free space
  • Resource planning: determining the amount of each item and resource needed
  • Scheduling: such as scheduling the shifts of workers
  • Location problems: placement of facilities to minimize transportation costs while ensuring the requested demand.

In this article, I will demonstrate solutions to some optimization problems, leveraging on linear programming, and using PuLP library in Python.

Linear programming deals with the problem of optimizing a linear objective function (such as maximum profit or minimum cost) subject to linear equality/inequality constraints on the decision variables.

There are three components of linear programming:

  1. Decision variable: variables that can be directly controlled by the decision-maker.
  2. Objective function: the linear function that mathematically expresses the quantities to be maximized or minimized.
  3. Constraints: mathematical expression of equalities or inequalities representing the restrictions on the decision variables.

Let’s look at some examples. You can check out the Jupyter Notebook on my GitHub for the full analysis.

1. Resource Planning Example

A bakery makes cakes and pies every day. It can make a total of 30 items in one day, which at least must be 5 cakes and 10 pies for its planned customers. The profit on each cake is $1.5 and the profit on each pie is $2.00. How many should be made to maximize the profit?


This example is simple, meaning it doesn’t require us to use PuLP or any functionalities of Python, yet it is a good exercise to understand the concepts.

The objective function of the question is to maximize the profit:

max Profit = 1.5*C + 2*P

Subject to following constraints:

C >= 5, P >= 10,

C + P <= 30

As can be seen in the graph below, if we draw the constraints: a vertical line at (5, 0), a horizontal line at (0,10) and a slanting line passes from (0,30) and (30,0) to fulfill the requirements, we can find the vertices of the polygonal region that are (5,10), (20,10) and (5,25).

Illustration of the vertices and polygonal region

According to our objective function Profit = 1.5*C + 2*P we can get the maximum profit of $57.5 with (5,25) vertex, which means the bakery should make 5 cakes and 25 pies.

But not all the problems are as easy as this one and if we want to calculate each constraint of each decision variable for complex problems we might have serious timing problems. 😅

Let’s look at another example:

A bakery makes cakes and pies every day of a month. There is: 1 oven, 2 bakers, 1 packaging packer that works only 22 days of the month. The cake requires to use the oven for 1 day and the pie requires 0.5 day. Each baker needs to work for cake 0.5 days and pie 2 days. Packer needs to work for cake 1 day and pie 0.5 days. The profit on each cake is $15 and the profit on each pie is $12. How many should be made to maximize the profit under given conditions?

To solve this problem using PuLP, we will follow the common modeling process.

1. Import PuLP and Initialize Model:

Inside LpProblem() method we define the problem name and sense of objective function which can either ‘LpMaximize’ or ‘LpMinimize’.

import pulp
from pulp import *model = LpProblem('Maximize Bakery Profits', sense= LpMaximize)

2. Define Decision Variables:

Inside LpVariable() method we define a name for the variable, values for lower and upper bound, and category type which can be ‘Integer’, ‘Binary’, or ‘Continuous’. Since we want an integer value for the number of cakes and pies, we choose integer.

C = LpVariable('C', lowBound=0, upBound=None, cat='Integer')
P = LpVariable('P', lowBound=0, upBound=None, cat='Integer')

3. Define the Objective Function:

Add the objective function to the initialized model using+=

model += 15 * C + 12 * P

4. Define the Constraints:

Add the constraints to the initialized model using += Notice that constraints are different from the objective function because they have (in)equalities on the right-hand side.

model += 1 * C + 0.5 * P <= 30
model += 0.5 * C + 2 * P <= 60
model += 1 * C + 0.5 * P <= 22

5. Solve Model

  • Call solve method model.solve()
  • Check the status of the solution LpStatus[model.status]
  • Print optimized decision variables C.varValue P.varValue
  • Print optimized objective function value(model.objective))
Produce 8.0 Cake
Produce 28.0 Pie

2. Scheduling Example

A post office is looking to hire postman, with the requirements to work 5 consecutive days and then 2 days off. The objective is to hire the minimum number of workers and the estimated number of postmen needed for each day is Monday: 25, Tuesday: 32, Wednesday: 22, Thursday: 18, Friday: 24, Saturday: 12, Sunday: 14. What would be the minimum number of postmen to hire?

To solve this problem, we need to write down the constraints in terms of the number of workers we need to start working on each day such as; x_0 is the number of workers starting to work on Monday, x_1 is the number of workers starting to work on Tuesday, etc. By doing so, we can store x_0 from Monday to Friday and x_1 from Tuesday to Saturday since they need to work 5 consecutive days.

Illustration of the model constraints
#Initialize model
model = LpProblem("Minimize Number of Workers", LpMinimize)#Define decision variables
days = list(range(7))
x = LpVariable.dicts('workers_', days, lowBound=0, upbound=None, cat='Integer')#Define model
model += lpSum([x[i] for i in days])# Define constraints
model += x[0] + x[3] + x[4] + x[5] + x[6] >= 25
model += x[0] + x[1] + x[4] + x[5] + x[6] >= 32
model += x[0] + x[1] + x[2] + x[5] + x[6] >= 22
model += x[0] + x[1] + x[2] + x[3] + x[6] >= 18
model += x[0] + x[1] + x[2] + x[3] + x[4] >= 24
model += x[1] + x[2] + x[3] + x[4] + x[5] >= 12
model += x[2] + x[3] + x[4] + x[5] + x[6] >= 14# Solve model
model.solve()#Print model status
print('Status:', LpStatus[model.status])#Print solution variables
for variable in model.variables():
print ('{} = {}'.format(, variable.varValue))

Status: Optimal
workers__0 = 7.0
workers__1 = 7.0
workers__2 = 0.0
workers__3 = 0.0
workers__4 = 10.0
workers__5 = 0.0
workers__6 = 8.0

As the solution suggests, we need to hire 32 postmen in total, 7 of them start working on Monday, the other 7 on Tuesday, 10 on Friday, and 8 on Sunday.

3. Location Problem Example

Assume that you need to optimize a manufacturing company’s supply chain network across 5 selling locations to meet with demand by location at the lowest cost. You can decide the plant size in each location where the options are low capacity and high capacity. One possibility is to set up a facility in each region with an advantage of low transportation costs and a disadvantage of having production plans sized to meet local demands and not exploiting the economies of scale. Another possibility is to set up a few manufacturing plants with an advantage of economies of scale but requiring higher transportation costs. Given, you have the estimated demand for each location, variable costs of transportation from one plant to another, fixed costs of having a plant-based on its size and the production capacity based on plant size is 500 for low capacity and 1500 for high capacity. How would you solve this problem with the minimum cost?

  • Demand column shows the estimated demand for each location
  • Columns A to E shows the transportation cost from locations in the index to each of the locations in columns (i.e. transportation cost from plant B to plant D is 14)
  • ‘High_C’ and ‘Low_C’ columns show the fixed cost of having high and low capacity plants in each location (i.e. having a low capacity plant at location E has $6500 fixed cost)

To solve this problem, we first need to initialize our model and decision variables. There are two decision variables;

  1. Production quantity that produced in plant i and shipped to plant j (continuous variable)
  2. Capacity of the production plant (binary variable: 1 if the plant at location i of capacity s is open, 0 if closed)
#Initialize model
model = LpProblem('Facility Location Problem', LpMinimize)#Define decision variables
loc = ['A', 'B', 'C', 'D', 'E']
size = ['Low_C','High_C']x = LpVariable.dicts('production_', [(i,j) for i in loc for j in loc], lowBound=0, upBound=None, cat='Continuous')
y = LpVariable.dicts('plant_', [(i,s) for s in size for i in loc], lowBound=None, upBound=None, cat='Binary')#Define model
model += (lpSum([fix_cost.loc[i,s] * y[(i,s)] for s in size for i in loc]) + lpSum([var_cost.loc[i,j] * x[(i,j)] for i in loc for j in loc]))

Constraints are:

  1. Total production needs to be equal to total demand
  2. Total production can be smaller or equal to total production capacity
# Define constraints
for j in loc:
model += lpSum([x[(i, j)] for i in loc]) == demand.loc[j,'Demand']for i in loc:
model += lpSum([x[i, j] for j in loc]) <= lpSum([capacity.loc[i,s] * y[i,s] for s in size])

# Solve

Since solve() method returns 1, the solution is optimal and we can print the results.

# Results for production quantities
[{'prod':'{} to {}'.format(i,j), 'quantity':x[(i,j)].varValue}
for i in loc for j in loc]# Results for plant capacities based on location
[{'lowCap':y[(i,size[0])].varValue, 'highCap':y[(i,size[1])].varValue}
for i in loc]# Objective Value
print('Objective = ', value(model.objective))


prod    quantity
A to A 145.4
A to B 0.0
A to C 0.0
A to D 0.0
A to E 1219.6
B to A 0.0
B to B 0.0
B to C 0.0
B to D 0.0
B to E 0.0
C to A 0.0
C to B 84.1
C to C 156.4
C to D 176.8
C to E 0.0
D to A 0.0
D to B 0.0
D to C 0.0
D to D 1500.0
D to E 0.0
E to A 0.0
E to B 0.0
E to C 0.0
E to D 0.0
E to E 1500.0 lowCap highCap
A 0.0 1.0
B 0.0 0.0
C 1.0 0.0
D 0.0 1.0
E 0.0 1.0Objective = 58850.9

As can be seen from results, the model suggests to open a low capacity plant in location C, and high capacity plants in locations A, D, and E. The demand at location B is suggested to be provided by the location C. Doing so, we can achieve the minimum cost of $58850.9.

Sanity Check the Model and Solution

Illustration of the workflow

Once you call the solve method model.solve() and print the status of the solution LpStatus[model.status] you can receive following outcomes.

  1. Optimal: Model and constraints worked well and the solution is feasible. You can continue by checking whether variables are in expected ranges.
  2. Infeasible: There are no feasible solutions. You can turn back and review the constraints.
  3. Unbounded: Meaning that the objective function is not bounded and it tends towards infinity. You can turn back and review the objective function.
  4. Undefined: The optimal solution may exist but may not have been found. You can consider reviewing the model.
  5. Not Solved: Status prior to solving the problem. You can consider reviewing the model.

If the model status is optimal, you can continue and check if the results are in expected ranges. If yes, then the model is completed.

Advanced Topic: Sensitivity Analysis

You can check the dual value if you are seeking to answer some business questions such as what would be the cost increase if the demand increases one unit, in other words, how much profit is needed to cover such a cost increase. You can also check the slack of each constraint if you are seeking to understand whether the constraints are binding, and if not which of them has room in production capacity for the future demand increase.

# Print Dual Value and Slack[{'name':name, 'dual value':c.pi, 'slack': c.slack}
for name, c in model.constraints.items()]name dual value slack
C1 8.0 -0.0
C2 13.0 -0.0
C3 8.0 -0.0
C4 10.0 -0.0
C5 12.0 -0.0
C6 0.0 135.0
C7 -7.0 -0.0
C8 0.0 82.7
C9 -7.0 -0.0
C10 -6.0 -0.0

The dual value represents the amount of the objective function changes for each unit change in the right-hand side of the constraints. Therefore if the demand in location A increases 1 unit, this will increase the overall cost $8, if the demand in location B increases 1 unit, this will increase the overall cost $13, etc.

The slack of each constraint represents the amount of resource that is unused for each constraint. If the slack is equal to 0, then the constraint is binding, if the slack value is greater than 0 then it means its not binding. In the example, we can see that constraints C6 and C8 are not binding, meaning that they do not use their production capacities in full. In location A there is room to increase production by 135, and in location C there is room to increase production by 82.7.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s