Improved Benders Decomposition and Feasibility Validation for Two

Subscriber access provided by UNIV OF TEXAS DALLAS

Process Systems Engineering

Improved Benders decomposition and feasibility validation for two-stage chance-constrained programs in process optimization Yu Yang Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b04777 • Publication Date (Web): 06 Mar 2019 Downloaded from http://pubs.acs.org on March 8, 2019

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.

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 38 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

Improved Benders decomposition and feasibility validation for two-stage chance-constrained programs in process optimization Yu Yang∗ Department of Chemical Engineering, California State University Long Beach, Long Beach, CA 90840 USA E-mail: [email protected]

Abstract The two-stage stochastic program with probabilistic constraints is studied by assuming normally distributed uncertain parameters. At stage I, the decision, represented by discrete variables, is made such that normal operations can be applied at stage II with high probability. At stage II, the uncertainty is realized and the optimal decision, represented by continuous variables, is determined accordingly. Different from the conventional two-stage decision process, recovery operations are introduced at stage II if normal operations are infeasible. This two-stage decision process can be approximately modeled by a scenario-based mixed-integer linear program (MILP). Two issues of this scheme are concerned. First, with more samples, the solution of the scenario-based MILP is closer to the original problem, but the computational demands increase significantly. Second, the feasibility of such a solution under unseen scenarios cannot be guaranteed. To overcome these two limitations, an improved Benders decomposition method is first proposed to calculate the globally optimal solution of the scenario-based MILP more efficiently than existing algorithms. Second, the probabilistic feasibility of a stage I solution from the improved Benders decomposition is rigorously evaluated through a chance-constrained program with linear decision-rule. A refinery ACS Paragon Plus Environment

1

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 2 of 38

optimization problem is solved through the proposed scheme to verify the computational efficiency and probabilistic feasibility.

Introduction Finding a globally optimal solution for the process design and operations under parametric uncertainties may enhance the profitability and safety of the plant. For example, a case study for the Sarawak Gas Production System shows 3.2 billion dollars improvement on the expected net present value by using the global optimization and stochastic program. 1 A risk management scheme for the supply chain planning, also based on the global optimization and stochastic program, shows a 6% reduction on the probabilistic financial risk. 2 Several methodologies have been proposed to search a global optimum for systems under uncertainties. 3–9 A very comprehensive review of this research has been done by Sahinidis. 10 Generally speaking, approaches for the optimization under uncertainty can be divided into two categories: robust optimization and stochastic programming. The early work of robust optimization dates back to Soyster in 1973, 11 who studied the columnwise uncertainty in the linear programming and obtained a feasible solution under all possible scenarios. However, such a method may yield a very conservative solution. 12 Thus, the “price of robustness” is proposed to adjust the level of conservatism by introducing the probabilistic bounds of constraint violations. 13 This method builds an analytical approximation of the problem and achieves a near optimal solution under bounded uncertainties. Motivated by this work, a series of papers have been published to study the robust counterpart optimization with different types of uncertainty set. 14–16 Efficient algorithms are proposed to find the smallest size of uncertainty set to meet the probabilistic bound of constraints violation. 17,18 The aforementioned works mainly focus on the robust optimization for linear program whereas more complex formulations, for example, quadratic program and linear matrix inequality, are discussed by Ben-Tal et al. 19 The second category, namely, stochastic programming, explicitly takes the distribution of uncertain parameters into account, and thus results in a less conservative solution. It is particularly applicable to the cases, in which sufficient data is available to characterize the probabilistic properties of uncertain parameters. Two methodologies can be found in the literature. The first ACS Paragon Plus Environment

2

Page 3 of 38 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


one uses the probabilistic constraint to quantify the reliability of a solution. That is the so-called chance-constrained program (CCP). If the uncertain parameters satisfy multivariate normal distribution, which is the case considered in this paper, the linear CCP with individual chance constraint (ICC) can be converted into a second-order cone program and solved to the global optimum. 20 If there are joint chance constraints (JCC), a conservative approximation can be built to yield a near optimal solution. If the distribution of uncertain parameter is not normal, the sampling method has to be used and the probabilistic feasibility is guaranteed if the number of scenarios is large enough. 21 However, such a number can be very conservative and impractically large. 22 The second methodology is the two-stage stochastic program which optimizes the sequential decision process with recourse actions. At stage I, the decision is made by taking all possible scenarios into account. At stage II, the uncertainty is realized and the recourse operations are made accordingly. The two-stage stochastic program also can be approximately solved by the sampling-based method. Namely, the samples are drawn from the distribution of uncertain parameters to generate a number of scenarios. Then, the original stochastic optimization problem is approximated by a deterministic formulation consisting of all generated scenarios. This approach suffers at least two issues: feasibility and computations. The unsampled scenarios are not considered, and the resulting large-scale formulation is difficult to solve. Directly using an optimization solver to handle all scenarios simultaneously is not efficient because its computational time grows exponentially as the number of scenarios increases. Several decomposition strategies thus have been developed to solve each scenario separately. For example, Lagrangian decomposition, 23–25 Benders decomposition 26,27 and cross decomposition, 29 their computational time grows linearly as the number of scenarios increases. A more complicated case, which combines the chance constraints and two-stage stochastic program, is proposed by Liu et al. 30 for the two-stage decision process with dual mode operations. At stage I, the decision is made to maximize profits and limit the probability of undesired outcomes at stage II. Once uncertainty is realized at stage II, either normal or recovery operations, are adopted to ensure the feasibility and maximize profits. The recovery operations sometimes may even yield more profits but are physically undesired or restricted by regulations. Thus, the prob-

ACS Paragon Plus Environment

3


Page 4 of 38

abilistic constraints must be incorporated into this scheme to avoid abusing recovery operations. A scenario-based approximation can be used to address the chance constraints, and the resulting large-scale optimization can be solved through a variation of the Benders decomposition. 30 In order to obtain the tightest cutting planes, the algorithm proposed by Liu et al. 30 solves an additional optimization problem for each scenario, which may result in significant computational demands. Moreover, how to handle the unsampled scenarios is not discussed. 30 Compared with the existing method, two innovations are presented in this paper to handle two-stage chance-constrained programs: (1) We develop an improved strategy to generate cutting planes without solving optimization for each scenario. The total computational time of the Benders decomposition thereby is reduced significantly. (2) To verify the probabilistic feasibility of the stage I solution, we develop a CCP formula through the linear decision-rule. 31 If such a CCP is feasible, then the probabilistic feasibility of the stage I solution is guaranteed. Even though the probabilistic feasibility of the stage I solution can be empirically verified by sampling and solving a large number of scenarios, the proposed verification method is sampling-free, and thus more rigorous. This paper is organized as follows. The problem formulation and its sampling-based approximation are presented in Section 2. The mathematical background of the Benders decomposition is shown in Section 3. In Section 4, the proposed methods including the improved Benders decomposition and probabilistic feasibility validation are presented. A refinery optimization problem with uncertain crude oil properties and operation parameters is solved to demonstrate the effectiveness of proposed algorithms in Section 5. In the final Section, the conclusion is drawn.

Problem Statement Two-Stage Chance-Constrained Program In this paper, we solve the two-stage chance-constrained program. Both stages I and II variables are determined to minimize the expected cost (maximize expected profits). Moreover, either normal ACS Paragon Plus Environment

4

Page 5 of 38 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


or recovery operations can be selected at stage II for each possible scenario to ensure the feasibility. However, the recovery operations are not preferred and thereby limited by chance constraints. A general mathematical formulation is introduced by Liu et al.: 30

min qxT x + P(Iω = 0)Eω (f (x, ω)|Iω = 0) + P(Iω = 1)Eω (fˆ(x, ω)|Iω = 1) x,Iω

(P)

s.t. P{x ∈ Fx (ω)} > 1 − , Iω = 0 ⇒ x ∈ Fx (ω), ˆ x (ω), Iω = 1 ⇒ x ∈ F x ∈ {0, 1}Dx ∩ Π, ω ∈ Ω, Iω ∈ {0, 1},

where Π is a polyhedron representing deterministic constraints on stage I variable x; Ω is the set of random parameter: ω; qx is the coefficient of objective function; indicator decision variable Iω ˆ x (ω); Eω is the expectation over ω ∈ Ω; is the threshold for ω ∈ Ω shows if x ∈ Fx (ω) or x ∈ F of constraint violation; P is the probability; f and fˆ are the nested cost functions for normal and recovery operations at stage II, defined as follows:

T T T f (x, ω) = min{qyT y + qsT s : aT i x + ω Bi y + ci s + di y 6 0, ∀i ∈ {1, 2, . . . , Dr }, y ∈ [y, y], s ∈ [s, s]} y,s

T ˆ ˆT fˆ(x, ω) = min{qˆyT yˆ + qˆsT sˆ : a ˆ + cˆT ˆ + dˆT ˆ 6 0, ∀i ∈ {1, 2, . . . , Dr }, yˆ ∈ [ˆ y , yˆ], sˆ ∈ [ˆ s, sˆ]} i x + ω Bi y i s i y y ˆ,ˆ s

