Data-Driven Chance Constrained and Robust Optimization under

Subscriber access provided by Nanyang Technological Univ

Article

Data-driven Chance Constrained and Robust Optimization under Matrix Uncertainty Yi Zhang, Yiping Feng, and Gang Rong Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.5b04973 • Publication Date (Web): 06 May 2016 Downloaded from http://pubs.acs.org on May 10, 2016

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 free 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 accessible to all readers and 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 61

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

Data-driven Chance Constrained and Robust

2

Optimization under Matrix Uncertainty

3

Yi ZHANG1, Yiping FENG1, Gang RONG1*

4 5

1

State Key Laboratory of Industrial Control Technology, Department of Control Science and Engineering, Zhejiang University, Hangzhou, 310027, China

6 7

ABSTRACT

8

To solve optimization problems with matrix uncertainty, a novel optimization

9

approach is proposed based on chance-constrained and robust optimization, which

10

focuses on constraints with continuous uncertainty, especially with matrix uncertainty.

11

In chance-constrained approach, constraints with matrix uncertainty are always

12

regarded as joint chance constraints (JCCs), which can be simplified into individual

13

chance constraints (ICCs) and can be further reformulated into algebraic constraints

14

by robust methods. Motivated by reformulation of chance constraints with right-hand

15

side (RHS) uncertainty, a novel formulation of constraints with LHS uncertainty is

16

proposed, where the uncertainty is described as intervals related to the confidence

17

level of chance constraints. Through using Kernel Density Estimation (KDE),

18

confidence sets of uncertain parameters are built to approximate unknown true

1

ACS Paragon Plus Environment


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

1

probability density functions. The approach is illustrated with a motivating and a

2

process industry scheduling example with energy consumption uncertainties.

3

KEYWORDS: Matrix uncertainty; Data-driven; Chance-constrained

4

1. INTRODUCTION

5

Optimization under uncertainty is an important and challenging topic in the process

6

systems engineering (PSE) community. There are many unpredictable and stochastic

7

factors which could cause uncertainties of production process, from which continuous

8

uncertainty mostly can be reflected by uncertain parameters with unknown probability

9

density functions. In general, according to the source of uncertainty in processes,

10

uncertainty can be classified into four categories: model-inherent, process-inherent,

11

external and discrete.1 While in optimization problems, uncertainty has two common

12

forms which are called right-hand side and matrix uncertainty.2 Similarly, for a

13

constraint with uncertain parameters, according to the position of uncertain

14

parameters, uncertainty is mostly classified into right-hand side (RHS) uncertainty

15

and left-hand side (LHS) uncertainty. To some extent, for dealing with uncertainties,

16

optimization approaches can be divided into two levels: passive optimization and

17

reactive optimization. The former considers uncertainties into the initial mathematical

18

model, while the latter prefers real-time optimization based on simulation or

19

re-optimization after uncertain events occurred. Current solutions of optimization

20

under continuous uncertainty mainly belong to passive approaches, such as robust

21

optimization, stochastic programming3, which introduce uncertain factors into

22

optimization model by random scenarios, chance constraints or other formulations.

2


Page 2 of 61

Page 3 of 61

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


1

One of the most classic methods of solving problems with uncertainties is to replace

2

them by intervals or bounded parameters. Meanwhile, scenario-based approaches are

3

also acceptable which will come to a bounded result finally. But when we focus on

4

continuous uncertainties, results from random scenarios cannot reflect the real

5

influence on the objective, as a result of the absence of probability density/distribution

6

of uncertainties. Thus, the approach of chance-constrained optimization was

7

introduced by Charnes (1958)4 and further discussed by many researchers5-7, and

8

recently Prékopa (2013)8 take the probability density and distribution of uncertain

9

parameters into consideration and define confidence intervals of the objective. While

10

the chance-constrained approach also has to face difficulties in the reformulation and

11

calculation of the original model, such as:

12

13

reformulation of chance-constrained problems with different kinds of uncertain parameters;

14

estimation of probability density/distribution of uncertainties;

15

choosing confidence coefficients of chance constraints;

16

solving individual and joint chance-constrained problems.

17

To deal with optimization problem which contains more than one uncertain

18

parameters, Li (2008) proposed a solution that the problem can be relaxed into an

19

equivalent nonlinear optimization problem9, based on the assumption that the

20

probability density of the uncertain parameter is known. Generally, with the help of

21

chance-constrained approach, when constraints with uncertainties are considered

22

individually, the optimization problem can be regarded as an individual

23

chance-constrained problem. While if the constraints must satisfy all possible values

24

of the uncertainties at the same time, the problem will become a joint 3



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

1

chance-constrained problem. Nemirovski (2006) has proved that joint constraints can

2

be decomposed into individual constraints10, then, the feasibility and the least

3

conservative solution were further discussed by Chen (2010)11. At the same time,

4

Ben-Tal (2000)12 and Bertsimas and Sim (2004)13 proposed a robust approach to

5

describe linear constraints with uncertainties, while Chen (2007)14 used the approach

6

to formulate linear chance constraints with probability bounds, and the boundary

7

depends on the amount of uncertainties of each chance constraint. Thus, the robust

8

approach can be used to solve chance-constrained problems.15

9

Another problem for optimization under uncertainty is to estimate probability

10

density/distribution of uncertainties. The most common method is driven by historical

11

data of all related uncertain factors. A data-driven chance-constrained approach was

12

proposed by Jiang (2015)16, which defined data-driven chance constraints (DCCs)

13

with density-based confidence set and illustrated the proposed approach by numerical

14

examples under RHS uncertainties. In different conditions, non-parametric estimation

15

methods, such as kernel density estimation (KDE), can be used to estimate the

16

probability density of the uncertainty. Through kernel smoothing, Calfa (2015) solved

17

individual and joint chance-constrained optimization problem under right-hand side

18

uncertainty17, which motivates us to research the solution of optimization problems

19

under matrix uncertainty. According to current researches, continuous uncertainties

20

are mostly discussed, which always arise from dynamic changes in production

21

operations and some are inherent, such as field uncertainty, product demand

22

uncertainty and capacity uncertainty of production equipments. In this paper, we take

23

energy consumption uncertainty as an example.

4


Page 4 of 61

Page 5 of 61

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


1

The content of this paper is organized as follows. With the help of

2

chance-constrained approach, uncertainties are reclassified into four categories, and a

3

brief statement of optimization problems under different kinds of uncertainties is

4

presented in Section 2. In Section 3, we introduce a data-driven optimization approach

5

on the basis of online kernel density estimation. In Section 4, the reformulation of

6

DCCs with different kinds of uncertainties is discussed in detail, where the robust

7

optimization approach is introduced to reformulate ICCs into algebraic constraints.

8

The proposed approach is illustrated with a motivating example and a case study from

9

the real world process industry in Section 5. At last, Section 6 draws conclusions of

10

the paper and indicates prospective for future researches.

11

2. OPTIMIZATION PROBLEM UNDER MATRIX UNCERTAINTY

12

In this paper, we are dealing with linear optimization problem under matrix

13

uncertainty in the form of Eq.(1), which contains n decision variables and m uncertain

14

constraints. min

15

s.t.

z = c′x Αx ≤ b%

(1)

x∈ X 16

In the uncertain model, Α represents matrix uncertainty, containing uncertain

17

constraint coefficients a% ij , for i = 1, K , m , j = 1, K , n , while b% is an uncertain vector

18

on the right side, which is the right-hand side uncertainty. And X represents a

19

computable bounded convex set in

n

.

20

Generally, feasibility is a fundamental issue of problem (1), which means that every

21

constraint should be satisfied for any possible value of Α and b% . Thus we can

5



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 6 of 61

1

research reformulation approach on individual constraints with LHS uncertainty in the

2

form of Eq.(2), not considering uncertainty on the right-hand side at first. n

3

∑ a%

ij

x j ≤ bi , i = 1,K , m

(2)

j =1

4

Similarly, constraints with RHS uncertainties can be described by Eq.(3). n

5

∑a x ij

j

≤ b% i , i = 1,K , m

(3)

j =1

6

With the help of chance-constrained approach, confidence level can be introduced

7

to the constraints in Eq.(2) and Eq.(3), thus the constraints can be represented as

8

individual chance constraints in Eq.(4), which we define as individual case.

9

 n  Ρ  ∑ a% ij x j ≤ bi  ≥ 1 − α i , i = 1,K , m  j =1 

(4-a)

10

 n  Ρ  ∑ aij x j ≤ b% i  ≥ 1 − α i , i = 1,K , m  j =1 

(4-b)

11 12

If all the chance constraints need to satisfy a certain confidence level 1 − α at the same time, then we get joint case in Eq.(5).

13

 n  Ρ  ∑ a% ij x j ≤ bi , i = 1,K , m  ≥ 1 − α  j =1 

(5-a)

14

 n  Ρ  ∑ aij x j ≤ b% i , i = 1,K , m  ≥ 1 − α  j =1 

(5-b)

15

According to the amount of uncertain parameters and the coupling degree of

16

uncertainties in the constraints, optimization problems are reclassified into four

17

categories (as shown in Fig. 1), which are Individual Case with RHS uncertainty,

18

Joint Case with RHS uncertainty, Individual Case with LHS uncertainty and Joint

19

Case with LHS uncertainty. As proved by Nemirovski(2006), joint chance 6


Page 7 of 61

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


1

constraints(JCCs) can be relaxed and decomposed into individual chance

2

constraints(ICCs). In general, compared to constraints with RHS uncertainty, LHS

3

uncertainty contains more uncertain parameters in a constraint, which cannot be easily

4

reformulated into algebraic form. Normally, Joint Case with LHS uncertainty is the

5

most complex in the four categories, and the reformulation approach of the four kinds

6

of cases will be discussed in Section 4 in detail. Optimization problem under uncertainty Optimization under Matrix Uncertainty Joint Case Constraints with LHS uncertainty Decomposed

Si milar decomposition strategy Complexity increased

Complexi ty increased

Optimization under Right-side Uncertainty Joint Case Constraints with RHS uncertainty Complexity increased

