How to Handle Optimization Problems?

👁 Photo by CHUTTERSNAP on Unsplash

Variables, Constraints and Objective

In order to define an optimization problem, you need three things: variables, constraints and an objective. The variables can take different values, the solver will try to find the best values for the variables. Constraints are things that are not allowed or boundaries, by setting these correctly you are sure that you will find a solution you can actually use in real life. The objective is the goal you have in the optimization problem, this is what you want to maximize or minimize. If it’s not completely clear by now, here is a more thorough introduction. Don’t hesitate to continue, because the examples will guide you and explain the concepts in more detail.

Code Examples

Imagine you are the CEO of a small package delivery company. The company discovered that the delivery process can be improved. They want to deliver the combination of most valuable packages first and determine the shortest route for every trip. The last thing they want to improve is the division of tasks. Some employees complain, because they want to deliver more often instead of sorting the packages.

If you want to try the examples yourself, you need a working python environment with the following packages installed: pyomo, pandas and numpy. You should also download a solver like cbc or glpk and use it in the executable path. Some knowledge about list comprehensions and python programming will be helpful and is not covered in this post.

Example 1: Choosing the packages to deliver

The higher the value of a package, the more important this package is to a customer. The delivery company has a table with packages and their value and weight:

👁 Package number with values and weights for every package.

The company wants to select the packages with the highest total value and a total weight that is less than 600 (that’s the highest amount the delivery vans can handle). Which packages should be selected? The objective is clear: we want to maximize the total value of the selected packages. The constraint is that the total weight of the selected packages should be less than or equal to 600. We can define the variables as a binary variable for every package that will be equal to zero if a package is not selected and one if a package is selected.

Let’s start programming.

First we import the packages, load the table and extract data from the table (this makes the code more readable in the next part):

import numpy as np
import pandas as pd
import pyomo.environ as pyo
from pyomo.opt import SolverFactory
import time

# load dataframe with package details
data = pd.read_excel('package_details.xlsx')

# extracting the indici, weights and values from the dataframe
indici = list(data.index.values)
weights = data['Weight'].values
values = data['Value'].values

Now we can start using pyomo:

# create a concrete model
model = pyo.ConcreteModel()

# define the VARIABLES (in the end 0=not selected, 1=selected)
model.x = pyo.Var(indici, within=pyo.Binary)
x = model.x

# define the CONSTRAINT, the total weight should be less than 600
model.weight_constraint = pyo.Constraint(expr= sum([x[p]*weights[p] for p in indici]) <= 600)

# define the OBJECTIVE, we want to maximize the value of the selected packages
model.objective = pyo.Objective(expr= sum([x[p]*values[p] for p in indici]), sense=pyo.maximize)

# print the complete model
model.pprint()

Recall that the variable is a binary decision variable that will be equal to 0 (not selected) or 1 (selected). The constraint makes sure the selected packages have a total weight of 600 or less. To accomplish that, we should multiply the x variable with the weights of the packages and sum these values. The objective calculates the sum of the values of the selected packages (multiplying x and the values of a package and summing) and we want to maximize this objective (sense=pyo.maximize). model.pprint() shows the defined model:

👁 Image

That’s all! Now it’s time to call the solver that will give us the packages we should select:

# replace the executable path by your own
opt = SolverFactory('cbc', executable='C:cbcbincbc.exe')
results = opt.solve(model)

Within a few seconds, the solver has found the optimal solution! We can add a solution column to the original dataframe.

solution = [int(pyo.value(model.x[p])) for p in indici]
data['solution'] = solution

👁 Solution column added.

Can you change this problem to select the highest number of packages instead of the highest total value of selected packages?

Problem solved, let’s continue to the next one!

👁 Photo by Stephen Monroe on Unsplash

Example 2: Determining the order to deliver the selected packages

In order to deliver the selected packages as fast as possible, we want to take the shortest route. We start at the pick up point and we want to end there to pick up new packages. How can we model this problem and solve it using optimization? This one is a bit more difficult than the previous problem. Let’s define the objective, constraints and variables.

The objective is minimization of the total distance of the route. The constraints are: we should start at the pick up point, end at the pick up point, and every package should be delivered, so we need to visit all the addresses at least once. We also need to make sure that we create one complete route instead of multiple small routes. Last but not least, the variables. We create a from, to matrix with binary variables, where a 1 means this way is selected and a 0 means not selected. If x[6, 2] == 1 we go from address 6 to address 2.

For this problem we need a distance matrix. Check in the matrix below that the distance from address 6 to address 2 is equal to 16. Address 0 is our starting (and ending) point. We need to deliver 8 packages, represented by the 8 other addresses.

👁 Distance matrix.

Code to create a random distance matrix:

import random
import pandas as pd

indici = range(9)
distances = pd.DataFrame(index=indici, columns=indici)

for from_ in indici:
 for to in indici:
 if from_ == to:
 distances.loc[from_, to] = 0
 elif pd.isna(distances.loc[from_, to]):
 dist = random.choice(range(1, 30))
 distances.loc[from_, to] = dist
 distances.loc[to, from_] = dist

