Unveiling the Power of Linear Optimization in a Delivery Network

Antonello Rabuffi

Published in

BIP xTech

15 min readDec 20, 2023

This article was written in collaboration with Gabriele Biagioni (Data Scientist @ BIP xTech)

The Problem

Imagine a distributor that wants to improve its delivery system. Their delivery network has some issues, they cannot find the best way to manage transportation to handle the required monthly volume of packages. They, on the other hand, struggle to observe the network constraints.

Solving this problem was part of a project, which will be described below, with the aim of optimizing the delivery routes to minimize costs and maximize efficiency.

The hypothesized solution is based on the concept of Linear Optimization and will be easy to understand for those with notions of operations research.

The choice of this approach is due to the nature of the problem, which is based on allocating product quantities to distribution centers as efficiently as possible and within the constraints imposed.

In section Unveiling the Puzzle the theoretical approach is explained.
In section Building the Model the process to obtain the model underneath the solution is depicted.
In section Real World Use Case details about the model are provided (overview, data, implementation, benefits).
The section Beyond the Problem: Real World Applications contains findings and closings.

Unveiling the Puzzle

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. More formally, linear programming is a technique for the optimization of a linear objective function, subject to either linear equality and/or linear inequality constraints.

A linear optimization problem consists of the following elements:

A linear Objective function f(X), in particular, in the case of minimization, one must solve min f(X) ;
A vector of real variables X,
Linear constraints on X.

This is an example where you want to maximize z=f(x,y) with some constraints on the variables x and y:

The constraints of the optimization problem can be of two types:

“Strong” constraints of the optimization problem, such as those in the above figure
“Weak” constraints, imposed as a penalty of the objective function

The solution to the problem must satisfy all the “strong” constraints. However, it may happen, it is not possible to find a feasible solution, the one which satisfies all of them.

To overcome the problem, some of the “strong” constraints, can be transformed into penalties of the objective function, which move the objective function away from the maximum or minimum, but allow the problem to be solved.

It is always possible to write an optimization problem in vector form in the following way:

Figure 2: LP problem formulation, vector form

Where c is the vector of the coefficients of the linear function, A is the matrix with the coefficients of the linear combinations defining the constraints, and b is the vector with the terms on the right-hand side of the inequalities defining the constraints.

With Python as the programming language of choice, there is a plethora of libraries that support linear programming and some of these have been listed below:

SciPy: it provides functions for minimizing (or maximizing) objective functions, possibly subject to constraints. It includes solvers for nonlinear problems (with support for both local and global optimization algorithms), linear programming, constrained and nonlinear least-squares, root finding, and curve fitting.

PuLP: it is a Python linear programming API for defining problems and invoking external solvers. Besides offering flexibility when defining problems and the ability to run various solvers, PuLP is less complicated to use than alternatives like Pyomo or CVXOPT, which require more time and effort to master.

Pyomo: it is a Python-based open-source software package that supports a diverse set of optimization capabilities for formulating, solving, and analyzing optimization models. It supports a wide range of problem types, including Linear programming, Quadratic programming, Nonlinear programming and so on.

Google OR-Tools: it is an open-source software suite for optimization, tuned for tackling the world toughest problems in vehicle routing, flows, integer and linear programming, and constraint programming. It allows the use of different solvers, each of which uses different algorithms for solving linear optimization (LP) problems.

Building the Model

As the project deepened, we defined the variables, such as package assignments to delivery routes and time windows for each stop. We carefully crafted the objective function in order to minimize the total travel distance and time. On top of that, the client could incorporate different constraints, like network capacity and delivery deadlines. Through clever formulation, we could ensure that the model produced feasible and efficient solutions.

The following classes of algorithms for LP are accessible via OR-Tools:

Simplex algorithm: The algorithm walks along the vertices (corner points) of the feasible region, iteratively improving the value of the objective function until reaching an optimal solution.

Barrier or interior-point methods: Barrier methods pair theoretical guarantees of efficient (polynomial time) convergence with reliable performance in practice. They complement the simplex algorithm when it performs poorly.

First-order methods: They are a family of algorithms that use exclusively gradient information (that is, first-order derivatives) to guide the iterations.

The table below lists the LP solvers available in OR-Tools, and indicates which of the three families of algorithms is implemented in each solver:

Figure 3: Table of Google OR-Tools LP solvers and related algorithms

Where ‘G’ indicates that the solver is developed by Google and ‘L’ indicates that the solver requires a license issued by the corresponding third-party developer.