Individual Case Constraints with LHS uncertainty

Decomposed

Individual Case Constraints with RHS uncertainty

7 8

Figure 1. Classification of optimization problems under uncertainty

9

(based on chance-constrained approach)

10

In optimization problems under matrix uncertainty, we should deal with the

11

difficulty in calculating and the coupling between uncertainties. For optimization

12

problem under RHS uncertainty, Calfa (2015) proposed a solution based on individual

13

and joint chance-constrained optimization. In this paper, we improve this method and

14

propose a series of reformulation strategy in following sections for problems under

15

matrix uncertainty. Before presenting the proposed approach, we introduce the initial

16

chance-constrained formulation with the above uncertainties based on the review of

17

modeling data-driven chance constraints (DCCs).

18

3. OPTIMIZATION APPROACH WITH DCCs 7



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 8 of 61

1

As discussed above, one of the key problems of optimization under uncertainty is to

2

estimate probability density or distribution functions of uncertainties. Based on the

3

modeling of data-driven chance constraints proposed by Ben-Tal (2011)18 and Jiang

4

(2015)16, we introduce an improved data-driven optimization approach, not only for

5

reformulation of DCCs with RHS uncertainty (Calfa, 2015)17.

6

3.1 Confidence Set Based on KDE

7

Assuming that every continuous uncertainty can be described by any probability

8

density, but with a distance that can measure the similarity between the real density

9

and the estimated result. Jiang(2015)16 built a confidence set for the uncertainty,

10

which is dependent on the specified tolerance. When the tolerance of the confidence

11

set gets smaller, we can get estimated result with higher accuracy. The density-based

12

confidence set is defined as follows:  Θ = Ρ ∈ Μ + : DKL f 

(

13

14

where ξ% ∈

K

)

f ≤ d, f =

dΡ   d ξ% 

(6)

represents K-dimensional random vectors within Α and b% ; Ρ K

induced by ξ% ; Μ + represents the set

15

represents the probability distribution on

16

of all probability distributions; d represents the pre-specified tolerance, which can

17

be approximated by χ M2 −1,1− β 2 N with large N, where χ M2 −1,1− β represents the

18

100 (1 − β ) % (e.g., β =0.05,0.10,0.15) percentile of the χM2 −1 distribution16. Here, M

19

represents the amount of intervals for constructing the whole interval of uncertain

20

parameter, while in KDE, M equals to the span of sampling data divided by the

21

bandwidth h ; N is the amount of effective sampling points, which means the bigger

22

the sample size is, the tolerance d is closer to zero. 8


Page 9 of 61

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


1 2

For modeling the distance between density functions, the Kullback-Leibler (KL) divergence was introduced to build the confidence set, which is defined as Eq.(7).

(

3

DKL f

)

f =∫

( )

f ξ% log

K

( ) ( )

f ξ% d ξ% f ξ%

(7)

4

where f and f denote the true density function and its estimation, respectively.

5

Obviously, when the distance is closer to zero, the estimated result is more accurate.

6 7

Here, we choose KDE to estimate the unknown density function, which is defined as f h ( x) =

1 n 1 n  x − Xi  Kh ( x − X i ) = ∑ ∑K  n i =1 nh i =1  h 

(8)

8

where Kh ( x ) = 1 K  x  , K ( ⋅) is the kernel, which is a non-negative function that

9

integrates to one and has mean zero, and h > 0 is a smoothing parameter called the

10

bandwidth. If Gaussian is chosen as the kernel, the bandwidth h can be calculated by

11

the following formulas (Silverman, 1998)19:

h

h

1

12

13 14

1  4σ 5  5 −  ≈ 1.06σ n 5 h=  3n   

(9)

Thus, we propose an online KDE approach in Fig. 2, which can help deal with dynamic estimation problem of production uncertainties.

9



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

Select the kernel function (Gaussian, Logistic, etc.) Initialize the sampling parameters (Period T, Frequency f ) Start sampling (t=0, k=0)

YES

Is the current sample point outlier?

t = t +1

NO

k = k + 1; t = t + 1

β k = β tsample

NO

t > Tf ? YES

Calculate the bandwidth h N = k , h ≈ 1.06σ N

( )

−

1 5

f N β = K ernel ( β 1 , β 2 , L , β N

1 2

( )

)

( )

f β = fN β

Figure 2. Flow chart of online sampling and kernel density estimation

3

For data-driven problems in the real world, sometimes the uncertainty is subject to a

4

definite distribution just for a short period, which means only the estimation within

5

this period is valuable. For example, when the production processing scheme is

6

changed, the dynamic data related to the machines would not reflect the uncertainty

7

we care about, such as field uncertainty. Thus, we propose the online estimation

8

approach in Fig. 2, which is configurable and contains a preprocessing step in the

9

process of sampling, and all the sampling points are online measurements taken from

10

the real-world production process. In this approach, T represents the total sampling

10


Page 10 of 61

Page 11 of 61

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


1

period, which is consistent with the scheduling period; f represents the sampling

2

frequency which is related to the period of batch cycle in the production process; t is

3

an integer parameter which is used for counting sampling points βtsample ; similarly, k

4

is another integer parameter for counting effective sampling points β k ; At last, the

5

density functions of the uncertainties will be estimated when t > Tf is satisfied.

6

Then we can get the kernel-based probability density function f N β

7

regarded as the final estimated result f β .

8 9 10

( )

, which is

( )

3.2 DCCs with Density-based Confidence Set With the help of chance-constrained optimization, the uncertain problem in Eq.(1) can be described as a model as follows: min

11

s.t.

z = c′x Ρ Αx ≤ b% ≥ 1 − α

{

}

(10)

x∈ X 12

which means the inequality constraints need to be satisfied by at the least

α value represents the risk level (or tolerance of constraint

13

probability of 1−α , and

14

violation) allowed by decision makers. The chance constraints in Eq.(10) are required

15

to be satisfied under every probability distribution in the confidence set Θ . Thus, the

16

constraints can be represented as follows:

{

}

inf Ρ Α x ≤ b% ≥ 1 − α

17

P∈Θ

(11)

18

⋅ implies the worst distribution in the confidence set Θ . where the operator inf {}

19

Through estimation, assuming that we can get a Ρ$ which satisfy the critical

20

P∈Θ

condition of the confidence set, then the constraints in Eq.(11) can be reformulated as

11



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

{

5

(12)

where α ′ is a reduced risk level, calculated by formula16

 e− d x1−α − 1    x −1 

α ′ = 1 − inf  x∈( 0,1)

3

4

}

P Αx ≤ b% ≥ 1 − α ′

1 2

Page 12 of 61

(13).

Obviously, the formula shows that when d is getting closer to zero, the reduced risk level α ′ is closer to the nominal risk level of original chance constraints.

6

In general, in individual cases with RHS uncertainty, every constraint can be

7

reformulated separately and for each constraint, only one uncertain parameter needs to

8

be dealt with and estimated. Usually, uncertainty in this kind of problems has limited

9

influence on the objective. However, for joint cases, when several constraints contain

10

RHS uncertainty at the same time, all the uncertain constraints should be satisfied for

11

every possible value of the RHS parameters. With data-driven chance-constrained

12

approach, this joint case with RHS uncertainty can be formulated as

  P ∑ aij x j ≤ b% i , i = 1,K , m  ≥ 1 − α ′  j 

13

14 15

Similarly, for left-hand side uncertainty, every single chance constraint is with an

(

uncertain row vector a% k = a% k1 ,K , a% kn

)

like

  P  ∑ a% ij x j ≤ bi  ≥ 1 − α i′ , i = 1,K , m  j 

16

17

(14).

Thus,

problem

under

matrix

uncertainty

usually

(15).

refers

to

the

joint

18

chance-constrained problem with LHS uncertainty, and the constraints can be

19

formulated into Eq.(16).

20

  P ∑ a% ij x j ≤ bi , i = 1,K , m  ≥ 1 − α ′  j  12


(16)

Page 13 of 61

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


1

which will be reformulated to algebraic constraints in the next section by robust and

2

joint chance-constrained approach.

3

4. ROBUST

4

REFORMULATION

OF

DCCs

UNDER

MATRIX

UNCERTAINTY

5

The data-driven chance-constrained formulations in Eq.(14-16) help us understand

6

problems under uncertainty more clearly, thus we can research on reformulation of

7

DCCs in three phases, where individual case and joint case with RHS uncertainty are

8

merged into one kind of problem. Many researches have been put forward to solve

9

problems under RHS uncertainty, which will be introduced briefly at first. Next, we

10

discuss data-driven robust reformulation of individual DCCs with LHS uncertainty.

11

At last, problem under matrix uncertainty will be taken as joint case with LHS

12

uncertainty, and the constraints will be reformulated as JCCs of every single

13

data-driven chance constraint.

14

4.1 Category 1&2: Optimization under Right-hand side Uncertainty

15

4.1.1 Individual Case with RHS Uncertainty

16

According to the formulation of joint case with RHS uncertainty, we can extract

17

every individual chance constraint from Eq.(14) as constraints with RHS uncertainty

18

in the form of

19

20 21

  P ∑ aij x j ≤ b% i  ≥ 1 − α i′ , i = 1,K , m  j 

(17).

The constraints represent that the probability of the random variable (r.v.) b% i to achieve a value greater than or equal to

∑a x ij

j

must be at least 1 − α i′ for

j

22

i = 1, K , m . If we use the estimated CDF F b% i ( ⋅) and its inverse CDF (quantile

13



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 61

−1

1

function) F b% i ( ⋅) to describe the constraints in Eq.(17), they will be transformed into

2

the equivalent constraints:   F b% i  ∑ aij x j  ≤ α i′ , i = 1,K , m  j 

3

4

(18)

or

∑a x

5

ij

−1

j

≤ F b% i (α i′ ) , i = 1,K , m

(19).

j

6

Likely, for constraints in the form of   P  ∑ aij x j ≥ b% i  ≥ 1 − α i′ , i = 1,K , m  j 

7

8

(20),

it can also be transformed into

∑a x

9

ij

−1

j

≥ F b% i (1 − α i′ ) , i = 1,K , m

