Global sensitivity analysis and uncertainty quantification of crude

Le Quang Minh1, Pham Luu Trung Duong1, Moonyong Lee*. School of Chemical Engineering, Yeungnam University, Gyeongsan, Rep. of Korea. Page 1 of 40...
0 downloads 0 Views 2MB Size
Subscriber access provided by Queen Mary, University of London

Article

Global sensitivity analysis and uncertainty quantification of crude distillation unit using surrogate model based on Gaussian process regression Le Quang Minh, Pham Luu Trung Duong, and Moonyong Lee Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b05173 • Publication Date (Web): 23 Mar 2018 Downloaded from http://pubs.acs.org on March 26, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

1

Global sensitivity analysis and uncertainty quantification of crude

2

distillation unit using surrogate model based on Gaussian process

3

regression Le Quang Minh1, Pham Luu Trung Duong1, Moonyong Lee* School of Chemical Engineering, Yeungnam University, Gyeongsan, Rep. of Korea

4

Submitted to Industrial & Engineering Chemistry Research

5 6 7 8

1

9

*

These authors contributed equally to this work. Correspondence concerning this article should be addressed to:

10

Prof. Moonyong Lee

11

Email: [email protected]

12

Telephone: +82-53-810-2512

13 14 15 16 17

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

18 19

Abstract

20

This study presents a global sensitivity analysis to simplify a surrogate-model-based

21

uncertainty quantification of a crude distillation unit with a large number of uncertainties. To

22

overcome the computational limitation of a conventional surrogate model based approach where

23

the number of simulations required grows exponentially as the input dimension increases, a

24

novel two-stage approach was proposed in this study: in the first stage, a multiplicative

25

dimensional reduction method is applied to identify factors that exert the highest influence on the

26

model outputs. In the second stage, the Gaussian process regression is exploited for uncertainty

27

quantification from the simplified model derived in the first stage. As a result, the computational

28

efforts for uncertainty quantification were significantly reduced (approximately more than 95 %)

29

compared to the conventional Quasi Monte Carlo, while the predicted density functions by the

30

proposed method closely matched with those from the Quasi Monte Carlo. The proposed two-

31

stage approach was executed for sensitivity analysis and uncertainty quantification of a crude

32

distillation unit by an interface between MATLAB and HYSYS. The economic revenue and the

33

operating cost per unit of crude oil processed were selected as the output of interests for the

34

crude distillation unit. The global sensitivity analysis result showed that the flow rates of crude

35

oil and naphtha products are critical for both of the economic revenue and the operating cost.

36

Keywords: Gaussian process regression; Global sensitivity analysis; Uncertainty

37

quantification; Crude distillation unit; Operating cost reduction; Maximum economic revenue.

38

2

ACS Paragon Plus Environment

Page 2 of 40

Page 3 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

39

Abbreviations

40

AGO: atmospheric gas oil

41

ASTMD86: American Society for Testing and Materials (standard test method for

42

distillation of petroleum products and liquid fuels at atmospheric pressure)

43

CDU: crude distillation unit

44

ER: economic revenue per unit of crude oil processed

45

GP: Gaussian process

46

GPR: Gaussian process regression

47

HP: high pressure

48

HSR: heavy straight run

49

LSR: light straight run

50

MC: Monte Carlo

51

MDRM: multiplicative dimensional reduction method

52

OP: operating cost per unit of crude oil processed

53

PA: pump around

54

PAQ: outlet temperature of the pump around from the cooler

55

QMC: Quasi Monte Carlo

56

SA: sensitivity analysis

57

UQ: uncertainty quantification 3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

58

1. Introduction

59

The use of fossil resources for the production of fuels and commodity chemicals has been

60

a focal point of discussion for the past few decades. It is anticipated that the total amount of

61

fossil fuels (coal, natural gas, and oil) available on the earth can support various needs for at least

62

another hundred years1. Identifying a profitable/efficient production process is a crucial task for

63

engineers which facilitates higher cost savings and is safer and cleaner. Most industries would

64

rather maximize available resources for maximum profitability than opt for large-scale

65

expansion2.

66

optimization, which is a major quantitative tool in decision-making for process availability3.

67

Mathematical models, which elucidate the relationship among input variables, play a vital role in

68

the final design to comply with the safety and serviceability constraints4.

A common approach that has gained widespread application in this regard is

69

A crude distillation unit (CDU) is among the largest energy-consuming processes owing

70

to the large amount of energy consumption and high operating temperature5. It is reported that

71

the operation of the CDU and vacuum distillation consumes over 35% of the energy required in

72

the refinery plants6. Therefore, research on crude distillation systems has received considerable

73

attention. In optimizing the CDU, the objective is to maximize the profits of a CDU while

74

minimizing the energy consumption (resulting in reduced CO2 emissions), subject to constraints

75

on the quality of the products7. Previous studies recommended that the separation and energy

76

performance of distillation processes can be represented by different type of models such as

77

rigorous, simplified, and statistical/empirical models8-10. Nevertheless, the development of such

78

models generally requires significant amount of effort because this involves various differential

79

and algebraic equations associated with a large number of input variables. Simultaneously,

80

uncertainties are sometimes correlated in the process design because of the connected process 4

ACS Paragon Plus Environment

Page 4 of 40

Page 5 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

81

and various random factors. The complex nature of CDUs, including their interactions with the

82

associated heat recovery network, the large number of degrees of freedom, and a large number of

83

uncertainties in feed mixture and process conditions, renders their optimization a highly

84

challenging task.

85

The goal of uncertainty quantification11, and in particular, a priority in the design and

86

modeling of an efficient process, is to assess the effect of the input uncertainties on the model

87

output and consequently, on the system performance12. Certain parameters such as the crude oil

88

property and tray efficiency, which depend on physical properties (e.g., tray hydraulics, liquid

89

mixing and entrainment), are likely to illustrate a natural disturbance, which is not revealed in

90

advance (this is referred to as an aleatory uncertainty). However, certain other factors such as

91

pressure, temperature, and flow rate, which are likely to exhibit a unique value, are directly

92

measurable and susceptible to lack of knowledge (epistemic uncertainty). The former are

93

typically observed to be irreducible, while the latter can be quantified and minimized13. For

94

considering the unknown factors, a global sensitivity analysis (SA), which is a widely used

95

technique for selecting representative parameters from mathematical models14, 15, is required to

96

reduce their quantity and select the ones to be incorporated.

97

Monte Carlo (MC)16 and Quasi Monte Carlo (QMC) methods are representative

98

probabilistic approaches for the propagation of uncertainties in the model inputs to its

99

output17. However, notwithstanding the simplicity of their implementation, the mean

100

convergence rate was estimated to be in the order of O(1/ M ) , where M is the number of

101

samples; this renders MC/QMC-based approaches computationally expensive, and they are used

102

only as the final option. To overcome the problem, a convenient-to-evaluate function may be

