Chance-Constrained Optimization for Refinery Blend Planning under

Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX-XXX

pubs.acs.org/IECR

Chance-Constrained Optimization for Refinery Blend Planning under Uncertainty Yu Yang,†,‡ Phebe Vayanos,†,¶ and Paul I. Barton*,† †

Process Systems Engineering Laboratory, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States Department of Chemical Engineering, California State University Long Beach, Long Beach, California 90840, United States ¶ Epstein Department of Industrial & Systems Engineering, University of Southern California, Los Angeles, California 90089, United States ‡

S Supporting Information *

ABSTRACT: A near-global optimization approach is proposed to design blending recipes for refinery products under uncertainties. In the refining industry, the most valuable products, such as gasoline and diesel, are produced by blending several intermediate feedstocks to maximize profit and ensure that all qualities are on specification. Because of the presence of property uncertainties, optimal blending recipes under linear mixing laws should be designed by solving a linear program with joint chance constraints. However, joint chance-constrained programming is generally intractable even with Gaussian distributions, and thereby, it is usually converted to an individual chanceconstrained (ICC) program to achieve a conservative approximation. To reduce this conservatism, we find a global optimal solution for the ICC program. In case studies, a multiproduct blending problem with 12 chance constraints and a crude oil procurement example with 14 chance constraints are studied to test and compare the proposed scheme with state-of-the-art optimization software to demonstrate its superior performance in terms of computational time.

■

INTRODUCTION Maximizing profit while meeting product quality specifications is always the primary goal of the refining industry in a highly competitive global market with increasingly stringent environmental regulations. For a refinery, the liquid products, such as gasoline and diesel, are the most important profit generators, usually yielding 60−70% of revenues,1 which presents opportunities for process improvement. Generally, each final product is produced by mixing different feedstocks according to certain recipes in a well-developed blender system. Those feedstocks usually include intermediate products such as light and heavy naphtha, catalytic cracked spirit (CN), light cycle gas oil (CGO), and so on. The blending recipe is crucial to determine the quality, cost, and sale price of the produced gasoline. Thus, the key task is to design blending recipes that minimize operational cost while keeping all final products on specification. In this context, a major challenge is that blendstock qualities are often unknown at the time of blend planning. In this study, we are concerned with the effects of uncertainties on the quality of blending products and developing a robust framework for optimal blend planning under uncertainty. This framework is highly practically relevant because blend plans are often made by oil and gas companies before crude oil qualities can be measured to inform © XXXX American Chemical Society

operational decisions (e.g., truck capacities, shipping schedules, etc.). Even after the blendstocks arrive at the refinery, and although advanced calibration equipment and online analyzers are available to measure the properties of crude oil and intermediate products to a certain extent, it is still very difficult to measure accurate properties in real time. In addition, noting that operating conditions and other external factors, including temperature, pressure, and catalysts in the distillation, cracker, and reformer change or drift even in the presence of process controllers, and so forth, the real properties of each feedstock fluctuate and are difficult to track in a timely fashion. As a result, without considering these uncertainties explicitly, blending recipes derived from an optimization formulation based on nominal parameter values can be infeasible (off specification) in reality, and the resulting economic loss could be large. In practice, and to avoid being off specification, refineries can add an offset to the target qualities, thereby requiring that the blended products meet more stringent target qualities or specifications. If the offset values are chosen appropriately, this approach can yield products that are on Received: Revised: Accepted: Published: A

June 13, 2017 September 17, 2017 October 4, 2017 October 4, 2017 DOI: 10.1021/acs.iecr.7b02434 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

difficult to solve even when the uncertain parameters are Gaussian distributed. For this reason, joint chance constraints are typically decomposed into individuals, which are equivalent to tractable second-order cone constraints when the uncertain parameters are Gaussian. The Bonferroni inequality then guarantees a feasible (yet often suboptimal) solution to the joint chance-constrained program that is efficiently computable.7 Nemirovski and Shapiro23 propose an approximate solution based upon the Bernstein approximation. Calafiore and Ghaoui26 extend the aforementioned result from the Gaussian to radial distribution case. In the work of Li et al.,27 the joint chance constraints are handled by relaxing the problem into a nonlinear program, and the major challenge is the calculation of the probability and its derivatives to satisfy inequality constraints. Two-stage optimization problems with chance constraints and Gaussian-distributed uncertainties for chemical process design are also considered.12,13 More recently, Zymler et al.28 derived a tractable formulation for CCP by assuming that only first- and second-order moments are given. Yang and Xu29 further considered nonlinear constraints in CCP. Both methods deal with joint chance constraints and introduce auxiliary parameters to control the accuracy of approximation. However, how to optimize these auxiliary parameters to achieve the best accuracy has not yet been discussed. The iterative risk allocation algorithm is proposed to achieve a suboptimal solution for joint chance-constrained programs.30 To see more approaches and applications of chance-constrained programming, we recommend a comprehensive review by Geletu et al.31 and references therein. In the present paper, we employ the popular Bonferroni approximation and directly optimize over the parameters of this approximation. This gives rise to a nonconvex individual chance-constrained problem that we solve to global optimality by using techniques from global optimization. This yields a near-global solution to the original joint chance-constrained blend planning problem. In particular, this is achieved by forming successive outer approximations for the inverse cumulative distribution function (CDF) of the Gaussian distribution (in the spirit of work by Cheng et al.32) combined with McCormick relaxations33 to obtain tighter bounds. The rest of this paper is organized as follows. The deterministic and chance-constrained formulations for optimal blending design are stated in section 2. A preliminary discussion of conservative approximation to the joint-chance constrained problem based on the Bonferroni approximation is introduced in section 3. The main method, a global optimization framework for ICC, is proposed in section 4. Case studies are presented in section 5 to show comparative results and highlight the effectiveness of the proposed method. Finally, conclusions are drawn in section 6. Notation. Throughout this paper, vectors (matrices) are denoted by lowercase (uppercase) boldface letters. The space of symmetric positive semidefinite matrices of dimension n is denoted by +n . We let e denote a vector of all ones. Its dimension will be clear from the context. For μ ∈ n and ϕ ̅ ( ·; μ , Σ): n →  Σ ∈ +n , we let and n Φ̅ ( ·; μ , Σ):  →  denote the probability density function and cumulative distribution function of the n-variate normal distribution with mean vector μ and covariance matrix Σ, respectively. Moreover, let Φ−1( ·):  →  and ϕ( ·):  →  represent the standard inverse cumulative function and probability distribution function for univariate normal distribu-