(21).

j

10

Noted: Because of the original constraints are linear, the reformulated constraints in

11

Eq.(19) and Eq.(21) remain linear, which means the reformulated optimization model

12

can be solved by the original algorithm when the uncertainty is on the right-hand side.

13

4.1.2 Joint Case with RHS Uncertainty

14

In the joint case, the difficulty mainly lies in the reformulation of joint chance

15

constraints. As defined in Eq.(14), the joint chance constraint has to be reformulated

16

into algebraic forms for calculating. Commonly, there are two solutions:

17 18 19

transform JCCs to approximate ICCs with reduced risk levels, which need to satisfy the following precondition:

∑α

i

≤α

i

14


(22)

Page 15 of 61

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


1 2 3 4

where

α

is the nominal risk level of JCC, and α i is the reduced risk levels

for ICCs before the reformulation.

calculate the equivalent integrated constraint Eq.(23) on the basis of estimated joint CDF of the random variables (r.vs.) b% i for i = 1, K , m ;   F b%1 ,K,b% m  ∑ a1 j x j ,K , ∑ amj x j  ≥ 1 − α ′ j  j 

5

(23)

6

In the former solution, the reformulated constraints remain linear, but can only get

7

the approximate results. However, although the latter solution is an equivalent

8

transformation, the nonlinearity arises from the integrated form on the left side. Thus,

9

before solving the joint chance-constrained problem, we should choose a proper

10

solution for reformulation of DCCs which will be referred and discussed in detail for

11

Category 3.

12

4.2 Category 3: Individual Case under Matrix Uncertainty

13

4.2.1 Decoupling of LHS Uncertain Parameters

14

For analyzing matrix uncertainty, we need to discuss the reformulation of individual

15

DCCs with LHS uncertainty defined as Eq.(15) first, where we should mainly focus

16

on the decoupling of the LHS uncertainties. Obviously, when the reformulation of

17

constraints in Eq.(15) is based on the joint probability density function (PDF), the

18

reformulated will certainly become nonlinear constraints, which is difficult for

19

calculating. In fact, in real-world problems, it is more important to measure the

20

influence of uncertainties on the objective value by simple optimization method,

21

especially focusing on the boundaries. As a result, we tend to replace the uncertain

22

parameters on the left-hand side by interval numbers. However, it is worthwhile to

23

discuss how to determine the intervals. 15



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 16 of 61

1

In some researches, people tend to find the worst case by scenario-based method,

2

which is to generate a number of scenarios with a series of random values for

3

uncertain parameters. For the problem with matrix uncertainty in Eq.(2), the most

4

conservative case is when all the uncertain parameters are treated as the upper bound

5

a ij , such as

∑a

6

j ∈J i

7

ij

x j + ∑ aij x j ≤ bi

, i = 1, K , m

(24),

j∉ J i

where J i represents the set of coefficients in row i that are subject to uncertainty,

8

and this formulation is more likely to become infeasible in some extreme cases. The

9

constraints will also be influenced by outlier points of uncertain parameters

10

significantly, because the upper bound is measured from sampling data.

11

However, robust method focuses on the possible objective values between the worst

12

case and the original deterministic case. At first, if we only consider the influence of

13

one uncertain parameter a% ij* on the left side in Eq.(15), the other uncertain

14

parameters are all regarded as deterministic parameters, here we replace them by their

15

median value aijmed , then we get Eq.(25).

16

 P  a% ij* x j* + 

17 18

19

 aijmed x j + ∑ aij x j ≤ bi  ≥ 1 − α i′ , i = 1,K , m j∉J i j ≠ j * , j ∈J i 

∑

(25)

Note that only when the decision variables x j are all non-zero positive, we can get an equivalent algebraic formulation of Eq.(25) as

Fa%−ij1* (1 − α i′ ) ⋅ x j* +

∑

j ≠ j* , j∈J i

aijmed x j + ∑ aij x j ≤ bi

, i = 1,K , m

(26).

j∉ J i

20

Compared with the case in Eq.(24), here we’ve generated another scenario for the

21

constraints in Eq.(15), which is more relaxed and the influence of outlier points 16


Page 17 of 61

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


1

decreased obviously. Motivated by this formulation, if the confidence level is greater

2

than 0.5, we will get a more conservative formulation as Eq.(27), where the uncertain

3

parameter is replaced by its quantile value.

∑ F (1 − α ′ ) ⋅ x + ∑ a

4

j∈ J i

−1 a% ij

j

ij

⋅ x j ≤ bi

, i = 1, K , m

j∉ J i

(27)

5

Here, we should note that this formulation is not an equivalent formulation of

6

chance constraints in Eq.(15), which is only a particular scenario (called extreme case

7

in the following sections) and it motivates us to consider intervals which are related to

8

the quantile values.

9

4.2.2 Robust Reformulation of Constraints with LHS Uncertainty

10

Motivated by the formulation in Eq.(27), we propose a robust method as follows,

11

which is based on bounded uncertain parameters and all the parameters are described

12

as intervals related to their quantile values.

13 14 15

STEP 1: Describe the uncertainty a% ij as a symmetric interval  aijc − aijw , aijc + aijw  with risk level α i and confidence level β i ; The risk level

αi

has been defined in the data-driven chance constraint which

16

represents the probability tolerance of constraint violation, but when the DCC is

17

reformulated into estimated chance constraint as Eq.(12), the reduced risk level

18

relies on the tolerance d

19

confidence level β i is specified, the tolerance d

20

determined. In a word, the risk level α i represents how much we consider the

21

uncertainty, while the confidence level

22

is to the real density function of the uncertainty.

αi′

of the confidence set. As illustrated in 3.1, once the

βi

of the confidence set is

represents how close the estimated result

17



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 18 of 61

1

Thus, through kernel density estimation, the probability density functions of

2

uncertainties can be obtained. Here, we connect the confidence level 1 − α i′ to the

3

span of uncertain intervals, which is to make sure the probability of random value

4

lying in the span is nearly 1 − α i′ . Thus when the confidence level becomes bigger,

5

the constraint becomes more conservative, because the probability limit is stricter and

6

the constraint need to satisfy more possible values on the boundaries. In other words,

7

this formulation magnifies the change in conservatism when the confidence level

8

changes, which can reflect the influence of uncertainties more clearly.

9

For researching optimization problems under uncertainty, we usually choose the

10

confidence level greater than 0.8. As a result, the uncertainty a% ij can be described as

11

a two side estimation  Fa%−1 (α i′ 2 ) , Fa%−1 (1 − α i′ 2 ) , which is better than one side   ij

12

ij

estimation with lower and upper bounds like  a ij , Fa%−1 (1 − α i′ )  or  Fa%−1 (αi′) , aij  , because 

ij





ij



13

the two side estimation contains the most possible values and make sure the

14

reformulated model is not affected by extreme values on the boundary. If the PDF of

15

uncertain parameters are centered around the median value, we can determine

16

intervals by the above method. Otherwise, the form of intervals should be

17

reconsidered, but the span should also cover most possible values and the probability

18

should be greater than or equal to the confidence level.

19

Then we can transform the original interval to a symmetric interval like

20

c w  aijc − aijw , aijc + aijw  , where aij represents the center of the interval, aij represents the

21

width of the interval, which can be calculated by

22

aijc =

Fa%−ij1 (αi′ 2 ) + Fa%−ij1 (1 − α i′ 2 ) 2 18


(28)

Page 19 of 61

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


a = w ij

1

Fa%−ij1 (1 − α i′ 2 ) − Fa%−ij1 (α i′ 2 ) 2

2

STEP 2: Transform the DCCs into constraints with uncertain intervals.

3

We can use

4

(29).

a% ij =  aijc − aijw , aijc + aijw  to represent uncertain intervals on the left

side of the ICCs in Eq.(16), thus, the ICCs can also be reformulated into

∑

5

a% ij ⋅ x j ≤ bi

, i = 1, K , m

(30)

j

6

which can also be described as

7

∑( a

c ij

+ aijw ⋅ eI ) ⋅ x j ≤ bi , i = 1,K, m

j

∑a x + ∑a e x c ij j

8

w I ij

j

9

j

≤ bi , i = 1,K, m

(31)

j

where e I represents an identity interval [-1,1] in the interval number theory.

10

STEP 3: Reformulate constraints into algebraic constraints by the robust approach.

11

Strictly, constraints in Eq.(31) are still not algebraic constraints, but are close to the

12

standard forms in robust optimization. Thus, when the original problem has LHS

13

uncertainties, the optimization model can be further transformed into Eq.(32) by the

14

robust approach13. min s.t.

c′x ∑ a x j + ∑ pij + zi Γi ≤ bi

i = 1,K , m

zi + pij ≥ aijw x j −yj ≤ xj ≤ yj

∀i, j ∈ J i ∀j

lj ≤ xj ≤ uj pij ≥ 0 yj ≥ 0 zi ≥ 0

∀j ∀i, j ∈ J i ∀j ∀i

c ij

j

15

j∈J i

19


(32)


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

1

where

Page 20 of 61

Γi represents the upper limit of the amount of considered uncertain

2

parameters. In other words, the reformulated model is protected against all cases that

3

up to Γi  of the uncertain coefficients are allowed to change.

4

As proved by Bertsimas (2004)13, the reformulated chance constraints has a lower

5

probability limit related to Γi . The probability that the ith constraint is not violated

6

satisfies

7

  Pr  ∑ a% ij x j ≤ bi  ≥ 1 − B ( n , Γ i )  j 

8

where B ( n, Γ i ) is the best possible bound; n = J i

9 10

,which is the amount of

uncertain parameters in the ith constraint. Here, the boundary can be approximated by the binomial distribution:  Γ −1  1 − B ( n, Γ i ) ≈ Φ  i  = f (Γi )  n 

11

12

(33)

When n is determined, we can just regard the boundary as a function of

(34)

Γi . At

13

last, we need to make sure the reformulated constraints are stricter than individual

14

constraints in Eq.(15), which means when the reformulated constraints are satisfied,

15

the original model cannot be violated. Thus, when choose the parameter Γ i , the