where qy ∈ RDy , qs ∈ RDs , ai ∈ RDx , Bi ∈ RDω ×Dy , ω ∈ RDω , ci ∈ RDs , di ∈ RDy , qˆy ∈ RDy , qˆs ∈ ˆ i ∈ RDω ×Dy , cî ∈ RDs and dî ∈ RDy ; Dr is the number of constraints; vectors ˆ i ∈ RDx , B RD s , a [y, s] and [ˆ y , sˆ] are the decision variables of normal and recovery modes, respectively. Throughout this paper, the overline and underline represent upper and lower bounds on a variable, respectively. The set of x, in which f (x, ω) > −∞, denoted by Fx (ω), is defined as: T T T Fx (ω) ={x|∃{y, s}, aT i x + ω Bi y + ci s + di y 6 0, ∀i ∈ {1, 2, . . . , Dr }, y ∈ [y, y], s ∈ [s, s]}.


5


Page 6 of 38

Similarly, the set of x in which fˆ(x, ω) > −∞ can be defined as:

T ˆ ˆ x (ω) ={x|∃{ˆ ˆT ˆ 6 0, ∀i ∈ {1, 2, . . . , Dr }, yˆ ∈ [ˆ y , yˆ], sˆ ∈ [ˆ ˆ + dˆT ˆ + cˆT s, sˆ]}. F y , sˆ}, a i y i s i x + ω Bi y

Here we do not assume complete recourse even with recovery operations. If fˆ(x, ω) > f (x, ω), ∀ω ∈ Ω and ∀x ∈ Π, then Iω = 0 will lead to an optimal solution and P{x ∈ Fx (ω)} = 1. In this study, we focus on a more complex case: ∀ω ∈ Ω and ∀x ∈ Π, fˆ(x, ω) 6 f (x, ω), which implies that the recovery operation is less costly. However, we can only allow the recovery operations to be used in a few scenarios. How to optimally choose these scenarios is the major challenge. Throughout this paper, we make the following assumptions: A1: The random variables satisfy multivariate normal distribution, namely ω ∈ N (µ, Σ). ˆ x (ω) ⊇ Fx (ω) and fˆ(x, ω) 6 f (x, ω), ∀ω ∈ Ω, ∀x ∈ Π. A2: F The first assumption enables us to use second-order cone programming to check the probabilistic feasibility of (P) given a solution x. The second assumption implies that recovery operation is more profitable and can handle more scenarios.

Scenario-Based Approximation The problem (P) is non-trivial since the support of random vector ω is not finite. A samplingbased method is utilized to obtain a deterministic approximation. Suppose that K scenarios are generated according to the distribution of ω with equal importance. Then the formula of a


6

Page 7 of 38 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


sampling-based approximation of (P) can be introduced from Liu et al.: 30

min x,z

s.t.

qxT x K X

K K 1 Xˆ 1 X f (x, ωk )(1 − zk ) + f (x, ωk )zk + K k=1 K k=1

(Pf)

zk 6 K,

k=1

zk = 0 ⇒ x ∈ Fx (ωk ), ˆ x (ωk ), zk = 1 ⇒ x ∈ F x ∈ {0, 1}Dx ∩ Π, ωk ∈ Ω, zk ∈ {0, 1}, ∀k ∈ {1, 2, . . . , K},

where zk replaces the Iω to indicate whether using recovery or normal operations for scenario k. A graphical illustration of (Pf) is shown in Fig. 1. The number of K should be large enough to yield a highly reliable solution for (P). However, the solving time of (Pf) thus will be very long. An improved Benders decomposition method is proposed in this paper to reduce the computational demands.

Since (Pf) is an approximation of (P), its optimal solution may not be even feasible

Figure 1: The graphical illustration of the sampling-based two-stage decision process. The master problem is solved to obtain stage I solution x. Different scenarios are generated by sampling uncertain parameters at stage II. Given x, if the number of scenarios that have to use recovery operations is less than K, then this x is feasible to (Pf). to (P). A straightforward approach to validate this solution is substituting it into a large number of newly generated scenarios and calculating the chance of feasibility. However, this empirical validation is not rigorous because the sufficient number of scenarios cannot be determined theoACS Paragon Plus Environment

7


Page 8 of 38

retically. An alternative method is proposed in this paper to verify the probabilistic feasibility of the solution more rigorously.

Mathematical Background Benders Decomposition The classical Benders decomposition 26 can efficiently solve two-stage stochastic program constructed by finite scenarios. To handle probabilistic constraints and dual operations, a variation of Benders decomposition is proposed by Liu et al. 30 In Benders decomposition, the original formula (Pf) is decomposed into a master problem and subproblems. The master problem is solved to find a lower bound of (Pf). The associated stage I solution is substituted into each subproblem/scenario and solving all subproblems can obtain an upper bound solution of (Pf) and dual cutting planes. Attaching dual cutting planes to the master problem will reduce the relaxation gap. The optimal solution will be reached by repeatedly solving master and subproblems several times. The convergence of Benders decomposition has been wellestablished, 27 and thus will not be further discussed. In the following subsections, details of the Benders decomposition will be proposed.


8

Page 9 of 38 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


Master Problem In this subsection, we describe the detail of master problem. Its formulation is shown as below:

min

x,z,η

s.t.

qxT x K X

K 1 X + ηk K k=1

(MP)

zk 6 K,

k=1 (n)

(n)

(n)

ηk > (αk )T x + βk + zk Mk , ∀k ∈ Ξ(x(n) ), ∀n ∈ {1, 2, . . . , N },

(1)

(n) (n) ˆ (n) , ∀k ∈ Ξ(x ˆ (n) ), ∀n ∈ {1, 2, . . . , N }, ηk > (α ˆ k )T x + βˆk + (1 − zk )M k

(2)

(n)

(n)

0 > (γk )T x + θk + zk Hk , ∀k ∈ Γ(x(n) ), ∀n ∈ {1, 2, . . . , N },

(3)

(n) (n) ˆ k , ∀k ∈ Γ(x ˆ (n) ), ∀n ∈ {1, 2, . . . , N }, 0 > (ˆ γk )T x + θˆk + (1 − zk )H

(4)

x ∈ {0, 1}h ∩ Π, zk ∈ {0, 1}, ∀k ∈ {1, 2, . . . , K}, (n) (n) (n) (n) (n) (n) (n) (n) where αk , α ˆ k , γk , γ ˆk , βk , βˆk , θk , θˆk are parameters derived from k th subproblem at (n) ˆ (n) (n) ˆ (n) nth iteration; x(n) is the solution of (MP) at nth iteration; Mk , M are non-positive k , Hk , Hk

