Trade-off table, blocking, and the first move toward optimization#

Source worksheets: yint.org/w8 and yint.org/w9 - weeks 8 and 9 of the applied DoE short course.

Module 5 ended with the defining relation and the resolution of a fractional design. This module zooms out to the trade-off table

  • the standard reference for picking k factors and p generators when budget is tight - then works through a small case study that combines a fractional factorial with blocking on the experimenter, and finishes by introducing the sequential experimentation mindset that drives Modules 7 and 8.

Tip

The trade-off table is the most frequently consulted quick reference in fractional-factorial design. The full table is at yint.org/tradeoff.

Q1-Q8 - reading the trade-off table#

A few quick lookups:

  • Five factors, sixteen runs: \(k = 5\), \(p = 1\). One generator; conventionally \(E = ABCD\). Defining relation \(I = ABCDE\). Resolution V.

  • Six factors, twenty runs: 20 is not a power of 2. The nearest workable designs are 16 or 32 runs. At \(k = 6\), \(p = 2\) (so \(2^{6-2} = 16\)) the standard generators are \(E = ABC\), \(F = BCD\); the defining relation has \(2^{2} - 1 = 3\) words.

The general pattern:

Each generator costs one half of the runs. \(p\) generators \(\rightarrow\) \(2^{k-p}\) runs \(\rightarrow\) \(2^{p} - 1\) words in the defining relation \(\rightarrow\) \(2^{p} - 1\) alias chains.

In Module 5 we saw \(p = 1\) (half-fraction, one word \(ABCD\)). With \(p = 2\) (quarter-fraction) the defining relation has three words: the two generators and their product.

[1]:

# 6 factors, 16 runs (k=6, p=2), generators E=ABC and F=BCD

generators = ["ABCE", "BCDF"]

# defining relation = I = ABCE = BCDF = ABCE * BCDF = (BC)^2 ADEF = ADEF

print(f"Generators (as words):   {generators}")
print(f"Defining relation: I = {' = '.join(generators)} = ADEF (= ABCE * BCDF, after cancellation)")

Generators (as words):   ['ABCE', 'BCDF']
Defining relation: I = ABCE = BCDF = ADEF (= ABCE * BCDF, after cancellation)

Solution

  • Q3. E = ABCD in the 2^(5-1) design, so the alias of E is the full four-factor interaction ABCD (all but invisible in a real chemistry system).

  • Q4. k = 6, p = 2, generators E = ABC and F = BCD give a resolution-IV quarter-fraction in 16 runs.

  • Q7-Q8. With p = 2 generators the defining relation has 2^p - 1 = 3 words: each generator plus their product.

Q9-Q13 - half-fraction in four factors, plus blocking on the experimenter#

The set-up: four factors A, B, C, D to study; only 8 runs in the budget; the work must be split between you and a colleague because the schedule does not allow one person to do all 8.

A natural design is a half-fraction with generator D = ABC (resolution IV). Splitting the work between two experimenters is a blocking problem: the experimenter is a nuisance factor - we do not care whether the value is “you” or “colleague”, but if we do not control for it any drift in technique gets blamed on the real factors. Block on it by treating “Person” as an extra column E whose generator is some interaction we are willing to lose.

The standard choice is E = AB (block confounded with the AB two-factor interaction). Two design “words” gives us a resolution-III design overall: every main effect is now confounded with a two-factor interaction.

Worksheet response values (Q13): y = [120, 76, 106, 90, 72, 74, 90, 55] in standard order.

[2]:

from process_improve.experiments import c, gather, lm

A = c(-1, +1, -1, +1, -1, +1, -1, +1, name="A")
B = c(-1, -1, +1, +1, -1, -1, +1, +1, name="B")
C = c(-1, -1, -1, -1, +1, +1, +1, +1, name="C")

# Generator D = A*B*C

D = c(-1, +1, +1, -1, +1, -1, -1, +1, name="D")

# Block factor E = "Person", with E = A*B

E = c(+1, -1, -1, +1, +1, -1, -1, +1, name="E")
y = c(120, 76, 106, 90, 72, 74, 90, 55, name="y")
design = gather(A=A, B=B, C=C, D=D, E=E, y=y)
design

