Last year, I was approached by a friend who works in a small, family-owned steel and metal business. He wanted to know if it was possible to create something that would help him solve the problem of minimising waste when cutting steel beams. Sounds like a problem for linear programming!
When I started out, there was not a huge amount of beginners articles on how to use Linear Programming in R that made sense for somebody not that versed in math. Linear programming with R was also an area where ChatGPT did not shine in early 2023, so I found myself wishing for a guide.
This series is my attempt at making such a guide. Hopefully it will be of use to someone.
This is part II of the series, if you are looking for a introduction to linear programming in R, have a look at part I.
But I don’t get the math..!?
If you read the theory behind linear programming, or linear optimisation, you end up with a lot of math. This can be off-putting if you don’t have a math background (I don’t). For me, it’s mostly because I never took enough math classes to understand a lot of the symbols. Initially, this made understanding the tutorials surrounding linear programming harder than it should have been. However, you don’t need to understand the math behind the theory to apply the principles of code in this article.

Defining the problem
For the rest of the article, I will have two subsections for each part; one where I try to explain the mathematical expressions, and how it relates to the code. The other will focus only on the code.
Hopefully, this will help – and if you only want the code, skip the math sections.
In this case, my friend was responsible for production planning, and had noticed that there was often a lot of waste after cuts. Planning production manually was time consuming, and together we came up with some problem definitions:
- Each work item in the production plan must be cut exactly once (to avoid having to weld together smaller cuts)
- The sum of lengths to be cut must not exceed the length available in the inventory
- Selected beams must match the dimension in the work order
- Selected beams must match the surface treatment in the work order
For the purpose of the article, I will simply some of the data structures, and provide mock code at the end so that you can create the model. The work items were kept in a data frame called workTable
, and similarly the inventory was kept in a data frame called inventoryTable
.
Some additional helper variables were defined to make the code a bit clearer, such as numberOfWorkItems
and numberOfInventoryItems
. Another useful variable is called the Big M number. This is a technique enabling logical conditions into MILP models. Choosing the Big M was done by the following code:
bigMNumber <- sum(workTable$Length)
The role of the Big M will be more clear once we get to the code, but we use it in this model to "activate" a steel beam, so that the model knows when to take a new inventory item.
The parameters (math)
The parameters of the model are as follows (again, an image as Medium isn’t very good at formatting mathematical expressions):

The decision variables (math)

The constraints (math)
And the constraints can be expressed as follows:

The parameters (code)
Above, we defined the problem formulation with math. Now, let’s do it with code. Below is a simplified version of the data I used, including mock work and inventory.
# Create mock data
inventoryTable <- data.frame(
ID = sapply(1:10, function(x) paste(sample(0:9, 5, replace = TRUE), collapse = "")),
Dimension = c(rep("HUP 120x120x5", 5), rep("HUP 150x150x5", 5)),
Length = c(1000, 2000, 3000, 4000, 5000, 1000, 2000, 3000, 4000, 5000),
Surface = c(rep("Black", 5), rep("Primed", 5))
)
workTable <- data.frame(
JobName = c("Job 1", "Job 2", "Job 3", "Job 4", "Job 5"),
Dimension = c("HUP 120x120x5", "HUP 120x120x5", "HUP 150x150x5", "HUP 150x150x5", "HUP 150x150x5"),
Surface = c("Black", "Black", "Primed", "Primed", "Primed"),
Length = c(500, 1000, 1500, 750, 2000)
)
The code above creates an inventory table with 10 beams of varying length, while the work table contains five jobs that we need to plan.
# Define model parameters
numberOfWorkItems <- nrow(workTable)
numberOfInventoryItems <- nrow(inventoryTable)
bigMNumber <- sum(as.numeric(workTable$Length)) * 1.1
Decision variables (code)
Next, we start defining the model with our decision variables and the constraints, as well as the objective function. The entire code for the model will be available at the end of the article, below is the code step by step:
library(ompr)
# Define the model
MIPModel() |>
# Add binary decision variable modeling which work item
# is cut from which inventory item
add_variable(cutSteel[work, inventory],
work = 1:numberOfWorkItems,
inventory = 1:numberOfInventoryItems,
type = "binary"
) |>
# Add binary decision variable: Take new item from inventory?
add_variable(takeItem[inventory],
inventory = 1:numberOfInventoryItems,
type = "binary"
) |>
# Add binary variable:
# Does the type of work item and inventory item match?
add_variable(typeMatch[work, inventory],
work = 1:numberOfWorkItems,
inventory = 1:numberOfInventoryItems,
type = "binary"
) |>
# Add binary variable:
# Does the surface of work item and inventory item match?
add_variable(surfaceMatch[work, inventory],
work = 1:numberOfWorkItems,
inventory = 1:numberOfInventoryItems,
type = "binary"
)
Our first decision variable models which work item is cut from which inventory item. The indices work represent the work items to be cut, and inventory the available beams to cut from. The decision variable is binary, which means that cutSteel[1,1] = 1 if work item 1 is cut from inventory item 1. This variable allows us to track which work item is assigned to which inventory item, and will be cruicial later on.
The second variable keeps track of whether or not we need to take another beam from the inventory. It is a one-dimensional, and the index is only inventory.
The third and fourth decision variables are there to check if the dimension and surface treatment of the selected work item matches that of the selected inventory item.
Constraints (code)
In total, we have six constraints:
# Constraint 1: Each item in workTable must be cut exactly once
add_constraint(
sum_over(cutSteel[work, inventory],
inventory = 1:numberOfInventoryItems
) == 1,
work = 1:numberOfWorkItems
) |>
# Constraint 2: The sum of each item used to cut from needs to be equal to or smaller than the work item #nolint
add_constraint(
sum_over(cutSteel[work, inventory] * workTable$Length[work],
work = 1:numberOfWorkItems
) <= inventoryTable$Length[inventory],
inventory = 1:numberOfInventoryItems
) |>
# Constraint 3: bigMNumber constraint to activate takeItem whenever a length is cut #nolint
add_constraint(
sum_over(cutSteel[work, inventory],
work = 1:numberOfWorkItems
) <= bigMNumber * takeItem[inventory],
inventory = 1:numberOfInventoryItems
) |>
# Constraint 4: typeMatch can only be 1 if types match
add_constraint(
typeMatch[work, inventory] == (workTable$Dimension[work] == inventoryTable$Dimension[inventory]), # nolint
work = 1:numberOfWorkItems, inventory = 1:numberOfInventoryItems
) |>
# Constraint 5: surfaceMatch can only be 1 if surfaces match
add_constraint(
surfaceMatch[work, inventory] == (workTable$Surface[work] == inventoryTable$Surface[inventory]), # nolint
work = 1:numberOfWorkItems, inventory = 1:numberOfInventoryItems
) |>
# Constraint 6: cutSteel can only be 1 if typeMatch and surfaceMatch are 1 #nolint
add_constraint(
cutSteel[work, inventory] <= typeMatch[work, inventory],
work = 1:numberOfWorkItems, inventory = 1:numberOfInventoryItems
) |>
add_constraint(
cutSteel[work, inventory] <= surfaceMatch[work, inventory],
work = 1:numberOfWorkItems, inventory = 1:numberOfInventoryItems
)
The first constraint expresses that each item must be cut from a single piece. I.e., if we have an order for a 4000mm beam, this needs to come from a single beam in the inventory (this was a customer requirement, since they did not want to cut it and then weld it together).
Important: Note the difference in the placement of the definition of the indices inventory and work. In the decision variables, we defined both inside the same parenthesis, but where we define inventory inside the sum_over()
function, while work is defined outside. Why? It matters because inventory is used as a summation index, that is to say sum over all items in inventory. Work however, is used to specify the applicability of the constraint, meaning that the constraint must hold true for all items in work. If you mess up the placement of the definitions, you will spend time debugging and tearing out your hair (or so I’ve heard…).
The second constraint is fairly obvious, and is used to ensure that the selected item in inventory is long enough (hence the <= inventoryTable$Length[inventory]
). Note that since this constraint applies to all items in inventory, that is the index that is defined outside the sum_over()
function.
Our third constraint involves the "Big M" mentioned early in the article. This constraint ensures that if a particular item in inventory is used, that item is ‘taken’ or ‘activated’. The role of the "Big M" number is to have a constant that effectively turns on or off a condition. This is a technique enabling logical conditions into MILP models. The way it works is that as long as takeItem[inventory]
is 0, the right-hand side is 0, and therefore the sum of cutSteel
for that item must be 0 (i.e., it is not activated). When takeItem[inventory]
is 1, the right-hand side becomes bigNumber
which is sufficiently big so that sum of cutSteel
for that item is practically unbounded. This allows for the constraint to act as a logical on/off switch.
The fourth, fifth and sixth constraints serve the same purpose – making sure that there is a match between the dimension and surface treatment for the item in the work queue, and the selected beam from the inventory. On the surface of it, the sixth constraint might seem unnecessary, but in reality it acts as a logical AND condition.
The objective function (math)

