Teaching Factory Physics Flow Benchmarking with R and Many-Objective Visuals

[This article was first published on R Blogs by Pedro N. de Lima, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Teaching to seasoned managers in MBE classes is challenging. While it’s important to bring new thoughts and ideas and not sound repetitive, it is necessary to provide a theoretical basis for experienced people with diverse backgrounds. One of the strategies I found to overcome these obstacles this week was to use a new analysis framework (in my case, Factory Physics concepts) to challenge their views about existing frames they already master. Using a combination of concepts, simulation sodels, many-objective tradeoffs visuals (like the gif below) and tasks was great to challenge their intuition about manufacturing systems.

This post shares some of the R code I developed while putting together course materials. This post is also an example of how to use simulation metamodeling and Arena to Run Factory Physics’ Flow Benchmarking.

Flow Benchmarking

Flow Benchmarking is an absolute benchmarking technique useful to determine how close a production flow is to its best possible performance. The technique has been introduced in the award-winning Factory Physics (FP) Book (Hopp and Spearman 2008), and is a key component of the science-based manufacturing management approach described in (Pound, Bell, and Spearman 2014).

Defining Factory Physics Laws

The Flow Benchmarking analysis is grounded on Little’s Law (WIP = TH * CT), and utilizes three general cases as absolute benchmarks for any real manufacturing system: The Best Case, the Worst Case and the Practical Worst Case .Please refer to (Hopp and Spearman 2008) and (Pound, Bell, and Spearman 2014) for the rationale for these laws and equations.

I’ll define these equations as R functions:

calc_w0 = function(rb, t0) {rb * t0}

ct_best = function(t0, w, w0, rb) {ifelse(w<=w0,t0,w/rb)}

th_best = function(t0, w, w0, rb) {ifelse(w<=w0,w/t0,rb)}

ct_worst = function(w,t0){w*t0}

th_worst = function(t0){1/t0}

ct_marginal = function(t0,w,rb){t0+(w-1)/rb}

th_marginal = function(w0,w,rb){rb*w/(w0+w-1)}

In summary, these equations provide a starting point to discuss how well a manufacturing system is doing in terms of converting inventory to Throughput. The initial analysis requires two inputs. The first input is the Bottleneck rate (rb), which is the production rate (parts, orders / time) of the bottleneck (defined as the process center with the highest long-term utilization). The second parameter is the Total Raw Processing Time (t0), which is the sum of the long-term average process times of each processing center. Based on these two parameters, it’s possible to draw benchmarking curves for the System’s Throughput and Cycle Time as a function of its Work in Process, assuming a CONWIP control system (SPEARMAN, WOODRUFF, and HOPP 1990).

Drawing Absolute Benchmarking Curves

Once I have the basic laws of manufacturing dynamics as R functions, I’ll create a benchmarck_flow function to execute the analysis. This function accepts the rb and t0 parameters and will calculate the system’s Throughput and Cycle time as a function of the wip under different scenarios for benchmarking purposes.

## Defining Cycle time and Throughput functions

benchmark_flow = function(rb, t0, step = 1, wip_mult = 5) {
  # First, calculate wip_crit
  w0 = calc_w0(rb = rb, t0 = t0)
  # Then, define WIP Range to consider:
  wip = seq.int(from = 1, to = w0 * wip_mult, by = step)
  # Then, calculate The Best Case Variables
  Best_Cycle_Time = ct_best(t0 = t0, w = wip, w0 = w0, rb = rb)
  Best_Throughput = th_best(t0 = t0, w = wip, w0 = w0, rb = rb)
  best_data = data.frame(WIP = wip,
                    Throughput = Best_Throughput,
                    CycleTime = Best_Cycle_Time,
                    Scenario = "Best Case")
  # Calculate the Marginal Cases:
  Marginal_Cycle_Time = ct_marginal(t0=t0,w=wip,rb=rb)
  Marginal_Throughput = th_marginal(w0=w0,w=wip,rb=rb)
  marginal_data = data.frame(WIP = wip,
                    Throughput = Marginal_Throughput,
                    CycleTime = Marginal_Cycle_Time,
                    Scenario = "Marginal")
  # Calculate Worst Case
  worst_data = data.frame(
    WIP = wip,
    Throughput = th_worst(t0 = t0),
    CycleTime = ct_worst(w = wip, t0 = t0),
    Scenario = "Worst Case"

  # Output A DataFrame with results:
  # I'm not including the Worst Case because it's unrealistic (and messes up my cycle time plot).
  rbind(best_data, marginal_data, worst_data)

# The First Penny Fab Example:
data_benchmark = benchmark_flow(rb = 0.5, t0 = 8)

WIP Throughput CycleTime Scenario
1 0.125 8 Best Case
2 0.250 8 Best Case
3 0.375 8 Best Case
4 0.500 8 Best Case
5 0.500 10 Best Case
6 0.500 12 Best Case

How would the Actual System Behave if…

Ok, now I have a table with all the basic benchmarking results. What if I have a better model of the system? We can accomplish this by building a discrete event simulation model of the actual system, and using a metamodel of this model to approximate its results (you can find the data from my penny fab model here). During my course, I used several Arena Simulation models to illustrate that adding variability to the system always degrades performance (as the variability law predicts!). Doing so allowed the students to build confidence into the model and the theory, which was great to see!

## Warning: package 'arena2r' was built under R version 3.5.2

arena_data = arena2r::get_simulation_results("2019-data/penny-fab")

# Filtering only Statistics of our Interest:

filtered_data = subset(arena_data, Statistic %in% c("w", "LeadTime", "Throughput"))

# Spreading and Data Wrangling

final_data = filtered_data %>% 
  tidyr::spread(Statistic, Value) %>%
  dplyr::select(LeadTime, Throughput, w)

colnames(final_data) = c("CycleTime", "Throughput", "WIP")

# Now, build a spline metamodel for CycleTime and Throughput as a function of WIP.

th_model = lm(Throughput ~ splines::bs(WIP), data = final_data)

ct_model = lm(CycleTime ~ WIP, data = final_data)

# Put Together a Final DataFrame like the Benchmarking:

model_data = data.frame(
  WIP = unique(data_benchmark$WIP),
  Throughput = predict(th_model, subset(data_benchmark, Scenario == "Best Case")),
  CycleTime = predict(ct_model, subset(data_benchmark, Scenario == "Best Case")),
  Scenario = "DES Model"
## Warning in splines::bs(WIP, degree = 3L, knots = numeric(0), Boundary.knots
## = c(1, : some 'x' values beyond boundary knots may cause ill-conditioned
## bases
# Adding our Model's Data to the DataFrame:

data_benchmark = rbind(

Once we have data from the basic FP laws and from our model, let’s plot it!

## Carregando pacotes exigidos: viridisLite
# Lets define a wrapper function for our plot:

plot_benchmarking = function(data) {
  data %>%
    gather(-WIP, -Scenario, key = "var", value = "Value") %>%
  ggplot(aes(x = WIP, y = Value, color = Scenario)) +
    geom_line(size = 1) +
    facet_wrap(~ var, scales = "free", nrow = 2, ncol = 1) +
    labs(title = "Flow Benchmarking Plot") +
    scale_color_viridis(discrete = TRUE, option = "D") + 

# Then let's just benchmark and plot!

plot_benchmarking(data = data_benchmark)

Cool! My simulation metamodel is still quite equivalent to the marginal case.

However, it has one advantage. I can now build a model that can simulate arbitrarily complex scenarios (e,g.: I can include different product routings, change product mix, include non-stationary demand, simulate setup time reduction, even maybe use a multi-method model, etc.) and my model will actually be a better approximation of the actual system than any simple queueing network model. Also, my model can simulate detailed improvement what-if scenarios, which queueing network models won’t be able to simulate.

Wrapping up with Tradeoffs and Many-Objective Visuals

I also used simulation models to illustrate tradeoffs implied by two simple decisions: How much WIP a manufacturing flow should have and what should be the reorder level of a part. Unfortunetly, trying to use R to this task wasn’t productive. I ended up using DiscoveryDV, which is a great tool for many-objective visualization.

For instance, plotting WIP, Throughput, Cycle Time and Utilization of the Practical Worse Case yields this:

And visualizing the tradeoffs implied by different reorder levels in a (Q,r) inventory system yields this:

At this point, many of the participants were excited to get to grips with models that illuminate tradeoffs they have been facing for years. Hopefully, their intuition was sharpened by these exercises and they will be better equiped to use these frontiers to promote productive and tradeoff-aware discussions.


Hopp, W.J., and M.L. Spearman. 2008. Factory Physics. Irwin/Mcgraw-Hill Series in Operations and Decision Sciences. McGraw-Hill. https://books.google.com.br/books?id=tEjkAAAACAAJ.

Pound, E.S., J.H. Bell, and M.L. Spearman. 2014. Factory Physics for Managers: How Leaders Improve Performance in a Post-Lean Six Sigma World. McGraw-Hill Education. https://books.google.com.br/books?id=B5sXAwAAQBAJ.

SPEARMAN, MARK L., DAVID L. WOODRUFF, and WALLACE J. HOPP. 1990. “CONWIP: A Pull Alternative to Kanban.” International Journal of Production Research 28 (5). Taylor & Francis: 879–94. https://doi.org/10.1080/00207549008942761.

To leave a comment for the author, please follow the link and comment on their blog: R Blogs by Pedro N. de Lima.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)