First we declare the model, the addresses (list with the index values of the distance matrix) and the variables:

model = pyo.ConcreteModel()

addresses = distances.index

# declaring the VARIABLE matrix, again with only binary values
model.x = pyo.Var(addresses, addresses, within=pyo.Binary)
x = model.x

In this example, the variable is a matrix: we create it by looping over the addresses twice. So x[2,3] corresponds to the route from address 2 to address 3.

Now it’s time to declare the constraints. We will use constraint lists this time, because we want to add multiple constraints at the same time.

The first constraint list will hold the values where we travel from and to the same address. We don’t want to activate those so we set all the variable values where the from and to point is the same equal to zero:

model.diagonal = pyo.ConstraintList()
for a in addresses:
 model.diagonal.add(expr= x[a, a] == 0)

The next constraint is about visiting every address once. This means every address should exactly be a ‘to’ address once, and a ‘from’ address once.

model.visit_once = pyo.ConstraintList()
for a in addresses:
 model.visit_once.add(expr=sum([x[a, to] for to in addresses])==1)
 model.visit_once.add(expr=sum([x[fr, a] for fr in addresses])==1)

And the final constraint: we don’t want to travel in circles, we want a complete route! So if we travel from address 0 to address 5, we don’t want to travel from address 5 to address 0! None or one of these two routes should be activated, so x[0, 5] + x[5, 0] <= 1 . This constraint makes sure we get one complete route, instead of multiple smaller ones. Look at the Miller–Tucker–Zemlin formulation to see an explanation of this constraint. We introduce a dummy variable called u here.

n = len(addresses)
model.u = pyo.Var(addresses, within=pyo.NonNegativeIntegers, bounds=(1, n))
u = model.u

model.no_circles = pyo.ConstraintList()
for a1 in range(1, n):
 for a2 in range(1, n):
 if a1 != a2:
 model.no_circles.add(expr= u[a1]-u[a2]+x[a1,a2]*n <= n-1)

Constraints, check! For the objective, we want to minimize the sum of the distances of the activated routes. We loop over the address combinations and multiply the value of x with the distance of the route. The sum of these values gives us the total distance. This time we need to call pyo.minimize, because we want to minimize the total distance.

# declaring the OBJECTIVE

sum_of_distances = sum([x[a1, a2]*distances.loc[a1, a2] for a1 in addresses for a2 in addresses])

model.obj = pyo.Objective(expr=sum_of_distances, sense=pyo.minimize)

Now it’s time to solve the problem. Call the solver and print the results:

opt = SolverFactory('cbc', executable='C:cbcbincbc.exe')
results = opt.solve(model)

print(pyo.value(model.obj)) # gives 62 (for this distance matrix)

solution = []
for a1 in addresses:
 solution.append([pyo.value(model.x[a1, a2]) for a2 in addresses])

The value of the objective is 62, so the total distance of this route is equal to 62. In your case it might be another number, because you probably have another distance matrix. The solution matrix looks like this:

👁 Image

Every row has one 1 and every column has one 1 too. So we visit every address once. To make this solution more readable, use this code:

def print_route(solution):
 sol_ind = list(np.where(np.array(solution[i]) == 1)[0][0] for i in range(len(solution)))
 one_dict = {k:sol_ind[k] for k, _ in enumerate(sol_ind)}

 sol_str = '0 -> '
 nxt = one_dict[0]
 sol_str += str(nxt)
 while True: 
 nxt = one_dict[nxt]
 sol_str += ' -> '
 sol_str += str(nxt)
 if nxt == 0:
 break

 print(sol_str)

print_route(solution)

Result:

👁 Solution pretty printed. Start and end at the pick up point (address 0).

So from the starting point (pick up point, address 0) we go to address 2, address 8, address 5, address 4, and so on until we are back at the pick up point. We delivered all packages and we can be sure this is (one of) the shortest route(s) available!

Note: This problem is a variation of the travelling salesman problem.

You already solved two problems, nice! Let’s continue to the last one!

👁 Photo by Annie Spratt on Unsplash

Example 3: Division of tasks

The last problem is about task division. Currently, you have six employees and you want to keep them happy and energetic. But at the same time, work needs to be done! So how can you accomplish both, finish work and keep your employees happy? Simple, you create another optimization model! You define the tasks that need to be done, like cleaning, sorting packages, delivering packages and being at the service center. For the upcoming month you define how many shifts you need for every task, and how long every shift takes. This gives you the total amount of hours you need employees to work, and you want to give every employee the same amount of working hours, performing the tasks they prefer.

# no of shifts needed for the upcoming month
DELIVERIES = 45
SORTING = 20
SERVICE_CENTER = 20
CLEANING = 10

task = {'DELIVERY': DELIVERIES, 'SORT': SORTING, 'SERVICE': SERVICE_CENTER, 'CLEAN': CLEANING}

# time needed for every shift
time_p_shift = {'DELIVERY': 8, 'SORT': 4, 'SERVICE': 6, 'CLEAN': 3}

