Introductory design of experiments

-6-
Introductory Design of Experiments
6.1 Introduction
Consider a chemical gas phase reaction where both pressure and temperature affect
equilibrium conversion (α) so that α = α (P, T). To study the effect of both parameters on
conversion, it is customary to fix one of the two variables and vary the other until a
maximum conversion is obtained then the second parameter is varied with fixing the first.
This method of performing experiments suffers from a serious drawback: that is, an
optimum value of one of the parameters obtained at a fixed value of the other does not
necessarily represent the optimum conditions of the experiment. To illustrate this point,
consider the following equation that was developed to relate equilibrium conversion of a
certain reaction to P (atm) and T (K):
  0.02P  .0011T  0.95 106 P 0.3 .T 2 (i)
If pressure is kept constant, say at 1 atm., then the maximum conversion is obtained by
d
differentiating the following equation and setting =0
dT
d
  0.02  .0011T  0.95 106 .T 2 , = 0.0011  2  0.95 10 6  T = 0
dT
This yields: T = 579 K
Now, substituting T = 579 in (i), we get:
d
  0.02P  .637  0.318P 0.3 and  0.01  0.3  0.318  P 0.7
dP
d
Setting = 0 , we get: P = 4.88 atm.
dP
Actually the optimum values of P and T should have been obtained by solving the system:
 
 0,  0 . We get:
T P

 0.0011  0.95  10 6  2  T .P 0.3  0 (ii)
T

 0.02  0.95 10 6  0.3.T 2 .P 0.7  0 (iii)
P
Solving (ii) and (iii), we get: P = 3.33 atm, T = 404 K
These values are different from those obtained previously.
That is why; in the present chapter is presented a method that takes into consideration the
simultaneous effect of all independent variables.

45
 Introductory design of experiments

6.2 Factorial design

6.2.1 The full factorial 2n design
A full factorial design is one in which all possible combinations of the n independent
variables (or factors) involved in the experiment are used. Each factor can be varied along k
levels. The number N of all possible combinations is:
N  kn (6.1)
For example, let in the example discussed in section 6.1, temperature be varied in steps of
50 K from 300 to 600, and pressure from 2 to 8 atm in steps of 1, then k = 7.
The number of experiments to be conducted to include all possible combinations of
temperature and pressure = 72 = 49.
When the number of factors is elevated, the total number of experiments increases
considerably and it is common in that case to use a two – level design in which each factor
is varied only twice and equation (6.1) then becomes:
N  2n (6.2)
Let the factors involved in the experiment be x1, x2, x3, …xn. To each of them is assigned a
lower limit xi min and an upper limit xi max. The mean value of xi over the interval (xi min, xi max)
is termed the center of the interval:
x  xi max
x 0i  i min (6.3)
2
And the deviation of either limit from the mean is:
x x
xi  i max i min (6.4)
2
The point with coordinates ( x 01 , x 02 , x 03 , …, x 0 n ) is called the center point of design.
These values are changed to coded dimensionless variables zi using the following
transformation:
xi  x 0i
zi  (6.5)
xi

These transformations are best explained by considering the following example.

Example 6.1
The yield of a chemical reaction is affected by three factors: Temperature (T oC), pressure
(P MPa) and residence time t (min.) Their upper and lower limits are respectively:100 –
200oC, 0.2 – 0.6 MPa and 10 – 30 min. The following table shows the % conversion
corresponding to each set of factors:

46
 Introductory design of experiments

Exp. no Temp (x1) Pressure (x2) Time (x3) % Conversion

1 100 0.2 10 2
2 200 0.2 10 6
3 100 0.6 10 4
4 200 0.6 10 8
5 100 0.2 30 10
6 200 0.2 30 18
7 100 0.6 30 8
8 200 0.6 30 12

Set the dimensionless factors table describing this experiment.

Solution:
The necessary transformations described by equations (6.3) to (6.5) have been made and the
results are presented in the following table.
100  200 0.2  0.6 10  30
x1   150 x2   0.4 x3   20
2 2 2
200  100 0.6  0.2 30  10
x1   50 x1   0.2 x3   10
2 2 2

Exp. no Temp (z1) Pressure (z2) Time (z3) % Conversion (y)

1 –1 –1 –1 2
2 +1 –1 –1 6
3 –1 +1 –1 4
4 +1 +1 –1 8
5 –1 –1 +1 10
6 +1 –1 +1 18
7 –1 +1 +1 8
8 +1 +1 +1 12