specification with high probability. However, they often result in suboptimal and overly conservative solutions. To mitigate infeasibility and suboptimality of the aforementioned approaches that are routinely employed in the refining industry, we propose to explicitly express in the refinery’s multiproduct blending problem the requirement that each product should be on specification with high probability. This gives rise to a so-called joint chance-constrained (JCC) optimization problem in which the number of joint chance constraints equals the number of end-products and each joint chance constraint ensures that all of the tracked qualities be jointly on specification with high probability, thereby guaranteeing that each end-product is on specification. We assume, as is commonplace among practitioners, that the uncertain product qualities are jointly Gaussian distributed. This assumption is very natural because, in this context, uncertainty stems from the addition of a large number of independent effects. Moreover, the variance of Gaussian distribution can be adjusted to approximate bounded distributions according to the three-sigma rule.2 We also assume that the mean and covariance of the distribution is perfectly known. This assumption is reasonable because refineries have access to large amounts of data recorded over long periods of time as a byproduct of their daily operations. The idea of explicitly accounting for uncertainty in process system optimization is not new. Indeed, the chemical engineering community has proposed numerous optimization approaches that explicitly account for uncertainty in problem parameters. These include design flexibility approaches,3−6 stochastic programming approaches,7−9 chance-constrained optimization approaches,10−13 and, more recently, robust optimization approaches.14−16 Stochastic programming assumes that uncertain parameters satisfy some probability distributions. Different sampling methods, for example, scenario approximation and sample average, are employed to generate a number of scenarios and approximate the true feasible region. Then, the resulting largescale formulation is solved to the global optimum even for nonconvex problems.17 However, this scheme suffers from at least two limitations. First, the computational demand is extremely large, which is usually addressed by Benders method and its generalizations.18−21 Second, how to decide the number of sampled scenarios is unclear. A threshold of sample complexity to guarantee the probabilistic feasibility for convex program is developed by Calafiore and Campi.22 However, the required number of samples grows quickly and makes it impractical for even a medium-size problem.23 Robust optimization tries to find a worst-case-oriented solution given a bounded uncertainty set without any assumption on the probability distributions of the uncertain parameters. This method, although providing an absolutely safe solution subject to bounded uncertainties, is usually very conservative. Recently, how to reduce the conservativeness of robust optimization is intensively studied.16,24,25 The key idea is to decide the type and size of the uncertainty set by using prior or posteriori probability bounds such that chance constraints instead of hard ones are satisfied. Chance-constrained programming (CCP) is the method mainly considered in this paper. The merit of CCP is that it introduces the probabilistic constraint into the formulation and is particularly suitable to handle uncertainty with known statistical information. However, if there are multiple probabilistic constraints, then the problem becomes notoriously B

DOI: 10.1021/acs.iecr.7b02434 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


multiplying left- and right-hand sides by eTxp. The linear programming formulation DB is generically referred to as the blending or mixing problem in operations research, which can be solved very efficiently with off-the-shelf solvers. In practice, refinery operators solve the blending problem DB to determine plant-wide operations (e.g., quantities of blendstocks to produce) and schedule enterprise-wide transportation (e.g., quantities of blendstocks and/or crude oil to ship to the refinery, quantities of products to ship from the refinery, frequency of shipping). Thus, at the time when DB is solved, the true qualities of blendstocks are actually unknown. This implies that, if the operator solves DB by fixing ξq to a nominal value (typically chosen to be the expected value), the proposed blending scheme will often yield off-spec products and thus cannot be sold in the market, leading to a significant loss. The issue mentioned above implies that the refinery operator should explicitly account for uncertainties when a blending plan is developed. Namely, cases when the blending plan leads to off-spec products should be infrequent. In the following section, we formulate the central problem of this paper, which precisely achieves this goal. Blend Planning under Uncertainty. Because qualities are uncertain, they henceforth are modeled as random variables. We denote the quality q for all blendstocks by ξq̃ ∈ B,

tion, respectively. For a random vector ξ ̃ ∈ n, we write ξ ̃ ∼ 5(μ , Σ) to express that ξ̃ is normally distributed with mean μ and covariance matrix Σ.

■

PROBLEM FORMULATION Deterministic Blend Planning. Optimal blending lies at the core of refinery operations. It is a major cost and profit driver for every oil company that operates refineries. To obtain product p ∈ 7 ≔ {1 ,..., P} (e.g., regular or premium gasoline) that meets certain quality specifications (or targets), it is concerned with determining the optimal quantities xb , p ∈ + of each component (or blendstock) b ∈ ) ≔ {1, ..., B} to blend or mix together. We denote the target (maximum acceptable value) for quality q ∈ 8 ≔ {1, ..., Q } in product p by tq,p and let ξq,b denote the value of quality q in blendstock b ∈ ) . Quality targets are imposed by law to prevent oil companies from selling substandard products, which may damage engines and/or pollute the environment. Typical qualities monitored include octane number, Reid vapor pressure (RVP), and sulfur and benzene contents. When all quality targets for a given product are met (i.e., when the product quality does not exceed the target for any monitored quality), the product is said to be “on-spec”. Conversely, it is said to be “off-spec” if any one of the quality values exceeds the target. Throughout this paper, we assume the linear mixing rule, implying that quality q of a product p is a linear function of the fraction of each blendstock used to produce it. Thus, the quality constraint is expressed as ξqTxp e Txp

∀ q ∈ 8 . We model uncertainty by the probability space (B × Q , - , ), which consists of the sample space B × Q , the σalgebra - and the probability measure . In practice, it is required that the product is on-spec with probability greater than 1 − ϵ, where ϵ ∈ (0, 1). Usually, ϵ ≪ 1 is a small positive constant chosen by the refinery operator and captures the tolerance for off-spec blends. This yields the chance-constrained blend planning problem

≤ tq , p , ∀ q ∈ 8, ∀ p ∈ 7 (1)

where vector ξq ≔ {ξq , b} ∀ b ∈ ) and xp ≔ {xb , p} ∀ b ∈ ) . Thus, for

P

each q ∈ 8 , the vector ξq represents the qualities q of all available blendstocks, and the vector xp determines the quantities of each blendstock mixed together to produce the end product p. Then, the numerator in (1) represents the value of quality q in product p, and the denominator expresses the quantity (units) of product p produced. The refinery operation incurs a cost cb ∈ + for acquiring (or producing) one unit of blendstock b. On the other hand, selling the end product on the market yields a revenue f p ∈ + , ∀ p ∈ 7 per unit. The objective of the refinery

min

X∈?

s.t.

s.t.

∑ r pTxp p=1

ξqTxp ≤ tq , pe Txp , ∀ q ∈ 8 , ∀ p ∈ 7

T (ξq̃ xp ≤ tq , pe Txp , ∀ q ∈ 8) ≥ 1 − ϵ, ∀ p ∈ 7

Thus, CB optimizes over all blends in ? for which each product will be on-spec with probability greater than 1 − ϵ. Note that a product is considered to be on-spec only when all qualities jointly meet their target specifications, giving rise to a so-called joint chance constraint. Several comments are in order. First, problem CB assumes perfect knowledge of the distribution of the uncertain qualities. This assumption is realistic in the present setting because uncertainty in Ξ ≔ {ξq}q ∈ 8 typically originates from either measurement error, with the test repeatability and reproducibility of the measuring equipment employed being generally well understood, or from, e.g., uncertainty in the suppliers, with the probability of a given supplier being chosen being wellknown. Second, the choice of ϵ in formulation CB is difficult to make. On the one hand, a low value of ϵ will result in qualified products with high probability, thus lowering the probability of and costs associated with reblending. On the other hand, a high value of ϵ will generate a more profitable blending plan when the resulting blend is feasible to certain scenarios. In this paper, we will assume that the oil company (or other refinery operator) has chosen ϵ in an optimal way to mitigate the tradeoff between these two costs. From experience, we have seen

P X∈?

p=1

(CB)

optimization is to maximize profit defined as the difference between revenues made from selling the end product and the cost of acquiring all blendstocks. The deterministic blend planning problem faced by the refinery can thus be formulated compactly as min

∑ r pTxp

(DB)

