Optimizing Multiple Responses

Tracking several criteria simultaneously.

Kris Sankaran true
12-02-2021

Readings 11.3, Rmarkdown

  1. In more complicated systems, we may want to optimize several objectives simultaneously. More often than not, the goals will be at odds with one another. For example, in our running chemical process example, we want to maximize yield while maintaining a target viscosity and minimizing molecular weight. How can we use response surface methods when we have several competing objectives?

Overlaying Contours

  1. The most direct approach is to simply fit several response surfaces. Then, we can visually inspect results to find factor configurations with desirable values across each response.

  2. To illustrate, we will work with a chemical process dataset, which varied temperature and time in some chemical process and observed changes in three responses – yield, viscosity, and molecular weight. The block below reads in the data and defines a version, mchem that puts all three responses into one column, using pivot_longer. We standardize all the responses in mchem so that when we plot them below, the colors are all on the same scale.

library(readr)
library(dplyr)
library(tidyr)

chem <- read_csv("https://uwmadison.box.com/shared/static/nbaj1m8j7tuaqmznjlrsgbzyhp9k61i8.csv")

mchem <- chem %>%
  mutate_at(vars(yield:molecular_weight), scale) %>% # standardize
  pivot_longer(yield:molecular_weight, names_to = "response")
  1. The plot below shows how the three responses vary across the input temperature and time. It seems that molecular weight is highest at high times and temperatures. Viscosity is highest at average temperatures, and doesn’t seem to depend on time. Yield is highest at the center point and at high temperature / time.
ggplot(mchem) +
  geom_point(
    aes(time, temp, col = value),
    position = position_jitter(w = 0.3, h = 0.3)
  ) +
  facet_grid(. ~ response) +
  scale_color_viridis_c()

  1. We can improve this view by fitting a response model surface model to each response. Then, we can view the estimated values across all times and temperatures, not just the ones that were sampled. Below, we code the original times and temperatures and fit a second order model associated with each run.
library(rsm)

chem_coded <- coded.data(chem, time_coded ~ (time - 35) / 5, temp_coded ~ (temp - 155) / 5)
fits <- list(
  "yield" = rsm(yield ~ SO(temp_coded, time_coded), data = chem_coded),
  "viscosity" = rsm(viscosity ~ SO(temp_coded, time_coded), data = chem_coded),
  "molecular_weight" = rsm(molecular_weight ~ SO(temp_coded, time_coded), data = chem_coded)
)
  1. Below, we visualize the contours associated with each of the three responses. The conclusions are consistent with the earlier discussion. Moreover, they highlight specific temperatures and times that would be optimal for each of the responses. However, it’s unclear how we would reconcile these three surfaces to decide on a single temperature or time for the process.
contour(fits[[1]], ~ time_coded + temp_coded, image = TRUE, asp = 1)
Three separate response surfaces, fit to yield, viscosity, and molecular weight, respectively.

Figure 1: Three separate response surfaces, fit to yield, viscosity, and molecular weight, respectively. 

contour(fits[[2]], ~ time_coded + temp_coded, image = TRUE, asp = 1)
Three separate response surfaces, fit to yield, viscosity, and molecular weight, respectively.

Figure 2: Three separate response surfaces, fit to yield, viscosity, and molecular weight, respectively. 

contour(fits[[3]], ~ time_coded + temp_coded, image = TRUE, asp = 1)
Three separate response surfaces, fit to yield, viscosity, and molecular weight, respectively.

Figure 3: Three separate response surfaces, fit to yield, viscosity, and molecular weight, respectively.

Constrained Optimization

  1. Whenever visual inspection is challenging, mathematical formalizations can offer support. One idea is to frame the multiple response surface problem as a constrained optimization.
  1. Formally, we look for a configuration of factors \(x_{\ast}\) that solves the optimization \[\begin{align*} \underset{x}{\text{maximize}}\medspace &y_{1}\left(x\right) \\ \text{subject to }\medspace &\left(y_{2}\left(x\right), \dots, y_{R}\left(x\right)\right) \in \mathcal{C} \end{align*}\] where \(C\) is the predefined acceptable region for the secondary responses.

Desirability Functions

  1. The main downside of the constrained optimization approach is that it forces us to choose one response to prioritize over all others. What if we care about each response more or less equally?

  2. One idea is to optimize a sort of (geometric) averaged response, \[\begin{align*} \underset{x}{\text{maximize}}\medspace \left[\prod_{r = 1}^{R} y_{r}\left(x\right)\right]^{\frac{1}{R}} \end{align*}\]

  3. The issue with this idea is that it treats all responses exactly equally. What if we want to maximize some, but minimize others? What if we want some to be near some target value?

  4. The solution is to use desirability functions. A few are plotted below. You can adjust their shape so that the \(r^{th}\) desirability function is large for the values of \(y_{r}\left(x\right)\) which are good (sloping down when you want to minimize, sloping up when you want to maximize).

Example desirability functions, for maximizing, minimizing, and achieving a target response.

Figure 4: Example desirability functions, for maximizing, minimizing, and achieving a target response.

  1. Then, instead of optimizing the raw averaged response, we optimize the averaged response after first passing through the desirability functions, \[\begin{align*} \underset{x}{\text{maximize}} \medspace \left[\prod_{r = 1}^{R} d_{r}\left(y_{r}\left(x\right)\right)\right]^{\frac{1}{R}} \end{align*}\]

Data Example

  1. For the chemical problem, we can define a desirability function per response, and then combine them into an overall objective. This is done using the desirability package.
library(desirability)
d_yield <- dMax(70, 95) # min / max / scale
d_viscosity <- dTarget(55, 65, 75) # min / target / max
d_weight <- dMin(2750, 4000) # min / max
objective <- dOverall(d_yield, d_viscosity, d_weight)
  1. Now, we can apply this objective function to predictions made by the three response surface fits we made above, in the section on overlaying contours. We’re evaluating the predictions over a grid of values in the range of the coded time and temperature.
coded <- as.data.frame(chem_coded)
x_grid <- expand.grid(
  time_coded = seq(min(coded$time_coded), max(coded$time_coded), .1),
  temp_coded = seq(min(coded$temp_coded), max(coded$temp_coded), 0.1)
)

desirabilities <- cbind(
  x_grid,
  y1 = predict(fits[[1]], x_grid),
  y2 = predict(fits[[2]], x_grid),
  y3 = predict(fits[[3]], x_grid)
)
y_hat <- desirabilities %>% 
  select(starts_with("y"))
desirabilities$objective <- predict(objective, y_hat)
  1. At this point, we can plot the desirability of each point in the time / temperature configuration space. Note that a second mode emerges along the low-temperature region – these are regions which have good viscosity and molecular weight properties, even though their yield isn’t as high. Try increasing the scale for yield to see what happens when you make the requirement for high yield more stringent.
ggplot(desirabilities) +
  geom_tile(
    aes(time_coded, temp_coded, fill = objective)
  ) +
  coord_fixed() +
  scale_fill_viridis_c()
Overall desirability, considering yield, viscosity, and molecular weight.

Figure 5: Overall desirability, considering yield, viscosity, and molecular weight.