6.2.2 The design matrix

The general from of the design matrix for a 23 experiment takes the following form. A
supplementary column is added consisting of a dummy variable z0 = +1.

Exp. no z0 z1 z2 z3 y
1 +1 –1 –1 –1 y1
2 +1 +1 –1 –1 y2
3 +1 –1 +1 –1 y3
4 +1 +1 +1 –1 y4
5 +1 –1 –1 +1 y5
6 +1 +1 –1 +1 y6
7 +1 –1 +1 +1 y7
8 +1 +1 +1 +1 y8

47
 Introductory design of experiments

This matrix satisfies the following two properties:

 The scalar product of any two column vectors = 0 (6.6)
N
 The sum of any column = z
i 1
ji 0 (6.7)

6.2.3 The regression equation

Once the design matrix is set, a regression equation is assumed in the general form:
yc  a0  a1 .z1  a2 .z2  a3 .z3  a12 .z1 .z2  a23 .z2 .z3  a31.z3 .z1  a123.z1 .z2 .z3
(6.8)
Where: y c is the calculated value of y from the regression equation.
The first four terms in the above equation represent a linear model while the remaining
terms are known as interaction terms. The determination of the 8 coefficients requires
extending the design matrix to include combinations of zi  zj (i ≠ j) as well as z1.z2.z3.
For example, the value of z1.z2 for the first set of conditions =  1 1  1 and for the
second set =  1 1  1 , so that using the data of example (6.1), the matrix shows as
follows:

no z0 z1 z2 z3 z1. z2 z2. z3 z3. z1 z1 z2 z3 y

1 +1 –1 –1 –1 +1 +1 +1 –1 2
2 +1 +1 –1 –1 –1 +1 –1 +1 6
3 +1 –1 +1 –1 –1 –1 +1 +1 4
4 +1 +1 +1 –1 +1 –1 –1 –1 8
5 +1 –1 –1 +1 +1 –1 –1 +1 10
6 +1 +1 –1 +1 –1 –1 +1 –1 18
7 +1 –1 +1 +1 –1 +1 –1 –1 8
8 +1 +1 +1 +1 +1 +1 +1 +1 12

The 8 coefficients can be written in form of a column vector as follows:

a0 2
a1 6
a2 4
a3 8
A The dependent variable vector is: Y =
a12 10
a 23 18
a 31 8
a123 12

Whereas the square design matrix is:

48
 Introductory design of experiments

+1 –1 –1 –1 +1 +1 +1 –1
M=
+1 +1 –1 –1 –1 +1 –1 +1
+1 –1 +1 –1 –1 –1 +1 +1
+1 +1 +1 –1 +1 –1 –1 –1
+1 –1 –1 +1 +1 –1 –1 +1
+1 +1 –1 +1 –1 –1 +1 –1
+1 –1 +1 +1 –1 +1 –1 –1
+1 +1 +1 +1 +1 +1 +1 +1

In matrix form equation (6.8) reads:

Y = M.A (6.9)
Hence A = M-1.Y (6.10)

Example 6.2
Find the values of the regression coefficients for the data of example (6.1)

Solution:
Using the EXCEL function MINVERSE, M-1 shows as follows:

0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

-0.125 0.125 -0.125 0.125 -0.125 0.125 -0.125 0.125
-0.125 -0.125 0.125 0.125 -0.125 -0.125 0.125 0.125
-0.125 -0.125 -0.125 -0.125 0.125 0.125 0.125 0.125
0.125 -0.125 -0.125 0.125 0.125 -0.125 -0.125 0.125
0.125 0.125 -0.125 -0.125 -0.125 -0.125 0.125 0.125
0.125 -0.125 0.125 -0.125 -0.125 0.125 -0.125 0.125
-0.125 0.125 0.125 -0.125 0.125 -0.125 -0.125 0.125

2
6
4
8
Where as the column vector Y =
10
18
8
12

49
 Introductory design of experiments

8 .5
2 .5
 0 .5
3 .5
Hence the column vector A is obtained from equation (6.10): A =
 0 .5
 1 .5
0 .5
 0 .5

Hence the coded regression equation is:

yc  8.5  2.5z1  0.5z2  3.5z3  0.5z1 .z2  1.5z2 .z3  0.5z3 .z1  0.5z1 .z2 .z3

6.2.4 Testing the coefficients of the regression equation