where vector rp ≔ (c − f pe) is the negative of the resulting profit of product p using each blendstock. Minimizing the negative of profit is equivalent to maximizing the profit. X = [x1 x2... xp] represents the blendstock quantities for each product. ? ⊂ +B × P imposes resource constraints (linear) on the blendstocks and any other requirements on blendstock quantities. The constraints in DB are obtained from (1) by C


Article

Industrial & Engineering Chemistry Research that oil companies typically aim for values of ϵ in the range from 1 to 5%. Finally, it is worth noting that problem CB is typically nonconvex and thus difficult to solve to global optimality. The difficulty in solving CB stems from two sources: first, the mathematical complexity of evaluating multivariate integrals, and second, the joint nature of the events whose probability is difficult to calculate. A common approach for solving stochastic problems of the form CB is to decompose the joint chance constraint into Q individual chance constraints and to enforce each constraint to hold with probability greater than ϵq,p, where ϵq,p are fixed (preselected) constants satisfying ∑q ∈ 8 ϵq , p = ϵ, ∀ p ∈ 7 . This yields problem

cone program)34 that can be solved very efficiently with off-theshelf solvers. Note in particular that this is the case if ϵ < 0.5, which is typical for the application considered. For any choice of ϵq,p such that ∑q ∈ 8 ϵq , p = ϵ: (a) any feasible solution to ICB′ is feasible in CB, and (b) the optimal value of problem ICB′ yields an upper bound to the optimal value of CB. It is thus natural to seek the best (i.e., smallest) upper bound by optimizing over ϵq,p, q ∈ 8 , p ∈ 7 , giving rise to the following nonconvex optimization problem min

X , ϵp , vp , ∀ p ∈ 7

P

min

X∈?

s.t.

s.t.

∑ r pTxp

p∈7

X ∈ ?, ϵp ∈ Q+ , vp ∈ Q+ , ∀ p ∈ 7 ⎫ μqT xp + vq , p xpTΣqxp ≤ tq , pe Txp , ⎪ ∀ q ∈ 8 ⎬ ⎪∀ p ∈ 7 vq , p = Φ−1(1 − ϵq , p), ⎭

p=1 T (ξq̃ xp

∑ r pTxp

T

≤ tq , pe xp) ≥ 1 − ϵq , p , ∀ q ∈ 8 , ∀ p ∈ 7 (ICB)

which, from the Bonferroni inequality, corresponds to a conservative (inner) approximation to problem CB. In the following section, we will detail the reformulation of problems ICB as tractable convex problems in the case when  is multivariate normal.7 Subsequently, we propose a method for jointly optimizing over X and ϵq , p , ∀ q ∈ 8, ∀ p ∈ 7 to obtain a global optimal solution of ICB even though it is nonconvex. This will enable us to obtain much less conservative approximations to CB.

∑ ϵq,p = ϵ, ∀ p ∈ 7 q∈8

(ICBG)

where vp = [v1,p v2,p...vq,p]T and ϵp = [ϵ1,p ϵ2,p...ϵq,p]T. The nonconvexity of ICBG stems from the term vq , p xpTΣqxp .

■

■

GLOBAL OPTIMIZATION OF ICBG

Although given a certain ϵq,p ICBG becomes a convex program and easy to solve, how to determine the ϵq,p is not well-studied systematically. In this section, we propose a global optimization-based method to determine ϵ and x simultaneously. The first concern is that the inverse CDF: Φ−1 does not have closed-form suitable for optimization solver. It is worthwhile to note that ϵq,p → Φ−1(1 − ϵq,p) is a convex function over the interval [0, 0.5], where any ϵq,p feasible in ICBG satisfies ϵq,p ∈ [0, 0.5] provided ϵ ∈ [0, 0.5]. Motivated by Cheng et al.,32 we note that Φ−1(1 − ϵq,p) can be outer approximated by a number of cutting planes, as seen in Figure 1. Thus, it is not difficult to develop the following relaxation of ICBG

THE BONFERRONI APPROXIMATION AND THE CASE OF NORMALLY DISTRIBUTED QUALITIES In this section, we briefly introduce the case when the ξq̃ , q ∈ 8 , are normally distributed and the exact reformulation of ICB (conservative approximation to CB) is in the form of a second-order cone program (SOCP). We formalize our present assumption as follows: 1. (AG) For each q ∈ 8 , it holds that ξq̃ ∼ 5(μq , Σq) for some known mean vector μq ∈ B and covariance matrix Σq ∈ +B . The reformulation of ICB as a second-order cone program under Assumption 1 relies on the following proposition. Proposition 17. Suppose ξ ̃ ∼ 5(μ , Σ), where μ ∈ B and Σ ∈ +B . Then, for any fixed x , g ∈ B and ϵ ∈ (0, 1), the following statements are equivalent: T (a) (ξ ̃ x ≤ g T x) ≥ 1 − ϵ, and

(b) μT x + Φ−1(1 − ϵ) xT Σx ≤ g T x . This result is well-known to the stochastic programming community.7 From Proposition 1, we obtain the following reformulation of ICB under Assumption 1 as a deterministic optimization problem minX ∈ ?

∑ r pTxp p∈7

s.t.

μqT xp + Φ−1(1 − ϵq , p) xpTΣqxp ≤ tq , pe Txp , ∀ q ∈ 8, ∀ p ∈ 7

(ICB′) −1