The primary OR-Tools linear optimization solver is GLOP, which is Google in-house linear programming solver. It is fast, memory efficient, and numerically stable. It is written in C++ and was released to the public as part of Google OR-Tools software suite in 2014.

GLOP uses a revised simplex algorithm optimized for sparse matrices. It uses Markowitz pivoting to reduce matrix fill-in, steepest-edge pricing to avoid degenerate pivots, and a LU decomposition tailored for sparse matrices.

Inside Google, GLOP is used to stabilize YouTube videos, and outside Google, it has been used to perform fast linear relaxations for reinforcement learning.

Suppose we want to solve the problem in Figure 1 using the GLOP solver and Python’s OR-Tools library.

These are the steps:

Import OR-Tools solver linear wrapper

# Import OR-Tools linear solver wrapper 
from ortools.linear_solver import pywraplp

2) Declare GLOP Solver

# Declare GLOP solver 
solver = pywraplp.Solver.CreateSolver('GLOP')

3) Create variables x and y, with x, y >= 0

# Create variables x and y, with x, y >= 0 
x = solver.NumVar(0, solver.infinity(), 'x') 
y = solver.NumVar(0, solver.infinity(), 'y')

4) Define linear optimization problem constraints

# Constraint 1: x + y <= 6. 
solver.Add(x + y <= 6.0) 

# Constraint 2: 3x + y <= \5. 
solver.Add(3 * x + y <= 15.0) 

# Constraint 3: x + 3y <= 15. 
solver.Add(x + 3 * y <= 15.0)

5) Define an objective function

# Objective function: 0.5x + 1.5y. 
solver.Maximize(0.5 * x + 1.5 * y)

6) Solve the problem invoking the solver

# Solver invocation 
status = solver.Solve()

7) Display the solution with elapsed time and number of iterations

if status == pywraplp.Solver.OPTIMAL:
  print('Solution:')
  print('Objective value =', solver.Objective().Value())
  print('x =', x.solution_value())
  print('y =', y.solution_value())
  print(f"Problem solved in {solver.wall_time():d} milliseconds")
  print("Iterations:", solver.iterations())
else:
  print('The problem does not have an optimal solution.')

8) Output screen

Solution:
Objective value = 7.5
x = 0.0
y = 5.0
Problem solved in 0 milliseconds
Iterations: 3

Real World Use Case

Overview of the Model

The Network Attribution Model (NAM) is a model that aims to determine the proper distribution of products (and their variants) to different networks and distribution centers. It takes as input an estimate of daily product quantities to be distributed to 4,668 ZIP codes across the national territory by four networks. The model then assigns 10 products (52 if variants are included), between mail and packages, to 1,112 distribution centers, taking into account network preferences, performance, costs, and capacities while considering organizational, operational, logistical and functional constraints. The model is implemented in Python using the Google OR-Tools library and runs on a virtual machine. The optimization process is efficient, reducing processing times from 8 hours to 1 hour. The model outputs two files, one in .csv format with approximately 4 million rows and one in .xlsx format, containing aggregated results and different business-related sheets for evaluating the model performance.

Efficiency in product allocation was measured by certain metrics (such as cost reduction) and by adherence to appropriate constraints (such as distribution center capacity limits). Therefore, the decision was made to solve this problem using a linear optimization model, implemented in Python using the OR-Tools library described earlier.

Model Data

The input files available were three Excel files, the first containing information about the products to be delivered, the second with details of the distribution centers, and the third containing the values of the input parameters.

Specifically, the former contained:

Total quantities of product expected to be distributed daily.
Percentages of daily product expected at the sorting center (SC) level.
Expected daily product percentages at the zip code level.
Product variant percentages at the zip code level.
Constraints on delivery of products only by certain networks.
Product delivery costs for each network at the zip code level.
Network preferences for each product variant.
Product variant preferences for each network.
Network preferences for each product variant in the allocation of remaining quantities.

The second contained:

Capacities of individual distribution centers.
Distribution center coverage at the zip code level.

The third and last contained:

Parameter values used as weights in the cost function of the optimization model.

Model Implementation

The input data just described are processed to obtain a dataframe suitable to:

Write the cost function of the optimization model.
Set the strong and weak constraints of the model.
Aggregate the data at the level of zip code, distribution center, product and product variant.
Easily generate the product quantity column assigned by the model (one value per row).

The following figure shows a screen of a part of the dataframe just described of which we will describe the individual fields shortly.