The coefficients of the regression equation obtained are not necessary significant; that is,
some of them may be eliminated without affecting the strength of correlation. This is done
by performing a set of replicate tests at the design center of the experiment which in
example (6.1) is (150, 0.4, 20) and determining the variance of the obtained values of y. Let
the values of y obtained for three such tests be 8, 9 and 8.8. Their mean value is y = 8.6
Their standard deviation s = 0.53
The standard deviation of the 8 coefficients is related to the standard deviation of replicates
by:
s
sa  (6.11)
N
0.53
So that in the present example: sa   0.187
8
The significance of each coefficient is determined using the t – test by calculating each time
the statistic:
aijk
t ijk  (6.12)
sa
These are compared to critical values of t obtained from TINV function at a suitable value
of α and r – 1degrees of freedom. (Where, r is the number of replications at center of
design). The tested hypothesis is:
H0: aijk = 0
Rejecting H0 means that the coefficient is significant. This will occur if
tijk > tcrit

Example 6.3
Estimate the significance of the regression coefficients for example (6.2)

50
 Introductory design of experiments

Solution:
The following table shows the steps undertaken to test the 8 coefficients. At α = 0.05, d.f. =
3 – 1 =2, tcrit = 4.3

aijk 8.5 2.5 -0.5 3.5 -0.5 -1.5 0.5 -0.5

tijk 46.434 13.363 2.67 18.71 2.67 8.02 2.673 -2.67
Result reject reject accept reject accept reject accept accept
H0 H0 H0 H0 H0 H0 H0 H0

Hence four coefficients are considered significant: 8.5, 2.5, 3.5 and – 1.5. The coded form
of the regression equation becomes y  8.5  2.5z1  3.5z3  1.5z 2 z3
Replacing the dimensionless variables (z) by the original variables using equation (6.5), we
x  x 0i
get: zi  i
xi
x1  150 x  0.4 x  20
So that z1  , z2  2 , z3  3 .
50 0.2 10
The regression equation becomes: yc  12  0.05x1  15x2  0.65x3  0.75x2 .x3

4.2.5 Testing the validity of the regression equation

A preliminary evaluation of the validity of the equation obtained can be made by calculating
the determination coefficient. To this aim, equation (9.11) is used:
n

(y
i 1
c  y)2
166
R2 = n = = 0.954
 ( y  y)
i 1
i
2 174

Calculations are shown next:

yobs ycalc (yobs - y )2 (ycalc - y )2 (ycalc - yobs)2

2 1 42.25 56.25 1
6 6 6.25 6.25 0
4 4 20.25 20.25 0
8 9 0.25 0.25 1
10 11 2.25 6.25 1
18 16 90.25 56.25 4
8 8 0.25 0.25 0
12 13 12.25 20.25 1
Mean 8.5 174 166 8

A more decisive criterion is to calculate an F – ratio defined by:

51
 Introductory design of experiments

(y  y
i 1
i ci )
2

Fcalc  / s2 (6.13)
N c

Where, s is the variance of replicate readings at center of design (0.28) and c the number of
eliminated coefficients in the regression equation.
N

(y  y
i 1
i ci )
2
8
In the present case, = = 2, so that Fcalc = 2/0.532 = 7.12
N c 84
This is compared to the critical F – value obtained from FINV function at degrees of
freedom: d.f.n = N – c and d.f.d = r – 1 (Where, r is the number of replications at center of
design). In the present case: d.f.n = 8 – 4 = 4 and d.f.d = 3 – 1 = 2. At a significance level α =
0.05, the critical F – value is 19.24
Since Fcalc. < Fcrit., then, the obtained regression equation fits the experimental data
adequately.

6.3 Optimization by steepest ascent method

The regression equation obtained in the
previous section has been derived over a
limited experimentation range. It is usually
required to use this equation to seek an
extremum value to optimize the parameter y; M
one method often use is the steepest ascent
method. K
Consider the simple case where y = f(x1, x2).
The surface representing this relation is called C
the response surface. Fig.(6.1) shows a typical
response surface exhibiting a maximum point.
In the first figure, point represents the center of Fig 6.1: Response surface
design point. Starting from that point, there is a with maximum value
particular path that would lead to the maximum
point M. This path is
known as the path of steepest ascent. Any other
path like CK will not lead to this maximum.

6.3.1 The linear model