The problem formulation has two parts, so let’s go through them one by one:

In short, this means calculate the total length of all the inventory items that are used. We’ll go through it step by step, but let’s start with an example.
Say that we have three beams that need to be cut, of lengths 1000mm, 1500mm and 2000mm (this is our production plan, or W). We have beams that are of length 3000mm, 4000mm and 2000mm (this is our inventory, or I). This rather simple example (not taking into account dimensions, surface treatment or the size of the actual cut) is easy to solve, and we can see that the optimal solution would be to use the 3000mm beam to cut the beams that are 1000mm and 2000mm, and the 1500mm order can be fulfilled by the beam that is 2000mm.
We start by translating some of the symbols that might be unfamiliar:
- The large strange looking E, means "Sum"
- Below the sum symbol, there are three symbols, i, ∈ and I. ∈ indicates a set membership, and can be interpreted as "is an element of". The small i means an item (a beam) and the I represents the set of inventory items.
- takeItem is a binary decision variable, and if an item (i) is taken out of the inventory to be used to cut a work item, it is equal to 1, otherwise it is 0.
- The lowercase l means lenght, and the subscript lowercase i means the lenght of item i. In the example above, this would for instance be 2000mm for one of the beams that was taken out of the inventory.
All in all, this part of the formula can be interpreted as:
Calculate the total length of all the inventory items that are used. In our example above, if the inventory beams are denoted A (3000mm), B (4000mm) and C (2000mm), we would use A and C, which means that the result of the summation would be 6000mm (4000 x 1 + 2000 x1).