16

following conditions must be satisfied:

f ( Γi ) ≥ 1− αi′ ∀i

17 18 19 20

(35)

Generally, we choose f ( Γi ) = 1 − α ′ for the reformulated model, which remains a linear optimization model.

4.3 Category 4: Joint Case under matrix uncertainty 20


Page 21 of 61

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


1

As a joint chance-constrained case, constraints with matrix uncertainty can be

2

reformulated into ICCs like Eq.(15) when Eq.(22) is satisfied. Then the ICCs can be

3

transformed into algebraic forms in Eq.(32), which is linear and easy to be solved.

4

In the following case, constraints with LHS uncertainty can be equally transformed

5

to the formulation in Eq.(27). When the uncertain coefficients a% ij in a constraint are

6

% i with different weights w like subject to the same distribution a ij a% ij = wij a% i

7 8

(36).

Then the DCCs in Eq.(15) will be simplified to  P  a% i ∑ wij x j +  j∈ J i

9

10

∀i, j ∈ Ji

When

∑w x ij

j

∑a j∉ J i

ij

 x j ≤ bi  ≥ 1 − α i′ , 

i = 1, K , m

(37).

remains positive, the constraint can be reformulated into algebraic

j∈J i

11

forms: Fa%−i 1 (1 − α i′ ) ⋅ ∑ wij x j +

12

j∈ J i

∑a

ij

x j ≤ bi

∀i

(38).

j∉ J i

13

After introducing the reformulation of DCCs under different kinds of uncertainties,

14

we will further discuss its conservatism and application in production scheduling

15

problems. Also, the accuracy of estimated results of the uncertainties will certainly

16

affect optimization results, which will be analyzed through the following numerical

17

examples.

18

5. NUMERICAL EXAMPLES

19

For illustrating the proposed approaches, two examples are introduced in this

20

section. The traditional State-Task Network (STN) case20 is cited as a motivating

21

example, and its continuous-time formulation has been further researched by 21



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 22 of 61

1

Ierapetritou (1998)21. Here, the discrete formulation will be reconsidered with

2

uncertainties and the influence of estimated results on the objective value is analyzed

3

at first. The second example is an industrial case with uncertainties in energy

4

consumption from the real world ethylene plant. On the basis of sampling data in

5

three

6

chance-constrained model are also illustrated.

7

months,

the

data-driven

estimation

approach

and

reformulation

of

5.1 Motivating Example

8

In this example, three feedstocks are used to produce two different products through

9

five processing stages: heating, reactions1, 2, 3 and separation. Four intermediate

10

products (Hot A, Int BC, Int AB, and Impure E) are produced during the process. In

11

the state-task network, the nine kinds of materials are represented as State 1~9, and

12

the processing stages are regarded as Task 1~5. Meanwhile, the capacity of the

13

equipment (Unit 1~4) has been defined by Kondili (1993)20 .Finally, the network is

14

shown in Fig. 3. Product 1 S8 40% 2h S1 Feed A

Heating (Task 1)

1h

S4

40%

Reaction 2 (Task 3)

Int AB 60% 2h

10% 2h

S5

Hot A

Int BC

60%

Impure E

S6

S7 80%

2h S2

50%

1h

Reaction 1 (Task 2)

Separation (Task 5)

90% 1h

S9

Product 2

Reaction 3 Feed C

(Task 4)

Feed B

15 16 17 18

50%

S3

20%

Figure 3. State-task network structure for the motivating example To optimize the scheduling of this process in 12 hours, the deterministic multi-period optimization model is given as follows:

22


Page 23 of 61

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


1 2 3

The objective is to maximize the total profit as Eq.(39a), which contains profit from final products and costs for inventory, raw material and utility. H

s

4 5

H

t =1

7 8

t =1

s

s

t =1

u

(39a)

t

The material balance constraints: State s , t = State s , t −1 +

∑ρ ∑B i,s

j∈ K i

i∈T s

6

H

max OBJ = ∑ ∑ pricesproduct ⋅ Ds ,t − ∑ ∑ Cost ssto,t ⋅ States ,t − ∑ ∑ Cost sfeedstock ⋅ Rs ,t − ∑ ∑ Costuutility ⋅ PCut ,t ,t ,t

i , j , t − pis

− ∑ ρ i,s i∈Ts

∑B

i , j ,t

− D s , t + R s ,t

∀s, t

(39b)

j∈ K i

Besides the difference of the amount produced and used in the state, product deliveries and raw material receipts are also considered. Allocation constraints:

9

∑W

i , j ,t

(39c)

≤ 1 ∀j , t

i∈ I j

t + pti −1

10

∑ ∑W t ′=t

i ′∈I j

i ′ , j ,t ′

− 1 ≤ M (1 − Wi , j ,t ) ∀i, j ∈ K i , t

(39d)

11

where M is a sufficiently large positive number. These constraints illustrate that

12

any unit can only start at most one task at any given time, and once the task is started,

13

any other tasks cannot be started until the current task is finished.

14 15

Capacity constraints: W i , j , tV i ,min ≤ B i , j , t ≤ W i , j , tV i ,mj ax j

16 17 18 19 20

0 ≤ States ,t ≤ Cs

∀ i , j ∈ K i , t = 1, 2, L , H

∀s , t

(39e) (39f)

The constraints illustrate capacity limitations of the amount of material processed and stored, and if Wi, j,t = 0 , the Bi, j,t are forced to zero. Utilities constraints:

∑ ∑ (γ% i

)

Wi , j ,t + β u ,i, j Bi, j ,t ≤ PCu ,t ∀u, t

u ,i , j

j∈Ki

23


(39g)


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

1

which is based on the assumption that the amount of utility u required by unit j

2

when processing task i in period t has a linear relationship with Bi, j,t . Here, we

3

introduce uncertainties to the utilities constraints, which means the amount of utility u

4

required by any unit is random. As a result, the optimized result from the original

5

deterministic model may be unable to satisfy the real demand of utilities. First, we

6

assume the uncertainties of the four units are subject to several specified distributions

7

as shown in Table 1.

8

Table 1

9

Specified distributions of uncertain parameters in utilities constraints

Unit

γ% u,i, j

β u ,i , j

Unit 1

N(2, 0.02)

N(0.1, 0.02)

Unit 2

N(5, 0.01)

N(0.1, 0.01)

Unit 3

N(3, 0.01)

N(0.06, 0.01)

Unit 4

N(3, 0.02)

N(0.1, 0.02)

10

In the reformulated model, the continuous uncertainties are replaced by intervals

11

with confidence levels, meanwhile, uncertain parameters are described by estimated

12

result of the sampling points from the exact distributions. In Fig. 4, we compare the

13

true density function and kernel-based estimations of γ% u,i, j for unit 1. Obviously,

14

estimation of more sampling points gets closer to the exact density function of the

15

uncertain parameter, but it will also be more complex for estimation.

24


Page 24 of 61

Page 25 of 61

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


1 2

Figure 4. Comparison of the exact density curve and estimations of sample points

3

(sample sizes: 1000 points, 800 points, 400 points).

4

For simplifying the calculation of the model, we assume that only one kind of the

5

utility is considered with a given negative value of -4 $/kg and a value of -2 $/kg is

6

given to the three feedstocks. Another important assumption is that the distributions

7

for the same unit are the same when different tasks are processed. In reality, the exact

8

distributions are usually unknown, which need to be estimated by the proposed

9

approach in Section 2. We will illustrate this approach by the real-world industrial

10

example.

11

Based on sampling points of different sizes from the exact distributions, the

12

estimated results of the uncertainties can be obtained to reformulate the constraints in

13

Eq.(39g) as follows:   Ρ$  ∑∑ γ% u,i, jWi, j ,t + β u ,i, j Bi, j ,t ≤ PCu,t ∀u, t  ≥ 1 − α ′  i j∈Ki 

(

14

15

)

or 25


(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

Page 26 of 61

  Ρ$  ∑ ∑ γ% u ,i, jWi, j ,t + β u ,i, j Bi, j ,t ≤ PCu ,t  ≥ 1 − αu′,t ∀u, t  i j∈Ki 

(

1

)

(41).

2

As discussed in previous sections, constraints with matrix uncertainties can be

3

regarded as a joint chance-constrained case in the form of Eq.(40), and can also be

4

simplified into individual chance-constrained case in Eq.(41). According to the

5

classification of the uncertainties, Eq.(41) belongs to Category 3, and Eq.(40) is the

6

case in Category 4, which also is the most complex formulation in the four categories.

7

Thus, the problem will be discussed in two parts: individual case and joint case.

8 9

According to the reformulation of DCCs in Section 4.2, the individual case can be solved by the following steps:

10

STEP 1: For any given tolerance d of the confidence set, the reduced risk level can

11

be obtained by Eq.(13), and the results are presented in Table 2.

12

Table 2

13

The reduced risk level of different confidence set.

14

d

0.075

0.085

0.1

α u ,t

0.01

0.05

0.10

0.15

0.01

0.05

0.10

0.15

0.01

0.05

0.10

0.15

α u′ ,t

0.00928

0.0464

0.0928

0.139

0.00919

0.0459

0.0919

0.136

0.00904

0.0452

0.0904

0.136

According to the estimated results from 400 sampling points, we choose the

15

tolerance of 0.1, and then the uncertainties γ% u,i, j and

16

intervals. (The result of intervals are shown in supporting information.)

17 18 19

Then the center

c

γ% u , i , j

,

β

c u ,i , j

and the width

w

γ% u , i , j

,

β

β u ,i , j

w u ,i , j

can be described as

of the intervals can be

calculated easily. STEP 2: DCCs in Eq.(41) can be transformed into the following constraints:

26


Page 27 of 61

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


∑ ∑ (γ%

1

j∈ K i

i

2

)

W i , j , t + β u , i , j Bi , j ,t + e I ∑ c

c u ,i , j

i

∑ (γ%

w u ,i , j

w

)

W i , j ,t + β u , i , j Bi , j ,t ≤ PC u ,t

j∈ K i

∀u, t

(42)

STEP 3: Then the robust formulation of the constraints in Eq.(42) is as follows:

)