coefficients of zk derived at nth iteration, which will be discussed later; Ξ(x(n) ) and Γ(x(n) ) are ˆ (n) ) and Γ(x ˆ (n) ) are also mutually exclusive; ηk is the lower bound value mutually exclusive; Ξ(x of operating cost at k th scenario and it will finally converge to f (x, ωk ) or fˆ(x, ωk ); (1) and (2) are the normal and recovery optimality cuts for k th scenario, derived at nth iteration; (3) and (4) are the normal and recovery feasibility cuts for k th scenario, derived at nth iteration to preclude infeasible {x(n) , z (n) }. Note that zk exists in all cutting planes to indicate if that constraint is associated with normal or recovery operations. Solving (MP) is time-consuming when Dx + K is too large. To address this issue, Luedtke 28 proposes a branch-and-cut method to solve (MP). The integer restriction on zk is relaxed, and thus only the binary variables x have to be considered. The non-integer zk in the intermediate solution will be branched at each iteration to reduce the relaxation gap. However, how to determine the branching sequence is not discussed in the reference.


9


Page 10 of 38

Subproblems and Cutting Planes At nth iteration, solving (MP) or (RMP) will obtain a solution of stage I variable: x(n) . We substitute x(n) into each scenario k and determine the optimal stage II solution under normal and recovery operations, respectively.

min

yk ,sk

qyT yk + qsT sk

(SPNk )

(n) T + ωkT Bi yk + cT s.t. aT i x i sk + di yk 6 0, ∀i ∈ {1, 2, . . . , Dr },

yk ∈ [y k , y k ], sk ∈ [sk , sk ].

min

y ˆk ,ˆ sk

qˆyT yˆk + qˆsT sˆk

(SPRk )

(n) ˆ i yˆk + cˆT sˆk + dˆT yk 6 0, ∀i ∈ {1, 2, . . . , Dr }, ˆT s.t. a + ωkT B i x i i

yˆk ∈ [ˆ y k , yˆk ], sˆk ∈ [ˆ sk , sˆk ]. For compact expression, following notations are used:

A = [a1 , a2 , . . . , aDr ]T , Gk = [ωkT B1 ; ωkT B2 ; . . . ; ωkT BDr ], C = [c1 , c2 , . . . , cDr ]T , D = [d1 , d2 , . . . , dDr ]T , ˆ k = [ω T B ˆ 1; ωTB ˆ 2; . . . ; ωTB ˆ Dr ], ˆ = [ˆ ˆ 2, . . . , a ˆ Dr ]T , G A a1 , a k k k ˆ = [ˆ ˆ = [dˆ1 , dˆ2 , . . . , dˆDr ]T . C c1 , cˆ2 , . . . , cˆDr ]T , D If subproblems (SPNk ) and (SPRk ) are feasible at nth iteration, then k belongs to Ξ(x(n) ) ˆ (n) ). Objective values f (x(n) , ωk ) = objSPN (x(n) ), fˆ(x(n) , ωk ) = objSPR (x(n) ) and dual and Ξ(x k k ˆ ∈ RDr to be dual multipliers cutting planes in (1) and (2) can be obtained. Let λ ∈ RDr and λ of constraints containing the x term in (SPNk ) and (SPRk ), respectively. The coefficients of x in (n) (n) (n) T ˆ (n) T ˆ ˆ (n) cutting planes are (αk )T = (λk )T A and (α k ) = (λk ) A, respectively. The constants βk


10

Page 11 of 38 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


(n) and βˆk are calculated as below:

(n)

(n)

βk = objSPNk (x(n) ) − (αk )T x(n) , ∀k ∈ Ξ(x(n) ) (n) (n) ˆ (n) ) ˆ k )T x(n) , ∀k ∈ Ξ(x βˆk = objSPRk (x(n) ) − (α

(n)

The optimality cuts shown in (1) and (2) have parameters Mk

ˆ (n) . If zk = 1, the recovery and M k

T ˆ(n) ˆ (n) operations are employed and thus the cost of scenario k is determined by ηk > (α k ) x + βk (n)

(n)

only. If zk = 0, then the cost at scenario k is determined by ηk > (αk )T x + βk

only. To model

the impact of zk , Liu et al. 30 introduce the following set Ψk for each scenario: Ψk = {x ∈ {0, 1}Dx ∩ Π, zk ∈ {0, 1}, ηk : ηk > (1 − zk )f (x, ωk ) + zk fˆ(x, ωk )} (n)

Mk

(5)

ˆ (n) should be designed to make inequalities (1) and (2) valid for the set Ψk . According and M k

T ˆ (n) to A2: fˆ(x, ω) 6 f (x, ω), the optimality cutting plane generated by (SPRk ) satisfies (α k ) x+ (n) ˆ (n) = 0 makes (2) valid βˆk 6 fˆ(x, ωk ) 6 f (x, ωk ), ∀x ∈ {0, 1}Dx ∩ Π, ∀ωk ∈ Ω. It implies that M k

for Ψk because when zk = 0, the right-hand side of (2) will not exceed f (x, ωk ). A highly negative (n)

value on Mk

can make (1) valid for Ψk , but it may increase the computational time of (MP) or

lower the objective value of (RMP). Therefore, the following Lemma is introduced 30 to determine (n)

Mk . Lemma 1. Define an optimization problem and its objective value: (n)

Vk

(n)

= min{fˆ(x, ωk ) − (αk )T x, x ∈ {0, 1}Dx ∩ Π ∩ Fx (ωk )}. x

(n)

Then the constraint (1) with Mk (n)

(n)

= Vk (n)

(n)

One can show that (αk )T x+βk +Mk

(n)

− βk

(BM)

is a valid cut of the set Ψk .

6 fˆ(x, ωk ) based on Lemma 1. Namely, if zk = 1, then

the right-hand side of optimality cut (1) is not greater than fˆ(x, ωk ), ∀x ∈ {0, 1}Dx ∩ Π ∩ Fx (ωk ). In such a situation, ηk will be only determined by the optimality cut of the recovery operation and finally converge to fˆ(x, ωk ).


11


Page 12 of 38

Proposed Methods Relaxed Master Problem The global optimum of (MP) is a lower bound solution of (Pf). However, it is not easy to solve (MP) since there are Dx +K binary variables. As the number of scenario increases, solving (MP) becomes more challenging. In the Benders decomposition, new optimality and feasibility cuts are attached to (MP) at each iteration, and the resulting (MP) has to be solved repeatedly. How to reduce computational burden in (MP) thereby is essential. In the proposed scheme, we solve the relaxed master problem but allow the integer constraints on several zk to be recovered simultaneously at each iteration. The proposed scheme is different from the work by Luedtke 28 because we rely on the CPLEX to assign the optimal values for several integer zk instead of branching them. The formula of proposed approach at (N + 1)th iteration is:

min

x,z,η

s.t.

qxT x K X

K 1 X ηk + K k=1

(RMP)

zk 6 K,

k=1 (n)

(n)

(n)

ηk > (αk )T x + βk + zk Mk , ∀k ∈ Ξ(x(n) ), ∀n ∈ {1, 2, . . . , N }, (n) (n) ˆ (n) , ∀k ∈ Ξ(x ˆ (n) ), ∀n ∈ {1, 2, . . . , N }, ηk > (α ˆ k )T x + βˆk + (1 − zk )M k (n)

(n)

0 > (γk )T x + θk + zk Hk , ∀k ∈ Γ(x(n) ), ∀n ∈ {1, 2, . . . , N }, (n) (n) ˆ k , ∀k ∈ Γ(x ˆ (n) ), ∀n ∈ {1, 2, . . . , N }, 0 > (ˆ γk )T x + θˆk + (1 − zk )H X X (n) xj − xj 6 |{i : xi = 1}| − 1, ∀n ∈ {1, 2, . . . , N }, (n) j∈{i:xi =1}

(6)

(n) j∈{i:xi =0}

x ∈ {0, 1}Dx ∩ Π, zk ∈ {0, 1}, ∀k ∈ ∆(N +1) , zk0 ∈ [0, 1], ∀k 0 ∈ {1, 2, . . . , K} \ ∆(N +1) ,

where ∆(N +1) is a scenario index set at (N + 1)th iteration, in which associated zk are binary. The main difference between (RMP) and (MP) lies in the relaxed variable zk0 , ∀k 0 ∈ {1, 2, . . . , K} \ ∆(N +1) . Due to the relaxation on binary variables zk0 , the solution of (RMP) is a lower bound of ACS Paragon Plus Environment

12

Page 13 of 38 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


(MP). Initially, we set ∆(1) = ∅, which means all binary variables are relaxed. Then, a heuristic strategy is employed to update ∆ at each iteration to continuously reduce the relaxation gap. (N )

(N )

Given the solution zk , ∀k ∈ {1, 2, . . . , K} \ ∆(N ) , if |zk

− 0.5| 0 α ˆ ki − αki 6 0

n o ⊆ x|xi ∈ {xi , xi }, ∀i and xi ∈ {0, 1}, xi ∈ {0, 1}. (n)

(n)

(n)

ˆ k − αk )T r˜ < Vk Proof. Owing to Lemma 1, it is sufficient to show that (α


13

(n)

− βk . Let us study


Page 14 of 38

the optimization problem for the scenario k in Lemma 1: (n)

Vk

(n) = min qˆyT yˆk + qˆsT sˆk − (αk )T x

(BM)

x,yˆk ,ˆ sk

ˆ +G ˆ k yˆk + C ˆ sˆk + D ˆ yˆk 6 0, s.t. Ax sk , sˆk ]. x ∈ {0, 1}Dx ∩ Π, yˆk ∈ [yˆk , yˆk ], sˆk ∈ [ˆ The problem (BM) is nonconvex due to the integer variable x. However, we can still make use of its dual formulation: (n) (n) ˆ (n) )T (Ax ˆ +G ˆ k yˆ + C ˆ sˆ + D ˆ y) ˆ Ok = min qˆyT yˆk + qˆsT sˆk − (αk )T x + (λ k x,yˆk ,ˆ sk

(BM1)

s.t. x ∈ {0, 1}Dx ∩ Π, yˆk ∈ [yˆk , yˆk ], sˆk ∈ [ˆ sk , sˆk ]. (n)

Here Ok

(n)

< Vk

ˆ (n) . The due to the weak duality and using the non-optimal dual multiplier λ k

objective function of (BM1) can be divided into two parts for x and {yˆk , sˆk }: (n) ˆ (n) )T (G ˆ k yˆk + C ˆ sˆk + D ˆ yˆk ) Ok = min qˆyT yˆk + qˆsT sˆk + (λ k yˆk ,ˆ sk

(n)

(BM1’)

(n)

ˆ k )T x + min − (αk )T x + (α x

s.t. x ∈ {0, 1}Dx ∩ Π, yˆk ∈ [yˆk , yˆk ], sˆk ∈ [ˆ sk , sˆk ]. The first minimization part of (BM1’) is optimized over {yˆk , sˆk }. Comparing it with the dual formulation of (SPRk ), there is

min y ˆk ∈[ˆ y k ,ˆ y k ],ˆ sk ∈[ˆ sk ,ˆ sk ]

(n)

(n)

(n) ˆ k yˆ + C ˆ )T (G ˆ sˆk + D ˆ yˆk ) = objSPR (x(n) ) − α ˆT = βˆk . qˆyT yˆk + qˆsT sˆk + (λ kx k k

(n)

(n)

(n)

ˆk − For the second minimization part of (BM1’), since (α ˆ k − αk )T r˜ is the minimal value of (α n o (n) αk )T x subject to x ∈ x|xi ∈ {xi , xi }, ∀i , which is a relaxation of the set x ∈ Π ∩ {0, 1}Dx , it (n) (n) (n) (n) (n) implies that βˆk + (α ˆ k − αk )T r˜ < Ok < Vk . (n)

The value of Mk

(n)

derived from Proposition 1 is always less than Vk ACS Paragon Plus Environment

14

− β (n) used in Lemma 1.

Page 15 of 38 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


Thus, the resulting optimality cut is weaker compared with the one in the reference. 30 However, it does not solve an additional optimization problem for each scenario and thereby the overall computational time is saved. Several types of cutting plane are developed by Liu et al. 30 based on (n)

the value of Vk

(n)

(n)

− βk . The Mk

derived in this proposition is also applicable in those cuts.

If a subproblem is infeasible for a given x(n) under normal or recovery mode, then k belongs ˆ (n) ). The objective value of its dual formula is unbounded and associated to set Γ(x(n) ) or Γ(x (n)

extreme ray can be obtained to determine γk

(n)

or γ ˆk . Luedtke proposes a new type of feasibility

cut while takes zk into account. 28 However, to develop such cut from a single infeasible subproblem, it needs solving many optimization problems to generate a pool of cutting planes. Then, the strongest cut is selected among them by finding the longest path in an acyclic graph. 33 This method is computationally expensive and thus may not be practical for a large-scale problem with many scenarios. Furthermore, based on our experience, the multi-cut method that uses a number of different extreme rays is more efficient to characterize the feasible set. 34 Thus, a simple, optimization-free approach of multi-cut generation for feasibility problem should be developed. In (n)

(3), γk

(n)

and θk

are obtained through the dual formulation of subproblems. According to A2:

ˆ x (ω) ⊇ Fx (ω), ∀ω ∈ Ω, we can set H ˆ k = 0 because the feasibility cut of recovery operations at F scenario k will not cut off any part of Fx (ωk ). Then it is only necessary to determine the value of Hk . To this end, the following proposition is presented: ˜ such that Proposition 2. Let us define vectors p˜ and t,    −xi if p˜i =   xi if

(n)

γki > 0 (n)

γki < 0

, t˜i =

   x if i

γˆki − γki > 0

  xi if

γˆki − γki < 0

(n0 )

(n)

(n0 )

(n)

n o where i is the index of vector, Π ∩ {0, 1}Dx ⊆ x|xi ∈ {xi , xi }, ∀i and xi ∈ {0, 1}, xi ∈ {0, 1}. If (n0 ) (n0 ) there exists a feasibility cut for recovery operations: 0 > (ˆ γk )T x + θˆk at previous iteration n0 , (n) (n) (n) (n0 ) (n) (n0 ) (n) makes (3) valid for the feasible then Hk = max − θk + (γk )T p, ˜ (ˆ γk − γk )T t˜ + θˆk − θk


15


Page 16 of 38

set of (Pf). (n)

Proof. If zk = 1 and Hk

(n)

(n)

= −θk + (γk )T p, ˜ then inequality (3) for scenario k does not impact (n)

(n)

the master problem because its right-hand side is always non-positive: (γk )T x + (γk )T p˜ = n o P P (n) (n) 0 0 , x ) + γ (x + x ) 6 0, ∀x ∈ x|x ∈ {x x }, ∀i , where h+ and h− are γ (x − 0 0 − + i i i i i i i ki ki i ∈h i∈h (n)

(n)

the index set for γki > 0 and γki < 0, respectively. 0 (n0 ) (n) (n0 ) (n0 ) ˜ θˆk(n ) −θk(n) , then we compare (ˆ γk )T x+θˆk with the right= (ˆ γk −γk )T t+ (n0 ) (n0 ) (n) (n) (n0 ) (n) (n0 ) (n) hand side of (3). There is (ˆ γk )T x + θˆk − (γk )T x + θk + (ˆ γk − γk )T t˜ + θˆk − θk = P P (n0 ) (n) (n0 ) (n) (n0 ) (n) ˜ = (ˆ γk − γk )T (x − t) γki − γki )(xki − xki ) + i0 ∈u− (ˆ γki0 − γki0 )(xki0 − xki0 ) > 0, i∈u+ (ˆ

(n)

If zk = 1 and Hk

where u+ and u− are the index set for γˆki − γki > 0 and γˆki − γki < 0, respectively. It implies (n0 ) (n) (n) (n) (n) (n0 ) (n0 ) (n) (n) (n0 ) (n) 0 > (ˆ γk )T x + θˆk > (γk )T x + θk + (ˆ γk − γk )T t˜ + θˆk − θk = (γk )T x + θk + Hk

for every admissible x. This further means that cutting plane (3) is dominated everywhere by an existing feasibility cut and thus redundant when zk = 1. (n)

Proposition 2 shows how to derive Hk by using single feasibility cut from the previous iteration (n)

n0 . The value of Hk

(negative) may increase if more feasibility cuts of (SPRk ) from previous

iterations are considered. This way will further tighten the cutting plane. At the end of this subsection, we present an efficient approach to derive tighter bounds: (n)

xi , xi , ∀i, which lead to a less conservative Mk , and thus reduce the computational time of (RMP). A simple scheme is to use the software of interval arithmetics, such as MC++. 35 Even though it is fairly fast, the resulting bounds are usually very loose. If Dx (dimension of x) is not very large, we recommend the optimization-based method because the time saving on (RMP) is


16

Page 17 of 38 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


more significant. The formulation of bound tightening at (N + 1)th iteration is shown below:

max / min xi x,z,η

(B)

x,z,η

s.t.

K X

zk 6 K,

k=1

qxT x

K 1 X ηk 6 UBD, + K k=1 (n)

(n)

(n)

(n)

(7) (n)

ηk > (αk )T x + βk + zk Mk , ∀k ∈ Ξ(x(n) ), ∀n ∈ {1, 2, . . . , N }, (n)

ˆ , ∀k ∈ Ξ(x ˆ (n) ), ∀n ∈ {1, 2, . . . , N }, ηk > (α ˆ k )T x + βˆk + (1 − zk )M k (n)

(n)

(n)

0 > (γk )T x + θk + zk Hk , ∀k ∈ Γ(x(n) ), ∀n ∈ {1, 2, . . . , N }, (n) (n) ˆ (n) , ∀k ∈ Γ(x ˆ (n) ), ∀n ∈ {1, 2, . . . , N }, 0 > (ˆ γk )T x + θk + (1 − zk )H k

x ∈ [0, 1]Dx ∩ Π, zk ∈ [0, 1], ∀k ∈ {1, 2, . . . , K},

where UBD is the best objective value of (Pf) found so far. The (B) is similar to (MP) but relaxes integer variables. The inequality (7) requires that the associated objective value (left-hand side of (7)) for any considered x should be less than the current optimal value: UBD. This inequality eliminates many feasible x of (Pf) but still preserves the global optimum. Hence the lower/upper bounds on x are tightened. In the minimization, if the solution is strictly positive: x∗i > 0, then xi can be fixed at 1. In the maximization, if the solution is less than 1: x∗i < 1, then xi can be fixed as 0.

The Algorithm of Improved Benders Decomposition In this subsection, we show the flowchart of the improved Benders Decomposition for a scenariobased two-stage stochastic program. Discussions about the implementation and comparisons with existing works are also presented.


17


Page 18 of 38

output: Optimal solution: x, z and objective function value Initialize: n = 1, x(1) , , τ , UBD, LB= −∞, ∆(1) = ∅ ; while |UBD − LB|/|UBD| 6 τ do ˆ (n) ) = ∅, Γ(x(n) ) = ∅, Γ(x ˆ (n) ) = ∅ ; o = 0; Ξ(x(n) ) = ∅, Ξ(x for k ← 1 to K do Substitute x(n) into (SPNk ) and solve it to obtain objSPNk (x(n) ); if (SPNk ) is feasible then Generate an optimality cut; o ← o + 1; Ξ(x(n) ) ← Ξ(x(n) ) ∪ {k}; else Generate a feasibility cut; Γ(x(n) ) ← Γ(x(n) ) ∪ {k}; end Substitute x(n) into (SPRk ) and solve it to obtain objSPRk (x(n) ); if (SPRk ) is feasible then ˆ (n) ) ← Ξ(x ˆ (n) ) ∪ {k}; Generate an optimality cut; Ξ(x else ˆ (n) ) ← Γ(x ˆ (n) ) ∪ {k}; Generate a feasibility cut; Γ(x Break; end end if o > K(1 − ) and P P (n) (n) obj (x ) + obj (x ) < UBD then qxT x + K1 (n) (n) ˆ SPN SPR k∈Ξ(x ) k∈Ξ(x ) k k P P (n) (n) UBD= qxT x + K1 obj (x ) + obj (x ) ; (n) (n) ˆ SPN SPR k∈Ξ(x ) k∈Ξ(x ) k k end Solve (B) and update xi , xi , ∀i; n ← n + 1; Add Balas cut; Solve (RMP) and generate x(n) ; Update LB and ∆(n) ; end Algorithm 1: Improved Benders decomposition In the flowchart, |UBD−LB|/|LB| is the formula of the relative gap. o is used to count the number of feasible scenarios under normal operations. If o is greater than K(1−), then it implies that x(n) ACS Paragon Plus Environment

18

Page 19 of 38 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


is a feasible solution of (Pf). Moreover, if the resulting objective value is lower than the current optimal solution, then UBD can be updated. LB represents the lower bound solution, obtained from solving (RMP). Given a solution x(n) from (RMP), this algorithm solves two linear programs for each scenario and then different dual multipliers will be obtained for optimality or feasibility cuts generation. The algorithm continuously solves subproblem (SPNk ) for all scenarios even if an infeasible scenario is found. If (SPRk ) is infeasible, then only the feasibility cut of the recovery operation for that scenario is attached to (RMP), and the scenario loop will be terminated. Based on these new cuts, a new x will be generated by solving (RMP). This procedure is repeated until the relative gap is small enough. Different from conventional Benders decomposition, since both z and x coexist in the cutting planes, the Balas cut should be used to exclude the binary values of x, which have been checked in previous iterations. Hence, even though the integer restriction on some zk in (RMP) is relaxed, the convergence still can be guaranteed. Here parameter τ is the threshold of the relative gap between lower and upper bounds. For the global optimization with a vast number of scenarios, it is practical to set τ as 0.1%. Compared with the Benders decomposition developed by Liu et al., 30 several improvements can be found in the proposed approach. First, the M (n) and H (n) are determined without using optimization. Considering that solving the MILP: (BM) for each scenario at each iteration could be very time-consuming, the proposed approach indeed saves the total computational time. Second, bound tightening is used to find x and x. Third, no branching in the proposed Benders decomposition scheme. Fourth, the Balas cut is adopted in the (MP). Fifth, the multi-cut generation is used in this paper whereas a single-cut scheme is adopted by Liu et al. 30

Validate the Probabilistic Feasibility A feasible solution of (P) requires x ∈ Fx (ω) with probability 1 − . Note that an optimal solution x∗ of (Pf) by using Algorithm 1 may not be feasible in unsampled scenarios. Thus, it is of importance to develop another algorithm to validate the probabilistic feasibility of solution x∗ in (P). Even though we may further generate more scenarios to verify such a solution ACS Paragon Plus Environment

19


Page 20 of 38

empirically, the required number of scenarios has not been rigorously determined for two-stage chance-constrained programs. Instead, a sampling-free method by combining linear decision-rule and chance-constrained program (CCP) is proposed in this section.

Revisit of Problem Formulation Among all constraints describing Fx (ω), the equalities should be processed first in order to develop a linear decision-rule. The feasible set of x for normal operation: Fx (ω) thus is rewritten as: n T T T Fx (ω) = x|∃{y, s},aT i x + ci s + (di + ω Bi )y 6 0, ∀i ∈ I, aT j x

+

cT j s

+

(dT j

o + ω Bj )y = 0, ∀j ∈ E, y ∈ [y, y], s ∈ [s, s] , T

(8)

where I ∩ E = ∅ and I ∪ E = {1, 2, . . . , Dr }. Design a Linear Decision-Rule The linear decision-rule assumes that some recourse variables can be represented by affine functions of uncertain parameters ω. 36 To obtain a linear formulation and enable the CCP later, we only let s to be dependent on ω. Then an assumption is made to handle equalities in (8): A3: |E| < Ds where |E| is the number of equalities in (8). Then the vector s is divided into two sub-vectors: Step 1: Choose an index subset S1 for vector s such that |S1 | = Ds − |E| and associated subvector is denoted by sS1 . The linear decision-rule is sS1 = U ω + e, where U and e are design parameters. This step is always feasible due to A3. Step 2: Define an index subset S2 such that S1 ∪ S2 = {1, 2, . . . , Ds } and S1 ∩ S2 = ∅. Here two sets S1 and S2 are mutually exclusive and collectively exhaustive. The vector s thus can be rearranged as [sS1 , sS2 ]T . Let AE , CE , DE + GE to denote the coefficient matrices of variables


20

Page 21 of 38 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


x, s and y in equalities, respectively. Then the equality constraints in (8) can be expressed as:

AE x + [CS1 , CS2 ]s + (DE + GE )y = 0,

(9)

where [CS1 , CS2 ] = CE . By assuming that: A4: CS2 is invertible. Then sS2 can be solved through the equality constraints (9): sS2 = −(CS2 )−1 (AE x + (DE + GE )y + CS1 sS1 ) = −(CS2 )−1 (AE x + (DE + GE )y + CS1 (U ω + e)) = −(CS2 )−1 AE x − (CS2 )−1 (DE + GE )y − (CS2 )−1 CS1 U ω − (CS2 )−1 CS1 e

(10)

Given a stage I solution x∗ , the s can be expressed as:

s = Tω ω + Ty y + Tc

(11)

where   Tω = 

 U −(CS2 )−1 CS1 U



   , Ty = 

 0 −(CS2 )−1 (DE + GE )



   , Tc = 

 e −CS−1 (AE x∗ + CS1 e) 2

 

Note that s is affine with respect to uncertain parameter w.

Chance-Constrained Programming Let us substitute (11) into inequalities and variable bounds in (8):

∗ T T T aT i x + (ω Bi + di )y + ci (Tω ω + Ty y + Tc ) 6 0, ∀i ∈ I,

y ∈ [y, y], Tω ω + Ty y + Tc ∈ [s, s].


21

(Ch1)


Page 22 of 38

The left-hand side of constraints (Ch1) are affine functions of ω with decision variables y, U and e. To meet the probabilistic constraint P{x∗ ∈ Fx (ω)} > 1 − , (Ch1) should be satisfied with probability 1 − : n ∗ T T T P aT i x + (ω Bi + di )y + ci (Tω ω + Ty y + Tc ) 6 0, ∀i ∈ I, o Tω,j ω + Ty,j y + Tc,j ∈ [sj , sj ], ∀j > 1 − .

(Ch2)

A deterministic equivalence for joint probabilistic constraints (Ch2) has not been developed. However, the theory of Bonferroni’s inequality can be used to build a conservative approximation for (Ch2):

∗ T T T 0 P{aT i x + (ω Bi + di )y + ci (Tω ω + Ty y + Tc ) 6 0} > 1 − i , ∀i ∈ I,

(12)

P{Tω,j ω + Ty,j y + Tc,j 6 sj } > 1 − 00j , ∀j,

(13)

P{Tω,j ω + Ty,j y + Tc,j > sj } > 1 − 00j , ∀j,

(14)

|I| X

0i + 2

i=1

Ds X

00j = ,

(15)

j=1

where 0i , ∀i ∈ I are risk levels; 00j is specified as 0.01% to require s bounded with high probability 99.99%. For more advanced linear decision-rule to handle variable bounds, readers are referred to Chen et al. 31 Given A1, chance constraints (12)-(14) can be equivalently converted to some deterministic inequalities. Namely, if ω ∈ N (µ, Σ), (12)-(14) are equivalent to:

p ∗ T T T T −1 0 (Bi y + TωT ci )T Σ(Bi y + TωT ci ) 6 0, aT x + c (T y + T ) + d y + µ (B y + T c ) + Φ (1 − ) y c i i i i i ω i ∀i ∈ I, Ty,j y + Tc,j + µT Tω,j Ty,j y + Tc,j + µT Tω,j

(16) q T + Φ−1 (1 − 00j ) Tω,j ΣTω,j 6 sj , ∀j, q T + Φ−1 (1 − 00j ) Tω,j ΣTω,j > sj , ∀j,

(17) (18)

where Φ−1 : R → R represents the standard inverse Normal cumulative distribution function. When 0i < 0.5, Φ−1 (1 − 0i ) is convex. Recently, the author has published an algorithm 37 for global ACS Paragon Plus Environment

22

Page 23 of 38 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


optimization of chance-constrained programs by using the outer approximation and McCormick relaxation. 39 In the current problem, given x∗ , the same algorithm can be adopted with minor revisions to find a feasible solution of y, U , e, and 0 , satisfying (16)-(18). The outer approximation is used to obtain a relaxation for Eq. (16): 38

∗ T T T T aT i x + ci (Ty y + Tc ) + di y + µ (Bi y + Tω ci ) + vi

p (Bi y + TωT ci )T Σ(Bi y + TωT ci ) 6 0 (19)

vi > κl,i + ρl,i 0i , ∀l ∈ {1, 2, . . . , L}, ∀i ∈ I,

(20)

where L is the number of cutting planes for Φ−1 (1 − 0i ). Eq. (20) represents the tangent of Φ−1 (1 − 0i ). Parameters κl,i and ρl,i , ∀l ∈ {1, 2, . . . , L} are defined at the point ∗l ∈ (0, 0.5) where cutting plane is built: dΦ−1 (1 − 0i ) −1 ρl,i = = , 0 −1 di φ(Φ (1 − ∗l )) 0 =∗ i

l

κl,i = Φ−1 (1 − ∗l ) − ρl,i ∗l ,

where φ is the density function of standard normal distribution. As more cutting planes are generated, the outer approximation will converge to the original function Φ−1 (1 − 0i ). A graphical illustration is shown in Fig. 2. The inequality (19) is still difficult to deal with due to its non-convexity. Some intermediate variables: wi,t = vi yt , Ri,j,h = vi Tω,j,h , ∀i ∈ I, ∀t ∈ {1, 2, . . . , Dy }, ∀j ∈ {1, 2, . . . , Ds } and ∀h ∈ {1, 2, . . . , Dω } have to be introduced to convexify bilinear terms. 39 Then, the resulting convex relaxation problem is:


23


Page 24 of 38

3.5 -1 (1- ) Cutting plane1 Cutting plane2 Cutting plane3 Cutting plane4

3

2.5

2

1.5 vi

l,i

+ l,i i

1 0

0.02

0.04

0.06

0.08

0.1

0.12

Figure 2: The outer approximation of Φ−1 (1 − i ). As i increases, cutting planes can approximate the inverse cumulative function more accurately.

min

0 ,v,w,R,y,Tω ,Tc

y1

(RC) 00j )

q

Ty,j y + Tc,j + µT Tω,j + Φ−1 (1 − 00j )

q

T

−1

s.t. Ty,j y + Tc,j + µ Tω,j + Φ (1 −

T Tω,j ΣTω,j 6 sj , ∀j ∈ {1, 2, . . . , Ds },

T Tω,j ΣTω,j > sj , ∀j ∈ {1, 2, . . . , Ds }, p ∗ T T T T aT x + c (T y + T ) + d y + µ (B y + T c ) + (Bi w + RT ci )T Σ(Bi w + RT ci ) 6 0 y c i i i i ω i

vi > κl,i + ρl,i 0i , ∀l ∈ {1, 2, . . . , L}, ∀i ∈ I, wi,t > vi y t + v i yt − v i y t , wi,t > vi y t + v i yt − v i y t ,

(21)

wi,t 6 vi y t + v i yt − v i y t , wi,t 6 vi y t + v i yt − v i y t ,

(22)

Ri,j,h > vi T ω,j,h + v i Tω,j,h − v i T ω,j,h , Ri,j,h > vi T ω,j,h + v i Tω,j,h − v i T ω,j,h ,

(23)

Ri,j,h 6 vi T ω,j,h + v i Tω,j,h − v i T ω,j,h , Ri,j,h 6 vi T ω,j,h + v i Tω,j,h − v i T ω,j,h ,

(24)

∀i ∈ I, ∀t ∈ {1, 2, . . . , Dy }, ∀j ∈ {1, 2, . . . , Ds }, ∀h ∈ {1, 2, . . . , Dω }, |I| X

0i

+2

Ds X

i=1

00j = ,

j=1

y ∈ [y, y], ACS Paragon Plus Environment

24

Page 25 of 38 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


where T ω,j,h and T ω,j,h are lower and upper bounds on Tω,j,h , respectively; v i and v i are lower and upper bounds on vi , respectively. Eqs.(21)-(24) are the McCormick relaxation. 39 Any decision variables can be chosen as the objective function in (RC) since we only need a feasible solution. Through the outer approximation and McCormick relaxation, (RC) is convex and can be solved by the second-order cone program. Its solution may not be feasible to (Ch2), but we can always substitute the resulting i , denoted by 0RC,i , ∀i ∈ I, into Eq. (16) and solve the following optimization problem:

min y1

(C)

y,Tω ,Tc

∗ T T T T s.t. aT i x + ci (Ty y + Tc ) + di y + µ (Bi y + Tω ci ) p + Φ−1 (1 − 0RC,i ) (Bi y + TωT ci )T Σ(Bi y + TωT ci ) 6 0, ∀i ∈ I, q T −1 00 T Ty,j y + Tc,j + µ Tω,j + Φ (1 − j ) Tω,j ΣTω,j 6 sj , ∀j ∈ {1, 2, . . . , Ds }, q T ΣTω,j > sj , ∀j ∈ {1, 2, . . . , Ds }, Ty,j y + Tc,j + µT Tω,j + Φ−1 (1 − 00j ) Tω,j

y ∈ [y, y].

Note that (C) is obtained from Eqs. (16)-(18) by fixing 0i as 0RC,i , ∀i ∈ I, without any convex relaxation. It is still a second-order cone program and also a conservative approximation of (Ch2). If it is feasible, then it means that the x∗ derived from (Pf) meets the probabilistic constraints in (P).

The Algorithm for Feasibility Validation The linear decision-rule enables us to validate the probabilistic feasibility of a stage I solution x∗ through chance-constrained programs (RC) and (C). Given 0RC from the solution of (RC), if (C) is infeasible, it implies that the relaxation gap between (RC) and (C) is still substantial. Therefore, we develop a branch-and-bound (B&B) method to continuously tighten the relaxation (RC) and search for 0RC . By branching v to reduce its upper or lower bound, the McCormick relaxation becomes tighter. In addition, the tangent lines are added adaptively to make (RC) more accurate. The procedure of this algorithm is similar to Yang et al., 37 but it can terminate when a feasible ACS Paragon Plus Environment

25


Page 26 of 38

solution of (C) is found. output: 0RC or problem (C) is infeasible for the given x∗ (t)

(t) (t) Initialize: l = 1, t = 1, [v (t) , v (t) ], [y (t) , y (t) ], [T (t) ω , T ω ], , ST = {[v , v ]} ;

while ST = 6 ∅ do Load a node {[v (t) , v (t) ]} from ST ; ST ← ST \ {[v (t) , v (t) ]} ; Solve (RC) ; if (RC) is feasible then Subsititute 0RC into (C) and solve it; if (C) is feasible then Return 0RC ; Break; else Calculate κl and ρl at = 0RC and add this cutting plane to (RC); L ← L + 1; end Let ST = ST ∪ {[v (t) , vRC ], [vRC , v (t) ]} ; end t←t+1 ; end Algorithm 2: B&B and Chance-constrained Programming where ST is a stack to store the lower and upper bounds on v; 0RC and vRC are solutions of (RC). Note that if (RC) is infeasible on the domain [v (t) , v (t) ], then the associated node can be abandoned. As variables in the vector v are gradually branched, the difference between (RC) and (C) will become smaller. Therefore, through continuously solving the (RC) on different ranges of v, we can either find a feasible solution of (C) or discard the entire range of v to terminate the algorithm. Finally, we need to remark that algorithm 2 only provides a sufficient condition to validate the probabilistic feasibility of x(∗) . The linear decision-rule and Bonferroni’s inequality result in a conservative approximation of joint chance constraints. Thus, if the proposed method cannot find a feasible solution in (C) for a given x∗ , it does not imply that x∗ is infeasible to the problem (P).


26

Page 27 of 38 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


Case Study In this section, the improved Benders decomposition and chance-constrained programming are implemented by using GAMS 24.8.3 with CPLEX 12.7.0. The hardware platform is an Intel Core i5-3230M CPU, 2.60GHz with 4GB memory. The profitability of a refinery is maximized by choosing crude oil and operations through optimization. The superstructure of this refinery is shown in Fig. 3 and its LP model is from Favennec 40 with modifications in parameters, specifications and recovery operations. Three types of crude oil can be purchased and processed with maximum capacity 1800000 bbl for each crude. The minimal procurement is 50000 bbl if such crude is selected. Thus, following constraints and stage I variables are considered: 41 9 X

2j−1 gCrudei ,j Q > 50000gpurchase , gCrudei ,j 6 gpurchasei , ∀i ∈ {1, 2, 3},

j=1

where gCrudei ,j ∈ {0, 1} is a binary variable to count purchase quantity, and gpurchasei ∈ {0, 1} represents the procurement decision on Crudei . Q is specified as 5000 bbl as one lot. Similar to the work of Yang and Barton, 41 the purchase quantity is an integer and substituted by binary variables to describe the real crude oil trade. The objective is to decide how many lots have to be purchased such that the profit is maximized and all quality specifications, as well as demands shown in the Support Information Table S3, are satisfied with 95% chance or more. The objective function used in the Benders decomposition is shown in Eq. (25) K 1 X (mCr,CGO,Mogas,k + mCr,CGO,AGO,k ) · 3 + mR,HN,Sev95,k · 2.7 + mR,HN,Sev100,k · 3.2 obj = − K k=1 X + mISO,LN,k · 0.6 + mDes,GOi ,k · Des-CostGOi + mDes,CGO,k · 1.4 i

+ mDIESEL,Import,k · 501.97 − mLG,Export,k · 203.1 − mLN,Export,k · 528 − mES95,Export,k · 960 − mPG98,Export,k · 1081.4 − mJet,Export,k · 428 − mDIESEL,Export,k · 497 − mHF,Export,k · 319 × 103 −

9 XX i

2j−1 gCrudei ,j QPriceCi ,

(25)

j=1


27


Page 28 of 38

where mass flow m is the stage II variable. Its subscript k is the scenario index. K is the total number of scenarios. The crude distillation unit (CDU) is used to separate different components in crude oil. The yields are shown in Support Information Table S1. After this processing, intermediate flows are mixed or further refined to produce gasoline 95, gasoline 98, jet fuel, diesel and heavy fuel oil (HF). Different from the original model, 40 there is no importation for any intermediate products. The ultra-low sulfur diesel is imported and mixed with the produced diesel from the refinery as a recovery operation to meet the specification of the final product. This operation is not favored even though it can be economically profitable, because such high-quality diesel may not be always available in the market. Thus assumption A2 satisfies and chance constraints should be adopted to limit the application of recovery operations. To this end, binary variable zk and extra constraints are introduced:

0 6 mDIESEL,Import,k 6 zk mDIESEL,Import,k , X zk 6 K, k

where mDIESEL,Import,k is the upper bound of importation. If zk = 1, then it implies that importation is allowed. If zk = 0, then the normal operation is feasible and no importation. Parameters, prices, specifications, demands, and uncertainties are listed in the Support Information, Tables S1-S4. Here three types of uncertain parameters in the model are considered, including sulfur content in gas oil (GO), viscosity in the vacuum residue (VR), the fraction of CN in cracker outcome. They are independent and normally distributed. The details of the original LP-formulation are shown in the Support Information. Scenarios are generated through the random sampling on these uncertain parameters. The problem (Pf) is constructed by 300, 1000, 2000, 3000 scenarios, respectively. The resulting MILPs are solved by the proposed improved Benders decomposition method. Since we specify the Big-M value for optimality cuts only using (B), which is different from Liu et al., 30 it is necessary to compare the computational time of these two schemes. For a fair comparison, both schemes solve (RMP) to avoid implementation issues on the branch-and-cut 28 . For feasibility cuts, we only ACS Paragon Plus Environment

28

Page 29 of 38 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


LG REFINERY FUEL RG

RG LG

Crude1~3

LN

RG

HN

LG R95 REFORMER R100

B

LN LG

PG98

LN ISOMERIZATION B

GASOLINE ISO R95 R100 CN

ES95

CDU F1 OR F2

KE

Importation

RG

GOi

VGO B

RG LG CN CRACKER CGO

JF

DES-GOi B

DESULFURIZATION

DES-CGO

DIESEL

HF

VRi

DIESEL

HF

Figure 3: The refinery flowchart. 40 specify the Big-M value according to Proposition 3 because it is very time-consuming to use the optimization scheme 28 to determine such a value. Moreover, we also directly use CPLEX to solve resulting MILPs as a benchmark solution. All algorithms terminate when the relative gap reaches 0.1%. The computational performances, including the clock and CPU time for different methods and cases, are shown in Table 1. For Benders decomposition, we further show the computational time for RMP, SP, and BM separately. Several discussions are provided in order. First, it is not difficult to see that CPLEX outperforms Benders decomposition when the number of scenarios is small. However, as the quantity of scenarios increases to 2000, the Benders decomposition becomes more efficient whereas the CPLEX cannot obtain a satisfactory solution within 4 hours. Second, since the proposed method relies on the bound tightening (B), Propositions 2 and 3 to calculate the Big-M value, thus it takes less time to obtain an optimal solution even though its Big-M value is weaker than the one from solving (BM). Third, the proposed method relaxes the integer constraint for certain zk when the relaxed master problem (RMP) is solved. If the Big-M parameter is too large such that the values of solution zk in (RMP) are far away from 0 and 1, then the integer constraint for such zk will be recovered to compensate the impact of large Big-M values. Fourth, ACS Paragon Plus Environment

29


Page 30 of 38

the clock time is important because loading different models usually is much more time-consuming than the computations, which is the major bottleneck of Benders decomposition implemented in GAMS. The proposed improved Benders decomposition does not load and solve (BM) for each scenario. Therefore, the total clock time is reduced and this advantage becomes more significant as more scenarios are considered in the formulation. Table 1: The CPU/Clock time (minute) of different methods for 300, 1000, 2000, 3000 scenarios. Here the method proposed by Liu et al. 30 solves (BM). The proposed method solves (B) for bound tightening. SP represents both (SPN) and (SPR). Scenarios in (Pf) CPLEX Method in Liu et al. 30 RMP SP BM Total time Proposed Method RMP SP B Total time

3000 > 240

2000 > 240

1000 15.77/15.78

300 0.23/0.24

0.95/1.14 2.31/2.58 0.67/0.90 0.08/0.12 7.29/120.41 5.14/77.80 4.76/79.98 0.68/ 7.16 2.64/46.38 1.84/32.14 2.01/32.95 0.22/2.84 10.88/167.93 9.29/112.52 7.44/104.83 0.98/10.12 1.91/2.10 6.42/87.31 0.14/0.33 8.47/89.74

1.68/2.01 5.02/77.31 0.00/0.02 6.70/79.34

0.60/0.80 4.64/79.84 0.06/0.15 5.30/80.79

0.08/0.12 0.73/7.53 0.00/0.01 0.81/7.66

Table 2: The optimal crude oil procurement (KT) and monthly profit (million$) calculated by solving (Pf) with 3000/2000/1000/300 scenarios, respectively. Scenarios in (Pf) 3000 2000 1000 300 Crude1 0 0 6.71 0 Crude2 205.993 208.006 198.612 220.083 Crude3 239.494 239.494 239.494 226.189 Profit 87.281 87.650 87.346 87.038

In the next step, the solutions in Table 2 are verified by the proposed linear decision-rule and chance-constrained programming. The set S1 , in which the variables are dependent on the uncertain parameters, is determined first. This refinery model satisfies assumption A3 and Ds − |E| = 3. Variables mES95,CN (CN to ES95), mDIESEL,CGO (CGO to Diesel) and mDES,CGO (CGO to Desulfurization unit) are selected to form set S1 . Since the mass flow of CN is only dependent on the cracker outcome fraction rather than the GO sulfur and VR viscosity, thus mES95,CN is an affine function of cracker outcome fraction only. mDIESEL,CGO and mDES,CGO are dependent on the cracker ACS Paragon Plus Environment

30

Page 31 of 38 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


outcome fraction and sulfurs in GO. Note that mHF,CGO is related to mDIESEL,CGO and mDES,CGO due to the mass balance. Moreover, mHF,CGO and mHF,VRi are also correlated because both of them determine the viscosity of the product HF. It implies that uncertainties on the viscosity of VR indirectly impact the solution of mHF,CGO , mDIESEL,CGO , and mDES,CGO . As a result, mDIESEL,CGO and mDES,CGO are affine functions of all uncertain parameters. Here m represents the mass flow. The first subscript is the product and the second subscript is the feedstock to make that product. Given the procurement, shown in Table 2, the chance-constrained program is solved to find a feasible solution for stage II by setting = 5% in (15). The probabilistic satisfactions of each quality and demand constraints are calculated and listed in Table 3. Such results are dependent on the parameters of algorithm 2. However, since only a feasible solution is needed, the conclusion will not change. One can see that the stage I solutions based on 300, 2000, 3000 sampled scenarios are able to meet quality and demand constraints with high probability whereas the solution based on 1000 scenarios cannot be verified by the proposed method. The main difference among all stage I solutions in Table 2 is the procurement of Crude1 . Without purchasing Crude1 , the total uncertain parameters used in the linear decision-rule is 6. On the contrary, with Crude1 , the number of uncertain parameters increases to 8, which either renders the problem infeasible or disables the linear decision-rule for stage II. We also show the probabilistic constraints satisfaction through the sampling method. This empirical validation indicates that the solution from 1000 scenarios indeed violates the chance constraint.

Conclusion In this paper, we first propose an improved Benders decomposition to solve the scenario-based two-stage stochastic program efficiently. Second, we combine the linear decision-rule and the chance-constrained program to validate the probabilistic feasibility of the solution derived from the improved Benders decomposition. This scheme inherits the advantages from both the scenariobased method and chance-constrained programming to find a near-optimal and reliable solution. The plan of crude oil procurement and plant operations based on a simplified refinery model is optimized to show the effectiveness of proposed algorithms. ACS Paragon Plus Environment

31


Page 32 of 38

Table 3: The probabilistic constraints satisfaction for the optimal solutions derived from 3000, 2000, 1000, 300 scenarios, respectively. The upper/lower bounds on variables and maximum capacity have to be satisfied with 99.99%. The probability of quality and demand satisfactions are determined by the proposed linear decision-rule and chance-constrained programming. Only the first feasible solution found by algorithm 2 is reported. Empirical validation is based on newly generated 10000 scenarios. Scenarios in (Pf) Constraints Butane 95 Butane 98 Gasoline-98 RVP (min) Gasoline-98 RVP (max) Gasoline-95 RVP (min) Gasoline-95 RVP (max) Gasoline-98 RON Gasoline-95 RON Gasoline-98 Sulfur (max, ppm) Gasoline-95 Sulfur (max, ppm) Gasoline-95 Sensitivity Gasoline-98 Sensitivity Diesel Sulfur Heavy fuel oil Viscosity (min) Heavy fuel oil Viscosity (max) Refinery Fuel (max) Gasoline98 (demand) Diesel (demand) HFO (demand) Empirical Validation

3000 99.974% 99.970% 99.970% 99.970% 99.974% 99.974% 99.970% 99.974% 99.963% 99.971% 99.974% 99.963% 99.965% 95.600% 99.965% 99.971% 99.967% 99.967% 99.965% 96.04%

2000 1000 Probabilistic Constraints Satisfaction 99.990% 99.990% 99.990% 99.990% 99.990% 99.990% 99.647% 99.990% 99.990% 99.990% 99.990% 99.990% 99.990% 95.584% 99.990% 99.990% 99.990% 99.990% 99.990% 96.06% 93.32%


32

300 99.969% 99.965% 99.965% 99.965% 99.969% 99.969% 99.965% 99.969% 99.955% 99.965% 99.969% 99.965% 99.959% 95.696% 99.959% 99.959% 99.966% 99.962% 99.959% 96.15%

Page 33 of 38 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


Supporting Information Crude yields and price, Price of the products, Specifications and demands, Mean/Variance of uncertain model parameters, process model. The Supporting information is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgment Author is grateful to the financial support of Faculty Startup Fund from California State University, Long Beach.

References (1) Li, X., & Barton, P. I. (2015). Optimal design and operation of energy systems under uncertainty. Journal of Process Control, 30, 1-9. (2) You, F., Wassick, J. M., & Grossmann, I. E. (2009). Risk management for a global supply chain planning under uncertainty: models and algorithms. AIChE Journal, 55, 931-946. (3) Grossmann, I. E., Halemane, K. P., & Swaney, R. E. (1983). Optimization strategies for flexible chemical processes. Computers & Chemical Engineering, 7, 439-462. (4) Grossmann, I. E., & Floudas, C. A. (1987). Active constraint strategy for flexibility analysis in chemical processes. Computers & Chemical Engineering, 11, 675-693. (5) Bansal, V., Perkins, J. D., & Pistikopoulos, E. N. (1998). Flexibility analysis and design of dynamic processes with stochastic parameters. Computers & Chemical Engineering, 22, 817-820. (6) Rooney, W. C., & Biegler, L. T. (2001). Design for model parameter uncertainty using nonlinear confidence regions. AIChE Journal, 47, 1794-1804.


33


Page 34 of 38

(7) Zhang, Y., Monder, D., & Forbes, J. F. (2002). Real-time optimization under parametric uncertainty: a probability constrained approach. Journal of Process Control, 12, 373-389. (8) Arellano-Garcia, H., & Wozny, G. (2009). Chance constrained optimization of process systems under uncertainty: I. Strict monotonicity. Computers & Chemical Engineering, 33, 1568-1583. (9) Al-Qahtani, K., & Elkamel, A. (2011). Robust planning of multisite refinery networks: Optimization under uncertainty. Computers & Chemical Engineering, 34, 985-995. (10) Sahinidis, N. V. (2004). Optimization under uncertainty: state-of-the-art and opportunities. Computers & Chemical Engineering, 28, 971-983. (11) Soyster, A. L. (1973). Convex programming with set-inclusive constraints and applications to inexact linear programming, Operational Research, 21, 1154-1157. (12) Bertsimas, D. & Thiele, A. (2006). Robust and data-driven optimization: modern decision making under uncertainty. INFORMS Tutorials in Operations Research: Models, Methods and Applications for Innovative Decision Making. (13) Bertsimas, D., & Sim, M. (2004). The price of robustness. Operations Research, 52, 35-53. (14) Li, Z., Ding, R., & Floudas, C. A. (2011). A comparative theoretical and computational study on robust counterpart optimization: I. Robust linear optimization and robust mixed integer linear optimization. Industrial& Engineering Chemistry Research, 50, 10567-10603. (15) Li, Z., Tang, Q., & Floudas, C. A. (2012). A comparative theoretical and computational study on robust counterpart optimization: II. Probabilistic guarantees on constraint satisfaction. Industrial & Engineering Chemistry Research, 51, 6769-6788. (16) Li, Z., & Floudas, C. A. (2014). A comparative theoretical and computational study on robust counterpart optimization: III. Improving the quality of robust solution. Industrial & Engineering Chemistry Research, 53, 13112-13124. (17) Li, Z., & Li, Z. (2015). Optimal robust optimization approximation for chance constrained optimization problem. Computers & Chemical Engineering, 74, 89-99. ACS Paragon Plus Environment

34

Page 35 of 38 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) Yuan, Y., Li, Z., & Huang, B. (2017). Robust optimization approximation for joint chance constrained optimization problem. Journal of Global Optimization, 67, 805-827. (19) Ben-Tal, A., Ghaoui, L. E., & Nemirovski, A. Robust optimization. Princeton, New Jersey: Princeton University Press, 2006. (20) Prékoba, A. Stochastic programming. Netherlands: Kluwer Academic Publishers, 1995. (21) Calafiore, G. C., Campi, M. C. (2006). The scenario approach to robust control design. IEEE Transaction on Automatic Control, 51, 742-753. (22) Nemirovski, A., & Shapiro, A. (2006). Convex approximation of chance constrained programs. SIAM Journal on Optimization, 17, 969-996. (23) Karuppiah, R., & Grossmann, I. E. (2008). A Lagrangean based branch-and-cut algorithm for global optimization of nonconvex mixed-integer nonlinear programs with decomposable structures. Journal of Global Optimization, 41, 163-186. (24) Mouret, S., Grossmann, I. E., & Pestiaux, P. (2011). A new Lagrangian decomposition approach applied to the integration of refinery planning and crude-oil scheduling. Computers & Chemical Engineering, 35, 2750-2766. (25) Shah, N. K., & Ierapetritou, M. G. (2015). Lagrangian decomposition approach to scheduling large-scale refinery operations. Computers & Chemical Engineering, 79, 1-29. (26) Benders, J. F. (1962). Partitioning procedures for solving mixed-variables programming problems, Numerische Mathematik, 4, 238-252. (27) Geoffrion, A. M. (1972). Generalized Benders decomposition. Journal of Optimization Theory and Applications, 10, 237-260. (28) Luedtke, J. (2014). A branch-and-cut decomposition algorithm for solving chance-constrained mathematical programs with finite support. Mathematical Programming, 146, 219-244.


35


Page 36 of 38

(29) Mitra, S., Garcia-Herreros, P., & Grossmann, I. E. (2016). A cross-decomposition scheme with integrated primaldual multi-cuts for two-stage stochastic programming investment planning problems. Mathematical Programming, 157, 95-119. (30) Liu, X., K¨ uçu ¨kyavuz, S., & Luedtke, J. (2016). Decomposition algorithms for two-stage chance-constrained programs. Mathematical Programming, 157, 219-243. (31) Chen, X., Sim, M., Sun, P., & Zhang, J. W. (2008). A linear decision-based approximation approach to stochastic programming. Operations Research, 56, 344-357. (32) Balas, E., & Jeroslow, R. (1972). Canonical cuts on the unit hypercube. SIAM Journal of Applied Mathematics, 23, 61-69. (33) G¨ unl¨ uk, O., & Pochet, Y. (2001). Mixing mixed-integer inequalities. Mathematical Programming, 90, 429-457. (34) Yang, Y., & Lee, J. M. (2012). A tighter cut generation strategy for acceleration of Benders decomposition. Computers & Chemical Engineering, 44, 84-93. (35) Chachuat, B. http://www.imperial.ac.uk/people/b.chachuat/research.html. (36) Ben-Tal, A., Goryashko, A., Guslitzer, E., & Nemirovski, A. (2004). Adjustable robust solutions of uncertain linear programs. Mathematical Programming, 99, 351-376. (37) Yang, Y., Vayanos, P., & Barton, P. I. (2017). Chance-constrained optimization for refinery blend planning under uncertainty. Ind. Eng. Chem. Res., 56, 12139-12150. (38) Cheng, J., Gicquel, C., & Lisser, A. (2012). A second-order cone programming approximation to joint chance-constrained linear programs. Lecture Notes in Computer Science, 7422, 71-80. (39) McCormick, G. P. (1976). Computation of global solutions to factorable nonconvex programs: Part I convex underestimating problems. Mathematical Programming, 10, 147-175. (40) Favennec, J. P. Refinery operation and management. Paris: Editions TECHNIP, 2001.


36

Page 37 of 38 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


(41) Yang, Y., & Barton, P. I. (2016). Integrated crude selection and refinery optimization under uncertainty. AIChE Journal, 62, 1038-1053.


37


Master Problem x

Scenario 1

Scenario 2

x

x

Scenario 3

Normal Operations

x

x

Scenario K-1

Scenario K

Recovery Operations Number

Improved Benders Decomposition and Feasibility Validation for Two

Recommend Documents