As a start an approximate linear model can be obtained to follow the path of steepest ascent.
To understand how this path can be followed, the vector Grad f is defined as:
f f ˆ f
Grad( f )  .î  . j  ... .k̂ (6.14)
z1 z 2 z N
The line of steepest ascent is the one defined by the direction of this vector.
In case of a linear regression model in the form:
yc  a0  a1 .x1  a2 .x2  a3 .x3  ...

52
 Introductory design of experiments

The partial derivatives are the coefficients of the different terms:

f f
 a1 ,  a2 ,...
z1 z 2
So that: to move along this path, one chooses an increment corresponding to each variable.
Let the increment corresponding to the variable zi be Δzi, then the corresponding increment
in any other variable should be proportional to the coefficient of this variable in the
regression equation. That is:
zi ai
 (6.15)
z j a j
We start at the design center (0, 0, 0, …) and choose the increment of the variable (i) which
has the highest coefficient Δzi, then we calculate the increments of the other variables using
equation (6.15).
For example, let the first dimensionless variable z1 that which has the highest coefficient
(a1). Let the chosen increment = Δz1
a2 a3
Then z 2  .z1 , z 3  .z1 , etc…
a1 a1
Hence, z11 = 0 + Δz1, z21 = 0 + Δz2, z31 = 0 + Δz3, … and
z12 = z11 + Δz1, z22 = z12 + Δz2, z32 = z13 + Δz3, …
Generally: z1, i+1 = z1i + Δz1, z2, i+1 = z1i + Δz2, z3,i+1 = z1i + Δz3, …
A series of experiments is then undergone at the obtained values of z1i, zi2, zi3, etc…for
value of i = 1, 2, 3, etc. until the value of y stabilizes which means that a maximum value is
obtained.

Example 6.4
Solid state sintering of magnesia is governed by soaking time ( h), temperature ( oC) and
compacting pressure ( MPa). A 23 factorial experiment gave the following results for bulk
density of pressed compacts (g/cm3).
Four replicate experiments performed at the center of design gave the following bulk
densities: 2.425, 2.415, 2.42, 2.4056.
Perform a full factorial design at significance level = 0.05 and use the steepest ascent
method to show how you would reach the optimum conditions

No Time Temp Pressure Density

1 1 1100 20 2.37
2 5 1100 20 2.44
3 1 1500 20 2.44
4 5 1500 20 2.39
5 1 1100 40 2.41
6 5 1100 40 2.45
7 1 1500 40 2.43
8 5 1500 40 2.41

53
 Introductory design of experiments

Solution:
Center of experiment: x1 = t = 3 h, x2 = T = 1300oC, x3 = P = 30 MPa
Deviations: Δx1 = 2, Δx2 = 200, Δx3 = 10
The dimensionless matrix is:
No Dummy Time Temp Pressure Density
1 1 -1 -1 -1 2.31
2 1 1 -1 -1 2.44
3 1 -1 1 -1 2.38
4 1 1 1 -1 2.43
5 1 -1 -1 1 2.38
6 1 1 -1 1 2.45
7 1 -1 1 1 2.42
8 1 1 1 1 2.41

The design dimensionless matrix M is:

x0 x1 x2 x3 x1.x2 x2.x3 x3.x1 x1.x2.x3

1 -1 -1 -1 1 1 1 -1
1 1 -1 -1 -1 1 -1 1
1 -1 1 -1 -1 -1 1 1
1 1 1 -1 1 -1 -1 -1
1 -1 -1 1 1 -1 -1 1
1 1 -1 1 -1 -1 1 -1
1 -1 1 1 -1 1 -1 -1
1 1 1 1 1 1 1 1

The inverse matrix is then obtained and multiplied by the column vector corresponding to
density values. This yields the coefficient column vector

A=

[ ]

The dimensionless regression equation is then:

54
 Introductory design of experiments

The next step is to eliminate insignificant coefficients by considering the four replicate
values at center of design: 2.425, 2.415, 2.41, 2.406.
Their standard deviation = 0.0082. The standard deviation of coefficients is therefore: sa =
0.0082
= 0.0029
8
The following table shows the steps undertaken to test the coefficients for the null
hypothesis:
H0: aijk = 0
At α = 0.05, d.f. = 4 – 1 = 3, tcrit = 3.18

aijk 2.4025 0.03 0.0075 0.0125 -0.02 -0.008 -0.015 0

tijk 828 10.34 2.58 4.31 6.89 2.58 5.17 0
reject reject accept reject reject accept reject accept
Result
H0 H0 H0 H0 H0 H0 H0 H0

The final equation is;