If ϵq , p < 0.5, ∀ q ∈ 8, ∀ p ∈ 7 , then Φ (1 − ϵq,p) > 0, and thereby, problem ICB′ is a convex optimization (second-order

Figure 1. Outer approximation of Φ−1(1 − ϵq,p). D


Article


directed to works by Moore36 and Neumaier.37 The McCormick relaxations (set of last four inequality constraints in LBP) construct the convex hull for a single bilinear term,38 and thereby, LBP is convex. Therefore, it is not difficult to solve it and derive a lower bound on the optimal objective function value. Given the solution ϵ*q , p ,LBP, we can easily calculate Φ−1(1 − ϵ*q , p ,LBP) and substitute it into ICBG. This may result

P

min

X ∈ ?, ϒ ∈ Q+ × P , V ∈ Q+ × P

s.t.

∑

r pTxp

p=1

μqT xp + vq , p xpTΣqxp ≤ tq , pe Txp , ∀ q ∈ 8, ∀ p ∈ 7

∑ ϵq,p = ϵ, ∀ p ∈ 7

in an upper bound on the optimal objective function value if it is feasible. To converge any gap between upper and lower bounds, we need to use an iterative procedure to solve lower and upper bounding problems several times. The main steps include adaptive outer approximation, branch-and-bound, and optimality-based bounds tightening. Remark: One may consider solving (3) directly by using global optimization software to obtain a tighter lower bound. However, we have tried BARON39 and ANTIGONE.40,41 Both of them cannot provide a satisfactory solution within 1 h. Another potential solver is SCIP,42 which can recognize and handle the SOCP structure. Its solution time is still much longer than our proposed method, and required memory space exceeds 1 GB in our case study. Adaptive Outer Approximation. The outer approximation is only accurate near the points where cutting planes are constructed. Thus, there will generally be a gap between the vq*, p ,LBP and Φ−1(1 − ϵ*q , p ,LBP). To avoid large gaps, one may consider adding cutting planes adaptively. For example, for large ϵq,p, Φ−1(1 − ϵq,p) is relatively smooth and requires a small number of cutting planes. For small ϵq,p, Φ−1(1 − ϵq,p) changes rapidly, and we need to place more cutting planes. Alternatively, the sandwich algorithm with different error rules may be used for building supporting planes.43 However, we employ another strategy to add cutting planes adaptively because we notice that the optimal ϵq,p usually concentrate in a small region. Instead of exhaustively searching the point with maximum outer approximation error, we let LBP tell us how to specify the locations of cutting planes. We first initialize a number of cutting planes from 1 − ϵ to 0.99999 with spacing 0.002. Then, given the ϵ*q , p ,LBP during the iterations, we add a new cutting plane to LBP if Φ−1(1 − ϵ*q , p ,LBP) ⩾ vq*, p ,LBP + θ , where θ is a prespecified parameter to balance accuracy and computational load. Let 2q , p ← 2q , p ∪ {l}, and the coefficients of the new cutting plane are

q∈8

vq , p ⩾ aq , p , k + hq , p , k ϵq , p , ∀ k ∈ 2q , p, ∀ q ∈ 8, ∀ p ∈ 7 (2)

where V = [v1 v2... vP], ϒ = [ϵ1 ϵ2... ϵP], 2q , p denotes the collection of cutting planes used to underestimate Φ−1(1 − ϵq,p) with parameters aq,p,k and hq,p,k, which are generated on the fly during optimization. By introducing wq , p = vq , pxp , ∀ p ∈ 7 , ∀ q ∈ 8 , we reformulate (2) as min

∑ r pTxp p∈7

s.t.

X ∈ ?, ϒ ∈ Q+ × P , V ∈ Q+ × P , Wq ∈ +B × P , ∀ q ∈ 8 ⎫ wqT, pΣqwq , p ≤ tq , pe Txp , ⎪ ⎪∀ q ∈ 8 ≥ aq , p , k + hq , p , k ϵq , p ∀ k ∈ 2q , p, ⎬ ⎪∀ p ∈ 7 ⎪ wq , p = vq , pxp , ⎭

μqT xp + vq , p

∑ ϵq ,p = ϵ ∀ p ∈ 7 q∈8

(3)

where Wq = [wq ,1wq ,2 ... wq , P], ∀ q ∈ 8 . Problem (3) is not convex due to the equality constraint involving bilinear terms. However, we can use McCormick relaxation to obtain a lower bounding problem in SOCP form min

∑ r pTxp p∈7

s.t.

X ∈ ?, ϒ ∈ Q+ × P , V ∈ Q+ × P , Wq ∈ +B × P , ∀ q ∈ 8 ⎫ wqT, pΣqwq , p ≤ tq , pe Txp ⎪ ⎪ vq , p ≥ aq , p , k + hq , p , k ϵq , p , ∀ k ∈ 2q , p ⎪ ⎪ wq , p ≥ vq̅ , pxp + vq , pxp̅ − vq̅ , pxp̅ ⎪∀ q ∈ 8 ⎬ wq , p ≥ v̲ q , pxp + vq , p x̲ p − v̲ q , p x̲ p ⎪ ∀ p ∈ 7 ⎪ wq , p ≤ vq̅ , pxp + vq , p x̲ p − vq̅ , p x̲ p ⎪ ⎪ wq , p ≤ v̲ q , pxp + vq , pxp̅ − v̲ q , pxp̅ ⎪ ⎭ μqT xp +

hq , p , l =

d Φ−1(1 − ϵq , p) d(ϵq , p)

|ϵq,p =ϵ*q,p,LBP =

−1 ϕ(Φ−1(1 − ϵ*q , p ,LBP))

aq , p , l = Φ−1(1 − ϵ*q , p ,LBP) − hq , p , l ϵ*q , p ,LBP

These equations in fact are tangent to Φ−1(1 − ϵq , p ,LBP) at the point ϵ*q , p ,LBP. As we mentioned before, the ϵ*q , p ,LBP are usually within a small region, and thereby, the outer approximation can be very accurate after just a few iterations. Note that adding planes at a specific ϵ will lift its corresponding objective value. Thus, it is unnecessary to locate the maximum outer approximation error point if it is not the optimal solution of LBP (minimization). Branch-and-Bound. Apart from the outer approximation of the inverse CDF, the relaxation of bilinear terms will also result in relaxation gap. The tightness of the McCormick

∑ ϵq , p = ϵ, ∀ p ∈ 7 q∈8

(LBP)

where the upper and lower bars represent the upper and lower bounds on variables, respectively. The bounds on x are initially determined according to the allowable quantity of feedstock. The bounds on v are initially determined by the properties of function Φ−1. Then, the forward/reverse interval analysis35 refines enclosures of x and v iteratively on the inequalities in LBP. For detailed discussion on interval analysis, readers are E


Article

Industrial & Engineering Chemistry Research relaxations is dependent on the variable bounds. Thus, we can use branch-and-bound (B&B) to partition the variable domain and obtain better lower bounds iteratively. Then, substituting ϵq,p from LBP into ICBG may yield a better upper bound solution. Here, we can partition either vq,p or xb,p according to the solution from LBP. For vq,p, there are in total QP variables that can be selected for partitioning. A conventional choice is the qualities and products corresponding to the maximum relaxation gap of the bilinear terms

s.t. μqT xp +

wqT, pΣqwq , p ≤ tq , pe Tx , ∀ q ∈ 8 , ∀ p ∈ 7

P

∑ r pTxp ⩽ UB p=1

∑ ϵq,p = ϵ, ∀ p ∈ 7 q∈8

⎫ ⎪ wq , p ⩾ v̲ q , pxp + vq , p x̲ p − v̲ q , p x̲ p ⎪ ⎪ wq , p ⩽ vq̅ , pxp + vq , p x̲ p − vq̅ , p x̲ p ⎪ ⎪∀ q ∈ 8 wq , p ⩽ v̲ q , pxp + vq , pxp̅ − v̲ q , pxp̅ ⎬ ⎪∀ p ∈ 7 vq , p ⩾ aq , p , k + hq , p , k ϵq , p , ∀ k ∈ 2q , p ⎪ ⎪ v̲ q , p ⩽ vq , p ⩽ vq̅ , p , ϵ̲ q , p ⩽ ϵq , p ⩽ ϵq̅ , p ⎪ ⎪ ⎭ wq , p ⩾ vq̅ , pxp + vq , pxp̅ − vq̅ , pxp̅

∑ |wb*,q ,p,LBP − vq*,p,LBPxb*,p,LBP| q ∈ 8, p ∈ 7

(q′, p′) ∈ arg max

b∈)

where wb*, q , p ,LBP is the bth element of vector wq*, p ,LBP. However, note that the objective function value of LBP is in fact highly impacted by the term wqT, pΣqwq , p . A more proper selection is to take Σq into account. Therefore, we use the following criterion to choose the partitioned vq,p T (q′, p′) ∈ arg max |wq*, pT,LBPΣqwq*, p ,LBP − (vq*, p ,LBPx*p ,LBP )

x̲ b , p ⩽ xb , p ⩽ xb̅ , p , ∀ p ∈ 7, ∀ b ∈ )

q ∈ 8, p ∈ 7

(OIRx)

Σq(x*p ,LBPvq*, p ,LBP)|