[2]:
A B C D E y
1 -1.0 -1.0 -1.0 -1.0 1.0 120.0
2 1.0 -1.0 -1.0 1.0 -1.0 76.0
3 -1.0 1.0 -1.0 1.0 -1.0 106.0
4 1.0 1.0 -1.0 -1.0 1.0 90.0
5 -1.0 -1.0 1.0 1.0 1.0 72.0
6 1.0 -1.0 1.0 -1.0 -1.0 74.0
7 -1.0 1.0 1.0 -1.0 -1.0 90.0
8 1.0 1.0 1.0 1.0 1.0 55.0 
[3]:

# Fit main effects plus the one two-factor interaction (A:C) that is
# not already absorbed by a generator.

m = lm("y ~ A + B + C + D + E + A:C", design)
print(m.get_parameters(drop_intercept=False).to_string())

Intercept    85.375
A           -11.625
B            -0.125
C           -12.625
D            -8.125
E            -1.125
A:C           3.375

Solution

With both generators D = ABC and E = AB the defining relation is

I = ABCD = ABE = ABCD * ABE = CDE.

Words: ABCD, ABE, CDE. Shortest word has length 3 (ABE, CDE), so the design is Resolution III.

The largest coefficients:

  • A = -11.6, C = -12.6, D = -8.1 - three large, all negative.

  • B and E are tiny: the blocking column E showing zero is exactly what you want from a nuisance factor.

  • A:C is modest but the only two-factor interaction we kept in the model.

Confounds to be aware of:

  • A is aliased with BE (multiply ABE by A) and with BCD (multiply ABCD by A). The 11.6 coefficient could in principle be b_A + b_BE + b_BCD.

  • E (“Person”) is aliased with AB and with CD. E coming out tiny tells you the operator effect and those two two-factor interactions together are negligible.

Guidance

Blocking is the practical answer to “the experimenter or batch probably affects the response.” You always know which block a run belongs to (which operator ran it, which day, which batch of reagent), so you can always model it. Picking which interaction to confound with the block is a deliberate choice: lose the one you are least likely to believe in.

Sequential experimentation - the bridge to optimization (w9)#

A response surface study is not one big design; it is a sequence of small designs, each pointing the next one in the right direction. The loop:

  1. Predict the outcome of the next experiment with the current model.

  2. Run the experiment.

  3. Compare the prediction with the measurement.

    • If they agree, the model is still useful. Decide:

      • extend the model with the new point and keep climbing, or

      • stop, because you are at the optimum.

    • If they disagree, the model has broken down. Switch to a higher-order model (add a center point, then axial points for a quadratic) or shift the design region.

  4. Plan where the next experiment will be.

  5. Repeat.

The mantra: the model is useful, but wrong - useful because it gives you a defensible direction, wrong because the true surface curves and the linear approximation will break eventually. Knowing when it breaks is the entire skill.

Check yourself

Before reading the next paragraph: in a 1-D problem (one factor), when is one-factor-at-a-time (OFAT) a good idea, and when is it a bad idea?

Solution

OFAT is fine when the system genuinely has one input - for example, finding the temperature that maximizes crystal size in a single-variable bioreactor process. Sequential OFAT is exactly what response surface optimization does in 1-D.

OFAT is bad the moment the system has two or more inputs that interact. Tweaking factor A while holding B fixed and then tweaking B while holding A at its OFAT optimum is guaranteed to miss the diagonal of the response surface, which is where interactions live. Module 1 already showed this on the yogurt example: the four corners of a 2x2 told you what no OFAT walk ever would.

A worked 1-D optimization#

To make the response-surface loop concrete, we will optimize a 1-D process by hand. Pick a single factor (think: temperature, mixer speed, or any other physical knob), take small steps, fit a linear model, swap to a quadratic when the surface starts to curve, predict the peak, and confirm.

The “true” system below is hidden from the optimizer; only the function call observe(t) is allowed - it returns the measured response at temperature t with a sprinkle of measurement noise.

[4]:

import numpy as np
import plotly.graph_objects as go
import plotly.io as pio

pio.renderers.default = "notebook_connected"

rng = np.random.default_rng(seed=42)


