Response surfaces

Why this matters

A two-level factorial tells you which factors matter and in which direction. It cannot tell you where the best setting is if the optimum lies inside the region, because two levels only see straight lines. Response surface methodology (RSM), developed by Box and Wilson in 1951 for chemical processes, is the standard toolkit for the next step: move toward better conditions, then fit a curved model and find the optimum. It is used to tune reaction yields, formulation recipes, machining parameters, fermentation conditions and many other processes.

The sequential strategy

RSM is a campaign, not a single experiment:

1. Screen many factors with a fractional factorial (module 2) and keep the few that matter.
2. Climb. Far from the optimum, the surface is roughly a tilted plane. Fit a first-order model and move along the path of steepest ascent (step sizes proportional to the coefficients) until the response stops improving.
3. Check for curvature with centre points. If the mean at the centre differs clearly from the mean of the factorial corners, you are near a peak (or valley).
4. Fit a second-order model with a design that supports it, and locate the optimum.
5. Confirm with new runs at the predicted best settings.

The second-order model

For two coded factors x1 and x2:

y = b0 + b1*x1 + b2*x2 + b11*x1^2 + b22*x2^2 + b12*x1*x2 + error

Six coefficients. The squared terms give curvature, and the cross term tilts and rotates the surface. You fit it by ordinary least squares exactly as in module 1, with the extra columns added to the model matrix.

Designs for a second-order model

A 2^k factorial cannot estimate the squared terms (every x² column is all ones). Two standard designs can:

• Central composite design (CCD): a 2^k factorial (the corners) plus 2k axial (star) points at distance α along each axis, plus several centre runs. For two factors with α = √2 ≈ 1.414, that is 4 + 4 + (3 to 5) runs. The axial points add the levels needed for curvature; the centre runs give pure error and stabilise the prediction variance.
• Box–Behnken design: three levels per factor, no corner runs. Useful when the extreme corners are unsafe or impossible (for example, maximum temperature with maximum pressure).

CCD, two factors, coded units
 x1      x2
-1      -1        factorial
+1      -1
-1      +1
+1      +1
-1.414   0        axial
+1.414   0
 0     -1.414
 0     +1.414
 0       0        centre (x3)

Coded units make the coefficients comparable and the design symmetric. Convert back to natural units at the end: temp = 75 + 10 * x1 if 65–85 °C is coded −1 to +1.

Finding the optimum

At the top of the surface both partial derivatives are zero:

∂y/∂x1 = b1 + 2*b11*x1 + b12*x2 = 0
∂y/∂x2 = b2 + 2*b22*x2 + b12*x1 = 0

That is a 2 × 2 linear system:

[ 2*b11   b12  ] [x1]     [ b1 ]
[  b12   2*b22 ] [x2] = − [ b2 ]

Its solution is the stationary point. Whether it is a maximum, minimum or saddle depends on the eigenvalues of the matrix B = [[b11, b12/2], [b12/2, b22]]:

• all negative: a maximum;
• all positive: a minimum;
• mixed signs: a saddle (the best point is on the edge of the region; use ridge analysis or constrained optimisation).

Worked example

Suppose the fitted model in coded units is

y = 80 + 4*x1 + 2*x2 - 3*x1^2 - 2*x2^2 + 1*x1*x2

The system is [[-6, 1], [1, -4]] x = -[4, 2]. Solving gives x1 = 18/23 ≈ 0.78, x2 = 16/23 ≈ 0.70. The eigenvalues of B = [[-3, 0.5], [0.5, -2]] are about −3.21 and −1.79, both negative, so it is a maximum. The predicted response there is 80 + 0.5 × (4 × 0.78 + 2 × 0.70) ≈ 82.3 (at a stationary point, ŷ = b0 + ½ bᵀx). The point lies inside the design region (|x| < 1.414), so the prediction is an interpolation and worth confirming with a few runs.

import numpy as np
b1, b2, b11, b22, b12 = 4, 2, -3, -2, 1
x_s = np.linalg.solve([[2 * b11, b12], [b12, 2 * b22]], [-b1, -b2])
print(x_s, np.linalg.eigvalsh([[b11, b12 / 2], [b12 / 2, b22]]))

For contour plots, evaluate the fitted model on a grid (np.meshgrid) and draw it with matplotlib.pyplot.contourf; for several responses at once (yield high, impurity low), overlay the contours or use a desirability function.

Common mistakes

• Fitting a quadratic to a 2^k design without axial points: the squared terms are not estimable.
• Trusting an optimum outside the design region. Extrapolated quadratics are unreliable; run another experiment centred nearer the new point.
• Skipping confirmation runs. The stationary point is a prediction with uncertainty.
• Optimising one response and ignoring others (cost, safety, impurity). Set constraints with the process owner first.
• Forgetting lack-of-fit. With replicated centre points, test whether the quadratic is adequate; if not, transform the response or narrow the region.

Quick recap

• RSM is sequential: screen, climb by steepest ascent, check curvature, fit a quadratic, confirm.
• CCD and Box–Behnken designs support second-order models.
• The stationary point solves a linear system; eigenvalues tell maximum, minimum or saddle.

Try it

In the Do tab you will fit a second-order model to an 11-run central composite design with NumPy and solve for its stationary point.