where UB is the current upper bound on the solution. The new bounds for xb,p are just the solution of minimization and max * * maximization problems: [xbmin , p ,OIR x , xb , p ,OIR x ], where the superscript max* and min* represent the solution values of the maximization and minimization problems, respectively. The constraint ∑p P= 1 rTp xp ⩽ UB is from the objective function (minimization). It is clear that any xb,p that results in a larger objective function value than the current upper bound is not necessary to be considered in the future search. Because we solve relaxations of ICBG for maximization and minimization with ∑pP= 1 rTp xp ⩽ UB embedded, the resulting solution is a valid bound for ICBG equipped with ∑pP= 1 rTp xp ⩽ UB and the same initial variable intervals. For another variable vq,p, we also use the similar formulation

Once the (q′, p′) is selected, we can divide the vq′,p′ domain into two sections, [ v̲ q ′ , p ′ , vq*′ , p ′ ,LBP] and [vq*′ , p ′ ,LBP , vq̅ ′ , p ′], and build the McCormick relaxations based on these two parts, respectively. When vq′,p′ is divided into two sections, then there are two nodes added to the search tree with objective function value objLBP. In addition, the variable bounds of these two nodes inherit the bounds in the current iteration except vq′,p′. In the next iteration, the node with lowest objective function value is selected to solve. Alternatively, we can also partition the xb,p. There are BP variables that can be selected for partitioning. The criterion to select xb,p is (b′, p′) ∈ arg max

∑ |wb*,q, p,LBPΣq, b ,bwb*,q, p,LBP

b ∈ ), p ∈ 7 q ∈ 8

min

X ∈ ?, ϒ ∈ Q+ × P , V ∈ Q+ × P , Wq ∈ +B × P , ∀ q ∈ 8

− (vq*, p ,LBPxb*, p ,LBP)Σq , b , b(xb*, p ,LBPvq*, p ,LBP)|

where Σq,b,b is the element in the bth row and bth column of matrix Σq. An empirical comparison of convergence rates for different branching schemes is conducted in the case study part. Optimality-Based Bounds Tightening. As mentioned in the previous section, the tightness of the McCormick relaxations is dependent on the bounds of variables in the bilinear terms. Therefore, bounds tightening is highly necessary to reduce the gap and speed up the algorithm. Integrated with the branch-and-bound framework, if a node includes the optimal solution, the bounds tightening approach will always generate a nonempty interval enclosing the optimal point; otherwise, the bounds tightening may be infeasible and fathom that node. Here, only two types of variable, x and v, are involved in the bilinear terms. Traditionally, interval analysis is used to obtain tight bounds with minimal calculation demands. However, here we employ a more aggressive way to obtain much tighter bounds. Namely, we directly solve minimization and maximization problems for variables subject to the constraints of LBP. For example, for xb,p, we have the formulation min

Q ×P ×P B×P X ∈ ?, ϒ ∈ Q + , V ∈ + , Wq ∈ + , ∀ q ∈ 8

s.t. μqT xp +

vq , p

wqT, pΣqwq , p ≤ tq , pe Tx , ∀ q ∈ 8 , ∀ p ∈ 7

P

∑ r pTxp ⩽ UB p=1

∑ ϵq,p = ϵ, ∀ p ∈ 7 q∈8

⎫ ⎪ wq , p ⩾ v̲ q , pxp + vq , p x̲ p − v̲ q , p x̲ p ⎪ ⎪ wq , p ⩽ vq̅ , pxp + vq , p x̲ p − vq̅ , p x̲ p ⎪ ⎪∀ q ∈ 8 wq , p ⩽ v̲ q , pxp + vq , pxp̅ − v̲ q , pxp̅ ⎬ ⎪∀ p ∈ 7 vq , p ⩾ aq , p , k + hq , p , k ϵq , p , ∀ k ∈ 2q , p ⎪ ⎪ v̲ q , p ⩽ vq , p ⩽ vq̅ , p , ϵ̲ q , p ⩽ ϵq , p ⩽ ϵq̅ , p ⎪ ⎪ ⎭ wq , p ⩾ vq̅ , pxp + vq , pxp̅ − vq̅ , pxp̅

x̲ b , p ⩽ xb , p ⩽ xb̅ , p , ∀ p ∈ 7, ∀ b ∈ )

xb , p

(OIRv) F


Article

Industrial & Engineering Chemistry Research * min * max * max * Let {ϵqmin , p ,OIR v , vq , p ,OIR v } and {ϵq , p ,OIR v , vq , p ,OIR v } be the solutions of OIRv by minimizing and maximizing vq,p, respectively. Then, straightforward new bounds for vq,p are max * * v̲ q , p ,OIR v = vqmin , p ,OIR v , vq̅ , p ,OIR v = vq , p ,OIR v

(4)

However, tighter bounds than (4) can be computed by (5) max * * v̲ q , p ,OIR v = Φ−1(1 − ϵqmin , p ,OIR v ), vq̅ , p ,OIR v = vq , p ,OIR v

(5)

This claim can be shown by the following proposition: Proposition 2. If OIR v is feasible, (5) generates valid and tighter bounds than (4) for ICBG. Proof. vq̅ , p ,OIR v in (5) is a valid upper bound of feasible vq,p for ICBG because the feasible region of OIR v is the same as LBP, which is a relaxation of ICBG. This implies that ∀ϵ′q,p, such that vq̅ , p ,OIR v ⩽ Φ−1(1 − ϵ′q , p), cannot be the solution of ICBG. Note that Φ−1(1 − ϵq,p) is a decreasing function of ϵq,p; thus, we can claim that ϵ̲ q , p ,OIR v = 1 − Φ( vq̅ , p ,OIR v) is a valid lower bound for ICBG. * We can also claim that for any ϵq,p ′ such that ϵ′q , p > ϵqmin , p ,OIR v , it cannot be feasible to OIR v . Otherwise, ∃v′q,p, feasible to the OIR v , such that

−1 min * * Figure 2. Difference between vqmin , p ,OIR v and Φ (1 − ϵq , p ,OIR v ).

{ϵ′q , p ,1, ϵ′q , p ,2 , ..., ϵ′q , p , |2q,p|}, ∀ q ∈ 8, ∀ p ∈ 7 , on which the outer approximation is constructed. Proposition 3. Given any ϵq,p ∈ [ϵ q,p, ϵ q,p], find i, such that ϵ′q,p,i ⩽ ϵq,p and ϵ′q,p,i+1 ⩾ ϵq,p, then aq , p , i + hq , p , i ϵq , p ⩾ aq , p , k + hq , p , k ϵq , p , ∀ k < i ,

min * * vqmin , p ,OIR v ⩾ max k ∈ 2q , p aq , p , k + hq , p , k ϵq , p ,OIR v

aq , p , i + 1 + hq , p , i + 1ϵq , p ⩾ aq , p , k + hq , p , k ϵq , p , ∀ k > i + 1

Proof. The cutting plane placed at point ϵq,p,k ′ can be rewritten as Φ−1(1−ϵq,p,k ′ ) + hq,p,k(ϵq,p−ϵq,p,k ′ ). Because Φ−1(1 − ϵq,p) is a convex function, then

> max k ∈ 2q aq , p , k + hq , p , k ϵ′q , p = vq′, p

The second inequality holds because the outer approximation max k ∈ 2q,p aq , p , k + hq , p , k ϵq , p is also a decreasing function of ϵq,p. However, the above inequality

* vqmin , p ,OIR v

Φ−1(1 − ϵ′q , p , i) ⩾ Φ−1(1 − ϵ′q , p , k ) + hq , p , k (ϵ′q , p , i − ϵ′q , p , k ), ∀ k < i

> vq′, p contradicts the

Note that h is a function of ϵ, thus