(

c c  γ% u ,i , jWi , j ,t + β u ,i , j Bi , j ,t + ∑ ∑ pu ,i , j + zu ,t Γu ,t ≤ PCu ,t ∑ ∑  i j∈K i j∈Ki i  w w  % zu ,t + pu ,i , j ≥ γ u ,i , jWi ′, j ,t + β u ,i , j Bi′, j ,t   −Wi ′, j ,t ≤ Wi , j ,t ≤ Wi ′, j ,t   − Bi′, j ,t ≤ Bi , j ,t ≤ Bi′, j ,t  pu ,i , j ≥ 0   Wi ′, j ,t = 0,1  Bi′, j ,t ≥ 0   zu ,t ≥ 0 

)

(

3

∀u, t ∀u, i, j , t ∀i, j , t ∀i, j , t

(43)

∀u, i, j ∀i, j , t ∀i, j , t ∀u, t

4

As only one kind of utility is considered (dim(u)=1), there are only t individual

5

chance constraints. Depending on the mapping relationship between units and tasks,

6

each constraint contains 16 uncertain parameters. Thus when the risk level changes,

7

the choice of Γu ,t can be seen in Table 3.

8

Table 3

9

Choice of

Γu , t

when the reduced risk level α u′,t changes.

α u′ ,t

0.00904

0.0452

0.0904

1−αu′,t

0.99096

0.9548

0.9096

Γu , t

38.3

27.3

21.5

10

The reformulated models were implemented in GAMS 24.1 which has 480 discrete

11

variables and 693 continuous variables. Then the mathematical model is solved on a

12

desktop computer with the following specifications: Dell Precision T5610 with Intel®

13

Xeon® CPU E5-2609 v2 at 2.350 GHz (total 4 threads), 16 GB of RAM, and running

14

Windows 7 Enterprise. The local MIP solver is CPLEX 12.5. The R ks package is

27



1

used for kernel smoothing estimation of the uncertain parameters. The scheduling

2

results are presented in Fig. 5. Equipment

Scheduling Results 1

1

1

52

20

52

Tasks

1.Heating

Heater

2.Reaction 1 Reactor 1

2

3

4

2

3

4

58

80

80

78

80

68.75 3.Reaction 2

2

3

3

4

3

4

50

50

50

50

50

50

Reactor 2

4.Reaction 3

5.Separation 5

5

5

80

50

118.75

Still

0

1

2

3

4

5

3

6 Time (h)

7

8

9

10

11

12

4

Figure 5. Scheduling results of the twelve hours’ production process when the utility

5

uncertainty is considered.

6

For discussing the influence of sample size and risk level, we repeated the

7

experiments on both the individual case and the joint case, containing twelve

8

scenarios, respectively. In Fig. 6, we compare the objective values for the

9

reformulated models with specified distributions and with estimated distributions of

10

different sample sizes.

Individual Case: Objective values Individual Case

Objective values ($)

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 28 of 61

1940 1930 1920 1910 1900 1890 1880 1870 1860 1850

Individual Case

Individual Case

1932.56 1932.01 1927.81

1925.98 1920.73 1917.37 1907.13 1903.93 1899.32 1889.91 1885.76 1880.01

90%

95% 99% Nominal Confidence Level

11 28


Extreme

Page 29 of 61

1

(a) Objective values for exact and kernel-based models with individual chance

2

constraints.

Joint Case: Objective values Joint Case


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


1940 1930 1920 1910 1900 1890 1880 1870 1860 1850 1840

Joint Case

Joint Case

1931.56 1927.16 1920.31 1916.21 1906.89 1902.6

1899.32 1899.01 1897.12 1885.47

1882.22 1875.21

90%


Extreme

3 4

(b) Objective values for exact and kernel-based models with joint chance constraints.

5

Figure 6. Comparison of the objective values of the motivating example for different

6

scenarios (exact vs. kernel-based estimations)

7

The comparison shows the objectives values over a span of confidence levels from

8

90% to 99%, meanwhile, the result of extreme cases with constraints in Eq.(27) is also

9

involved in the comparison. Obviously, the profits for the exact model are always the

10

highest in both the individual and the joint case. Compared with the cases with

11

N=400, the objective values of cases with N=800 are closer to the profit in exact case,

12

which reveals that the more sampling points are used for the data-driven estimation,

13

the closer is the optimization result to the real world performance. In the experiment,

14

the reduced risk levels of the exact and kernel-based models are getting closer when

15

the confidence level becomes larger. In the new formulation of the uncertain model,

16

the extreme points of uncertainties are ignored, and the intervals mainly focus on the 29



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 30 of 61

1

middle part of the possible values of uncertainties. Finally, the joint case is proved

2

more conservative than the individual case. Besides, the new formulation in this paper

3

is proved more reasonable than the formulation in Eq.(27).

4

5.2 Process Industry Example

5

In this section, we take the real-world ethylene plant as our industry example,

6

mainly focusing on the entire production process through cracking furnaces, cooling

7

section, compression section and separation section. Here, four kinds of crude oil

8

(ETHA, NAP, AGO, HVGO) are considered, and nine main products, which are

9

shown at the right of Fig. 7, can be obtained finally. And the objective of the

10

optimization problem in the plant is to obtain ethylene as much as it can, meanwhile,

11

penalties such as cost of energy and crude oil are all considered. Generally, five kinds

12

of energy materials are utilized to run the petro-chemical industry, which are water,

13

electricity, steam, wind and fuel. Natural gas

Ethane Methane Hydrogen

Propane BA101

Ethylene

Fuel gas BA102

ETHA

BA103

Demethanizer

Quench tower

BA104

Deethanizer

Ethylene fractionator

C4

BA105

Butadiene

BA106

NAP

Compression Section

Cooling Section

BA107 BA108

Depropanizer

C5 Propylene fractionator

BA109

Quench tower

Debutanizer

Benzene

BA110

AGO

Separation Section

BA111

Cracked gasoline

Depentanizer

BA112

Propylene

BA113

HVGO

MS

SS

BA114

MS BA115

Natural gas

Cracked fuel oil

MS

SS

MS

Fuel gas

Steam network Bf101

Outer plant

Cracking furnace

Rectification

Gas tank

Product tank

Compressor

Distillation column

Material tank

Cooling section

14 15

Figure 7. Structure of the integrated production and utility system of ethylene plant

30


Page 31 of 61

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


1

For discussing the scheduling problem under matrix uncertainty in the real chemical

2

plant, we introduce uncertainties in fuel gas consumption of several furnaces, which

3

can be represented by constraints in Eq.(44). The related deterministic scheduling

4

model has been summarized in Appendix, and the details can be seen in Supporting

5

Information.

FGf ,t = f (FCf ,r,t , DS f ,r,t , Cot f ,r,t ) ∀f , r, t

6 7

(44)

Here, FGf ,t represents fuel consumption of furnace f during period

t , which is

8

mainly related to the amount of raw material and dilution steam consumed by furnace

9

f during period

t when processing crude oil

r

, represented by FC f ,r ,t and

10

DS f ,r ,t , repectively. Meanwhile, the outlet temperature of furnace

11

when consuming material r Cot f ,r ,t can also help predict the fuel consumption.

12

Normally, the fuel consumption can be predicted by a linear model in Eq.(45).

13

FG f ,t = ∑ a f ,r FC f ,r ,t + ∑ b f ,r DS f ,r ,t + ∑ c f ,r COT f ,r ,t + ∑ d f ,r r

14

f of period

r

r

t

(45)

r

where a f ,r , b f ,r , c f ,r , d f ,r represent coefficients of the linear model of furnace

r . But for real-world problem, the model sometimes

15

f when processing material

16

needs to be modified, such as COTf ,r ,t is invariant when the furnace is running

17

steadily, which makes the model in Eq.(39) become a bivariate model in Eq.(46).

18

% f ,r is the uncertain parameter which we concern about, because it represents Here, a

19

the relationship between fuel consumption and the amount of processing material. In

20

general, the other parameters just fluctuate around it experimental value, which can be

21

regarded as deterministic parameter. 31



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

FG f ,t = ∑ a% f ,r FC f ,r ,t + ∑ b f ,r DS f , r ,t + ∑ d ′f , r

1

r

2

where

∑ d′

f ,r

r

r

∀f , t

Page 32 of 61

(46)

r

= ∑ c f ,r COT f ,r ,t + ∑ d f ,r is a constant in the model. According to r

r

3

the sampling and estimation approach in Fig. 2, we monitor several furnaces for 80

4

days, resulting in 11520 records of FC f ,r ,t , DS f ,r ,t and FGf ,t . After removing the

5

outliers and nonsense points, we finally got 10637 valuable records. For any furnace,

6

only one kind of crude oil can be processed at the same time, as a result, we should

7

identify the processing oil from the data at first. Fortunately, the ratio of FC f ,r ,t and

8

DS f ,r ,t can help us identify, and sample results of furnace BA105 are shown in Fig.

9

8.

10 11

Figure 8. Sample data of the ratio of dilution steam consumption and raw material

12

consumption. 32


Page 33 of 61

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


1

In fact, the higher the ratio is, the more dilution steam is utilized in the pyrolytic

2

process. The ratio represents different processing schemes. For example, when the

3

ratio is within [0.55, 0.65], the furnace is processing Naphtha (NAP), and the ratio for

4

processing Heavy-vacuum gas oil (HVGO) is within [0.70, 0.75], otherwise, the

5

furnace is decoking or shifting processing schemes.

6

Thus we can identify that during B→C and D→E, the furnace BA105 is processing

7

HVGO; the period of F→G is for decoking; and the crude oil is NAP for the other

8

periods. Thus we can obtain the uncertain values of consumption rate a% f ,r through

9

the corresponding experimental model. Then through kernel density estimation, the

10

% f ,r of furnace BA105 and BA106 density curves of the uncertain consumption rate a

11

are shown in Fig. 9. Similarly, sampling data of the other furnaces can also be used to

12

identify processing condition and to estimate the probability density function of the

13

fuel consumption rate.

14

33



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 34 of 61

1

Figure 9. Kernel density estimation of consumption rate when processing NAP and

2

HVGO：(a) Furnace-BA105 processing NAP, (b) Furnace-BA105 processing

3

HVGO, (c) Furnace-BA106 processing NAP, (d) Furnace-BA106 processing HVGO.

4

For the optimization model, the fuel consumption constraints with uncertainty are as

5 6

follows:

∑ ∑ ( a% f

7 8 9

10 11

f ,r

)

FC f , r ,t + b f , r DS f , r ,t + d ′f , r ≤ FPPfuel ,t + PC fuel ,t

∀t

(47)

r

Obviously, the uncertainty in Eq.(47) makes the optimization model become a problem with matrix uncertainty, and the data-driven chance-constrained form is:  inf Ρ  ∑ ∑ a% f , r FC f , r ,t + b f ,r DS f , r ,t + d ′f ,r ≤ FPPfuel ,t + PC fuel ,t P∈Θ  f r

(

)

 ∀t  ≥ 1 − α (48) 

The joint chance constraint can be simplified to individual chance constraints with reduced risk level α i′ .

12

  Ρ$  ∑∑ a% f ,r FC f ,r ,t + b f ,r DS f ,r ,t + d ′f ,r ≤ FPPfuel ,t + PC fuel ,t  ≥ 1 − α i′ i =1,K, T (49)  f r 

13

Thus, the uncertain parameters can be replaced by intervals according to the

14

estimated density functions in Fig. 9. Then the DCCs in Eq.(49) can be reformulated

15

to algebraic forms in Eq.(50) by robust approach in Section 4.2.

16

(

)

∑∑ ( a cf ,r FC f , r ,t + b f , r DS f ,r ,t + d ′f , r ) + ∑ p f + zt Γ t ≤ FPPfuel ,t + PC fuel ,t  f r f  w ′f ,r ,t z + p ≥ a FC t f f ,r   − FC ′f , r ,t ≤ FC f , r ,t ≤ FC ′f , r ,t   pf ≥ 0  FC ′f ,r ,t ≥ 0   zt ≥ 0 

∀t ∀f , t ∀f , r , t ∀f ∀f , r , t ∀t

(49)

17

For optimizing the distribution and purchasing of fuel gas in 10 days (T=10 days),

18

the reformulated models were implemented in GAMS 24.1 which has 650 discrete 34


Page 35 of 61

1

variables and 10773 continuous variables. Several experiments are implemented for

2

discussing the probable consumption amount of fuel gas each day. The nominal

3

confidence levels are chosen as 0.99, 0.95 and 0.90. After the model is solved on the

4

same computer, the optimization results of the consumption of fuel gas are presented

5

in Fig. 10. Consumption of Fuel Gas

1,250.00 1,200.00

Consumption of Fuel Gas Per Day (t)

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


1,150.00 1,100.00 1,050.00 1,000.00

Mean Value Nominal Confidence Level=0.99 Nominal Confidence Level=0.95 Nominal Confidence Level=0.90 Extreme Case

950.00 900.00 850.00 800.00

1

2

3

4

5

6

7

8

9

10

Period (day)

6 7

Figure 10. Optimization Results of the consumption of fuel gas with point-wise

8

binomial confidence intervals.

9

Comparing to Energy Daily of the real plant, the optimization results match the real

10

consumption amount of fuel gas well. The probable intervals in Fig. 10 shows:

11

a) The maximum fluctuation of consumption of fuel gas in a day is 75 tons;

