Ind. Eng. Chem. Res. 2009, 48, 1973–1979
1973
Convex Hull Discretization Approach to the Global Optimization of Pooling Problems Viet Pham, Carl Laird, and Mahmoud El-Halwagi* Department of Chemical Engineering, Texas A&M UniVersity, College Station, Texas 77843-3122
The pooling problem is an important optimization problem that is encountered in process operation and scheduling. Because of the presence of bilinear terms, the traditional formulation is nonconvex. Consequently, there is a need to develop computationally efficient and easy-to-implement global-optimization techniques. In this paper, a new approach is proposed based on three concepts: linearization by discretizing nonlinear variables, preprocessing using implicit enumeration of the discretization to form a convex-hull which limits the size of the search space, and application of integer cuts to ensure compatibility between the original problem and the discretized formulation. The continuous quality variables contributing to the bilinear terms are first discretized. The discretized problem is a mixed integer linear program (MILP) which is globally solvable in a computationally effective manner using the branch and bound method. The merits of the proposed approach are illustrated with case studies from literature and comparison with published results. 1. Introduction Market demands with particular qualities may require the blending of several process streams to meet the desired specifications. For instance, intermediate streams from various processing units of a petroleum refinery are typically blended to produce value-added products satisfying quality specifications and demands. As an example, intermediate streams from reforming, cracking, and naphtha treatment units are typically mixed to yield gasoline. Quality specifications such as octane number, vapor pressure, and sulfur and aromatic concentrations are among the stipulated constraints for gasoline. Prior to blending and final storage, intermediate streams are mixed and stored in intermediate tanks, called pools. Such pools enhance the operational flexibility of the process but complicate the factors involved in making optimum decisions. The problem to identify the least costly mixing recipe from intermediate streams to pools and then from the pools to final products is referred to as the pooling problem. Much interest has been given to the solution of pooling problems because of the potential savings in refineries and other process industries. One of the earliest solution procedures is the recursion approach used for the pooling problems introduced by Haverly.1 This procedure is sensitive to the starting point and may not converge or may converge to a local optimum. Lasdon et al.2 applied generalized reduced gradient and successive linear programming (SLP) algorithms to Haverly’s pooling problems. Griffith and Stewart3 were the first to introduce SLP under the name mathematical approximation program and apply in Shell Oil. Palacios-Gomez et al.4 proposed an efficient SLP algorithm for linearly constrained formulation and showed more successful computational results than the generalized reduced gradient algorithm did. Baker and Lasdon5 suggested a multiplicative formulation for the linearized subproblems to be solved by SLP, with nonnegative deviation variables to prevent the occurrence of infeasibility. Decomposition approaches have also been used to solve the pooling problem. Several extensions of the generalized Bender’s decomposition (GBD) method have been developed. Floudas et al.6 used a decomposition-based global optimization approach * To whom correspondence should be addressed. E-mail:
[email protected].
to solve nonlinear programs (NLPs) and mixed integer nonlinear programs (MINLPs) and, later,7 applied it specifically to the pooling problem. The subproblems were structured in such a manner that they are solved for respective global optima in every iteration. However, there is no guarantee that this procedure produces the global optimum of the original optimization problem because every local optimum of a nonlinear program gives rise to a local optimum in the projected problem of Bender.8 Androulakis et al.9 proposed a distributed implementation of a global optimization algorithm (GOP). Branch-and-bound methods have been the basis for many approaches to the pooling problem. Foulds et al.10 partitioned the rectangle of feasible region into four subrectangles and linearized these subproblems using the convex underestimators proposed by McCormick11 and Al-Khayyal and Falk.12 The subrectangle with the best solution was chosen to be the base point for the next partition step. The algorithm continued with smaller and smaller subrectangles and was proven to converge to the optimal solution. Sherali and Alameddine13,14 proposed a two-stage reformulation-linearization technique (RLT) to solve bilinear programming problems. In the first stage, the problem wass reformulated by generating additional valid nonlinear constraints from pairwise multiplications between the original constraints and nonnegative variable bounds. In the second stage, the resulting formulation was linearized by substituting the nonlinear terms with new defined variables. Quesada and Grossmann15 used a branch-and-bound search to solve models of general process networks using this reformulation-linearization technique. Ben-Tal et al.16 partitioned the flow rate proportion variables in order to reduce the duality gap between the primal and its Lagrangian dual problem until convergence was achieved. Piecewise approximation for global optimization of bilinear programs has been demonstrated as a promising technique to improve tightness of lower bounds. The basic idea of piecewise relaxation is to partition the variable range into many intervals and linearize the nonconvex terms in each interval. Meyer and Floudas17 proposed a piecewise linear formulation based on the RLT technique to solve the superstructure model of generalized pooling problems. Before the RLT is applied, the continuous space of each quality is partitioned into many subintervals. In
10.1021/ie8003573 CCC: $40.75 2009 American Chemical Society Published on Web 01/05/2009
1974 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009
Figure 1. General pooling problem.
a large scale industrial problem, the approach can reduce the gap between the upper bound (found by DICOPT) and lower bound to 1.2%. Bergamini et al.18 developed an outer approximation algorithm and applied it to bilinear MINLPs involving concave terms. The algorithm consists of two phases. The outer phase finds a feasible solution for the associated piecewise relaxed problem. Using the solution from the outer phase, the inner phase fixes values of the binary variables to find an upper bound and solve an incremental piecewise relaxed problem for the lower bound. Wicaksono and Karimi19 derived many piecewise relaxed formulations by combining variable partitions (differential or incremental), reformulations (big-M, convex combination, or incremental cost), and segmentation schemes (arbitrary or identical segment length). These formulations can be much more compact than those found in previous literature and may reduce calculation time without loosing the tightness of relaxation. A branch-and-bound procedure using selective Lagrangian relaxation proposed by Adhya et al.20 produced a tighter lower bound than McCormick estimator-based linearization relaxation. Audet et al.21 formulated pooling problems on three models which were based on flow variables, flow proportion variables, and their hybrid ones. This proposed formulation was shown to be suitable for the branch-and-cut quadratic algorithm introduced by Audet et al.22 Sahinidis23 reviewed the theory and algorithms of the branch and reduce approach for the global optimization of NLPs and MINLPs, which has evolved from the traditional branch and bound approach. The branch-andreduced approach consists of two steps. The preprocessing step, applied before relaxation, reduces the range of all variables. After the relaxed problem is solved, the postprocessing step utilizes the solution to further reduce the variable ranges. As a result, this implemented branch and bound converges faster. This approach has been integrated in the computational system BARON. This paper proposes a discretization approach that produces global or near global optimum results in a reasonable computational time. To linearize the bilinear terms (flow variables multiplied by quality variables), the quality variables are finitely discretized. Implicit enumeration of the discretization is used to obtain a convex hull which limits the size of the search space. Next, integer cuts are added to ensure correspondence with the original formulation. The result is a mixed integer linear programming (MILP) formulation which can be easily described and solved by commercial optimization software, e.g. LINGO,24 for its global optimum. 2. Problem Statement The following generalized pooling problem (Figure 1) is investigated: Given is a set SOURCES ) {i|i ) 1, 2,..., l) of intermediate streams. Each source has a given available capacity, Si, one unit cost, Ci, and known values of Nq characterizing
qualities, ai,q, where q is an index for properties (e.g., octane number, Reid vapor pressure, sulfur concentration). The amount (or flow rate) from source i to pool j is denoted by a variable Xij. The sources must be blended because there are not enough pools (m < n) or there are insufficient pool capacities to store them separately. Sources can be sent to some or all of pools. As a result of mixing the sources, each pool, j, has unknown values for the qualities bj,q. The flow from pool j to sale product k is Yjk and is to be determined. The n products with price Pk have constraints on known demand Dk and minimum values of desired quality specifications ckq. The pooling problem is formulated as follows: Objective function: margin return is maximized )
(
n
m
∑ P ∑Y k)1
k
jk
j)1
) ( ) l
-
m
∑ C∑X i
i)1
ij
i)1
This is subject to the following constraints: m
∑X
Available supply:
for i ) 1, ... , l
ij e Si
j)1 l
Mass balance on pools:
n
∑X )∑Y
for j ) 1, ... , m
∑Y
for k ) 1, ... , n
ij
jk
j)1 m
Product demand:
k)1
jk e Dk
j)1 l
Mixing rule for pools: bjq
∑ j)1
l
Xij )
∑X a
ij iq
i)1
for j ) 1, ... , m and q ) 1, ... , Nq Assuming that the desired quality is an upper bound, the following mixing rules can be written for products (similar rules can be written when the desired quality is a lower bound): m
ckq
∑ j)1
m
Yjk g
∑Y b
jk jq
for k ) 1, ... , n and q ) 1, ... , Nq
j)1
Non-negativity constraints: Xij g 0,
Yjk g 0 for i ) 1, ... , l j ) 1, ... , m k ) 1, ... , n
The foregoing formulation is nonconvex because of the bilinear l terms appearing in the mixing rule constraints bjq∑i)1 Xij and m Yjkbjq and may have many local solutions. ∑j)1 Here, we reformulate this problem and develop a global solution approach. 3. Description of Proposed Approach The proposed global optimization approach is based on three concepts: (a) Discretization of qualities for each pool: The characterizing qualities for each pool, bjq, are unknown. Let us discretize the search space of the qualities of the pools into a set of Np vectors of known values: {(bj1 bj2, ..., bjNq)|j ) 1,..., Np). The rationale for the selection of these values will be discussed later. The discretization of the unknown pool qualities into known values transforms the formulation into a linear program. (b) Application of integer cuts for the pools: Since the actual number of pools must be limited to m, Np - m should not be selected in the final solution. Therefore, to ensure consistency between the original problem and the discretized formulation, an integer cut is used to select m pools from among the Np discretized pools. Consequently, the formulation becomes a
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 1975 l
Ljfj e
∑X
for j ) 1, .... , (t + 1)Nq
ij e Ujfj
(2)
i)1
Figure 2. Pools are enumerated on their assigned quality values.
mixed integer linear program (MILP) that can be globally solved using the branch-and-bound method. (c) Convex hull search: A potential problem with the large number of discretizations for several qualities is the large size of the resulting MILP. In order to reduce the problem dimensionality, a convex hull is constructed by invoking physical limits on the possible combinations of pool qualities. The following sections provide more details on the concepts, rationale, and implementation of the above-mentioned steps. 3.1. Discretization by Exclusively Enumerating bj Values. This discretization approach limits the search space to pools whose values are bound between aq,min and aq,max. The quality of any pool is bound by qualities of the blended sources. Using the following notation, aq,min ) argi{aiq} and aq,max ) argi{aiq}, then the domain of bjq is aq,min e bjq e aq,max. When the quality range is discretized into known values, the bilinear terms become linear which is conducive to the global solution of the problem. The quality range may be discretized in various ways (e.g., random, structured). One way of discretizing the range of quality q is to divide it into tq equal intervals. In such cases, the values of discretized quality q are calculated through the following expression: bqr ) aq,min + (rq - 1)
aq,max - aq,min tq
where Lj and Uj are given positive real numbers that correspond to the minimum and maximum capacities of the pools. If the original problem (P) has no constraint on pool capacity, Uj should be at least the sum of available supply of total sources. l Xij From these two inequalities, fj is forced to be 0 when ∑i)1 l ) 0, i.e. no flow rate to discretized pool j, or 1 when ∑i)1 Xij l > 0. As a characteristic of pooling problems, ∑i)1 Xij tends to be positive to contribute to the objective function. Hence, the l Xij can be removed without affecting the inequality Ljfj e ∑i)1 result. Next, the restriction on the number of pools is included in (t+1)Nq fj e m. For simplicity in the formulation simply as ∑j)1 notation, let us denote (t + 1)Nq by M. Therefore, reformulation is given by:
k
M
∑X
ij e Si for
l
∑ P ∑Y k)1
Available supply :
( ) ( ) M
n
Objective function is maximized :
-
jk
j)1
M
∑ C∑X i
i)1
ij
j)1
i ) 1, ... , l
j)1
l
∑
Mass balance on pools :
n
Xij )
i)1
∑Y
for j ) 1, ... , M
jk
k)1
M
Product demand :
∑ eD
for k ) 1, ... , n
k
j)1
l
Mixing rule for pools :
bjq
∑
l
Xij )
j)1
∑X a
ij iq
for j ) 1, ... , M
i)1
Number of pools: fi ) {0, 1}
for j ) 1, ... , M
l
∑X
ij e Uj, fj
for j ) 1, ... , M
i)1 M
for rq ) 1,...,(tq + 1)
∑f em j
(1) When tq f ∞, the solution of the discretized and the original formulations become equivalent. However, for all practical purposes, the number of discretizations is selected to strike the right balance between computational time and proximity of the discretized solution to the original solution. For Nq qualities, let us denote the number of interval discretizations for qualities 1, 2,..., Nq by t1, t2,..., tNq, respectively. Therefore, the total number of discretized pools is (t1 + 1), (t2 + 1),... (tNq + 1). If t1 ) t2 )... ) tNq ) t, there are (t + 1)Nq pools. Figure 2 shows (t + 1)Nq pools, where bj ) {bq,r|q ) 1,..., Nq and r ) 1,..., (t + 1)}. This number of pools (t + 1)Nq is independent of the number of pools in the original problem m. Instead, it depends only on the quality increment t is selected and the number of qualities. 3.2. Integer Cuts for Limitation on Number of Pools. At this point, the problem is to blend l sources into (t + 1)Nq pools with known qualities. It is a linear optimization formulation and a global result can be obtained. However, according to the original problem statement, the number of pools which can be physically used is m. To limit the number of pools to ensure correspondence between the original problem and the discretized version, we add an integer cut. A binary variable fj is introduced. If the pool j is used in the solution of (P′), then fj is assigned a value of 1, otherwise, fj is 0. The number of pools can be restricted using the following linear constraint:
j)1
l
Mixing rule for products : ckq
l
∑Y g∑Y b jk
j)1
jk jq
for k ) 1, ... , n
j)1
Positive variables: Xij g 0, Yjk g 0 for i ) 1, ... , l j ) 1, ... , m k ) 1, ... , n
where bjq values are known parameters defined by (1). This formulation is a mixed integer linear program (MILP). While this reformulation has resulted in a globally solvable problem, it has a limitation in its current form. The number of pools has increased from m to M, resulting in (M - m) more parameters for bj, M variables of fj and l(M - m) more variables for each set of Xij and Yjk. It also increases rapidly when the increment t decreases. Consequently, the selection of the increment size is critically important for the problem dimensionality and computational efficiency of the devised discretization approach. To overcome this limitation, we introduce a convex hull formulation as described in the following section. 3.3. Rationale for a Convex Hull Search To Reduce the Problem Size. Figure 3 depicts a quality diagram of an example adapted from Greenberg,25 in which three sources are presented by S1, S2 and S3, two qualties are investigated (research octane number and sulfur content), and discretized pools are shown as dots. When using (1) to enumerate bj values, the search space of discretized pools spreads the whole boundary rectangle in Figure 3. However, the search space can be reduced.
1976 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009
Figure 3. Quality diagram of an example adapted from the work of Greenberg.25
Every pool’s qualities should be a linear blending resulting from the three sources. For each pool inside the triangle S1S2S3 called a pool region, there is one and only one blending recipe. The positions of these pools can be specified using the levelarm rule. The pool region is, in fact, the convex hull bounded by the quality values of the sources It is not possible to obtain qualities outside the convex hull. Therefore, these pools can be eliminated from the problem formulation. Next, we introduce a new discretization approach to enumerate bj values only in the attainable region. This approach sufficiently reduces the problem size and is described in the next section. 3.4. Implicit Enumeration of Discretized Qualities bj Values Using Flow Rate Proportion. This section introduces a new approach that implicitly enumerates the discretized pool qualities bj values within the convex hull of the attainable region of the pools. For those pools inside the convex hulls, the pool qualities relate to the source qualities by the mixing rule l l Xij ) ∑i)1 Xijaiq. constraints: bjq∑i)1 l Xij and let Divide both sides of this equation by ∑i)1
Figure 4. Discretization-based algorithm.
The discretization of bjq values is a preprocessing step, and the values are given as known inputs to the problem. Using this discretization strategy, xij can also be used as known parameters to the problem and the problem size can be further reduced. If we let Zj be the total flow through the pool, consistent with the flow propotion discretization, then Xij can be replaced by Xij ) xijZj. The formulation is ultimately given by n
k
k)1
xij )
Xij l
∑
M
l
jk
j)1
i)1
M
(3)
Available supply :
∑x Z eS ij j
[
M
]
∑ P ∑ Y - ∑ C ∑ (x Z )
Objective function is maximized :
i
ij j
j)1
for i ) 1, ... , l
i
j)1
Xij
n
i)1
Zj )
Mass balance on pools :
we have:
∑Y
jk
for j ) 1, ... , M
i)1
M
l
bjq )
∑x a
ij iq
(4)
i)1
where xij values are the flow rate proportion from source i to pool j with respect to the total flow through the pool and the convex hull condition:
Product demand :
∑Y
jk e Dk
for k ) 1, ... , n
j)1
Number of pools: fi ) {0, 1} for j ) 1, ... , M Zj e Uj, fj for j ) 1, ... , M M
∑f em j
j)1
l
l
∑
xij ) 1
(5)
i)1
The bjq points are indirectly enumerated within the convex hull by discretizing the values of xij within the interval [0, 1] while l xij ) 1. Therefore, the values of satisfying the condition ∑i)1 the discretized xij values can not be randomly selected. In each set of discretization, they must add up to 1. If there are l sources and t intervals chosen to discretize xij values, then Total number of discretized pools ) (t + 1)(t + 2)...(t + l - 1) (t + l - 1)! ) (6) 1 × 2 × 3...(l - 1) (l - 1) ! t! It is worth noting that the above formula resulting from the emplicit enumeration concept is that the number of discretized pools independent from the number of qualities Nq. This ia very attractive attribute for problems with large number of qualities.
Mixing rule for products : ckq
∑ j)1
l
Yjk g
∑Y b
jk jq for k ) 1,
... , n
j)1
Positive variables: Zij g 0, Yjk g 0 for i ) 1, ... , l j ) 1, ... , m k ) 1, ... , n
where bjq and xij values are known parameters defined in (2). 3.5. Discretization-Based Algorithm To Solve Large-Scale Pooling Problems. For multiple-quality pooling problems which are commonly encountered in practice, the implicit enumeration approach is preferred. The size of the MILP formulation (P′′) mainly depends on the number of sources l and the number of discretization intervals (t). A desired value of t may result in such a large MILP formulation that solvers can not handle or take too much calculation time. When it comes to this issue, the following algorithm is recommended (Figure 4). At first, the problem is solved with a “coarse” discretization scheme, in which t is low enough for the size of (P′′) not to exceed
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 1977
Figure 5. Results for problem Ben-Tal 5 using exhaustive and implicit enumeration approaches. Table 1. Descriptions and Calculation for the Pooling Problems (a) Pooling Problems with One Attribute problem name
Haverly 1
Haverly 2
Haverly 3
Foulds 2
Foulds 3
Foulds 4
Foulds 5
Ben-Tal 4
number of sources pools products attributes discretized pools variables constraints global optimum found optimum % error from optimum calculation time (discretization), s calculation time (global solver), s
3 1 2 1 4 22 16 400 400 0
