Integer vs. Linear Programming in Python

Mixed Integer Programming for optimization with Google OR-Tools

linear programming
Author

Maxime Lbonne

Published

March 5, 2022

Why is linear programming called that way?

Both terms are confusing:

In summary, it has nothing to do with code: linear or not. It’s about optimizing variables with various constraints.

In this article, we’re gonna talk about another type of optimization: integer programming. We’ll see why a good understanding of the problem we face is necessary to choose the right solver. Finally, we will write a model that can take on a bigger challenge and actually solve a whole class of optimization problems.

You can run the code from this tutorial with the following Google Colab notebook.

📊 I. Optimization problem types

In the introduction to linear programming, we optimized an army composition. Here was the result:

================= Solution =================
Solved in 87.00 milliseconds in 2 iterations

Optimal power = 1800.0 💪power
Army:
 - 🗡️Swordsmen = 6.0000000000000036
 - 🏹Bowmen = 0.0
 - 🐎Horsemen = 5.999999999999999

How can we have 5.999… horsemen? We specified that our variables should be integers with VarInt. What was wrong with our code?

The problem is not the model but the choice of the solver.

GLOP is a pure linear programming solver. This means that it cannot understand the concept of integers. It is limited to continuous parameters with a linear relationship.

This is the difference between linear programming (LP) and integer linear programming (ILP). In summary, LP solvers can only use real numbers and not integers as variables. So why did we declare our variables as integers if it doesn’t take them into account?

GLOP cannot solve ILP problems, but other solvers can. Actually, a lot of them are mixed integer linear programming (MILP, commonly called MIP) solvers. This means that they can consider both continuous (real numbers) and discrete (integers) variables. A particular case of discrete values is Boolean variables to represent decisions with 0–1 values.

Other solvers like SCIP or CBC can solve both MILP and MINLP (mixed integer nonlinear programming) problems. Thanks to OR-Tools, we can use the same model and just change the solver to SCIP or CBC.