103

used to replace the original computational model18. Shen and Braatz19 developed a new 5

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

104

polynomial chaos based algorithm on the design of batch and continuous-flow chemical reactors

105

with probability uncertainties. Duong et al.20 studied the UQ and SA of chemical processes using

106

the standard polynomial chaos method for systems with a small number of random inputs. Using

107

an efficient-to-evaluate surrogate model, Celse et al.21 investigated the influence of each input on

108

the construction of the model of the hydrodesulfurization/hydrodenitrogenation reaction in

109

vacuum gas oil hydrotreatment. Duong et al.22 and Minh et al.23 recently developed a two-stage

110

approach using a compressive polynomial chaos method to overcome the computational

111

limitation in a system with a moderate/large number of uncertain parameters.

112

Surrogate models, also known as meta-models, can approximate detailed mechanistic

113

models, which are capable of simplifying highly nonlinear and computationally expensive

114

problems. The Kriging technique, which originated from interpolating geographical data in

115

mining, is also known as Gaussian process (GP) modeling24, 25. In general, the GP modeling is a

116

Bayesian approach, which provides a complete posterior distribution over feasible functions

117

from the observed data. It is completely determined by its mean and covariance function. The

118

mean is used to encode the prior knowledge of the output. The GP prediction is a distribution.

119

Use of the predictive GP can facilitate the selection of the new data. A high predictive variance

120

implies an absence of the training data, whereas, a low predictive variance implies a high

121

confidence in the model. Note that the higher confidence in the model results in a more reliable

122

design. The covariance function is used to capture a straightforward intuition: if the inputs are

123

close, the function outputs are likely to be closed. It can also be used to encode the periodic

124

property of the output. They are widely used in diverse fields of engineering, such as structural

125

reliability analysis and design optimization26-28.

6

ACS Paragon Plus Environment

Page 6 of 40

Page 7 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

126

For the Gaussian process regression (GPR), however, the number of hyper-parameters

127

increases with the dimension of the input space so that their statistical inference can be

128

problematic in high-dimensional situations29. The number of simulations required to learn a

129

generic multivariate response grows exponentially as the input dimension increases. This is

130

because the GPR defines an input space correlation with Euclidean (2-norm) distance. Since the

131

Euclidean distance becomes noninformative as the dimensionality of uncertain inputs increases30,

132

more evaluations of the expensive-to-evaluate models are required to construct an accurate GPR.

133

Therefore, it is desirable to reduce the dimension of the problem from the perspectives of both

134

accuracy and efficiency. In this study, a two-stage approach was presented to tackle this issue:

135

firstly, a multiplicative dimensional reduction method is used for detecting non-influential inputs,

136

then a GPR model is constructed with the influential inputs only. Although the GPR possesses

137

the automatic relevant determination mechanism wherein the kernel scales can be used to

138

determine the relevancy of inputs, this approach may be not reliable23.

139

This work aimed to overcome a current limitation of the conventional GPR where the

140

number of model evaluations grows exponentially and may not be applicable to systems with a

141

moderate/large number of uncertainties such as a CDU. One of main contributions of this work is

142

to propose a novel two-stage approach by combining the multiplicative dimensional reduction

143

with Gaussian process kriging regression methods. The proposed two-stage approach could

144

remarkably reduce the simulation time (one of the most burdens in handling UQ) so that the UQ

145

for complex systems with a large/moderate number of uncertainties can be popular in practice.

146

Another main contribution of this work is to carry out a global SA and UQ of the CDU which is

147

highly nonlinear and complex. In terms of our understanding, there is no report of the global SA

148

that includes a discussion on the uncertain inputs effects on the process performance for the

7

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

149

conventional CDU, particularly with a large number of random variables. Thus, this work is

150

expected to describe the uncertain effects of the process inputs on main process interests of the

151

CDU such as the economic revenue and operating costs. Also, this approach attracts considerable

152

attention in refinery modeling and optimization because it directly determines the quality of final

153

products and its profitability.

154

This paper is organized as follows: Section 2 addresses the global SA of the

155

multiplicative dimensional reduction method, and subsequently, the UQ through the GPR of the

156

identified influential factors from the previous step. Section 3 briefly describes the CDU and

157

analysis methods of the system including the economic revenue and operating costs. Section 4

158

discusses the results, and Section 5 concludes the study.

159 160

2. Methodologies

161

UQ requires a large number of original model evaluations, particularly for high-nonlinear

162

and high-dimensional problems. The GPR is a Bayesian modeling technique for UQ purposes

163

that has recently attracted engineers from different fields such as biological, chemical, and

164

control engineering24, 31, 32 owing to its computational efficiency. Although the GPR is flexible,

165

the optimizing approach by selecting hyper-parameters either got stuck in local solutions24, 31 or

166

achieved non-satisfactory results when it is applied to model systems with a large number of

167

random inputs. To address these problems, an efficient two-stage UQ method, which combines

168

the multiplicative dimensional reduction method (MDRM) with the GPR, is proposed in this

169

study. In the proposed method, the MDRM is first adopted by determining critical inputs for a

170

model-dimension reduction. The meta-model by the GPR is then consequently generated

8

ACS Paragon Plus Environment

Page 8 of 40

Page 9 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

171

corresponding to the critical inputs determined in the previous step. The work intends to explain

172

the efficient and practical framework, and the typical steps involved in SA and UQ in order to

173

reduce the computational efforts.

174

2.1. MDRM for SA

175

SA is a procedure aimed at determining the inputs that are relatively more critical than

176

the others. It is of two types: local SA and global SA. The first utilizes the information from the

177

partial derivative of the model output with respect to the model inputs. Thus, it is used only to

178

study the vicinity of a nominal value of the selected inputs. That is, the global SA involving

179

Sobol indices14 is addressed with the global output disturbance for the entire ranges of input

180

variation. Therefore, it is preferable over the local SA for an incomplete determination of the

181

sensitivities of a specific system. Several approaches for global SA with Sobol indices are

182

available, such as Fourier amplitude sensitivity test15, extended Fourier amplitude sensitivity

183

test33, MC, and MDRM. A multivariate function in MDRM34 is approximated as a product of

184

univariate ones, which demonstrates a reasonable tradeoff between the accuracy and

185

computational load. The main theory of global SA by the MDRM is briefly described in this

186

subsection.

187

A system defined by a computational model (1) is considered:

188

y ( ξ ) = Μ (ξ )

189

where M denotes the computational model, y is the output of interest, and

190

represents the vector of uncertain inputs.

(1)

9

ACS Paragon Plus Environment

ξ = ( ξ1,...,ξn )

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 40

191

The uncertain inputs are assumed to be randomly independent variables with density

192

functions pi (ξi ): Γi → R . If the output y exhibits a finite variance, it can be decomposed into

193

the sum of increasing dimension terms:

+