The second part of the problem formulation is similar to the first. It calculates the total length of all the beams that are actually cut. The binary decision variable cutSteel dictates if the length of of the item is included or not.
Bringing it together, the two parts first calculate the total length of the beams used from the inventory, and then subtract the length of steel cut. We are then left with some material, and we want to minimise this amount.
The objective function (code)
# Set objective function to minimize scrap / waste
set_objective(
sum_over(
takeItem[inventory] * inventoryTable$Length[inventory],
inventory = 1:numberOfInventoryItems
)
- sum_over(cutSteel[work, inventory] * workTable$Length[work],
work = 1:numberOfWorkItems,
inventory = 1:numberOfInventoryItems
),
sense = "min"
)
If you read the math-part of the objective function, the code will be familiar. If not, the code above has two summations:
First, it summarises the total length of the beams used from the inventory. For each item in the inventory where takeInventory
equals 1 (i.e. the item is used), we multiply it with it’s length. So if a beam that is 3000mm long is taken, we get 1 * 3000 = 3000
. The index 1:numberOfInventoryItems
indicates that this should be done for all items in the inventoryTable
. For items that are not taken, again using a 3000mm beam, the formula would be 0 * 3000 = 0
.
The second part does the same, but this time for the decision variable cutSteel
. This variable indicates which beam in the production plan is cut from which inventory item. Since it too is binary, if beam 1 from the production plan is cut from beam 4 in the inventory cutSteel[1,4]
, the formula would (if beam 1in the production is 4000mm) 1 * 4000 = 4000
. Note that while we define both the work
and inventory
indices, it is the item in the production plan whos length we are interested in.
We then subtract the second from the first, and we are left with the excess material from the selected beams. The parameter sense = "min"
indicates the direction of the optimalisation. In other words, we want to minimize the amount of leftover material.
Given two configurations that both meet the constraints:
Option 1: Cut beams of 1000mm (A), 2000mm (B) and 4000mm (C) from an inventory of beams that are 1000mm (1), 20000mm (2), 3000mm (3), 5000mm (4) and 6000mm (5). We cut A from 1, B from 2 and C from 5. This gives us the total sum of length of items selected (the first part of the objective function): 1000 + 2000 + 6000 = 9000
. The total sum of cut material (the second part of the objective function) is 1000+2000+4000 = 7000
. The sum of both will then be 9000 - 7000 = 2000
.
Option 2: Cut beams of 1000mm (A), 2000mm (B) and 4000mm © from an inventory of beams that are 1000mm (1), 20000mm (2), 3000mm (3), 5000mm (4) and 6000mm (5). We cut A from 1, B from 2 and C from 4. This gives us the total sum of length of items selected (the first part of the objective function): 1000 + 2000 + 5000 = 8000
. The total sum of cut material (the second part of the objective function) is 1000+2000+4000 = 7000
. The sum of both will then be 9000 - 7000 = 1000
.
Given that we want to minimise the leftover material, our optimalisation engine would select option 2 as the best.
Solving the model (code only)
Once we have our model, we can pass it on to a solver. As mentioned in the first part, there are various commercial solvers available, but for this we will still use the freely available "glpk" solver. Given a model called model
, the following gets us the result and the solution:
# Solve the model
result <- solve_model(model, with_ROI(solver = "glpk"))
# Get the solution
solution <- get_solution(result, cutSteel[work, inventory])
The resulting data.frame from get_solution()
can be a bit confusing to read, since it can become quite large. We define the cutSteel
variable as a matrix, so the resulting data.frame from a matrix of workTable
(5 rows) and inventoryTable
(10 rows) will be 50 rows. For each item in workTable
, the engine checks if it should be cut from any of the 10 items in the inventory. Since we are interested in which items got cut from which inventory items, we can filter the ones that were not used by adding the following code:
# Solve the model
result <- solve_model(model, with_ROI(solver = "glpk"))
# Get the solution
solution <- get_solution(result, cutSteel[work, inventory]) |>
filter(value > 0)
This gets us the following 5 row data.frame:
variable work inventory value
1 cutSteel 1 2 1
2 cutSteel 2 2 1
3 cutSteel 3 10 1
4 cutSteel 4 10 1
5 cutSteel 5 10 1
This is still somewhat clunky to read for a human, and we’ll clean it up momentarily. However, it can be read as follows:
- Work items 1 and 2 were cut from inventory item 2
- Work items 3, 4 and 5 were cut from inventory item 10
To create a somewhat nicer looking production plan that adds some useful information, we can do the following:
# Get the solution
solution <- get_solution(result, cutSteel[work, inventory]) |>
filter(value > 0) |>
mutate(cutLength = workTable$Length[work])
First we add a variable called cutLength
to the solution, so that we know how much was cut.
Next, we add back in some of the information from workTable
and inventoryTable
using dplyr
functions mutate()
and select()
:
productionPlan <- solution |>
select(-c(variable, value)) |>
arrange(work) |>
mutate(
JobName = workTable$JobName[work],
Dimension = workTable$Dimension[work],
Surface = workTable$Surface[work],
Length = workTable$Length[work],
InventoryID = inventoryTable$ID[inventory],
`Start length` = inventoryTable$Length[inventory],
`Cut length` = cutLength,
`Leftover` = inventoryTable$Length[inventory] - cutLength,
) |>
select(-c(work, inventory, cutLength))
This gives us the following:
JobName Dimension Surface Length InventoryID Start length Cut length Leftover
1 Job 1 HUP 120x120x5 Black 500 84508 2000 500 1500
2 Job 2 HUP 120x120x5 Black 1000 84508 2000 1000 1000
3 Job 3 HUP 150x150x5 Primed 1500 48222 5000 1500 3500
4 Job 4 HUP 150x150x5 Primed 750 48222 5000 750 4250
5 Job 5 HUP 150x150x5 Primed 2000 48222 5000 2000 3000
We can now check which job needs to be cut from which item (note that the inventoryID is randomly generated, so if you use the code supplied to genererate, the IDs will differ but the selected material is the same). What other information you want to supply is up to you, in this case I have given some examples. For instance, the production plan supplies information about how much length was available at the start, how much is cut, and how much was left after the cut. Notice that for beams that are cut multiple times, the "leftover" changes with each cut.
When adding back information like this, you use the same matrices that you have supplied in the model. This allows us to get the name of each job using the [work]
and [inventory]
indices respectively.
Summary
Part II of the series was intended to give a more detailed introduction to linear programming, using a real life example. I’ve touched upon the math used, but as I said initially, you don’t really need to understand or be able to read the math to code the solution.
That being said, if you plan on spending some time doing linear programming, I would recommend brushing up on it, at least enough to learn to read the formulas, as a lot of the theory and documentation will reference it.
The code used in this article was written as part of the project (with some adaptations to the article), but was used in a Shiny app. The app is a lot more complicated though, since it also deals with real life problems such as planning jobs, reservering material, releasing planned material back into the inventory, placement etc. Should there be enough interest, I might write up an article about the experience, but we will see.
Hopefully, this has been an interesting second installment. If you have any questions or comments, feel free to reach out.
As promised, the entire code can be found on my GitHub.