# Create the linear solver using the CBC backend
solver = pywraplp.Solver('Maximize army power', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

# 1. Create the variables we want to optimize
swordsmen = solver.IntVar(0, solver.infinity(), 'swordsmen')
bowmen = solver.IntVar(0, solver.infinity(), 'bowmen')
horsemen = solver.IntVar(0, solver.infinity(), 'horsemen')

# 2. Add constraints for each resource
solver.Add(swordsmen*60 + bowmen*80 + horsemen*140 <= 1200)
solver.Add(swordsmen*20 + bowmen*10 <= 800)
solver.Add(bowmen*40 + horsemen*100 <= 600)

# 3. Maximize the objective function
solver.Maximize(swordsmen*70 + bowmen*95 + horsemen*230)

# Solve problem
status = solver.Solve()

# If an optimal solution has been found, print results
if status == pywraplp.Solver.OPTIMAL:
  print('================= Solution =================')
  print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
  print()
  print(f'Optimal value = {solver.Objective().Value()} 💪power')
  print('Army:')
  print(f' - 🗡️Swordsmen = {swordsmen.solution_value()}')
  print(f' - 🏹Bowmen = {bowmen.solution_value()}')
  print(f' - 🐎Horsemen = {horsemen.solution_value()}')
else:
  print('The solver could not find an optimal solution.')
================= Solution =================
Solved in 3.00 milliseconds in 0 iterations

Optimal value = 1800.0 💪power
Army:
 - 🗡️Swordsmen = 6.0
 - 🏹Bowmen = 0.0
 - 🐎Horsemen = 6.0

Strictly speaking, our variables are still floats (type(swordsmen.solution_value()) = float) but we can see that they don’t have weird decimals anymore: the CBC solver really considered them as integers.

In this example, we would generally just round up these values since the error is insignificant. However, it is important to remember to choose the appropriate solver according to the studied problem:

  • LP for continuous variables;
  • MIP/MILP for a combination of continuous and discrete variables.

There are other types such as quadratic (QP) or nonlinear (NLP or MINLP, with an exponential objective function or constraints for instance) problems. They’re applied in different contexts, but follow the same principles as LP or MIP solvers.

🧱 II. Building a general model

But what if our resources change? Or if the cost of a unit evolved? What if we upgraded horsemen and their power increased?

One of the best perks of OR-Tools is that it uses a general-purpose programming language like Python. Instead of static numbers, we can store our parameters in objects like dictionaries or lists.

The code won’t be as readable, but it becomes much more flexible: actually, it can be so flexible that we can solve an entire class of optimization problems without changing the model (just the parameters).

Let’s transform our input parameters into Python lists and feed them to the solver through a function.

# Inputs
UNITS = ['🗡️Swordsmen', '🏹Bowmen', '🐎Horsemen']

DATA = [[60, 20, 0, 70],
        [80, 10, 40, 95],
        [140, 0, 100, 230]]

RESOURCES = [1200, 800, 600]


def solve_army(UNITS, DATA, RESOURCES):
  # Create the linear solver using the CBC backend
  solver = pywraplp.Solver('Maximize army power', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

  # 1. Create the variables we want to optimize
  units = [solver.IntVar(0, solver.infinity(), unit) for unit in UNITS]

  # 2. Add constraints for each resource
  for r, _ in enumerate(RESOURCES):
    solver.Add(sum(DATA[u][r] * units[u] for u, _ in enumerate(units)) <= RESOURCES[r])

  # 3. Maximize the objective function
  solver.Maximize(sum(DATA[u][-1] * units[u] for u, _ in enumerate(units)))

  # Solve problem
  status = solver.Solve()

  # If an optimal solution has been found, print results
  if status == pywraplp.Solver.OPTIMAL:
    print('================= Solution =================')
    print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
    print()
    print(f'Optimal value = {solver.Objective().Value()} 💪power')
    print('Army:')
    for u, _ in enumerate(units):
      print(f' - {units[u].name()} = {units[u].solution_value()}')
  else:
    print('The solver could not find an optimal solution.')

solve_army(UNITS, DATA, RESOURCES)
================= Solution =================
Solved in 2.00 milliseconds in 0 iterations

Optimal value = 1800.0 💪power
Army:
 - 🗡️Swordsmen = 6.0
 - 🏹Bowmen = 0.0
 - 🐎Horsemen = 6.0

We obtain the same results: our code seems to work. Now let’s change the parameters to tackle a slightly more complex problem.

Imagine we have a lot more resources: 🌾183000, 🪵90512, and 🪙80150, so we can also produce a lot more units! This is the new table:

Unit 🌾Food 🪵Wood 🪙Gold 💪Attack ❤️Health
🗡️Swordsman 60 20 0 6 70
🛡️Man-at-arms 100 0 20 12 155
🏹Bowman 30 50 0 5 70
❌Crossbowman 80 0 40 12 80
🔫Handcannoneer 120 0 120 35 150
🐎Horseman 100 20 0 9 125
♞Knight 140 0 100 24 230
🐏Battering ram 0 300 0 200 700
🎯Springald 0 250 250 30 200

Notice that we transformed the 💪power into two values: 💪attack and ❤️health, which is a little more detailed. Health values are higher than attack values, which is why we want to add a weighting factor to make them more comparable.

Let’s take 10 as an example, so power = 10 \times attack + health. Our objective function becomes:

\text{maximize} \ \sum_{u \in units} (10\times attack + health) \cdot u

Adapting our code to this new problem is actually quite simple: we just have to change the input parameters and update the objective function.

UNITS = [
    '🗡️Swordsmen',
    '🛡️Men-at-arms',
    '🏹Bowmen',
    '❌Crossbowmen',
    '🔫Handcannoneers',
    '🐎Horsemen',
    '♞Knights',
    '🐏Battering rams',
    '🎯Springalds',
    '🪨Mangonels',
]

DATA = [
    [60, 20, 0, 6, 70],
    [100, 0, 20, 12, 155],
    [30, 50, 0, 5, 70],
    [80, 0, 40, 12, 80],
    [120, 0, 120, 35, 150],
    [100, 20, 0, 9, 125],
    [140, 0, 100, 24, 230],
    [0, 300, 0, 200, 700],
    [0, 250, 250, 30, 200],
    [0, 400, 200, 12*3, 240]
]

RESOURCES = [183000, 90512, 80150]


def solve_army(UNITS, DATA, RESOURCES):
  # Create the linear solver using the CBC backend
  solver = pywraplp.Solver('Maximize army power', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

  # 1. Create the variables we want to optimize
  units = [solver.IntVar(0, solver.infinity(), unit) for unit in UNITS]

  # 2. Add constraints for each resource
  for r, _ in enumerate(RESOURCES):
    solver.Add(sum(DATA[u][r] * units[u] for u, _ in enumerate(units)) <= RESOURCES[r])

  # 3. Maximize the new objective function
  solver.Maximize(sum((10*DATA[u][-2] + DATA[u][-1]) * units[u] for u, _ in enumerate(units)))

  # Solve problem
  status = solver.Solve()

  # If an optimal solution has been found, print results
  if status == pywraplp.Solver.OPTIMAL:
    print('================= Solution =================')
    print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
    print()
    print(f'Optimal value = {solver.Objective().Value()} 💪power')
    print('Army:')
    for u, _ in enumerate(units):
      print(f' - {units[u].name()} = {units[u].solution_value()}')
  else:
    print('The solver could not find an optimal solution.')

solve_army(UNITS, DATA, RESOURCES)
================= Solution =================
Solved in 74.00 milliseconds in 412 iterations

Optimal value = 1393145.0 💪power
Army:
 - 🗡️Swordsmen = 2.0
 - 🛡️Men-at-arms = 1283.0
 - 🏹Bowmen = 3.0
 - ❌Crossbowmen = 0.0
 - 🔫Handcannoneers = 454.0
 - 🐎Horsemen = 0.0
 - ♞Knights = 0.0
 - 🐏Battering rams = 301.0
 - 🎯Springalds = 0.0
 - 🪨Mangonels = 0.0

This problem would take a long time for humans to address, but the ILP solver did it in the blink of an eye. Better than that: it also gives us the guarantee that our solution is optimal, which means that our enemy cannot find a better army composition for the same cost!

We could increase the number of units and give billions of resources but you get the picture: it would just take longer to obtain a solution, but it wouldn’t change the problem.

⚔️ III. Combining constraints

Now, let’s say we scouted our enemy and know that their army has a 💪power of 1,000,000. We could build a much better army, but our resources are precious and it wouldn’t be very efficient: all we have to do is to build an army with a 💪power higher than 1,000,000 (even 1,000,001 would be enough).

In other words, the total power is now a constraint (💪 > 1,000,000) instead of the objective to maximize. The new goal is to minimize the resources we need to produce this army. However, we can reuse our input parameters since they didn’t change.

The new constraint can be translated as “the sum of the power of the selected units must be strictly greater than 1,000,000”.

\sum_{u \in units} (10\times attack + health) \cdot u > 1\,000\,000

In code, we can loop through our units and resources to design this constraint.

The objective function also has to change. Our goal is to minimize the sum of resources spent to build the army.

\text{minimize} \ \sum_{u \in units} (food + wood + gold) \cdot u

Once again, we can loop through our resources to implement it in OR-Tools.

def solve_army(UNITS, DATA, RESOURCES):
  # Create the linear solver using the CBC backend
  solver = pywraplp.Solver('Minimize resource consumption', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

  # 1. Create the variables we want to optimize
  units = [solver.IntVar(0, solver.infinity(), unit) for unit in UNITS]

  # 2. Add power constraint
  solver.Add(sum((10 * DATA[u][-2] + DATA[u][-1]) * units[u] for u, _ in enumerate(units)) >= 1000001)

  # 3. Minimize the objective function
  solver.Minimize(sum((DATA[u][0] + DATA[u][1] + DATA[u][2]) * units[u] for u, _ in enumerate(units)))

  # Solve problem
  status = solver.Solve()

  # If an optimal solution has been found, print results
  if status == pywraplp.Solver.OPTIMAL:
    print('================= Solution =================')
    print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
    print()

    power = sum((10 * DATA[u][-2] + DATA[u][-1]) * units[u].solution_value() for u, _ in enumerate(units))
    print(f'Optimal value = {solver.Objective().Value()} 🌾🪵🪙resources')
    print(f'Power = 💪{power}')
    print('Army:')
    for u, _ in enumerate(units):
      print(f' - {units[u].name()} = {units[u].solution_value()}')
    print()

    food = sum((DATA[u][0]) * units[u].solution_value() for u, _ in enumerate(units))
    wood = sum((DATA[u][1]) * units[u].solution_value() for u, _ in enumerate(units))
    gold = sum((DATA[u][2]) * units[u].solution_value() for u, _ in enumerate(units))
    print('Resources:')
    print(f' - 🌾Food = {food}')
    print(f' - 🪵Wood = {wood}')
    print(f' - 🪙Gold = {gold}')
  else:
      print('The solver could not find an optimal solution.')

solve_army(UNITS, DATA, RESOURCES)
================= Solution =================
Solved in 4.00 milliseconds in 0 iterations

Optimal value = 111300.0 🌾🪵🪙resources
Power = 💪1001700.0
Army:
 - 🗡️Swordsmen = 0.0
 - 🛡️Men-at-arms = 0.0
 - 🏹Bowmen = 0.0
 - ❌Crossbowmen = 0.0
 - 🔫Handcannoneers = 0.0
 - 🐎Horsemen = 0.0
 - ♞Knights = 0.0
 - 🐏Battering rams = 371.0
 - 🎯Springalds = 0.0
 - 🪨Mangonels = 0.0

Resources:
 - 🌾Food = 0.0
 - 🪵Wood = 111300.0
 - 🪙Gold = 0.0

The solver found an optimal solution: we need to build 371 🐏battering rams for a total cost of 111,300 🪵wood. Wait, what if we don’t have that much wood? In the previous section, we only had 🪵90512: we cannot produce 371 🐏battering rams. 😱

So is it possible to take these limited resources into account and still try to build the best army? Actually, it’s super easy: we just have to copy/paste the constraints from the previous section.

In this version, we have two types of constraints:

  • The total power must be greater than 1,000,000;
  • We cannot spend more than our limited resources.
def solve_army(UNITS, DATA, RESOURCES):
  # Create the linear solver using the CBC backend
  solver = pywraplp.Solver('Minimize resource consumption', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

  # 1. Create the variables we want to optimize
  units = [solver.IntVar(0, solver.infinity(), unit) for unit in UNITS]

  # 2. Add constraints for each resource
  for r, _ in enumerate(RESOURCES):
    solver.Add(sum((10 * DATA[u][-2] + DATA[u][-1]) * units[u] for u, _ in enumerate(units)) >= 1000001)

  # Old constraints for limited resources
  for r, _ in enumerate(RESOURCES):
    solver.Add(sum(DATA[u][r] * units[u] for u, _ in enumerate(units)) <= RESOURCES[r])

  # 3. Minimize the objective function
  solver.Minimize(sum((DATA[u][0] + DATA[u][1] + DATA[u][2]) * units[u] for u, _ in enumerate(units)))

  # Solve problem
  status = solver.Solve()

  # If an optimal solution has been found, print results
  if status == pywraplp.Solver.OPTIMAL:
    print('================= Solution =================')
    print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
    print()

    power = sum((10 * DATA[u][-2] + DATA[u][-1]) * units[u].solution_value() for u, _ in enumerate(units))
    print(f'Optimal value = {solver.Objective().Value()} 🌾🪵🪙resources')
    print(f'Power = 💪{power}')
    print('Army:')
    for u, _ in enumerate(units):
      print(f' - {units[u].name()} = {units[u].solution_value()}')
    print()
    
    food = sum((DATA[u][0]) * units[u].solution_value() for u, _ in enumerate(units))
    wood = sum((DATA[u][1]) * units[u].solution_value() for u, _ in enumerate(units))
    gold = sum((DATA[u][2]) * units[u].solution_value() for u, _ in enumerate(units))
    print('Resources:')
    print(f' - 🌾Food = {food}')
    print(f' - 🪵Wood = {wood}')
    print(f' - 🪙Gold = {gold}')
  else:
      print('The solver could not find an optimal solution.')

solve_army(UNITS, DATA, RESOURCES)
================= Solution =================
Solved in 28.00 milliseconds in 1 iterations

Optimal value = 172100.0 🌾🪵🪙resources
Power = 💪1000105.0
Army:
 - 🗡️Swordsmen = 1.0
 - 🛡️Men-at-arms = 681.0
 - 🏹Bowmen = 0.0
 - ❌Crossbowmen = 0.0
 - 🔫Handcannoneers = 0.0
 - 🐎Horsemen = 0.0
 - ♞Knights = 0.0
 - 🐏Battering rams = 301.0
 - 🎯Springalds = 0.0
 - 🪨Mangonels = 0.0

Resources:
 - 🌾Food = 68160.0
 - 🪵Wood = 90320.0
 - 🪙Gold = 13620.0

Since we now have a limited resource of 🪵wood, the number of 🐏battering rams sadly dropped from 371 to 301. In exchange, we got 681 🛡️men-at-arms and 1 lost 🗡️swordsman (welcome to them).

The total cost of the army is 172,100, which is much higher than the 111,300 we previously found (+65% increase) but it truly is the optimal solution under these constraints. It shows that we should produce more wood because these 🐏 battering rams are extremely cost-efficient!

This example shows how modular LP models can be. It is possible to reuse parts of the code, like constraints, in another model to combine them and solve more complex problems.

🧠 IV. Linear Programming vs Machine Learning

Let’s talk about the elephant in the room. Why not use machine learning (in a broad sense) instead of linear programming? It’s not like this problem cannot be solved with a genetic algorithm for instance.

Mathematical optimization is often neglected in favor of machine learning techniques, but both have their merits:

  • Linear programming can produce an optimal solution in an undetermined amount of time (it can take years), while machine learning can approximate complex functions in no time.
  • There is no training in LP, but an expert is required to build a mathematical model. Machine learning needs data, but the models can be used as black boxes to solve a problem.
  • As a rule of thumb, problems that do not have a particular time constraint and/or are not extremely complex can be advantageously solved with linear programming.

Conclusion

In this tutorial, we dived deeper into our understanding of mathematical optimization.

  • We talked about solvers and types of optimization problems: LP, MIP, NLP;
  • We modeled and solved an extremely common optimization problem in an optimal way and generalized our model through a function;
  • We reframed this problem and merged two sets of constraints to obtain the best army composition for the lowest price;
  • We compared the pros and cons of linear programming and machine learning.

There are a lot more problems where optimization can be applied. For instance, how to create school timetables that satisfy everybody’s requirements? How to deliver 1,000 different orders in a minimum amount of time? Where to create a new metro line to maximize its usefulness?

In future articles, we’ll talk about new types of applications for these techniques, including satisfiability and nonlinear problems.

I hope you enjoyed this more advanced article. If you like machine learning and optimization, let’s connect on Twitter!

🥇 Linear Programming Course

🔎 Course overview

📝 Chapter 1: Introduction to Linear Programming

📝 Chapter 2: Integer vs. Linear Programming

📝 Chapter 3: Constraint Programming

📝 Chapter 4: Nonlinear Programming for Marketing Budget Allocation