* ϵqmin , p ,OIR v

is in fact a valid upper bound for result of OIRv. Thus, ϵq,p in OIR v and ICBG. In addition,

Φ−1(1 − ϵq , p) dϕ(Φ−1) dh 1 = 2 −1 >0 (1 − ϵq , p) = 2 −1 dϵq , p ϕ (Φ (1 − ϵq , p)) dϵq , p ϕ (Φ (1 − ϵq , p))

min * * vqmin , p ,OIR v = max k ∈ 2q , p aq , p , k + hq , p , k ϵq , p ,OIR v

This implies hq,p,i ⩾ hq,p,k because ϵq,p,k ′ < ϵp,q,i ′ , ∀k < i. Given condition ϵq,p ⩾ϵq,p,i ′ , it follows that

Then, it is not difficult to see that −1

Φ (1 −

* ϵqmin , p ,OIR v )

⩾ max k ∈ 2q,p aq , p , k +

* hq , p , k ϵqmin , p ,OIR v

=

Φ−1(1 − ϵ′q , p , i) + hq , p , i(ϵq , p − ϵ′q , p , i)

* vqmin , p ,OIR v

⩾ Φ−1(1 − ϵ′q , p , k ) + hq , p , k (ϵ′q , p , i − ϵ′q , p , k ) + hq , p , i(ϵq , p − ϵ′q , p , i) ⩾ Φ−1(1 − ϵ′q , p , k ) + hq , p , k(ϵq , p − ϵ′q , p , k )

. This shows that (5) is a valid and tighter lower bound for vq,p in ICBG. Given new bounds for vq,p, we can always find the corresponding new bounds for ϵq,p

This means that aq,p,i + hq,p,iϵq,p ⩾ aq,p,k + hq,p,kϵq,p. The second inequality in the proposition can be shown in the same way and is thus omitted. According to this proposition, given the interval [ϵ q,p, ϵ q,p], we need to identify the i and j, such that ϵq,p,i ′ ⩽ ϵ q,p ⩽ ϵq,p,i+1 ′ and ϵ′q,p,j−1 ⩽ ϵ q,p ⩽ ϵ′q,p,j. Then, only the cutting planes constructed based on points within [ϵ′q,p,i, ϵ′q,p,j] need to be considered. Finally, we show the flowchart of the algorithm in Figure 3, which is guaranteed to converge to the globally optimal solution of (ICBG).

ϵq̅ , p ,OIR v = 1 − Φ( v̲ q , p ,OIR v), ϵ̲ q , p ,OIR v = 1 − Φ( vq̅ , p ,OIR v) (6)

With the new bounds for x and v, the LBP will become tighter. One may consider to resolve OIR v and OIR x until the bounds are not changed. However, in our implementation, we only conduct the bounds tightening once for each node. The bounds in (5) are very useful especially when cutting planes are not placed well. Figure 2 shows an example of the difference * −1 min * between vqmin , p ,OIR v and Φ (1 − ϵq , p ,OIR v ). Dropping Constraints. For the size of LBP to be kept relatively small, a constraints dropping procedure is investigated for the outer approximation when the vq,p are branched. We claim the following proposition to show that only a small portion of the outer approximation cuts have to be considered in our scheme. Let us define points set with ascending order,

■

CASE STUDIES In this section, we apply our algorithm and compare it with ANTIGONE, BARON, and SCIP. The problems are solved on a Intel Core i7-7500U CPU, 2.70 GHz with 8 GB memory. The software is GAMS 24.8.3 with CPLEX 12.7.0, CONOPT, BARON 17.1.2, SCIP 3.2, and ANTIGONE 1.1. In the preprocessing step, we use MC++ 1.044 to calculate the G


Article


Figure 3. Algorithm flowchart.

variables bounds. The initial upper bound of ϵq , p , ∀ q ∈ 8, ∀ p ∈ 7 is set as 0.999999. Case 1. The first case is a pure blending problem. By blending 10 types of intermediate flows, three types of gasoline are produced. The parameters of blendstocks are shown in Supporting Information S1. The specifications and prices for three types of gasoline are shown in Supporting Information S2. The maximum production of each type of gasoline is 50. Among all quality specifications, RVP, average octane number, sulfur, and benzene are subject to uncertainties and thereby modeled by chance constraints. The allowed violation chance for each type gasoline is ϵ = 5%. We assume that uncertain parameters satisfy independent Gaussian distributions. This enables us to calculate the satisfaction probability of each product, namely, (ξqTxp ≤ tq , pe Txp , ∀ q ∈ 8) =

∏

Table 1. Profits and Violation Chances for Different Schemesa method

profit

joint violation chance

deterministic LP ICBG (equal ϵq,p) ICBG (global optimization)

744.568 532.863

(87.52%, 87.52%, 87.58%) (3.71%, 3.70%, 3.69%)

585.600

(4.99%, 4.99%, 4.98%)

a

computational time (s) 0.030 0.034 496.949

Relative gap = 0.1%, branching x, for global optimization.

and maximum flow rate of blendstock. To demonstrate the importance of ϵq,p, we compare two schemes: simply distributing the ϵ to probabilistic constraints of each product equally23 and using global optimization to calculate ϵq,p, ∀ q ∈ 8 , ∀ p ∈ 7 . All computational results are shown in Table 1. One can see that the deterministic LP without any protection has the highest profit if it is feasible. Unfortunately, its solution violates specifications with probabilities more than 87%, thereby being unacceptable. The global optimization indeed gains much more profit, 9.9%, than simply specifying the ϵq,p equally. To compare the risk allocation, we know that distributing violation chance equally leads to 1 − ϵq , p = 1 − ϵ/|8| = 1 − 0.05/4 = 0.9875, whereas the proposed scheme derives 1 − ϵq , p , ∀ q ∈ 8, ∀ p ∈ 7 , from the global optimization, as shown in Table 2. By reducing termination tolerance, the results in Tables 3 and 4 show that iterations and computational time increase moderately if

(ξqTxp ≤ tq , pe Txp)

∀ q∈8

The mean of distribution is just the nominal value of feedstock parameters. The standard deviation is shown in Supporting Information S3. We first use the conventional LP method to solve this problem. However, because of the uncertainty, the quality constraints are violated with high probability based on the blending recipes derived from this deterministic optimization scheme, see Table 1. Alternatively, an offset is usually employed to make the LP method more robust. Namely, one can add an offset term to each constraint in the LP formulation. However, the optimal offset is usually decided empirically and rarely discussed systematically. For small-scale problems, this offset method may provide a feasible solution. However, for largescale problems, its performance cannot be guaranteed. In fact, robust optimization or chance-constrained programming are the best ways to satisfy constraints in the face of uncertainty. Thus, the second method is chance-constrained programming. This model has three products and 12 chance constraints in total. There are also deterministic constraints for the density

Table 2. 1 − ϵq,p in the near-Global Optimal Solution

H

gasoline

RVP (max)

(RON + MON)/2 (min)

sulfur (max)

benzene (max)

type 1 type 2 type 3

0.998760 0.999093 0.998683

0.952148 0.951491 0.951806

0.999093 0.999417 0.999512

0.999999 0.999999 0.999999


Article


