Hydrocarbon Processing - January 2021 - 46

Process Optimization
point; since the center point of the design represents the current SOP, operators will have no discomfort returning to
this point multiple times in the investigation. TABLE 1 gives such a design with
sample data. The points s1 through s4 are
run in random order interspersed with
repeated interrogations of p1 to ensure
that the replicates are genuine.
In TABLE 1, the response varies for each
center-point replicate due to process
noise. If the noise is excessive, it may indicate that we have not taken account of
all factors, that our step size is too small,
or that other issues with equipment or
operation should be corrected. Statistical tests (described presently) will assess
this formally.
Our task is to construct an empirical
model to correlate the data over the factor ranges. In general, factorial designs
support models of the form (Eq. 7):

Here, α0 to αjk are coefficients to be determined, nf is the number of factors under consideration, and j and k are indices
for the summation operator. In our case,
nf = 2 and our model reduces to (Eq. 8):
ŷ = c 0 + c 1 x1 + c 2 x2 + c 12 x1 x2(8)
The term c 12 x1 x2 accounts for interactions (synergy or moderation) between
the factors, if they exist. For the small
steps taken in EVOP, interactions are unlikely. We test for them, nonetheless.
The method of least squares. According to Eq. 8 and TABLE 1, nine observations
must be accounted for (Eq. 9):
c 0 + c1x1,2 + c 2 x 2,2 + c12 x1,2 x 2,2 + ε 2

y2

!

!

(9)

c 0 + c1x1,9 + c 2 x 2,9 + c12 x1,9 x 2,9 + ε 9

y9

nf
(7)
ˆy = α 0 + ∑ k=1
α k xk +
n f −1

c 0 + c1x1,1 + c 2 x 2,1 + c12 x1,1x 2,1 + ε1

y1

where the second subscript indicates the
associated run number of each observation. The same system of equations may
be put in matrix form as shown in Eq. 10:

∑ j < k ∑ k=1 α jk x j xk
nf

TABLE 1. The run number gives the sequential order in which the points were run;
the next column shows the label of the point (see FIG. 6); the x1 and x2 columns
hold the factor coordinates for each point; the y column holds the response in the
form of a yield index (higher is better)
Run

Point

x1

x2

y

Run

x1

Point

x2

y

1

p1

0

0

165

6

s2

1

1

144

2

s3

1

-1

118

7

p1

0

0

152

3

p1

0

0

162

8

s1

-1

1

203

4

s4

-1

-1

172

9

p1

0

0

161

5

p1

0

0

161

A
Run

B
y

C

1

ŷ

D
yr

E
ȳ

2
3
4
5

1
2
3
4

165
118
162
144

160
119
160
145

160
118
160
144

6

5

161

160

7
8
9
10
11
12
13

6
7
8
9

203
152
172
161

204
160
173
160

X T X: = MMULT(TRANSPOSE(F2:I10),F2:I10)
(X T X )-1 : = MINVERSE(J2:M5)

14

ŷ = Xa: = MMULT(F2:I10,W7:W10)

F

or more simply (Eq. 11):
y = Xc + ε(11)
Here, the bolded lower-case letters refer
to the respective column vectors of Eq.
10, and X refers to the matrix of factors
that has both rows and columns (i.e., Eqs.
9, 10 and 11 are all equivalent expressions). The values for the coefficients (c)
must now be determined. If we presume
that the errors of Eq. 9 are normally distributed, the best method for fitting the
coefficients is to minimize the sum of the
squared errors (ε). For equations of the
form (Eq. 11), least squares reduces to
the following matrix equation (Eq. 12):
X T y = (X T X)c(12)
where T indicates the transpose operator
(which swaps the rows and columns of X).
Excel implements the transpose with
the function = TRANSPOSE (y) where y
is the column vector holding the responses in Excel and X is the matrix in Eq. 10.
In fact, X T X is an nc-by-nc square matrix,
where nc is the number of coefficients.
For Eq. 9, which has four coefficients, X
is a 4-by-4 matrix. To find c, we need only
pre-multiply by the inverse of X T X, denoted as (X T X)-1. This is implemented in
Excel using = MINVERSE (M) where M
is a square matrix. This provides Eq. 13:
c = (X T X)-1 (X T y)(13)

G

H

I

J

0

1

2

12

160
160
160
160

1
1
1
1

0
1
0
1

0
-1
0
1

0
-1
0
1

160

160

1

0

0

0

203
160
172
160

160
160
160
160

K

L

M

N
XTy

XTX
0
0
4
0

0 1438
0 113
0 -57
4
-5

1
-1
1
-1
0.111
0
0
1
0
0
0
0
0.25
0
1
-1
-1
1
0
0
0.25
1
0
0
0
0
0
0
CALCULATION LEGEND
X T y: = MMULT(TRANSPOSE(F2:I10),B2:B10)
c = (X T X )-1 X T y: = MMULT(S7:V10,W2:W5)

0 159.8
0 -28.3
0 14.3
0.25
-1.3

y r : = AVERAGE(B2,B4,B6,B8,B10)

9
0
0
0

0
4
0
0
(X T X )-1

y� : = AVERAGE

FIG. 7. Spreadsheet sample: the spreadsheet shows the values and results of the
Excel calculations

46 JANUARY 2021 | HydrocarbonProcessing.com

⎛ y1 ⎞ ⎛ 1 x1,1 x 2,1 x1,1x 2,1 ⎞ ⎛ c 0 ⎞ ⎛ ε 0 ⎞
⎜ y2 ⎟ ⎜ 1 x1,2 x 2,2 x1,2 x 2,2 ⎟ ⎜ c ⎟ ⎜ ε1 ⎟
⎟ ⎜ 1 ⎟ + ⎜ ⎟ (10)
⎜ ⎟ =⎜
!
! ⎟ ⎜ c2 ⎟ ⎜ ! ⎟
⎜ ! ⎟ ⎜! !
⎟⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜
⎝ y9 ⎠ ⎝ 1 x1,9 x2,9 x1,9 x2,9 ⎠ ⎝ c12 ⎠ ⎝ ε 9 ⎠

c

(B 2 .B 10 )

We may now calculate ŷ = Xα using the
Excel formula = MMULT (X,c). FIG. 7
shows the spreadsheet and all of the calculations up to this point, along with a calculation legend that details every formula.
The columns yr and ȳ will become
important when statistical significance is
TABLE 2. General analysis of variance
(ANOVA): the table provides a
convenient consolidation of statistical
properties for the data set
Term

SS

DF

MS

F

p

M

4,011

3

1,337

69.05

0.0002

R

97

5

19

T

4,108

8

s:

4.4

R : 0.9764
2
http://www.HydrocarbonProcessing.com
Hydrocarbon Processing - January 2021

Table of Contents for the Digital Edition of Hydrocarbon Processing - January 2021