n

n

n

y (ξ ) = y0 + ∑ yi (ξi ) + ∑∑ yij (ξi , ξ j ) + ... + yi1 ...in (ξi1 ,..., ξ n ) i =1

194



= y0 +

i =1 i < j

(2)

yA (ξ A )

A ⊂{1,..., n} A ≠∅

A ={i1,..., is} ⊂{1,.., n} denotes a subset of indices, and ξA is a subset of ξ that contains

195

where

196

only parameters indexed by A . The terms in Eq. (2) are computed recursively as

y0 = E [ y(ξ)] 197

yi (ξ) = E [ y(ξ ) | zi ] − y0 (t )

yij (ξi , ξ j ) = E  y (t , ξ ) | ξi , ξ j  − yi (ξi ) − y j (ξ j ) − y0 (t )

(3)

...

E(ξ | ξA ) is the conditional expectation. Because the

198

where E is the expectation operator, and

199

functions defined in Eq. (3) are orthogonal to each other, the variance of the model output can be

200

decomposed into

201

Vy =



Var ( yA (ξ A ))

(4)

A ⊂{1,..., n} A ≠∅

202

The first order Sobol indices (Si), which measure the contribution of a single input factor

203 204

ξi

205

by Eq. (5):

on the variance of the output without considering the joint effect with other inputs, is defined

10

ACS Paragon Plus Environment

Page 11 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Vi Vy

206

Si =

207

where Vi denotes the variance of

(5)

yi (ξi ) .

The Sobol indices for a subset of parameters A are defined by Eq. (6), which quantifies

208 209

the joint effect of the subset of parameters A .

210

SA =

VA Vy

(6)

Finally, the total Sobol indices (Ti), including all the interactions with the other inputs,

211 212

quantify the total impact of a given input ξi :

213

Ti = ∑ SA

(7)

i∈A

214

215

The main concept behind the MDRM is to approximate the model output in Eq. (1) as n

n

i =1

i =1

y(ξ) ≈ y(c)1−n ∏ y(c1 ,..., ci−1 , ξi , ci+1 ,..., cn ) = h01−n ∏ y(ξi , c−i )

(8)

216

where c = (c1,..., cn ) is an anchor point, and c−i = (c1,..., ci−1, ci+1,..., cn ) . Generally, the mean

217

vector of the uncertain input is used as the anchor point.

218

219

Let us denote the ith mean and the mean square of the dimensional function as µ i = E [ y (ξ i , c− i )]

(9)

σ i = E [ y (ξ i , c− i ) 2 ]

11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Utilizing a Gaussian quadrature with q nodes {wi , ξi }l =1 , the ith mean and mean (l )

220 221

Page 12 of 40

square of the dimensional function ( µi and

(l ) q

σi , respectively) can be computed as

q

µi = ∑ wi (l ) y(c1 ,..., ci −1 , ξi (l ) , ci +1 ,..., cn ) 222

l =1

(10)

q

σ i = ∑ wi (l ) ( y (c1 ,..., ci −1 , ξi (l ) , ci +1 ,..., cn ))2 l =1

223 224 225

226

Using Eq. (8), it can be demonstrated that the first order and total sensitivity functions can be approximated34 as

Si ≈

−1 + σ i / µi 2 n

−1 + (∏σ k / µk 2 )

;

Ti ≈

k =1

1 − µi 2 / σ i

(11)

n

1 − (∏ µk 2 / σ k ) k =1

227

where the ith mean and mean square of the dimensional function are computed as in Eq. (10).

228

Remarks: If Ti is approximately equal to zero, it implies that the input

229

can be fixed at any point in its distribution without affecting the variance of the output (4). The

230

difference is an indicator of the presence of interactions in the model.

ξi

is non-influential and

231 232

2.2. GPR method

233

We assume that m (< n) important inputs ( ξ m = (ξ i , ..., ξ i 1

234

m

) ) were detected in the previous

{

}

step. Let us make p simulations and collect the following data set ξm , y = M (ξm ) (i )

(i )

12

ACS Paragon Plus Environment

(i )

p

i=1

.

Page 13 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

235 236

237

Industrial & Engineering Chemistry Research

A GPR model is a generalization of the multivariate Gaussian random variables to an infinite dimension. It is defined as

y(ξm ) = h(ξm )T β + f (ξm )

(12)

238

where f (ξm ) is a zero mean Gaussian random process with the covariance function k ( ⋅, ⋅) ,

239

and h(ξm) is a set of basic functions or a trend which transforms the original input space into a

240

new feature vector h(ξm) in Rm, a vector of basis coefficients. That is, the GPR is a statistical

241

model that exhibits an instant response of the form

242

P( y(i ) | f (ξm(i ) ), ξm(i ) ) ≈ N ( y(i ) | h(ξm(i ) )T β + f (ξm(i ) ),σ 2 )

243

where N () is a normal distribution in which the polynomial functions of the inputs can be

244

normally used as the basis functions. It is noted that the latent function, f , is introduced for

245

each observation. Various covariance functions are described in the literature25, 32. In this study,

246

the automatic relevance determination squared exponential is used, which is defined as follows:

247

k (ξ m(i ) , ξ m( j ) | θ) = θ f exp(−

(13)

1 m (i ) ∑ (ξl − ξl ( j ) )2 / θl ) 2 l =1

(14)

248

m +1 where θ ∈ R = (θ f , θ1 ,..., θm ) is a vector of hyper-parameters. It is expected that similar values of

249

inputs should exhibit similar responses. Thus, the covariance function can be used to characterize

250

the similarity between two input vectors ξm , ξm . The length scale

251

extent to which one is required to move along a particular input space for its response to become

252

uncorrelated. Such a covariance function implements the automatic relevance determination

(i )

( j)

θ1,...,θm

13

ACS Paragon Plus Environment

quantifies the

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 40

253

mechanism24, which determines the relevance of an input is owing to the inverse of the length-

254

scale. Therefore, the GPR model is defined by the coefficient vectors β , hyper parameters ɵ,

255

σ . The values for these parameters will be updated and learnt

256

and measurement noise variance,

257

from the given data by maximizing the logarithm likelihood expressed in Eq. (15): log P ( y | ξ m , β , θ , σ )

258

=

1 2

259

T

'

2

p 2

log(2π ) −

1 2

'

log K ( ξ m , ξ m | θ ) + σ I p 2

(15)