approximate Φ−1 more accurately, which requires solution of (3) several times. Thus, we can claim that our global optimization approach is fairly more efficient in this problem than general solvers. The evolution of upper and lower bounds is shown in Figure 7. In the implementation, we employ CONOPT to solve (3) directly and substitute the resulting ϵ into ICBG. Solving this ICBG yields an upper bound solution (UB), which is very close to true optimality. Even though this is not necessary, a good UB will benefit the variable tightening, thereby speeding up the algorithm significantly. In addition, Figure 7 also shows that most computational time is for improving the lower bound in this global optimization scheme. This highlights the significance of optimality-based bound tightening method to speed up convergence. Case 2. In this case, a crude oil procurement problem with refinery operations is studied. The flowsheet, shown in Figure 8, and LP model are from the work of Favennec45 with modifications. In total, three types of crude oil can be purchased and processed with upper bound 200 KT. After the initial processing, the intermediate flows are mixed to produce gasoline 95, gasoline 98, jet fuel, diesel, and heavy fuel oil. Different from the original model in Favennec,45 we do not allow importation for any intermediate products, and one more pipeline is added to allow gas oil (GO) flow into the heavy fuel oil tank even though this is not a very economic operation. Some parameters, specifications, and uncertainties are listed in Supporting Information S4−S8. Other data and the LP formulation can be found in the Favennec.45 The objective is to maximize the profit and meet the specifications for each product with 95% chance. This can be modeled by chance-constrained programming. There are in total four products. For gasoline 98 and 95, each has four specifications, including maximum or minimum RVP, RON, sulfur, and sensitivity. For heavy fuel oil, it has three specifications. The diesel only has one specification. For each product, we also assume that the properties are independent such that their joint violation chance can be calculated to verify our results from chance-constrained programming. As in the previous case study, we compare the deterministic LP, SOCP with equal ϵq, and global optimization of ICB. In Table 5, we can see that simply specifying equal ϵq cannot

bounds tightening for both x and v are integrated into the scheme. Table 3. Convergence by Branching v with and without x, v Boundinga time to converge (s)

a

iterations to converge

relative gap %

with bounding

without bounding

with bounding

without bounding

5 1 0.5 0.1

14.662 51.194 185.593 774.164

0.837 218.374 ⩾4059.896 ⩾4059.896

1 7 27 123

1 1053 ⩾8000 ⩾8000

Maximum iteration: 8000.

Table 4. Convergence by Branching x with and without x, v Boundinga time to converge (s)

a

iterations to converge

relative gap %

with bounding

without bounding

with bounding

without bounding

5 1 0.5 0.1

11.589 52.580 107.581 496.949

0.523 11.238 285.050 ⩾3716.340

1 7 15 74

1 61 1545 ⩾8000

Maximum iteration: 8000.

The blending recipes for the three types of gasoline are shown in Figures 4−6. Clearly, the solution of the deterministic LP is far different from the chance-constrained programming. Moreover, even for the chance-constrained programming, the recipes, based on the different distributions of ϵq,p, are also distinct. We also use ANTIGONE, BARON, and SCIP to solve (3) directly. However, ANTIGONE and BARON cannot return a satisfactory solution (relative gap: BARON 2.8%, ANTIGONE 14.5%) of (3) within 1 h even integrated with the interval analysis method (MC++). SCIP takes the SOCP structure into account and thus is suitable for chance-constrained programming. It spends 1413.6 s to find a solution with 0.1% relative gap. Moreover, we may still need to add cutting planes to

Figure 4. Blending recipe for type 1 gasoline. I


Article


Figure 5. Blending recipe for type 2 gasoline.

Figure 6. Blending recipe for type 3 gasoline.

the crude oil purchase. Crude2 and Crude3 are equal or close to the upper bound. The deterministic solver does not prefer Crude1 because it is more expensive than others. However, notice that GO1 has the lowest sulfur content; thus, it can be used to reduce the sulfur fraction in the final products such as diesel and heavy fuel oil. This decision will reduce the profit but increase the reliability of refinery operations. We also show the 1 − ϵq,p by equally distributing violation chance and global optimization in Tables 6 and 7, respectively. One can see that the sulfur and RON constraints are allowed to be violated with more chances because they have more impact on the profitability. The convergence of the algorithm for this case is shown in Figure 9.

■

CONCLUSIONS This paper presents a near-global optimization method to design blending recipes under uncertainties to maximize the profit and guarantee the qualities with high probability. The original joint chance-constrained program is first decomposed into an individual chance-constrained formulation. With an assumption of Gaussian distributions for uncertain parameters,

Figure 7. Upper and lower bounds vs iterations (relative gap = 0.1%).

obtain a feasible solution. Deterministic LP still generates a solution with high risks. Table 5 also shows the procurement of crude oil based on the different schemes. One can see that the uncertainties in the blending part indeed impact the decision of J


Article


Figure 8. Refinery flowsheet.45

Table 5. Profit, Violation Chance, and Computational Time for Different Schemes method

profit (103$)

joint violation chance (98, 95, diesel, HF) %

crude (KT)

solving time (s)

LP ICBG (equal ϵq) ICBG

13323 infeasible 10948

(75.01, 50, 95.05, 49.91)

(43.511, 177.657, 200)

0.026

(4.98, 4.99,5, 4.97)

(50.833, 178.275, 200)

20.488

Table 6. 1 − ϵq,p by Distributing the Violation Chance for Each Product Equally specifications gasoline 98 rvp gasoline 95 rvp gasoline 98 ron gasoline 95 ron gasoline 98 sulfur gasoline 95 sulfur gasoline 95 sensitivity gasoline 98 sensitivity diesel sulfur heavy fuel oil sulfur heavy fuel oil viscosity

max 0.99 0.99

0.99 0.99 0.99 0.99 0.95 0.9833 0.9833

min 0.99 0.99 0.99 0.99

0.9833

Figure 9. Upper and lower bounds vs iterations (relative gap = 0.1%).

Table 7. 1 − ϵq,p in the Global Optimal Solution of ICB specifications

max

min

gasoline 98 rvp gasoline 95 rvp gasoline 98 ron gasoline 95 ron gasoline 98 sulfur gasoline 95 sulfur gasoline 95 sensitivity gasoline 98 sensitivity diesel sulfur heavy fuel oil sulfur heavy fuel oil viscosity

0.999991 0.999991

0.999991 0.999991 0.954902 0.951411

0.999991 0.999991 0.998621 0.995130 0.95 0.954988 0.995026

This convex program can be solved by CPLEX or GUROBI very efficiently to obtain a lower bound. The branch-and-bound approach enables us to further reduce the relaxation gap and finally achieve global optimality of the individual chanceconstrained program. Finally, two case studies are considered to demonstrate the effectiveness of the proposed method.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.7b02434. Blendstock parameters, gasoline specifications and prices, std dev of uncertain parameters, crude yields and prices, specifications, and mean/variance of parameters for gasoline blending, of sulfur for diesel blending, and of sulfur and viscosity for heavy fuel oil blending (PDF)

0.999989

using McCormick relaxations and outer approximation of the inverse cumulative function, we can convert the chanceconstrained formulation into a second-order cone program. K


Article


■

