First-order Model

Response Surface Designs

A response variable 𝑌 depends on several explanatory variables. The goal is to optimize the value of 𝑌 .

𝑌 = 𝑓(𝑥¹, 𝑥²) + 𝜖, function f is unknown.

𝑓(𝑥1, 𝑥2) is called response surface.

ETH – p. 1/16

Maximal yield of a chemical process

The yield 𝑌 of a chemical process depends on reaction time (A) and temperature (B). Current conditions are 35 min. and 155˚C.

2²-factorial:

- 6

r

r r

r

30 ^𝑥1

𝑥²

40 160

150

𝐵

𝐴

𝑥¹ = time ₋ 35 5

𝑥² = temp ₋ 155 5

First-order Model

run A B 𝑥1 𝑥2 𝑦

1 30 150 –1 –1 39.3 2 40 150 +1 –1 40.9 3 30 160 –1 +1 40.0 4 40 160 +1 +1 41.5

First-order model: 𝑌 = 𝛽⁰ + 𝛽¹𝑥¹ + 𝛽²𝑥² + 𝜖

ETH – p. 3/16

R Output

> summary(mod1)

Call: lm(formula = y ˜ x1 + x2)

Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 40.425 0.025 1617 0.000394 ***

x1 0.775 0.025 31 0.020529 *

x2 0.325 0.025 13 0.048875 *

Residual standard error: 0.05 on 1 df

Multiple R-squared: 0.9991, Adj R-squared: 0.9973 F-statistic: 565 on 2 and 1 DF, p-value: 0.02974

First order response

ˆ

𝑦 = 40.425 + 0.775𝑥¹ + 0.325𝑥²

ETH – p. 5/16

Model with interaction

𝑌 = 𝛽⁰ + 𝛽¹𝑥¹ + 𝛽²𝑥² + 𝛽¹²𝑥¹𝑥² + 𝜖 𝛽ˆ¹² = −0.025, perfect fit, no estimation of 𝜎.

Additional runs are necessary to estimate the

experimental error and test the interaction effect.

We add five observations at the center to get an independent estimation of 𝜎 and check the fit of the first-order model or a curvature.

Factorial with center points

r

r r

r

r

30 40

160

150

𝐵

𝐴

A B 𝑥¹ 𝑥² 𝑦

35 155 0 0 40.3 35 155 0 0 40.5 35 155 0 0 40.7 35 155 0 0 40.2 35 155 0 0 40.6

ETH – p. 7/16

Checking Curvature and Interaction

Second-order model

𝑌 = 𝛽⁰ + 𝛽¹𝑥¹ + 𝛽²𝑥² + 𝛽¹²𝑥¹𝑥² + 𝛽¹¹𝑥²₁ + 𝛽²²𝑥²₂ + 𝜖

Estimation of curvature: mean in factorial – mean at center = ¯𝑦_𝑓 − 𝑦¯_𝑐 = 40.425 − 40.46 = −0.035.

¯

𝑦_𝑓 − 𝑦¯_𝑐 is an estimate for 𝛽¹¹ + 𝛽²²

Estimation of 𝜎²:

𝑀 𝑆_𝑟𝑒𝑠 =

∑⁵

(𝑦_𝑖 − 𝑦¯_𝑐)²

4 = 0.043

Can construct confidence intervals for 𝛽¹² and

𝛽¹¹ + 𝛽²². First-order model is adequate.

Method of steepest ascent

ETH – p. 9/16

Example: 𝑦 ˆ = 40 . 44 + 0 . 775 𝑥

+ 0 . 325 𝑥

contour lines 𝑥² = −

0.775

0.325𝑥¹ + 𝑐,

direction of steepest ascent: ⁰₀^._.³²⁵₇₇₅ = 0.42

additional observations:

A B 𝑥₁ 𝑥₂ 𝑦

35 155 0 0

40 157.1 1 0.42 40.5 45 159.2 2 0.84 51.3 50 161.3 3 1.26 59.6 55 163.4 4 1.68 67.1 60 165.5 5 2.10 63.6

Second-order model around (55,163)

Central Composite Design

r r

r r

(-1,-1) (1,-1)

(-1,1) (1,1)

(0,0)

r r

r r r

ETH – p. 11/16

Data

Time Temperature 𝑥₁ 𝑥₂ Yield

50 160 -1 -1 65.3

60 160 1 -1 68.2

50 170 -1 1 66.0

60 170 1 1 69.8

48 165 -1.414 0 64.5

62 165 1.414 0 69.0

55 158 0 -1.414 64.0

55 172 0 1.414 68.5

55 165 0 0 68.9

55 165 0 0 69.7

55 165 0 0 68.5

55 165 0 0 69.4

55 165 0 0 69.0

R Output

> summary(mod1)

Call: lm(formula=yield˜x1+x2+I(x1ˆ2)+I(x2ˆ2)+x1*x2) Coefficients:

Value Std.Error t value Pr(>|t|) (Intercept) 69.0999 0.3506 197.0815 0.0000

x1 1.6331 0.2772 5.8913 0.0006 x2 1.0830 0.2772 3.9070 0.0058 I(x1ˆ2) -0.9688 0.2973 -3.2585 0.0139 I(x2ˆ2) -1.2189 0.2973 -4.0996 0.0046 x1:x2 0.2250 0.3920 0.5740 0.5839 Res. st.error: 0.784 on 7 df Mult R-Sq: 0.9143

F-statistic: 14.93 on 5 and 7 df, p-value 0.00129

ETH – p. 13/16

Response surface

Time

Temper ature Yield

Contour lines

Time

Temperature

65.5 66 66

66.5 66.5

67

67.5 68

68.5 69 69.5

70

55 60 65

160165170

ETH – p. 15/16

R code

mod1=lm(Yield˜x1+x2+I(x1ˆ2)+I(x2ˆ2)+x1*x2,data=surf) x=seq(-1,1.5,0.2); y=seq(-1,1.5,0.2)

f = function(x,y) {new.x=data.frame(x1=x,x2=y) predict(mod1,new.x)}

z = outer(x,y,f)

persp(x,y,z,theta=30,phi=30,expand=0.5,col="lightblue", xlab="Time",ylab="Temperature",zlab="Yield")

image(x,y,z,axes=F,xlab="Time",ylab="Temperature")

contour(x,y,z,levels=seq(65.5,80,by=0.5),add=T,col="peru") axis(1,at=seq(-1,1,by=1),labels=c(55,60,65))

axis(2,at=seq(-1,1,by=1),labels=c(160,165,170)) box()