def observe(temperature_C: float) -> float:
    """Return a noisy measurement: true peak at 60 degC, sigma = 1."""
    true = -0.05 * (temperature_C - 60) ** 2 + 75.0
    return float(true + rng.normal(scale=1.0))


# First three experiments: anchor a linear trend.

xs = [40.0, 50.0, 60.0]
ys = [observe(t) for t in xs]
for t, y in zip(xs, ys, strict=True):
    print(f"  t = {t:5.1f} degC -> measured y = {y:.2f}")

  t =  40.0 degC -> measured y = 55.30
  t =  50.0 degC -> measured y = 68.96
  t =  60.0 degC -> measured y = 75.75

[5]:

# Fit a simple linear model on the first three points - is the surface
# still going up?

slope, intercept = np.polyfit(xs, ys, deg=1)
print(f"Linear fit (3 points): y = {intercept:.2f} + {slope:.3f} * t")
print("Slope is positive -> step further uphill.")

# Step uphill by +10 degC.

next_t = max(xs) + 10
xs.append(next_t)
ys.append(observe(next_t))
print(f"\nNext run at t = {next_t} degC -> y = {ys[-1]:.2f}")

Linear fit (3 points): y = 15.56 + 1.022 * t
Slope is positive -> step further uphill.

Next run at t = 70.0 degC -> y = 70.94

[6]:

# Add a couple more steps; switch to a quadratic fit once we have 5
# points so curvature can show up.

for next_t in [70.0, 80.0]:
    xs.append(next_t)
    ys.append(observe(next_t))

coefs = np.polyfit(xs, ys, deg=2)
poly = np.poly1d(coefs)
predicted_peak = -coefs[1] / (2 * coefs[0])
print(f"Quadratic fit: y = {coefs[0]:.4f}*t^2 + {coefs[1]:.3f}*t + {coefs[2]:.2f}")
print(f"Predicted peak temperature: {predicted_peak:.1f} degC")
print(f"Predicted response at peak: {poly(predicted_peak):.2f}")

Quadratic fit: y = -0.0514*t^2 + 6.139*t + -108.45
Predicted peak temperature: 59.7 degC
Predicted response at peak: 74.92

[7]:

# Confirm the predicted optimum with a fresh run, then plot the journey.

xs.append(predicted_peak)
ys.append(observe(predicted_peak))
print(f"Confirmation run at t = {predicted_peak:.1f} degC -> y = {ys[-1]:.2f}")

grid = np.linspace(35, 85, 200)
fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=ys, mode="markers+text",
                         text=[f"{i+1}" for i in range(len(xs))],
                         textposition="top center",
                         name="Observations", marker={"size": 9}))
fig.add_trace(go.Scatter(x=grid, y=poly(grid), mode="lines",
                         name="Quadratic fit",
                         line={"dash": "dash"}))
fig.update_layout(width=720, height=420,
                  title="1-D sequential optimization",
                  xaxis_title="Temperature [degC]",
                  yaxis_title="Response y")
fig

Confirmation run at t = 59.7 degC -> y = 75.12

Solution

The sequence above is exactly the steepest-ascent strategy:

  1. Anchor with a small linear study (3 points).

  2. As long as the slope is positive, step further uphill.

  3. The moment the surface starts curving, switch to a quadratic model so the peak can be predicted.

  4. Confirm the predicted optimum with one more run. If it matches, you are at the optimum (within noise). If it does not, the quadratic was still wrong and the search continues.

For the simulated system above the true peak is at 60 degC. The quadratic fit on the first five points predicts a peak within a degree or two of that, depending on the noise.

Guidance

The general loop (“predict -> run -> compare -> decide”) is the single most important habit in DoE-driven optimization. Module 7 generalizes it to two factors. The same loop runs in 4 or 5 dimensions; the only thing that changes is the visualization.

Wrap-up#

Three transferable points from this module:

  • The trade-off table sets the budget conversation. Pick the resolution you can afford, name the generators, write the defining relation, and then run the experiment.

  • Blocking is how you keep nuisance factors honest. Confound them with the highest-order interaction you are willing to lose.

  • Optimization is sequential. No single design lands you on the peak; a string of small designs does.

Next: Module 7 takes the same loop into two factors with the classical response-surface trio of factorial + center + axial points.