# total time needed 
total_hours = sum([task[k]*timepershift[k] for k in task.keys()])

Next up you define your employees and you let them decide what shifts they prefer to do.

employees = {'John':['DELIVERY', 'CLEAN'], 'Jack':['DELIVERY'], 'Joshua':['SORT', 'CLEAN'], 'Janice':['DELIVERY', 'SERVICE'], 'Joyce':['DELIVERY', 'CLEAN', 'SERVICE', 'SORT'], 'Joy':['CLEAN', 'SERVICE', 'SORT']}

emp_list = employees.keys()

John prefers delivery and cleaning, while Joy likes to clean, stay at the service desk and sort packages (she probably doesn’t have a driver’s license). Now you can calculate the average number of hours everyone should work, you will use this value in the objective:

avg_hours_pp = total_hours/len(emp_list)

final_df = pd.DataFrame(index=emp_list, columns=task.keys())

The final_df will receive the optimal values in the end:

👁 The final_df will receive the best variable values later on.

If x['John', 'DELIVERY'] == 6 John should get six delivery shifts this month.

After the data processing, we can start modeling. We create a variable matrix using the employees and the tasks.

model = pyo.ConcreteModel()

model.Employees = pyo.Set(initialize=emp_list)
model.Tasks = pyo.Set(initialize=task.keys())

# define the VARIABLE matrix x using the employees and the tasks
model.x = pyo.Var(model.Employees, model.Tasks, within=pyo.Integers, bounds=(0, None))

# define the hours of work per task
model.hours = pyo.Param(model.Tasks, initialize=time_p_shift)

We want every shift to be filled, so we need to add constraints to make sure we have enough persons for every task:

model.shifts_needed = pyo.ConstraintList()
for t in model.Tasks:
 model.shifts_needed.add(expr= sum([model.x[employee, t] for employee in model.Employees]) == task[t])

We want to exclude certain tasks people don’t want to do:

model.excluded = pyo.ConstraintList()
for employee in model.Employees:
 incl = employees[employee]
 excl_tasks = list(task.keys()-incl)
 for t in excl_tasks:
 model.excluded.add(expr= model.x[employee, t] == 0)

Our goal is to give all the employees an equal amount of work and only tasks they want to do. How can we accomplish that? We calculated the average number of hours avg_hours_pp earlier. We want to minimize the absolute distance from the hours people receive to the avg_hours_pp. Because if that distance is equal to zero, everybody gets exactly the same amount of work.

Calculating an absolute distance is possible, but not as straightforward as using abs(). Pyomo can’t handle the absolute function, because the problem has become nonlinear. So we introduce two new variables here: model.posdelta for the positive distance from avg_hours_pp and model.negdelta for the negative distance. Using a constraint list we set these variable values:

model.posdelta = pyo.Var(model.Employees, bounds=(0, None))
model.negdelta = pyo.Var(model.Employees, bounds=(0, None))

# defining the variables posdelta en negdelta using CONSTRAINTLIST
model.distance_from_avg = pyo.ConstraintList()
for employee in model.Employees:
 model.distance_from_avg.add(expr= sum([model.x[employee, t]*model.hours[t] for t in model.Tasks]) + model.posdelta[employee] - model.negdelta[employee] == avg_hours_pp)

# defining the OBJECTIVE, minimize posdelta + negdelta
model.obj = pyo.Objective(expr= pyo.summation(model.posdelta) + pyo.summation(model.negdelta), sense=pyo.minimize)

Time to call the solver and show the result:

opt = SolverFactory('cbc', executable='C:cbcbincbc.exe')
results = opt.solve(model)

def create_result_frame(model, df):
 solution = pd.DataFrame(np.reshape(pyo.value(model.x[:, :]), (df.shape[0], df.shape[1])), index=emp_list, columns=task.keys())
 solution = solution.astype('int')
 return solution

solution = create_result_frame(model, final_df)

The solution DataFrame looks like this:

👁 Image

It shows us how many shifts everyone should perform. And the most important question: Are the tasks divided equally?

for emp in model.Employees:
 print(f'{emp}:', int(sum([pyo.value(model.x[emp, t])*pyo.value(model.hours[t]) for t in model.Tasks])), 'hours')

👁 Image

Yes, everyone works about the same amount!

If you want more information about handling absolute values, linearization, and other advanced concepts, take a look here.

---

Conclusion

The problems we handled in this post are examples of integer linear programming. It is a technique used in many different industries and has great capabilities. If your problems get bigger, an opensource solver like CBC will take a long time to run. You might want to switch to a commercial solver like Gurobi or CPLEX. If you want to take your optimization skills to the next level, I can recommend this Udemy course.

---

Related

Why Every Data Scientist Should Learn Mathematical Optimization

Taking Your Optimization Skills to the Next Level

---

Don’t forget to subscribe if you’d like to get an email whenever I publish a new article.

Share This Article

• Share on Facebook
• Share on LinkedIn
• Share on X