(19) Geoffrion, A. M. Generalized Benders decomposition. J. of Opt. Theory and Appl. 1972, 10, 237−260. (20) Li, X.; Tomasgard, A.; Barton, P. I. Nonconvex generalized Benders decomposition for stochastic separable mixed-integer nonlinear programs. J. of Opt. Theory and Appl. 2011, 151, 425−454. (21) Li, X.; Tomasgard, A.; Barton, P. I. Decomposition strategy for the stochastic pooling problem. J. Global. Optim. 2012, 54, 765−790. (22) Calafiore, C. C.; Campi, M. C. The scenario approach to robust control design. IEEE Trans. Autom. Control 2006, 51, 742−753. (23) Nemirovski, A.; Shapiro, A. Convex approximation of chance constrained programs. SIAM J. Optimiz. 2006, 17, 969−996. (24) Ben-Tal, A.; Ghaoui, L. E.; Nemirovski, A. Robust Optimization; Princeton University Press: New Jersey, 2006. (25) Li, Z.; Floudas, C. A. A comparative theoretical and computational study on robust counterpart optimization: III. Improving the quality of robust solutions. Ind. Eng. Chem. Res. 2014, 53, 13112−13124. (26) Calafiore, G. C.; Ghaoui, L. E. On distributionally robust chance-constrained linear programs. J. of Opt. Theory Appl. 2006, 130, 1−22. (27) Li, P.; Arellano-Garcia, H.; Wozny, G. Chance constrained programming approach to process optimization under uncertainty. Comput. Chem. Eng. 2008, 32, 25−45. (28) Zymler, S.; Kuhn, D.; Rustem, B. Distributionally robust joint chance constraints with second-order moment information. Math. Program. Ser. A 2013, 137, 167−198. (29) Yang, W.; Xu, H. Distributionally robust chance constraints for non-linear uncertainties. Math. Program. Ser. A 2016, 155, 231−265. (30) Ono, M.; Williams, B. C. Iterative risk allocation: a new approach to robust model predictive control with a joint chance constraint. In Proceedings of the 2008 IEEE Conference on Decision and Control, 2008. (31) Geletu, A.; Klöppel, M.; Zhang, H.; Li, P. Advances and applications of chance-constrained approaches to systems optimization under uncertainty. International J. of Systems Science 2013, 44, 1209− 1232. (32) Cheng, J.; Gicquel, C.; Lisser, A. A second-order cone programming approximation to joint chance-constrained linear programs. Lecture Notes in Computer Science 2012, 7422, 71−80. (33) McCormick, G. P. Computation of global solutions to factorable nonconvex programs: Part I convex underestimating problems. Math. Program. 1976, 10, 147−175. (34) Boyd, S.; Vandenberghe, L. Convex Optimization; Cambridge University Press, 2004. (35) Jaulin, L.; Kieffer, M.; Didrit, O.; Walter, E. Applied Interval Analysis; Springer-Verlag: London, U.K., 2001. (36) Moore, R. E. Methods and Applications of Interval Analysis; SIAM: Philadelphia, PA, 1979. (37) Neumaier, A. Interval Methods for Systems of Equations; Cambridge University Press: Cambridge, U.K., 1990. (38) Al-Khayyal, F.; Falk, J. Jointly constrained biconvex programming. Math. Oper. Res. 1983, 8, 273−286. (39) Sahinidis, N. V. BARON: A general purpose global optimization software package. J. of Global. Optim. 1996, 8, 201−205. (40) Misener, R.; Floudas, C. A. GloMIQO: Global mixed-integer quadratic optimizer. J. Global. Optim. 2013, 57, 3−50. (41) Misener, R.; Floudas, C. A. Antigone: algorithms for continuous/integer global optimization of nonlinear equations. J. Glob. Optim. 2014, 59, 503−526. (42) Achterberg, T. SCIP: solving constraint integer programs. Math. Prog. Comp. 2009, 1, 1−41. (43) Tawarmalani, M.; Sahinidis, N. V. Global optimization of mixedinteger nonlinear programs: A theoretical and computational study. Math. Program. Ser. A 2004, 99, 563−591. (44) Chachuat, B. http://www.imperial.ac.uk/people/b.chachuat/ research.html. (45) Favennec, J. P. Refinery operation and management; Editions TECHNIP: Paris, France, 2001.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Yu Yang: 0000-0002-7282-1452 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors are grateful to BP for funding support, and Y.Y. also thanks financial support of Faculty Startup Fund from California State University, Long Beach.

■

REFERENCES

(1) Mendéz, C. A.; Grossmann, I. E.; Harjunkoski, I.; Kaboré, P. A simultaneous optimization approach for off-line blending and scheduling of oil-refinery operations. Comput. Chem. Eng. 2006, 30, 614−634. (2) Pukelsheim, F. The Three Sigma Rule. Am. Stat. 1994, 48, 88− 91. (3) Grossmann, I. E.; Sargent, R. W. H. Optimum design of chemical plants with uncertain parameters. AIChE J. 1978, 24, 1021−1028. (4) Grossmann, I. E.; Halemane, K. P.; Swaney, R. E. Optimization strategies for flexible chemical processes. Comput. Chem. Eng. 1983, 7, 439−462. (5) Grossmann, I. E.; Floudas, C. A. Active constraint strategy for flexibility analysis in chemical processes. Comput. Chem. Eng. 1987, 11, 675−693. (6) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N. Flexibility analysis and design of dynamic processes with stochastic parameters. Comput. Chem. Eng. 1998, 22, 817−820. (7) Prékoba, A. Stochastic programming; Kluwer Academic Publishers: The Netherlands, 1995. (8) Li, X.; Armagan, E.; Tomasgard, A.; Barton, P. I. Stochastic pooling problem for natural gas production network design and operation under uncertainty. AIChE J. 2011, 57, 2120−2135. (9) Gebreslassie, B. H.; Yao, Y.; You, F. Q. Design under uncertainty of hydrocarbon biorefinery supply chains: multiobjective stochastic programming models, decomposition algorithm, and a comparison between CVaR and downside risk. AIChE J. 2012, 58, 2155−2179. (10) Zhang, Y.; Monder, D.; Forbes, J. F. Real-time optimization under parametric uncertainty: a probability constrained approach. J. Process Control 2002, 12, 373−389. (11) Arellano-Garcia, H.; Wozny, G. Chance constrained optimization of process systems under uncertainty: I. Strict monotonicity. Comput. Chem. Eng. 2009, 33, 1568−1583. (12) Ostrovsky, G. M.; Ziyatdinov, N. N.; Lapteva, T. V. Optimal design of chemical processes with chance constraints. Comput. Chem. Eng. 2013, 59, 74−88. (13) Ostrovsky, G. M.; Ziyatdinov, N. N.; Lapteva, T. V. Optimization problem with normal distributed uncertain parameters. AIChE J. 2013, 59, 2471−2484. (14) Rooney, W. C.; Biegler, L. T. Design for model parameter uncertainty using nonlinear confidence regions. AIChE J. 2001, 47, 1794−1804. (15) Al-Qahtani, K.; Elkamel, A. Robust planning of multisite refinery networks: Optimization under uncertainty. Comput. Chem. Eng. 2011, 34, 985−995. (16) Li, Z.; Tang, Q. H.; Floudas, C. A. A comparative theoretical and computational study on robust counterpart optimization: II. Probabilistic guarantees on constraint satisfaction. Ind. Eng. Chem. Res. 2012, 51, 6769−6788. (17) Yang, Y.; Barton, P. I. Integrated crude selection and refinery optimization under uncertainty. AIChE J. 2016, 62, 1038−1053. (18) Benders, J. F. Partitioning procedures for solving mixedvariables programming problems. Numerische Mathematik 1962, 4, 238−252. L


Chance-Constrained Optimization for Refinery Blend Planning under

Recommend Documents