To understand the structure of this dataframe, the main columns of one of its rows can be analyzed.

In particular, the following columns indicate the maximum quantity (column VOL_VAR_DC) of the fourth variant of the third product, which can be assigned to the distribution center 77586 (column DC) delivering to the zip code 7010 (column ZIP_CODE):

The columns ZIP_CODE, PRODUCT, PRODUCT_VARIANT, NETWORK and DC can be considered the keys of the dataframe.

In addition, the capacity of the distribution center for each product category is also indicated on each row:

There is also other information related to the delivery of the row product variant, such as delivery preferences and delivery cost:

In the end, the optimization model will assign a quantity of product variant for each row in this dataframe, respecting all the constraints of the problem, i.e., considering the costs and delivery preferences just shown, as well as the capacities of the various distribution centers.

The fields of the obtained dataframe can be divided into groups.

Columns associated with logistics

ZIP_CODE - Delivery zip code
NETWORK — Delivery network to which the DC belongs
DC — Distribution Center
AREA — Delivery area not connected with outside

Columns associated with products to be assigned to distribution centers

PRODUCT_CATEGORY — Category the product belongs to
PRODUCT — Product identifier
PRODUCT_VARIANT — Product variant identifier

The order of these columns respects the hierarchical order, each product belongs to a category and there may be one or more variants of the same product.

Columns associated with the estimated quantities of product to handle on a daily basis

VOL — Estimated daily quantity of row product distributed from the row sorting center to the row zip code
VOL_VAR — Estimated daily quantity of the row product variant distributed from the row sorting center to the row zip code
VOL_VAR_DC — Estimated daily quantity of the row product variant distributed from the row sorting center to the row zip code, handled by the row distribution center.

Columns associated with row distribution center capacities

MIN_DC_CATEGORY_1 — Minimum contractual quantity of the product belonging to the first category managed by the row distribution center
MIN_DC_CATEGORY_2 — Minimum contractual quantity of the product belonging to the second category managed by the row distribution center
MAX_DC_CATEGORY_1 — Maximum quantity of the product belonging to the first category managed by the row distribution center
MAX_DC_CATEGORY_2 — Maximum quantity of the product belonging to the second category managed by the row distribution center
MAX_DC_CAPACITY — Maximum capacity of row distribution center, regardless of product category

Columns associated with delivery preferences (network, product, etc.).

NETWORK_ORDER — Integer number from 1 to 4 indicating the delivery preference of the row product variant by the row distribution center
PRODUCT_ORDER — Integer number from 1 to 10 indicating for each row network, the delivery preference of the row product
NETWORK_ORDER_REVERSE — Negative number from -1 to -4 obtained by associating the highest preference 1 with the negative number -4, preference 2 with the number -3, and so on
PRODUCT_ORDER_REVERSE — Negative number from -1 to -10 obtained by associating the highest preference 1 with the negative number -10, preference 2 with the number -9, and so on.

Fields indicated with REVERSE are required to have a higher penalty in the objective function when, for example, products are assigned to networks with low preference.

Columns associated with delivery cost

COST — Cost of delivery of the row product by the row network
NORM_COST — Normalized cost of delivery of the row product by the row network

Using the dataframe just described, it is very easy to define the optimization problem to be solved, with its variables and constraints, using the Python OR-Tools library.

Below is the screen related to the definition of the problem variables:

The defined X variables are all bounded above by the row value taken by VOL_VAR_DC (defined earlier). Each variable is assigned a unique name that depends on DC, PRODUCT, PRODUCT_VARIANT and ZIP_CODE.

Each row in the dataframe described above will correspond to a variable that will have a unique name, defined as f-string in the code using all row references (i.e., DC, SC, ZIP_CODE, NETWORK, PRODUCT, and PRODUCT_VARIANT).

After defining the variables, it is possible to easily define the strong and weak constraints associated with the optimization problem.

Strong constraints are easily defined, as can be seen from the following code screens:

The first strong constraint is on the maximum capacity of the distribution center, so for each DC the sum of the Xs associated with that DC must be less than or equal to its maximum capacity.

The second and final strong constraint requires that the sum of the quantities of product variant X assigned to individual DCs at the zip code level do not exceed VOL_VAR.

So-called “weak” constraints have been used to upper and/or lower limit the capacity of packages or mail to be allocated to certain distribution centers; these constraints were applied as objective function penalties, as we can see in the code below:

where:

Penalty 1: penalties for complying with the contractual minimum product category for certain distribution centers.
Penalty 2: penalties for complying with the maximum product category for certain distribution centers.

The definition of the objective function requires the use of the parameters set in the third input file.

Their role is to be the coefficients (or weights) of the linear combination that defines the objective function of the problem; in particular, the greater the value associated with the parameter, the greater the significance of the function terms associated with the parameter.

An example of parameter enhancement is shown in the figure below:

Different weights can be assigned for each delivery network and the weights are associated with the “importance” of delivery cost, product and network orders, and penalties within the objective function.

The greater weight given to the first penalty is due to a contractual obligation that defines the minimum amount of product that the distribution center must handle daily.

This means that in the optimizer’s resolution of the problem, some emphasis will be placed on meeting this obligation.

Using the first three parameters, it is possible to define the part of the objective function of the optimization problem that does not contain the penalties of weak constraints, as shown in the following code screen:

It is important to note that each row of the initial dataframe and, therefore, each variable X in the problem will have these weights associated and this will influence the algorithm assignment of values to the variables X.

The last two parameters, associated with penalties, indicate the relevance of the penalty within the objective function.

In our case, since it is a minimization problem, they will indicate the order of magnitude of the terms of the objective function associated with noncompliance with the weak constraints of the problem.

The objective function defined up to this point will therefore have the following form:

with:

Note that, for example, a product order of 1 corresponds to a po_reverse of -10, so the larger the quantity

(always positive), the more the function will move towards the minimum.

The optimizer, therefore, will tend to assign as much product as possible to the networks with the highest priority, as its goal is to minimize the objective function.

As a last thing, you instruct the OR-Tools solver to solve the minimum problem in the following way:

Solving the optimization problem can take up to one and a half hours if it is performed on all areas.

At the end of solving the optimization problem, the status of the solver can be displayed, which can be of three types: OPTIMAL, FEASIBLE and INFEASIBLE.

We call solutions that satisfy all the strong constraints of the optimization problem feasible, and the feasible solution that has, in our case, minimum value an optimal solution.

The solutions obtained in our case are always of the OPTIMAL type, and in addition to this, OR-Tools allows us to display other details such as the number of iterations and the value taken by the objective function at its point of minimum (or objective value):

As a result of solving the problem, the n variables X are valued, with n equal to the number of rows in the dataframe described above.

Therefore, a column called VALUE is added to the dataframe and will contain all the values assigned to the Xs for the minimization of the objective function.

Solving the problem within all constraints does not allow the allocation of all the input products, so it was necessary to redistribute the remaining quantity according to certain allocation logic, defined in the first input file in the part containing data on network preferences for each product variant in the allocation of remaining quantities.

A final section of the code was therefore devoted to the allocation of these remaining quantities.

We give below an example in which remaining quantity is redistributed among the available distribution centers:

In this case, the quantity of product variant to be delivered to zip code 20037 is 618.732 units.

As a result of the optimization, a product quantity of 452.841 units was allocated to the two distribution centers 38693 and 38663.

Therefore, 165.891 units (column REMAINING_PROD_VAR) remain to be allocated and they are redistributed between the two DCs according to logic associated with network preferences but also with the remaining capacities of the DCs (column REMAINING_PROD_VAR_ASSIGNED).

Model Benefits

This model has brought many advantages. It automatically estimates, for each type of product, a quantity to be delivered daily to the individual distribution centers, and within all constraints. Even though it is not possible to comply with every constraint and allocate the expected quantities of products, the result is still useful, because it suggests a different workload distribution among the distribution centers. The model, however, manages to distribute more than 95% of the planned quantities.

Beyond the Problem: Real-World Applications

Linear optimization emerges as a powerful and flexible mathematical methodology that is finding increasing applications in distribution networks and a wide range of other sectors. Its ability to solve optimization problems efficiently and precisely has been extensively demonstrated, significantly contributing to resource optimization, cost reduction, and increased efficiency.

In the context of delivery networks, linear optimization offers robust solutions for route planning, fleet management, and production scheduling, enabling companies to enhance product delivery to customers and optimize logistic resources. However, the applications of linear optimization do not stop there. This technology is rapidly spreading into sectors such as finance, manufacturing, logistics, healthcare, and many others, where it has proven to be equally effective in optimizing complex processes.

In conclusion, linear optimization is a legacy, yet very versatile application that will continue to play a significant role in process optimization across various domains. Its widespread adoption and contribution to solving complex problems are the testimony of its growing importance in modern enterprises.