12

b) The probable consumption amount of fuel gas of all furnaces ranges from

13

925.04t to 1156.13t.

14

The actual statistics of energy cost from the real plant shows the consumption

15

amount of fuel gas in a day ranges from 963.32t to 1129.44t. However, the results of 35



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

1

the extreme case are more conservative, which represent more fuel gas will be cost for

2

10 days’ production. Furthermore, the optimization approach in this paper cannot

3

always be suitable and precise for any real-world problem. There are several causes

4

which can lead to deviations, for instance, new changes appeared in the real plant

5

which were not contained in historical production data; the optimization model is

6

much stricter for the actual scheduling of fuel gas resources.

7

In a word, the interval result can help schedulers and decision makers optimize the

8

distribution and purchasing of energy resources. And this model is used for

9

optimizing the profits in 10 days, which is much more sensible than the

10

average-production mode in real plants. For the other continuous uncertainties, the

11

kernel-based estimation approach can also be used to estimate their probability

12

density functions. In the real-world industry case, when matrix uncertainty is

13

introduced to the optimization model, the constraints can be reformulated and

14

simplified by the approach above. For future researches, the influence of the interval’s

15

form on the optimization result is worth for further discussion, and scenarios with

16

nonlinearity and discrete uncertainties should also be studied.

17

6. CONCLUSION

18

On the basis of classic optimization approaches for uncertain problems, such as

19

chance-constrained programming and robust optimization, we reclassify the problem

20

under uncertainty into several categories, where we mainly focuses on problems under

21

matrix uncertainty. Motivated by the reformulation of constraints with RHS

22

uncertainty, we propose a novel construction of uncertain parameters, which is to

23

replace them by intervals with confidence levels. The new formulation magnifies the

24

change in conservatism when the confidence of chance constraints changes, and the 36


Page 36 of 61

Page 37 of 61

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


1

reformulated model remains linear. Numerical examples mainly illustrate approaches

2

for solving two problems, which are kernel-based estimation of the continuous

3

uncertainty and reformulation of data-driven chance constraints (DCCs) with matrix

4

uncertainty.

5

The meaning of the study is to build optimization models closed to the actual

6

uncertain condition of real-world scenarios. Our data-driven approaches seize several

7

key characteristics of actual uncertain conditions, such as multivariate, dynamic and

8

continuous, and this is our motivation for studying problem with matrix uncertainty

9

and using sampling data to estimate continuous uncertainties. However, compared

10

with the deterministic problem, the proposed formulation is still very conservative,

11

especially when all the uncertainties are introduced in the model at the same time. For

12

the future work, the conservativeness will be discussed with an improved optimization

13

approach, moreover, improved algorithms for solving the optimization model have to

14

be researched in detail.

15

APPENDIX:

INTEGRATION

OPTIMIZATION

MODEL

16

PRODUCTION AND UTILITY SYSTEM FOR ETHYLENE PLANT

17

Mass Balance Constraints (Raw Material & Product)

18

MI c ,t = M I c ,t −1 + SM c ,t +

∑

19

MI c ,t = MI c ,t −1 +

∑

FPu , c ,t − SC c ,t −

∀ c ∈ CR , f , t

(B.1)

∑

FBu ,c ,t

∀ c ∈ CP , t

(B.2)

u∈UB

SCc,t ≥ DPc,t ∀c ∈CP, t

20

22

FC f ,c ,t

f ∈ FI u ,c

u ∈SOu ,c

21

∑

FPu , c ,t −

u ∈SOu ,c

OF

(B.3)

Mass Balance Constraints (Utilities)

PU c ,t +

∑

u∈UFB

SBu ,c ,t ≥

∑ STC

u , c ,t

+ LSc ,t

u∈UT

37


∀c ∈ CS , t

(B.4)


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

1

∑ STE

u , c ,t

∑ ED

∀c ∈ CS , t

+ ∑ EDu ,c ,t

∀c ∈ CW , t

+ LSc ,t ≥ SCc,t +

u∈UT

2

∑ WG

u ,t

u , c ,t

u∈US

∑W

+PU c ,t ≥

u∈UT

u ,t

u∈UC

Page 38 of 61

(B.5) (B.6)

u

3

WGut ,t ≥ Wuc,t ∀t,{ut , uc}∈ Mbuc,ut

(B.7)

4

EDu ,c,t = f ( FUu ,t , Tpu ,t ) ∀u ∈US , c ∈ CU , t

(B.8)

6

PMIc,t ≥ MIc,t − MISc,t ∀c, t

(B.9)

7

−PMIc,t ≤ MIc,t − MISc ∀c, t

(B.10)

8

MIc,t ≤ MIUc ∀c, t

(B.11)

5

Inventory Constraints

∑

M C c ,t =

9

FC f , c ,t ∀ c ∈ C R , f , t

(B.12)

f ∈ FI u ,c

10

Cracking Furnace Model (Part 1: Production System)

FPf ,c,t = αr , f ,c FC f ,r ,t ∀r, f , c ∈ FOf ,c , t

11

∑α

12

r , f ,c

≤ 1 ∀f , r , c ∈ FO f ,c

(B.13) (B.14)

c

13

FCu , r ,t ≥ xu , r ,t FULu , r

∀u ∈ UF , r , t

(B.15)

14

FCu ,r ,t ≤ xu ,r ,t FUU u ,r

∀u ∈UF , r , t

(B.16)

15

∑x

∀f , t

(B.17)

f , r ,t

= y f ,t

r

∑ FC

16

∀f , r ∈ FI f ,c , t

(B.18)

zu ,t +1 ≥ xu ,r ,t +1 − xu ,r ,t ∀u ∈UF , r, t

(B.19)

f , r ,t

= FU f ,t

r

17 18

Cracking Furnace Model (Part 2: Utility System)

19

FB f ,c,t = f ( FC f ,r ,t , Cot f ,r ,t , RD f ,r ,t , Fcd f ) ∀f , r , c ∈ CF , t

(B.20)

20

SB f ,c,t = f (FC f ,r ,t , Cot f ,r ,t , RDf ,r ,t , Tle f ,r ,t , Fcd f ) ∀f , r, c ∈ CS , t

(B.21)

21

Separation Unit Model 38


Page 39 of 61

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


∑

1

FC u ,c ,t = FU u ,t

∀ u ∈ US , t

(B.22)

c∈ SI u ,c

FPu ,c,t = f ( FUu ,t , Tpu ,t , Fsu,c ) ∀u ∈US , c ∈ SOu ,t , t

2 3

Boiler Model FBu , c ,t = ∑ ( Aq ,u ,c ,t XFBq −1,u , c + β q ,u , c ,t ( XFBq ,u ,c − XFBq −1,u ,c ))

4