y  2.4025  0.03z1  0.0125z3  0.02 z1 .z2  0.015z2 .z3 (6.16)
This is tested for validity using the R2 criterion. Following equation (9.11), we get: R2 =
0.999
Also, from equation (6.13), for N = 8 and c = 3, we get:
N

( y
i 1
i  yci )2
0.0009
= = 0.00018
N c 83
0.00018
The calculated value of F = = 2.67
0.00822
This is compared to the critical F – value obtained from FINV function at degrees of
freedom: d.f.n = n – c and d.f.d = r – 1. In the present case: d.f.n = 8 – 3 = 5 and d.f.d = 4 – 1
= 3. At a significance level α = 0.05, the critical F – value is 9.013 obtained regression
equation is adequate.
The actual equation can be obtained from:
x1  3 x  1300 x  30
z1  , z2  2 , z3  3
2 200 10
This way we get:
y  1.8325  0.08x1  0.000375x2  0.011x3  0.00005x1 .x2  0.0000075x2 .x3 (6.17)
To follow the steepest ascent direction, we arbitrarily choose an increment of time Δx1 = 0.5
h. This corresponds to a dimensionless increment Δz = 0.25
The increment in temperature, from the obtained dimensionless model is zero. So, we can
move along with a fixed temperature of 1300oC.

55
 Introductory design of experiments

0.0125
As for pressure, the increment is calculated from equation (6.15) to be:  0.25 =
0.03
0.104 (corresponding to a pressure increment = 0.104  10 = 1.04 MPa)
We have to note, however, that a linear approximation has been made as the regression
equation contains a non – linear term z2.z3
The suggested experimental conditions to be followed towards maximum density will then
yield the following values for density by substitution in equation (6.17):
Step no Time Temp. Pressure ρ [Link]-3
1 3 1300 30 2.373
2 3.5 1300 31.04 2.380
3 4 1300 32.08 2.388
4 4.5 1300 33.12 2.396
5 5 1300 34.16 2.404

The last entry for time is 5 hours, as this is its maximum value. Therefore, according to that
model, a maximum density would be obtained at the following conditions: t = 5 h, T =
1300oC and P = 34.16 atm. One gets a maximum value of density = 2.404 g/cm3.

6.4 Exercise problems

(1) The strength of polystyrene boards (kPa) depends on their porosity and temperature of
exposure. In a 22 experiment, the following design center is chosen: porosity = 0.6 with
step = 0.15 and temperature = 30oC and step = 10oC. The following results were
obtained:

Step Porosity Temp. Strength

1 0.45 20 450
2 0.75 20 180
3 0.45 40 300
4 0.75 40 125

Five replicate experiments performed at the center of design gave the following values
of strength: 256, 284, 270, 278, 250 kPa
Perform a full factorial design at significance level = 0.025 and use the steepest ascent
method to show how you would reach the optimum conditions.

(2) To study the effect of temperature, time and particle size on the extraction of a certain
liquid from a solid a full factorial two – level design was performed. The results are
listed in the following table:
The last three rows represent a replicate at the center of the design.
Derive a first order regression equation showing the dependence of yield on the three
factors then check the significance of the different coefficients.
Also check the validity of the system. (Take α = 0.05)

56
 Introductory design of experiments

Time Temp Particle size Yield

20 )
(min ( o40
C ) 0.1 )
( mm 16/ g )
( mg
20 40 0.3 14
20 80 0.1 21
20 80 0.3 17
40 40 0.1 23
40 40 0.3 21
40 80 0.1 26
40 80 0.3 24
30 60 0.2 21
30 60 0.2 22
30 60 0.2 23

(3) The following data have been obtained on studying the effect of three variables:
temperature of calcination (oC), time of leaching (h) and acid concentration (%) on the
concentration of aluminum sulfate in solution upon treatment of kaolin with sulfuric
acid. Four specimens investigated at the center of design (750oC, 3 h, 40%) yielded
concentrations of 6.82%, 7.27%, 6.77% and 7.8%. Derive a regression equation
showing the dependence of yield on the three factors then check the significance of the
different coefficients.
Also check the validity of the system. (Take α = 0.05)

# Temp. oC Time h Conc. % Yield %

1 600 2 20 8.27
2 600 2 60 7.08
3 600 4 20 7.36
4 600 4 60 9.31
5 900 2 20 4.95
6 900 2 60 4.45
7 900 4 20 6.22
8 900 4 60 2.57

57