where Ip is an identity matrix of size p, and

 (ξ m (1) )T  (2 T  (ξ ) X = m ...   (ξ ( P ) )T  m 260

−1

( y − H β )  K ( ξ m , ξ m | θ ) + σ I p  ( y − H β ) −

  h(ξ m(1) )T   (2) T   h(ξ m ) = , H   ...     h(ξ ( P ) )T   m

      

 y (1)    y =  ...   y( P)   

(16)

 k (ξ m (1) , ξ m (1) ) ... k (ξ m(1) , ξ m( P ) )    K (ξ m , ξ m ) = ... ... ...   k (ξ ( P ) , ξ (1) ) ... k (x( P ) , ξ ( P ) )  m m  m  261

where H p×r is a matrix of values of the basic vector at training points, r is the number of basic

262

function, and K is the covariance matrix at the training points.

14

ACS Paragon Plus Environment

Page 15 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

263

Regarding the estimated optimal parameters {βˆ , θˆ, σˆ } , which are obtained by the

264

maximization of the logarithm likelihood function, the density function of the response at a new

265

data point ξ m is

266

N (h(ξm* )βˆ + µ,σˆ 2 + Σ)

267

where

*

(17)

µ = K (ξ m*T , ξm | θˆ)( K (ξm , ξm | θˆ) + σˆ I p )-1 ( y - H β ) (18)

268

K (ξ

*T m

, ξm | θˆ) = ( k (ξm* , ξm(1) ),..., k (ξm* , ξm (n ) ) )

269

and

270

∑ = k (ξ m * , ξ m * ) - K (ξ m * T , ξ m | θˆ )( K (ξ m , ξ m | θˆ ) + σˆ I p ) -1 K (ξ m , ξ m * T | θˆ )

(19)

271

The mean value calculated in Eq. (19) is used as the prediction for the new data, whereas

272

the variance ∑ is an indicator of the uncertainty of the prediction. One can use Eq. (18) as a

273

surrogate to sample large numbers of output responses and to estimate density functions and

274

statistics of the output interest.

275 276

3. Crude oil distillation unit

277

Table 1 presents the compositions of the crude oil processed in this study35.

278

simulation is carried out in Aspen HYSYS v10. The property package is Peng–Robinson, which

279

utilizes the properties of equation of state in its enthalpy calculations. Owing to the assay

280

characterization, the feed contains 36 components, including 5 light components and 31 pseudo15

ACS Paragon Plus Environment

The

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

281

components. Figure 1 shows the oil cut production distribution. Table S1 lists the product

282

specifications (see Supporting Information). The temperature of each cut is referred at 95%

283

liquid fraction recovery (ASTMD86-95%). The boiling point curves of column products are

284

shown in Figure S1 (see Supporting Information); these products are the naphtha product, three

285

side products (kerosene, diesel, and atmospheric gas oil), and the residues at the bottom of the

286

column. In this study, it is assumed that the light and heavy straight run naphtha streams are

287

mixed to be the naphtha product stream. It is also assumed that the input crude oil is at a

288

temperature of 232.2 °C and pressure of 5.17 bar, operating at 100,000 barrels/day.

289

Table 1. Feed compositions Type

Light component

Hypothetical components

Components Mole fraction (%) Methane Ethane Propane i-Butane n-Butane Water NBP49 NBP79 NBP111 NBP144 NBP176 NBP208 NBP240 NBP272 NBP304 NBP336 NBP368 NBP400 NBP433 NBP464 NBP496 NBP528 NBP560

0.031 0.069 0.952 0.601 2.131 0.000 0.000 0.000 5.773 2.995 3.040 5.381 5.277 4.821 4.489 4.105 3.811 3.718 3.837 4.275 4.301 4.053 3.617 16

ACS Paragon Plus Environment

Page 16 of 40

Page 17 of 40

NBP592 NBP624 NBP656 NBP688 NBP720 NBP752 NBP784 NBP830 NBP888 NBP947 NBP1009 NBP1062 NBP1124

3.205 2.824 2.438 2.138 1.893 1.691 1.490 1.876 2.173 2.043 2.250 3.355 5.346

290

0.5 Light straight run naphtha Heavy straight run naphtha Kerosene Diesel Atmospheric gas oil Residues

0.4

Liquid volume fraction

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

0.3

0.2

0.1

0.0 0

200

400

600 o

291 292

ASTM D86 - 95% ( C) Figure 1. Cut distributions of crude oil.

17

ACS Paragon Plus Environment

800

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

293

Figure 2 demonstrates a rigorous CDU model as a substitute for the actual plant35. The

294

tower consists of 29 stages in the main shell and three side columns associated with three pump-

295

arounds for stripping the middle products. Furthermore, a pre-flash separator is used to enhance

296

the separation efficiency in separating the light components from the liquid phase36. After

297

leaving the pre-flash separator, the crude oil enters the CDU furnace (see Figure 2). The crude is

298

vaporized through the furnace in order to recover the products from the side strippers. The outlet

299

temperature of the furnace is normally set so that the vaporized fraction equals the sum of the

300

column products plus a small percentage. This small percentage excess is the “over-flash,” which

301

indicates the amount of the heavy residues recovered in the light fractions37. The two-phase

302

mixture enters the column at 333.4 °C and is flashed in the bottom section (flash zone).

303

Furthermore, a significant amount of the main stripping steam is included at the bottom, which

304

serves to strip any residue and prevent the excessive thermal cracking of crude oil because of the

305

high temperature at the flash zone. Other light fractions of the crude were retrieved by side

306

strippers at different locations using steam or reboilers. Both the atmospheric gas oil and diesel

307

side-strippers use stripping steams, while the kerosene side-stripper includes a reboiler. The

308

presence of side coolers or pump-arounds is an additional feature in the CDU. These units permit

309

the reduction of the condenser cooling requirements by reducing the vapor flow in the column

310

and allowing for heat recovery. The residue from the column contains valuable hydrocarbons,

311

which is typically further separated in a vacuum distillation column (not discussed in this study).

18

ACS Paragon Plus Environment

Page 18 of 40

Page 19 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

312 313

Figure 2. Process flow diagram of crude oil distillation unit.

314 315

3.1. Process decision variables

316

Rigorous models, which simulate the distillation column according to the material and

317

energy balances as well as equilibrium relations in each stage of the column38, are sufficiently

19

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

318

realistic to provide accurate predictions39. In this work, Aspen HYSYS v10 is used to simulate

319

the CDU. To facilitate the data transfer required for the global SA and the UQ, an interface

320

between MATLAB and HYSYS is established using the automation client–server application

321

provided by HYSYS; meanwhile, the algorithm is coded in MATLAB.

322

For UQ and SA, the inputs of the surrogate model are varied, including the flow rate of

323

the crude and that of the main stripping steam, product flow rates, outlet temperature of crude oil

324

from the furnace, pump-around rates and their temperature changes, stripping steam rates, and

325

reboiler duty of the side strippers. Sensitivity analysis indicates the parameters that are most

326

influential and are to be focused upon and those which can be discarded. Furthermore, the UQ

327

study can examine their relative influences on the outputs. All the inputs are assumed to be

328

distributed uniformly or normally and to independently affect the process performance

329

(economic revenue and operating costs). Table 2 lists all the uniformly distributed input variables

330

associated with the confidence boundaries. Note that the normally distributed input variances, i.e.,

331

the Gaussian uncertainties, were also investigated in the case study without fixing product flow

332

rates. In that case study, the same means and variances as in the uniform case were used.

333 334

Table 2. Boundaries of random inputs under uniform uncertainties Variables

Base

Lower

Upper

Unit

Crude oil

4166. 7

3958.3

4375.0

bbl/h

Naphtha

958.3

910.4

1006.2

bbl/h

Kerosene

387.5

368.1

406.9

bbl/h

Diesel

802.1

762.0

842.2

bbl/h

AGO

187.5

178.1

196.9

bbl/h

20

ACS Paragon Plus Environment

Page 20 of 40

Page 21 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

PA1

2083.3

1979.2

2187.5

kmol/h

PA1Q

95.0

90.2

99.8

°C

PA2

1250.0

1187.5

1312.5

kmol/h

PA2Q

80.0

76.0

84.0

°C

PA3

1250.0

1187.5

1312.5

kmol/h

PA3Q

72.0

68.4

75.6

°C

Diesel steam

75.5

71.8

79.3

kmol/h

AGO steam

62.9

59.8

66.1

kmol/h

Main steam

188.8

179.4

198.3

kmol/h

Kerosene HP steam

157.1

149.3

165.0

kmol/h

Hot crude temperature

343.3

326.1

360.5

°C

335 336

3.2. Objective functions in the CDU

337

The economic revenue (ER) and operating cost (OP) per unit crude oil processed are

338

selected as the performance indicator for UQ and SA. Fundamentally, the ER is estimated on the

339

basis of the operating costs in the column, as follows40: n

OP = 340

Fsteam C steam + ∑ Qx C x x =1

Fcrude n

n

∑ Fj C j − [ FcrudeCcrude + FsteamCsteam + ∑ QxCx ] ER =

j =1

(20)

x =1

Fcrude

341

where Fj and Cj are the flow rates and costs, respectively of the jth product; and Fsteam and Csteam

342

are the flow rates and costs, respectively, of the low pressure steam. Qx and Cx are the heat

343

requirement and its utility cost, respectively. The estimation is based on the assumption of 8600

344

hour per year. 21

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

345

The ER is based on the product prices, raw material cost, and utility costs to calculate the

346

operating costs within a refinery. Furthermore, the fixed costs such as capital modification, labor,

347

and maintenance were not considered in this study. Table 3 lists the prices of the feed, products,

348

and utility, respectively.

349 350

Table 3. Feed, product, and utility prices40,41 Items Crude oil Naphtha Kerosene Diesel Atmospheric gas oil Residue Fired heating Cooling water Kerosene HP steam Stripping steam

Cost 80.00 136.00 122.70 121.70 95.29 89.71 150.00 5.25 17.70 0.14

Unit $/bbl $/bbl $/bbl $/bbl $/bbl $/bbl $/kJ $/kJ $/GJ $/kmol

351 352

4. Results and discussions

353

Here, we address the UQ of a complex CDU by integrating a rigorous simulation with a

354

global SA. The column is simulated rigorously in HYSYS, while the independent input

355

variables are uniformly or normally distributed. All the variables are linked in a spreadsheet-

356

variables window in HYSYS and executed to MATLAB via the ActiveX/COM connection.

357

The ER and the OP were calculated as the outputs in HYSYS spreadsheet cells and used for

358

the global SA. The Sobol’s sensitivity indices were derived in Eq. (11).

359 360

4.1. Global SA and UQ (unfixed product flow rates with uncertain inputs of uniform distributions)

22

ACS Paragon Plus Environment

Page 22 of 40

Page 23 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

361

The independent input variables are uniformly varied along their upper and lower

362

boundaries. According to the MDRM, a set of 160 Gauss Legendre quadrature nodes was

363

generated from the orthogonal polynomial toolbox42 passed to HYSYS. Figure 3 illustrates the

364

percentage of the total SA indices obtained from the MDRM, which indicates the influential

365

factors for the two outputs. The values of total SA indices are listed in Table S2 in Supporting

366

Information. As a result, for the ER and OP, two random variables such as the flow rates of the

367

crude oil and naphtha products were detected as being critical. Other random variables become

368

non-influential factors that can be omitted in the next stage11.

369 370

Figure 3. Percentage of total SA indices obtained by the MDRM for (a) the ER and (b) the OP

371

under uniform uncertainties.

Comment [ML1]: I have changed the Figure.

372

The effective detection of the non-influential inputs from the SA step enables one to

373

simplify the model. Two simulation sets of size 500 and 300 were considered for training the

374

GPR24 for the ER and the OP, respectively. Because precise estimates of the process outputs are

375

not accessible, the results of the proposed method were compared with those of the QMC method

23

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

376

with a sufficient number of simulations. For an accurate estimation of the probability for the

377

QMC method, the number of samples was selected from Chernoff bound following Tempo et

378

al.43. Figure 4 compares the density functions of the ER obtained by the GPR (with 500 training

379

data points) and by the QMC (with 10,000 simulations) methods with two influential random

380

inputs (flow rates of crude oil and naphtha). Furthermore, these density functions were compared

381

with that by the QMC method with all 16 random inputs (using the 10,000 simulations from

382

Halton sequence). Similarly, as observed in Figure 5, the density functions of the OP were

383

obtained using the GPR and the QMC methods with two critical random variables, namely, the

384

flow rates of crude oil and naphtha product, which were then compared with those by the QMC

385

method with all 16 random inputs. Table 4 lists the statistical properties and computational time

386

achieved by the proposed (MDRM + GPR) and QMC methods for the ER and the OP. The

387

computational time of the proposed method includes the computational time for running 160

388

simulations for SA permitted by the MDRM and for constructing two surrogate GPR models

389

with 500 quadrature nodes for UQ. The QMC/MC methods require a large number of

390

simulations (approximately 10,000) to estimate the expected values, variances, and densities

391

accurately; hence, they are computationally expensive. As compared with the conventional QMC

392

method, the computation time was reduced by 95 % in the proposed method. Also, it was

393

observed by comparing Figures 4 and 5 that the proposed method can achieve acceptable results

394

with only a small computational cost in comparison with the conventional QMC method.

395 396

Table 4. Statistical data and simulation time obtained by the proposed and QMC methods

397

under uniform uncertainties Methods

No. of simulations

Runtime (sec.)

Mean µ (ER)

24

ACS Paragon Plus Environment

Mean µ (OP)

Page 24 of 40

Page 25 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

QMC – sixteen inputs for the ER and OP QMC – two inputs for the ER QMC – two inputs for the OP Proposed (MDRM+GPR) for the ER Proposed (MDRM+GPR) for the OP

Industrial & Engineering Chemistry Research

10,000

11,877.4

255,729.4

1,742.66

10,000

11,809.7

255,196.8

-

10,000

11 820.2

-

1730.84

660

786.7

255,182.2

-

460

555.6

-

1730.62

398 399 400

25

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

401 402

Figure 4. Density profiles of the ER obtained by the GPR and QMC methods under

403

uniform uncertainties.

404 405 406

26

ACS Paragon Plus Environment

Page 26 of 40

Page 27 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

407 408

Figure 5. Density profiles of the OP obtained by the GPR and QMC methods under

409

uniform uncertainties.

410

Figure 6 illustrates the accuracy of the surrogate GPR model with respect to the number

411

of fitting points by using the Nash-Sutcliffe model efficiency coefficient. The closer the Q value

412

is to 1, the more accurate the surrogate can be. As seen in Figure 6, the Q value seems to be

413

constant at 500, which is a reasonable value for the fitted points. The Nash-Sutcliffe model

414

efficiency coefficient is defined as44: 10000

415

Q =1−

∑ (µ

i

i =1 10000

∑ (f

− fi ) 2 ;

i

− f)

2

f =

1 10000 ∑ fi 10000 i =1

i =1

27

ACS Paragon Plus Environment

(21)

Industrial & Engineering Chemistry Research

416

where µi is the prediction of the trained GPR model at the ith test point of the test sample of size

417

of 10,000 nodes from Halton sequence, fi is the real model output from HYSYS with the ith test

418

point.

0.9447420

0.9447415

0.9447410

0.9447405

Q

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.9447400

0.9447395

0.9447390

0.9447385 0

200

300

400

500

600

n fitting points

419 420

100

Figure 6. Nash-Sutcliffe model efficiency coefficient vs. number of training points

421 422

4.2. Global SA and UQ (unfixed product flow rates and uncertain inputs of Gaussian

423

distributions)

424

In this case study, 16 input variables were assumed to be of Gaussian distribution with the same

425

means and variances as in the uniform distribution case. The ER may vary due to numerous

426

factors, and thus, the refinery engineers have to adjust its operations frequently to maintain its

427

profitability. Although the optimal solution based on the nominal parameter values were

428

assumed to achieve the highest profit, it possibly leads to serious losses in other cases. Therefore, 28

ACS Paragon Plus Environment

Page 28 of 40

Page 29 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

429

a reliable optimization approach should take uncertainties into account and guarantee that its

430

solution is robust enough to parameter perturbations. Table S3 in Supporting Information shows

431

the mean and standard deviation of the parameter distributions for the objective. After

432

identifying the influential variables, the proposed GPR method used 500 training simulations.

433

From the simulation result, the surrogate GPR model was obtained for ER using the theory

434

described in Section 2.2. Table 5 lists the statistical properties of the ER obtained from the GPR

435

and MC/QMC methods. To obtain accurate estimations, the QMC sampling requires 10,000

436

simulations on these 16 parameters independently according to their probability distribution.

437

Figure 7 compares the density functions of the ER obtained from the proposed and QMC

438

methods. The results from the proposed method closely match with those from the traditional

439

QMC method. Note that the number of simulations required by the proposed method is

440

significantly lower than those by the QMC method (about 20 times). Also, another unique

441

advantage of the proposed method is that Sobol’s sensitivity indices, which are used to identify

442

the influential inputs in the propagation of process uncertainty, can be obtained directly from the

443

surrogate analytical model. Thus, the number of simulations is significantly reduced from

444

10~100 times, which leads to a drastic decrease in the computational efforts needed for UQ.

445

Note that to obtain the Sobol indices by using the MC/QMC methods, a huge number of

446

simulations are required which hurdles its application in a simplifying model for UQ. In this

447

example, the SA indices for the process output with respect to 16 uncertain inputs (see Table S4

448

in Supporting Information) were obtained from the MDRM for the Gaussian uncertainties.

449

Interestingly, the SA indices identify three important inputs, the flow rates of crude oil feed,

450

naphtha, and diesel products. Others are non-influential and can be maintained at the median

451

values for the uncertainty propagation. It is worth to say that the Sobol indices obtained from the

29

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

452

Gaussian distribution were similar to those estimated by uniformly distributed uncertainties.

453

Note that in the Gaussian distribution case, independent variables have no limits of upper and

454

lower boundaries (the mean and variance of the input are selected to be equal to the case of

455

uniform distribution). For the case associated with independent variables as continuous variables

456

with upper and lower boundaries, the GPR method can make the prediction when the inputs

457

exceed the fixed boundaries, and using Eq. (19), one can quantify the reliability of the prediction

458

at the outlier point. The variance Σ is expected to be large in this case, which indicates the

459

unreliable prediction.

460 461

Table 5. Statistical data and simulation time obtained by the proposed and QMC methods

462

under Gaussian uncertainties

No. of

Runtime

simulations

(sec.)

QMC – sixteen inputs

10,000

11,777.7

256,208.6

QMC – two inputs

10,000

11,736.4

257,483.7

Proposed (MDRM+GPR)

660

775.8

257,416.2

Methods

Mean µ (ER)

463 464

30

ACS Paragon Plus Environment

Page 30 of 40

Page 31 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

465 466

Figure 7. Density profiles of the ER obtained by the GPR and QMC methods under Gaussian

467

uncertainties.

468 469

4.3. Global SA and UQ (fixed product flow rates with uncertain inputs of uniform

470

distributions)

471

In this study, the product flow rates of naphtha, diesel, kerosene, and AGO were specified

472

at the upper bounds when the maximum flow rates of each product were obtained. Twelve

473

random inputs were run at 120 simulations with 120 Gauss-Legendre quadrature nodes to derive

474

the Sobol indices for the outputs; Table S5 in Supporting Information lists the Sobol’s total SA

475

indices. As a result, the flow rate of the crude oil was the most influential factor for both ER and

476

OP when product flow rates were fixed. 31

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

477

To train the GPR for the simplified model, 100 training data points were used (from

478

Halton sequence). Figures 8 and 9 present the density function results for the ER and the OP,

479

respectively. The density functions of the ER obtained by the GPR and QMC methods (with

480

10,000 simulations) were compared with those by the QMC method with all 12 random inputs

481

(using the 10,000 simulations from Halton sequence). Similarly, the density functions of the OP

482

obtained by the GPR and QMC methods with two influential random variables (the flow rates of

483

crude oil and naphtha product), were compared with those by the QMC method with all 12

484

random inputs. A summary of the statistical results and computational time for this case is

485

presented in Table 6.

486

Table 6. Statistical data and simulation time obtained by the proposed and QMC methods

487

with fixed product flow rates Methods

No. of simulations

Runtime (sec.)

Mean µ (ER)

Mean µ (OP)

QMC - twelve inputs

10,000

11,830.4

254,381.6

1721.16

QMC - one input

10,000

11,597.7

254,385.3

-

QMC - one input

10,000

11,528.5

-

1720.64

proposed (MDRM+GPR) for the ER

220

264.6

254,382.9

-

proposed (MDRM+GPR) for the OP

220

267.3

-

1720.75

488

32

ACS Paragon Plus Environment

Page 32 of 40

Page 33 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

489 490

Figure 8. Density profiles of the ER obtained by the GPR and QMC methods with fixed product

491

flow rates under uniform uncertainties.

33

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

492 493

Figure 9. Density profiles of the OP obtained by the GPR and QMC methods with fixed product

494

flow rates under uniform uncertainties.

495 496

To summarize, in order to reduce the cost of simulations, one can replace the rigorous

497

simulation model with an inexpensive surrogate based on conventional GPR. However, the

498

number of simulations required to learn a generic multivariate response can increase significantly

499

as the input dimension increases. Furthermore, the optimization of hyper-parameters for the GPR

500

can be stuck straightforwardly in a local solution when a large number of random inputs are

501

considered. This paper addresses the dimensionality problem with a two-stage approach: Firstly,

502

the inputs are screened with the MDRM method; secondly, the GPR is constructed with respect 34

ACS Paragon Plus Environment

Page 34 of 40

Page 35 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

503

to the critical inputs to reduce the number of simulations required to learn the model. The

504

proposed approach is used for modeling the ER and the OP of CDU. Under this circumstance,

505

the result demonstrated that the flow rates of crude oil and naphtha product are critical for the ER

506

and the OP. Prediction results by the proposed method are verified by comparison with those of

507

the QMC method. It is demonstrated that there is precise agreement between the proposed and

508

the QMC methods. Moreover, the proposed method is computationally superior to the standard

509

QMC method.

510 511

5. Conclusions

512

This work focused on UQ and SA for testing the model robustness against uncertainties

513

in the CDU with the following objectives: obtaining a better understanding of the effects of

514

inputs on the output, determining the key inputs that influence the uncertainty of model outputs

515

(such as the EP and OP), and achieving model simplification. To address a large number of

516

uncertainties in the CDU with the computational limitation, an efficient two-stage SA approach

517

combining the MDRM and GPR was presented and successfully applied to UQ and SA of

518

complex CDU problem. Density functions of the outputs predicted by the Gaussian regression

519

from the simplified model showed a good agreement with those of the conventional QMC

520

methods. The proposed method is superior to the popular QMC approach primarily in terms of

521

the computational efforts (approximately more than 95 %). The SA and UQ results indicate that

522

the ER and the OP are mostly influenced by the fluctuation of the flow rates of crude oil and

523

light fraction (e.g., naphtha) in the CDU. The proposed approach is expected to offer an effective

524

method to address both SA and UQ for other complex chemical processes associated with a large

525

number of uncertainties. In future work, another study of SA is planned in order to focus on the 35

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

526

scheduling of the blending of crude oils to emulate an actual refinery optimization process and

527

therefore provide more realistic results regarding the practical application.

528 529

Acknowledgment

530

This study was supported by Engineering Development Research Center (EDRC) funded

531

by the Ministry of Trade, Industry & Energy (MOTIE) (No. N0000990), the Priority Research

532

Centers Program through the National Research Foundation of Korea (NRF) funded by the

533

Ministry of Education (2014R1A6A1031189), and the R&D Center for Reduction of Non-CO2

534

Greenhouse Gases (201700240008) funded by the Ministry of Environment as a ‘Global Top

535

Environment R&D Program’.

536 537 538 539

Supporting Information

The supporting information contains Sobol’s sensitivity analysis indices from all case studies.

540 541 542

Conflict of Interest

The authors declare no competing financial interest.

543

36

ACS Paragon Plus Environment

Page 36 of 40

Page 37 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

544

Literature cited

Industrial & Engineering Chemistry Research

545

(1) Kong, Q.; Shah, N. Development of an Optimization-Based Framework for Simultaneous Process

546

Synthesis and Heat Integration. Ind. Eng. Chem. Res. 2017, 56, 5000-5013.

547

(2) Saltelli, A.; Ratto, M.; Tarantola, S.; Campolongo, F. Update 1 of: Sensitivity Analysis for Chemical

548

Models. Chem. Rev. 2012, 112, PR1-PR21.

549

(3) Basak, K.; Abhilash, K. S.; Ganguly, S.; Saraf, D. N. On-Line Optimization of a Crude Distillation

550

Unit with Constraints on Product Properties. Ind. Eng. Chem. Res. 2002, 41, 1557-1568.

551

(4) Osuolale, F. N.; Zhang, J. Energy efficiency optimisation for distillation column using artificial

552

neural network models. Energy 2016, 106, 562-578.

553

(5) Kim, Y. H. An Energy-Efficient Crude Distillation Unit with a Prefractionator. Chem. Eng. Technol.

554

2017, 40, 588-597.

555

(6) Waheed, M. A.; Oni, A. O. Performance improvement of a crude oil distillation unit. Appl. Therm.

556

Eng. 2015, 75, 315-324.

557

(7) Mahalec, V.; Sanchez, Y. Inferential monitoring and optimization of crude separation units via

558

hybrid models. Comput. Chem. Eng. 2012, 45, 15-26.

559

(8) More, R. K.; Bulasara, V. K.; Uppaluri, R.; Banjara, V. R. Optimization of crude distillation system

560

using aspen plus: Effect of binary feed selection on grass-root design. Chem. Eng. Res. Des. 2010, 88,

561

121-134.

562

(9) Enríquez-Gutiérrez, V. M.; Jobson, M.; Ochoa-Estopier, L. M.; Smith, R. Retrofit of heat-integrated

563

crude oil distillation columns. Chem. Eng. Res. Des. 2015, 99, 185-198.

564

(10) Ochoa-Estopier, L. M.; Jobson, M.; Smith, R. The use of reduced models for design and

565

optimisation of heat-integrated crude oil distillation systems. Energy 2014, 75, 5-13.

566

(11) Gu, W.; Huang, Y.; Wang, K.; Zhang, B.; Chen, Q.; Hui, C.-W. Comparative analysis and

567

evaluation of three crude oil vacuum distillation processes for process selection. Energy 2014, 76, 559-

568

571.

569

(12) Doostan, A.; Owhadi, H. A non-adapted sparse approximation of PDEs with stochastic inputs. J,

570

Comput. Phys. 2011, 230, 3015-3034.

571

(13) Uusitalo, L.; Lehikoinen, A.; Helle, I.; Myrberg, K. An overview of methods to evaluate uncertainty

572

of deterministic models in decision support. Environ. Model. Softw. 2015, 63, 24-31.

573

(14) Sobol′, I. M. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo

574

estimates. Math. Comput. Simulation 2001, 55, 271-280.

575

(15) Saltelli, A. Making best use of model evaluations to compute sensitivity indices. Comput. Phys.

576

Commun. 2002, 145, 280-297.

37

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

577

(16) Webster, M. D.; Tatang, M. A.; McRae, G. J., Application of the Probabilistic Collocation Method

578

for an Uncertainty Analysis of a Simple Ocean Model. Report 4 (MIT Joint Program on the Science and

579

Policy of Global Change). In 1996; Vol. 32.

580

(17) Abubakar, U.; Sriramula, S.; Renton, N. C. Reliability of complex chemical engineering processes.

581

Comput. Chem. Eng. 2015, 74, 1-14.

582

(18) Schobi, R.; Sudret, B.; Wiart, J. Polynomial-chaos-based Kriging. Int. J. Unc. Quant. 2015, 5, 171-

583

193.

584

(19) Shen, D. E.; Braatz, R. D. Polynomial chaos-based robust design of systems with probabilistic

585

uncertainties. AlChE J. 2016, 62, 3310-3318.

586

(20) Duong, P. L. T.; Ali, W.; Kwok, E.; Lee, M. Uncertainty quantification and global sensitivity

587

analysis of complex chemical process using a generalized polynomial chaos approach. Comput. Chem.

588

Eng. 2016, 90, 23-30.

589

(21) Celse, B.; Costa, V.; Wahl, F.; Verstraete, J. J. Dealing with uncertainties: Sensitivity analysis of

590

vacuum gas oil hydrotreatment. Chem. Eng. J. 2015, 278, 469-478.

591

(22) Duong, P. L. T.; Minh, L. Q.; Pham, T. N.; Goncalves, J.; Kwok, E.; Lee, M. Uncertainty

592

quantification and global sensitivity analysis of complex chemical processes with a large number of input

593

parameters using compressive polynomial chaos. Chem. Eng. Res. Des. 2016, 115, Part A, 204-213.

594

(23) Minh, L. Q.; Duong, P. L. T.; Goncalves, J.; Kwok, E.; Lee, M. A two-stage approach of

595

multiplicative dimensional reduction and polynomial chaos for global sensitivity analysis and uncertainty

596

quantification with a large number of process uncertainties. J. Taiwan Inst. Chem. Eng. 2017, 78, 254-264.

597

(24) Clifton, L.; Clifton, D. A.; Pimentel, M. A. F.; Watkinson, P. J.; Tarassenko, L. In Gaussian

598

process regression in vital-sign early warning systems, 2012 Annual International Conference of the

599

IEEE Engineering in Medicine and Biology Society, Aug. 28 2012-Sept. 1 2012, 2012; 2012; pp 6161-

600

6164.

601

(25) Rasmussen, C. E.; Williams, C. K. I., Gaussian Processes for Machine Learning (Adaptive

602

Computation and Machine Learning). The MIT Press: 2005.

603

(26) Quirante, N.; Javaloyes, J.; Caballero, J. A. Rigorous design of distillation columns using surrogate

604

models based on Kriging interpolation. AlChE J. 2015, 61, 2169-2187.

605

(27) Kajero, O. T.; Chen, T.; Yao, Y.; Chuang, Y.-C.; Wong, D. S. H. Meta-modelling in chemical

606

process system engineering. J. Taiwan Inst. Chem. Eng. 2017, 73, 135-145.

607

(28) Kaymaz, I. Application of kriging method to structural reliability problems. Struct. Safety 2005, 27,

608

133-151.

38

ACS Paragon Plus Environment

Page 38 of 40

Page 39 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

609

Industrial & Engineering Chemistry Research

(29) Liu, X.; Guillas, S. Dimension Reduction for Gaussian Process Emulation: An Application to the

610

Influence of Bathymetry on Tsunami Heights. SIAM/ASA Journal on Uncertainty Quantification 2017, 5,

611

787-812.

612

(30) Bengio, Y.; Delalleau, O.; Roux, N. L. In The curse of highly variable functions for local kernel

613

machines, Advances in Neural Information Processing Systems 18 (NIPS'05), 2005; Weiss, Y.; Schölkopf,

614

B.; Platt, J., Eds. Cambridge, MA: MIT Press: 2005; pp 107- 114.

615

(31) Kalaitzis, A. A.; Lawrence, N. D. A Simple Approach to Ranking Differentially Expressed Gene

616

Expression Time Courses through Gaussian Process Regression. BMC Bioinformatics 2011, 12, 180.

617

(32) Kocijan, J., In Modelling and Control of Dynamic Systems Using Gaussian Process Models,

618

Springer Cham, 2016.

619

(33) Saltelli, A.; Chan, K.; Scott, E. M., Sensitivity analysis. Wiley New York: 2000; Vol. 1.

620

(34) Zhang, X.; Pandey, M. D. An effective approximation for variance-based global sensitivity analysis.

621

Reliab, Eng. Sys. Safety 2014, 121, 164-174.

622

(35) HYSYS, A. Aspen HYSYS customization guide. http://www.aspentech.com/core/aspen-hysys.aspx

623

(36) Errico, M.; Tola, G.; Mascia, M. Energy saving in a crude distillation unit by a preflash

624

implementation. Appl. Therm. Eng. 2009, 29, 1642-1647.

625

(37) Chang, A.-F.; Pashikanti, K.; Liu, Y. A., Atmospheric Distillation Unit. In Refinery Engineering,

626

Wiley-VCH Verlag GmbH & Co. KGaA: 2012; pp 57-116.

627

(38) Javaloyes-Antón, J.; Ruiz-Femenia, R.; Caballero, J. A. Rigorous Design of Complex Distillation

628

Columns Using Process Simulators and the Particle Swarm Optimization Algorithm. Ind. Eng. Chem. Res.

629

2013, 52, 15621-15634.

630

(39) Ibrahim, D.; Jobson, M.; Guillén-Gosálbez, G. Optimization-Based Design of Crude Oil Distillation

631

Units Using Rigorous Simulation Models. Ind. Eng. Chem. Res. 2017, 56, 6728-6740.

632

(40) Osuolale, F. N.; Zhang, J. Thermodynamic optimization of atmospheric distillation unit. Comput.

633

Chem. Eng. 2017, 103, 201-209.

634

(41) Turton, R., Analysis, Synthesis, and Design of Chemical Processes. Prentice Hall: 2012.

635

(42) Gautschi, W., Orthogonal Polynomials: Computation and Approximation. Clarendon Press: 2004.

636

(43) Tempo, R.; Calafiore, G.; Dabbene, F., Randomized algorithms for analysis and control of

637

uncertain systems: with applications. Springer Science & Business Media: 2012.

638

(44) N. Moriasi, D.; G. Arnold, J.; W. Van Liew, M.; L. Bingner, R.; D. Harmel, R.; L. Veith, T. Model

639

Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations. Transactions

640

of the ASABE 2007, 50, 885.

641 642

39

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

643

Table of Contents (TOC) Graphic

644

645

40

ACS Paragon Plus Environment

Page 40 of 40