q∈Q

(B.24)

∀c ∈ CF , u ∈ UB , t SBu , c ,t = ∑ ( Aq ,u ,c ,t XSBq −1,u ,c + β q ,u ,c ,t ( XSBq ,u ,c − XSBq −1,u ,c ))

5

q∈Q

(B.25)

∀c ∈ CF , u ∈ UB, t

∑∑A

6

≤ 1 ∀u ∈ UB , t

q ,u , c , t

(B.26)

q∈Q c∈CF

0 ≤ βq,u,c,t ≤ Aq,u,c,t ∀c ∈CF, u ∈UB, t

7

yu ,t = ∑

8

∑A

q ,u , c ,t

∀u ∈ UB, t

q∈Q c∈CF

9

(B.23)

(B.27) (B.28)

Cooling Unit Model FU u ,t = ∑ FU f ,t

10

∀f ∈ UF , u ∈ UCL, t

(B.29)

f

11 12 13

SCCu,c,t = f (FUu,t ,Timcu,t ,Tomcu,t ,Tiscu,c,t ,Toscu,c,t ) ∀u ∈UC, c ∈CS, t Turbine Model STC u , c ,t = f ( STE u , c ,t , STD u , c ,t , Tit u ,t , Tot u ,t , Tct u ,t , Pit u ,t , Pct u ,t , Po t u ,t )

16 17

W G u , t = f ( STE u , c ,t , STD u , c ,t , Tit u , t , Tot u , t , T ct u ,t , Pit u ,t , Pct u , t , P ot u , t )

19

(B.31) (B.32)

∀ u ∈ U T , c ∈ C S , t (B.33)

Compressor Model FU u ,t = ∑ f

18

∀ u ∈ UT , c ∈ C S , t

STCu,c,t = STEu,c,t + STDu,c,t ∀u ∈UT, c ∈CS, t

14 15

(B.30)

∑

FPf ,c ,t

∀ f ∈ UF , u ∈ UC ,{ f , u} ∈ M f ,u , t

(B.34)

c∈C PI u ,c

Wu,t = f (FUu,t ,Ticu,i,t ,Tocu,i,t , Picu,i,t ) ∀u ∈UC, c, t Supply and Demand Balance Constraints of Fuel 39


(B.35)


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

PCc ,t + FPPc ,t ≥

1

∑

FBu ,c ,t

Page 40 of 61

∀c ∈ CF , t

(B.36)

∀c ∈ CF , u ∈US , t

(B.37)

u∈UFB

FPPc ,t = ∑ FPu ,metha ,t + ∑ FPu , prop ,t

2

u

3

u

Objective Function Pr = ∑

Max

4

t

∑

c∈CP

pric ,t SC c ,t + ∑ t

∑

− ∑ ∑ PMI c ,t IC c − ∑ ∑ FU u ,tπ u − t

c

t

u∈us

pric ,t SC c ,t −

c∈CU

∑ ∑ pri

c ,t

c∈CR

∑ ∑ pri

c∈CU

t

SM c ,t

t

c ,t PU c , t − ∑ ∑ z u , t SEC f t

(B.38)

f

5 6

SUPPORTING INFORMATION

7

1. Calculating results of intervals of uncertain parameters in motivating example

8

2. Deterministic optimization model of process industry example

9

This information is available free of charge via the Internet at http://pubs.acs.org/.

10 11

12

Corresponding Author

13

*Tel: 086-87953145. Email: [email protected].

AUTHOR INFORMATION

14 15

ACKNOWLEDGEMENTS

16

The authors gratefully acknowledge financial support from the National Program on

17

Key Basic Research Project (973 Program) of China (2012CB720506) and Projects of

18

International Cooperation and Exchanges NSFC (61320106009).

19 20

21

1. Motivating Case

NOMENCLATURE

40


Page 41 of 61

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


1

Indices

2

i = processing tasks

3

j = equipment units

4

s

5

t = time periods

6

u = utilities

7

= material states

Sets

8

Ij = processing tasks which can be performed by unit j

9

Ki

= units capable of performing task i

10

Si

11

Si = output states of task i

12

Ts

13

T s = tasks producing material in state s

14

= input states of task i

= tasks requiring material from state s

Parameters

15

ρi ,s = proportion of input of task i from state

16

ρ i , s = proportion of output of task i from state

17

pi , s = processing time for the output of task i to state

18

pti

= processing time of task i

19

Cs

= maximum storage capacity dedicated to state s

20

Vi,min = minimum capacity of unit j when used for performing task i j

21

Vi,max = maximum capacity of unit j when used for performing task i j

22

pricesproduct = price of product in state s at time t ,t

23

Costuutility = unit cost of utility u at time t ,t

24

Costssto,t = running cost of keeping in storage a unit of material of state s at time t

s ∈ Si

s ∈ Si s ∈ Si

41



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

1

Costsfeedstock = unit cost of feedstock state s at time t ,t

2

γ% u,i, j = uncertainty factor of fixed demand for utility u by task I at time period t

3

β u ,i , j = uncertainty factor of variable demand for utility u by task I at time t,

4

related to Bi, j,t

5 6

c γ% u,i, j = center of the interval uncertainty factor γ% u,i, j when the confidence level is determined

7 8

w γ% u,i, j = width of the interval uncertainty factor γ% u,i, j when the confidence level is determined

9 10

β u , i , j = center of the interval uncertainty factor β u ,i , j when the confidence level is determined

11 12

β u , i , j = width of the interval uncertainty factor β u ,i , j when the confidence level is determined

13

Variables

14 15

c

w

Bi, j,t = amount of material which starts undergoing task i in unit j at the beginning of time period t

16

Ds,t = amount of material in product state s due for delivery at time period t

17 18

Rs ,t = amount of material of feed state s received from external sources at time period t

19

States,t = amount of material stored in state s at the beginning of time period t

20

PCu,t = amount of utility u purchased at time period t

21 22

Wi, j,t = binary variable designated the assignment of processing task i of unit j at

23

time period t

OBJ = objective value of the optimization model

24 25

2. Process Industry Case

26

Indices

27

u

= units (cracking furnaces/separation units/compressors/turbines/cooling unit) 42


Page 42 of 61

Page 43 of 61

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


1

f = cracking furnaces (subset of

2

c

= commodities

3

r

= raw material (subset of

4

t = time periods

5

i = segments of compressor

6

u)

c)

Sets

7

CF = fuel products

8

CR = crude oils

9

CRE = recycle material of ethane

10

CV = intermediate products

11

CP = final products

12

CU = utilities (steam/power)

13

CW = power

14

CS = different pressure steam

15

UB = boilers

16

UF = cracking furnaces

17

UFB = cracking furnaces and boilers

18

U S = separation units

19

UC = compressors

20

UT = turbines

21

UCL = cooling units

22

FI f ,c = input material c of furnace f

23

FOf ,c = output products c of furnace f

24

CPI u ,c = input commodities of compressor u

25

SIu,c = input material c of separation unit u

26

SOu ,c = output products c of separation unit u 43



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

1

M f ,u = units association between cracking furnace and compressor

2

Mbuc ,ut = units association between compressor and turbine

3

Parameters

4

DPc ,t = market demand of final product c of period t

5 6

Af f ,r = fuel consumption coefficient of cracking furnace f related to unit load

7 8 9 10 11 12 13 14 15 16

when cracking material r

Bf f = fuel consumption coefficient of cracking furnace f related to outlet temperature

Cf f = fuel consumption coefficient of cracking furnace f related to dilution ratio with steam

As f ,r = super steam generation coefficient of cracking furnace f related to unit load when cracking material r

Bs f = super steam generation coefficient of cracking furnace f related to outlet temperature

Cs f = super steam generation coefficient of cracking furnace f related to dilution ratio with steam

17

Awu

18 19

Bwu ,i = power consumption coefficient of compressor u related to inlet material

20 21

= power consumption coefficient of compressor u related to unit load

temperature of segment i

Cwu,i = power consumption coefficient of compressor u related to inlet pressure of segment i

22

SMc,t = supply of material c of period t

23

SPU c ,t = the amount of utility c purchased limitation of period t

24 25

α r , f ,c = the yield ratio of the product c of cracking furnace f consuming raw material r

26

FULu ,r = lower bound of capacity of unit u on processing material r

27

FUUu,r = higher bound of capacity of unit u on processing material r

28

M IU

u

= inventory capacity of commodities c 44


Page 44 of 61

Page 45 of 61

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


1

M IS c

2

Fcd f = constant factor of energy consumption of furnace f

3

IC c

4

SECf = material switching cost coefficient for cracking furnace f

5

pric,t = price of commodity c

6

XFBq ,u ,c = maximum amount of fuel c that can be consumed by boiler u in q piece

7 8

XSBq,u,c = maximum amount of steam that can be generated by boiler u by

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

= safety inventory level of commodity c

= inventory cost of commodity c

consuming fuel c in q piece

DRf ,r = dilution steam ratio of furnace f when consuming material r RDU f ,r ,t = upper bound of dilution steam consumption of furnace f when consuming material r

RDL f , r ,t = lower bound of dilution steam consumption of furnace f when consuming material r

CotU f ,r ,t = upper bound of outlet temperature of furnace f of period t when consuming material r

CotL f ,r ,t = lower bound of outlet temperature of furnace f of period t when consuming material r

TleU f ,r ,t = upper bound of outlet temperature of transfer line exchanger for material r of furnace f in period t

TleL f ,r ,t = lower bound of outlet temperature of transfer line exchanger for material r of furnace f in period t

a f ,r = coefficients which represent fuel gas consumption rate of furnace f when consuming material r

b f ,r = coefficients related to the dilution steam consumption rate of furnace f when consuming material r

c f ,r = coefficients related to the outlet temperature of furnace f when consuming material r

45



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

1 2

d f ,r , d ′f ,r = constant of equations for predicting fuel gas consumption of furnace f when consuming material r

a% f ,r = uncertainty parameter represent dynamic fuel gas consumption rate of

3 4

furnace f when consuming material r

5 6

c a% f ,r = center of the interval uncertainty parameter a% f ,r when the confidence level is determined

7 8

w a% f ,r = width of the interval uncertainty parameter a% f ,r when the confidence level is determined

9

Variables

10

FGf ,t = amount of fuel gas consumed by furnace f of time period t

11

FC f ,r ,t = amount of material r consumed by furnace f of time period t

12 13

DS f , r ,t = amount of dilution steam consumed by furnace f when processing

14 15

material r at time period t

COT f ,r ,t = outlet temperature of furnace f of time period t when consuming material r

16

FPPc,t = amount of commodity c produced at period t

17

PCc ,t = amount of commodity c purchased at period t

18

Aq,u,c,t = binary variable of piecewise segment of boiler

19

EDu , c ,t = amount of utility c consumed of period t of unit u

20

FBu , c ,t = fuel c consumed by furnace or boiler of period t

21

FCu ,c ,t = amount of commodities c consumed c of period t of unit u

22

FPu ,c ,t = amount of commodities c produced of period t of unit u

23

FUu,t = flow rate of unit u of period t

24

LSc,t

25

MCc,t = raw material consumed of period t

26

MI c ,t = material inventory of c of period t

= letdown valve flow rate of steam c

46


Page 46 of 61

Page 47 of 61

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


1 2

PMI c ,t = penalty difference between real inventory and expected inventory of commodity c

3

Pr = overall profit

4

PUc,t = amount of utility c purchased of t

5

SBu,c,t = super pressure steam generated by furnace or boiler of period t

6

SCc ,t = the amount of commodity c sold/delivered in period t

7

SCCu,c,t = steam c consumed by cooling unit u of period t

8

STCu,c,t = steam c consumed by turbine u of period t

9

STDu,c,t = steam c condensed by turbine u of period t

10

STEu,c,t = steam c extracted by turbine u of period t

11

Wu ,t = power consumption of compressor u f of period t

12

WGu,t = power generation of turbine u f of period t

13

βq,u,c,t = continuous piecewise variable of boiler u in segment q from 0 to 1

14

xu,c,t = binary variable denotes whether the unit u operates with consuming

15

material r in period t

16

yu,t = binary variable denotes whether the unit u operates in period t

17

z f ,t = binary variable denoting whether the material consumed of cracking

18

furnace is changed in the period t

19

Tpu,t = top pressure of separation unit

20

Fsu,c = separation factor of unit u

21

Ticu,i,t = inlet material temperature of segment i of compressor u in period t

22

Picu ,i ,t = inlet pressure of segment i of compressor u in period t

23

Titu,t = inlet steam temperature of back pressure turbine/condensing turbine u in

24

period t 47



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

1 2

Totu,t = extracted steam temperature of back pressure turbine/condensing turbine u in period t

3 4

Tctu,t = condensing steam temperature of back pressure turbine/condensing turbine u in period t

5 6

Pitu,t = inlet steam pressure of back pressure turbine/condensing turbine u in period t

7 8

Pctu,t = condensing steam pressure of back pressure turbine/condensing turbine u in period t

9 10

Potu,t = extracted steam pressure of back pressure turbine/condensing turbine u in period t

11

Timcu,t = inlet material temperature of cooling u in period t

12

Tomcu,t = outlet material temperature of cooling u in period t

13

Tiscu,c,t = inlet steam c temperature of cooling u in period t

14

Toscu,c,t = outlet steam c temperature of cooling u in period t

15

RDf ,r ,t = dilution steam consumption of furnace f when consuming material r

16

Tle f ,r ,t = outlet temperature of transfer line exchanger for material r of furnace f

17

in period t

18 19 20 21 22 23 24 25 26

REFERENCES (1) Pistikopoulos E N. Uncertainty in process design and operations[J]. Computers & Chemical

Engineering, 1995, 19: 553-563. (2) Bertsimas D, Thiele A. Robust and data-driven optimization: Modern decision-making under uncertainty[J]. INFORMS tutorials in operations research: models, methods, and applications for innovative decision making, 2006: 137. (3) Li Z, Ierapetritou M. Process scheduling under uncertainty: Review and challenges[J]. Computers & Chemical Engineering, 2008, 32(4): 715-727.

48


Page 48 of 61

Page 49 of 61

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


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

(4) Charnes A, Cooper W W, Symonds G H. Cost horizons and certainty equivalents: an approach to stochastic programming of heating oil[J]. Management Science, 1958, 4(3): 235-263. (5) Olson D L, Swenseth S R. A linear approximation for chance-constrained programming[J]. Journal of the Operational Research Society, 1987: 261-267. (6) Iwamura K, Liu B. A genetic algorithm for chance-constrained programming[J]. Journal of Information and Optimization sciences, 1996, 17(2): 409-422. (7) Nemirovski A, Shapiro A. Convex approximations of chance-constrained programs[J]. SIAM Journal on Optimization, 2006, 17(4): 969-996. (8) Prékopa A. Stochastic programming[M]. Springer Science & Business Media, 2013. (9) Li P, Arellano-Garcia H, Wozny G. Chance-constrained programming approach to process optimization under uncertainty[J]. Computers & Chemical Engineering, 2008, 32(1): 25-45. (10) Nemirovski A, Shapiro A. Convex approximations of chance-constrained programs[J]. SIAM Jour1nal on Optimization, 2006, 17(4): 969-996. (11) Chen W, Sim M, Sun J, et al. From CVaR to uncertainty set: Implications in joint chance-constrained optimization[J]. Operations research, 2010, 58(2): 470-485. (12) Ben-Tal A, Nemirovski A. Robust solutions of linear programming problems contaminated with uncertain data[J]. Mathematical programming, 2000, 88(3): 411-424.

18

(13) Bertsimas D, Sim M. The price of robustness[J]. Operations research, 2004, 52(1): 35-53.

19 20

(14) Chen X, Sim M, Sun P. A robust optimization perspective on stochastic programming[J].

21 22 23 24 25 26 27 28 29 30

Operations Research, 2007, 55(6): 1058-1071. (15) Margellos K, Goulart P, Lygeros J. On the road between robust optimization and the scenario approach for chance-constrained optimization problems[J]. Automatic Control, IEEE Transactions on, 2014, 59(8): 2258-2263. (16) Jiang R, Guan Y. Data-driven chance-constrained stochastic program[J]. Mathematical Programming, 2015, Series A: 1-37. (17) Calfa B A, Grossmann I E, Agarwal A, et al. Data-driven individual and joint chance-constrained optimization via kernel smoothing[J]. Computers & Chemical Engineering, 2015, 78: 51-69. (18) Ben-Tal A, Den Hertog D, De Waegenaere A, et al. Robust solutions of optimization problems affected by uncertain probabilities[J]. Management Science, 2013, 59(2): 341-357.

49



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

1 2 3 4 5 6 7 8 9

(19) Silverman, B.W. (1998). Density Estimation for Statistics and Data Analysis. London: Chapman & Hall/CRC. p. 48. ISBN 0-412-24620-1. (20) Kondili E, Pantelides C C, Sargent R W H. A general algorithm for short-term scheduling of batch operations—I. MILP formulation[J]. Computers & Chemical Engineering, 1993, 17(2): 211-227. (21) Ierapetritou M G, Floudas C A. Effective continuous-time formulation for short-term scheduling. 1. Multipurpose batch processes[J]. Industrial & engineering chemistry research, 1998, 37(11): 4341-4359. (22) Soyster A L. Technical note—convex programming with set-inclusive constraints and applications to inexact linear programming[J]. Operations research, 1973, 21(5): 1154-1157.

10

Table of Contents Graphic

11 12

For Table of Contents Only

50


Page 50 of 61

Page 51 of 61

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


Figure 1 Classification of uncertainty based on chance-constrained approach 81x59mm (96 x 96 DPI)



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

Figure 2 Flow chart of online sampling and kernel density estimation 88x165mm (96 x 96 DPI)


Page 52 of 61

Page 53 of 61

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


Figure 3 State-task network structure for the motivating example 201x97mm (96 x 96 DPI)



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

Figure 4 Comparison of the exact density curve and estimations of sample points 262x191mm (96 x 96 DPI)


Page 54 of 61

Page 55 of 61

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


Figure 5 Scheduling results of the twelve hours' production process when the utility uncertainty is considered 138x60mm (96 x 96 DPI)



Individual Case: Objective values


Individual Case EXACT 1940 1930 1920 1910 1900 1890 1880 1870 1860 1850

Individual Case N=800

Individual Case N=400

1932.56 1932.01 1927.81

1925.98 1920.73 1917.37 1907.13 1903.93 1899.32 1889.91 1885.76 1880.01

90%


Extreme

(a) Objective values for exact and kernel-based models with individual chance constraints.

Joint Case: Objective values Joint Case EXACT


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 56 of 61

1940 1930 1920 1910 1900 1890 1880 1870 1860 1850 1840

Joint Case N=800

Joint Case N=400

1931.56 1927.16 1920.31 1916.21 1906.89 1902.6

1899.32 1899.01 1897.12 1885.47

1882.22 1875.21

90%


Extreme

(b) Objective values for exact and kernel-based models with joint chance constraints.


Page 57 of 61

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


Figure 7 Structure of the integrated production and utility system of ethylene plant 469x282mm (96 x 96 DPI)



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

Figure 8 Sample data of the ratio of dilution steam consumption and raw material consumption 177x177mm (96 x 96 DPI)


Page 58 of 61

Page 59 of 61

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


Figure 9 Kernel density estimation of consumption rate when processing NAP and HVGO 268x177mm (96 x 96 DPI)



Consumption of Fuel Gas

1,250.00 1,200.00

Consumption of Fuel Gas Per Day (t)

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 60 of 61

1,150.00 1,100.00 1,050.00 1,000.00

Mean Value Nominal Confidence Level=0.99 Nominal Confidence Level=0.95 Nominal Confidence Level=0.90 Extreme Case

950.00 900.00 850.00 800.00

1

2

3

4

5

6

7

Period (day)


8

9

10

Page 61 of 61

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


86x64mm (300 x 300 DPI)


Data-Driven Chance Constrained and Robust Optimization under

Recommend Documents