A Two-Stage Method for the Approximate Solution of General

Optimization of Water Network Synthesis for Single-Site and Continuous Processes: ... Reactive Scheduling by a Multiparametric Programming Rolling Hor...
1 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/IECR

A Two-Stage Method for the Approximate Solution of General Multiparametric Mixed-Integer Linear Programming Problems Martina Wittmann-Hohlbein and Efstratios N. Pistikopoulos* Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College London, SW7 2BY London, U.K. ABSTRACT: In this work, we focus on the approximate solution of multiparametric mixed-integer linear programming (mpMILP) problems involving uncertainty in the objective function coefficients and in the entries of the constraint matrices and vectors. A two-stage algorithmic procedure is proposed. In the first stage, the model is partially immunized against uncertainty using the worst-case oriented approach which leads to a partially robust mp-MILP model, whereas in the second stage explicit solutions of the robust model are derived by applying a suitable multiparametric programming algorithm for mp-MILP problems. Computational studies are presented, demonstrating that the proposed two-stage robust optimization/multiparametric programming procedure is computationally efficient and that it provides an upper bound on the overall solution of the general mp-MILP problem.

■

INTRODUCTION We consider the general multiparametric mixed integer linear programming (mp-MILP) problem (P), given by

partitioned into subregions, so-called critical regions, in which a specific solution of (P) remains optimal. In the open literature, the classes of (P) most widely studied are the RHS-mp-MILP problem with uncertainty in the righthand-side constraint vector, the OFC-mp-MILP problem with uncertain objective function coefficients, and the RIM-mp-MILP problem, which involves both right-hand-side and objective function uncertainty. A pioneering algorithm to solve RHS-mp-MILP problems was presented by Acevedo and Pistikopoulos1 which employs a branch and bound method. At each leaf node of the search tree the integer variables are fixed to their assigned values and the resulting multiparametric linear (mp-LP) problem is solved. Appropriate comparison procedures are employed to update the search tree. This method benefits from early fathoming of leaf nodes if the corresponding objective value is found to be larger than the current best upper bound. A branch and bound algorithm for the solution of the general mp-MILP problem (P) can be found in the work by Li and Ierapetritou.10 Alternatively, a number of algorithms have been developed that decompose the mp-MILP problem (P). The algorithms iterate between a master MILP or MINLP problem and a corresponding mp-LP subproblem. After each iteration integer and parametric cuts are introduced that exclude previously visited integer nodes and aim to find solutions that yield a lower objective value in a specific region. The algorithm for RHS-mpMILP problems outlined by Dua and Pistikopoulos11 was among the first to apply the ideas of decomposition from the single parametric case12 to the multiparametric framework. The algorithms proposed by Faísca et al.13 and Li and Ierapetritou5 further extend the theory and address the more general problem of type RIM-mp-MILP.

⎧ z(θ) := min((c + Hθ)T x + (d + Lθ)T y) ⎪ x ,y ⎪ s.t. ⎪ ⎪ (P)⎨ A(θ)x + E(θ)y ≤ b + Fθ ⎪ ⎪ x ∈ n , y ∈ {0, 1} p ⎪ ⎪ θ ∈ Θ := {θ ∈ q |θ min ≤ θ ≤ θ max , l = 1, ..., q}, ⎩ l l l

where θ denotes the vector of parameters and c ∈ n, H ∈ n × q, d ∈ p , L ∈ p × q, A(θ ) ∈ m × n, E(θ ) ∈ m × p, b ∈ m, and F ∈ m × q. To comply with the linear structure of (P) we further impose that A(θ) and E(θ) are affine mappings with respect to θ, that is, q

A(θ ) := AN +

∑ θlAl l=1

with A ∈  and Al ∈ m × n for all l. The matrix AN denotes the nominal part of A(θ). A similar expression holds for E(θ). In the following, we will denote by the lower case letter with the subscript i, for instance [ai], the column vector of entries related to the ith row of the corresponding matrix. Problem (P) has widespread application in many areas of process systems engineering such as process synthesis, planning, and scheduling.1−8 When confronted with data that have not yet revealed their true value, the introduction of uncertainty in the underlying mathematical model as in (P), while common, poses an additional challenge on its solution. Multiparametric programming is a powerful tool to account for the presence of uncertainty in a mathematical model during the optimization stage if at the time of decision making all data are known.9 The objective of multiparametric programming is to determine the effect of the parameter variation on the optimal solution of (P). The optimal solution of (P) is expressed as a function of the varying parameters while the parameter space is N

m×n

© 2012 American Chemical Society

Received: Revised: Accepted: Published: 8095

July 1, 2011 April 21, 2012 May 2, 2012 May 2, 2012 dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

constraints. Incorporating the worst-case oriented approach eq 1 into (P) results in the partially robust counterpart of the original mp-MILP problem. Each component of θ ∈ Θ takes values within a given interval. Let rl denote the range,

Despite major advances for certain classes of (P), the mpMILP problem with left-hand side (LHS) uncertainty remains the least studied. For the single parametric case two algorithms have been outlined.14 However, an extension of this solution method to address LHS-mp-MILP problems is nontrivial. LHSuncertainty introduces bilinear terms into the constraints with respect to both the optimization variables and the parameters. The algorithm by Li and Ierapetritou10 that addresses the general mp-MILP problem (P) uses discretization of the feasible parameter set to solve LHS-mp-LP subproblems. On the other hand, regarding (P) as nonlinear and nonconvex optimization problem, global optimization based solution techniques may also be applied.15 The aim in this work is to find solutions of (P) that are good approximations of the optimal solution and can be efficiently obtained. We follow the lines of Ben-Tal and Nemirovski16 and Lin et al.17 who postulate a worst-case-oriented robust counterpart of the linear problem and the mixed integer linear problem, respectively, with bounded type of uncertainty. Our approach differs as we are foremost interested in a partial immunization against the uncertainty. We propose to immunize the mp-MILP problem (P) solely against the variation in the constraint matrices A and E. This yields a partially robust counterpart that closely resembles the parametric nature of the original mp-MILP problem. The resulting partially robust model is a problem of type RIM-mp-MILP. The first step of this approach, denoted as a two-stage method for the approximate solution of general mp-MILP problems, involves recasting (P) as a partially robust model, whereas in the second step, the robust counterpart is solved employing suitable multiparametric algorithmic procedures. The paper is organized as follows. First, we formulate the partially robust counterpart of the general mp-MILP problem (P). We then outline the steps of a decomposition algorithm to solve RIM-mp-MILP problems, and introduce the proposed two-stage method for the approximate solution of (P). We also briefly discuss the key issues and difficulties in solving LHS-mpMILP problems. Finally, we present a number of numerical examples and discuss the quality of the approximate solution obtained by our proposed approach with respect to the exact solution and in comparison with the results obtained by the conventional robust optimization approach.

rl := (θlmax − θlmin)/2

and θNl the nominal value of θl,

θl N := θlmax − rl for l ∈ {1, ..., q}. It follows from eq 1 that a partially robust feasible solution (x,y ̅ )̅ of (P) satisfies q

[aiN ]T x ̅ + [eiN ]T y ̅ + = [aiN ]T x ̅ + +

q

q

[aiN ]T x ̅ + [eiN ]T y ̅ +

l=1

∑ γlEly ̅

∑ θlN([ail]T x ̅ + [eil]T y ̅ ) l=1

q

+

∑ rl|[ail]T x ̅ + [eil]T y ̅ | l=1

≤ bi + [ fi ]T θ

(3)

for i = 1, ..., m, θ ∈ Θ. Reversely, every feasible solution of eq 3 is also feasible for eq 1. By introducing q auxiliary variables and additionally 2q linear constraints for each constraint, eq 3 can be transformed into a system of linear inequalities. Consequently, the partially robust counterpart (RC) of the general mp-MILP problem (P) is given by ⎧ r(θ ) := min((c + Hθ )T x + (d + Lθ)T y) ⎪ x , ui , y ⎪ ⎪ s.t. q ⎪ ⎪[aiN ]T x + [eiN ]T y + ∑ θl N ([ail]T x + [eil]T y) ⎪ l=1 q ⎪ i T ⎪ + ∑ ru i = 1, ..., m l l ≤ bi + [ fi ] θ , ⎪ l=1 ⎪ (RC)⎨ ⎪− uli ≤ [ail]T x + [eil]T y ≤ uli ⎪ ⎪ l = 1, ..., q , i = 1, ..., m ⎪ x ∈ n , y ∈ {0, 1} p ⎪ ⎪ ui ∈ q , i = 1, ..., m ⎪ ⎪ q N N ⎪ θ ∈ Θ = {θ ∈  |θl − rl ≤ θl ≤ θl + rl , ⎪ l = 1, ..., q} ⎩

q

+

(2)

for every constraint i, i = 1, ..., m, and any θ ∈ Θ. Equation 2 represents those constraints of (P) that are most difficult to maintain. Equation 2 is equivalent to

Every instance of the constraints needs to be satisfied at the point (x,y ̅ ). ̅ For our needs it is sufficient to use a weaker concept of robust feasibility than the one presented above. We make the following adjustments: The pair (x,y (θ)) is called a ̅ )̅ := (x(θ),y ̅ ̅ partially robust feasible solution of (P) if

∑ γlAl x ̅ + EN y ̅

max ρl ([ail]T x ̅ + [eil]T y ̅ )

≤ bi + [ fi ]T θ

A(θ )x ̅ + E(θ )y ̅ ≤ b + Fθ

AN x ̅ +

∑

l = 1 −rl ≤ ρl ≤ rl

THE PARTIALLY ROBUST COUNTERPART OF THE GENERAL MP-MILP PROBLEM (P) In the context of robust optimization, see Ben-Tal et al.18 for an excellent textbook, the pair (x,y ̅ )̅ is called a robust feasible solution of (P) if it is feasible for every θ ∈ Θ, that is,

∀ γ ∈ Θ:

∑

q

■

∀ θ ∈ Θ:

max max γl([ail]T x ̅ + [eil]T y ̅ ) min l = 1 θl ≤ γl ≤ θl q N T [ei ] y ̅ + θl N ([ail]T x ̅ + [eil]T y ̅ ) l=1

∑

≤ b + Fθ

l=1

Problem (RC) is a multiparametric mixed integer linear problem of type RIM-mp-MILP. Every feasible solution of (RC) is a partially robust feasible solution of (P). The set of partially robust feasible solutions of (P) is contained in the set of feasible

(1)

for any θ ∈ Θ. Solutions of (P) satisfying eq 1 are immune against data uncertainty in the left-hand side matrices of the 8096

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

the model (P).17,20 The conventional robust counterpart of (P), denoted as (cvRC), can be obtained in the same way as outlined above if a reformulation of (P) is used. In (P), the objective function is replaced by the auxiliary variable t ∈  and the constraint

solutions of (P). Therefore, the optimal objective value r(θ) of (RC) is an upper bound on the optimal objective value z(θ) of (P) for every θ ∈ Θ. Note that the formulation of (RC) accounts for dependencies between uncertain data entries in each constraint. This includes the special case that the uncertain entries of the constraint matrices are modeled as independent which is often assumed in the open literature.16,17,19 Problem (RC) is infeasible if the set of constraints

(c + Hθ )T x + (d + Lθ )T y − t ≤ 0

is introduced. Additionally, new variables vi, i = 1, ..., m with a fixed value of 1 are introduced, and the constraints

q

[aiN ]T x + [eiN ]T y + T

∑

max

l = 1 θl

≤ bi + max [ fi ] γ ,

min

≤ γl ≤ θl

max

[ai(θ )]T x + [ei(θ )]T y ≤ bi + [ fi ]T θ

γl([ail]T x + [eil]T y)

i = 1, ..., m

γ ∈Θ

(5)

are rewritten as (4)

which can be expressed as a set of linear constraints, has no feasible solution (x,y) ∈ n × {0,1}p. Example 1. Consider problem (P1) with LHS-, RHS-, and OFC-uncertainty,10 given by

[ai(θ )]T x + [ei(θ )]T y − ([ fi ]T θ )vi ≤ bi

(6)

1 ≤ vi ≤ 1

(7)

for i = 1, ..., m. The robust counterpart (cvRC) of (P) is a deterministic mixed integer linear problem. The optimal objective value v of (cvRC) is an upper bound on r(θ) of (RC) and, consequently, on z(θ) of (P) for every θ ∈ Θ. Note that problem (cvRC) is infeasible if (P) is infeasible for any θ ∈ Θ. Example 1 (continued). The conventional robust counterpart (cvRC1) of (P1) is formulated as follows,

⎧ z(θ ) := min(( −3 + θ1)x1 − 8x 2 + 4y + 2y ) 1 2 ⎪ x ,y ⎪ ⎪ s.t. ⎪ x1 + x 2 ≤ 13 + θ2 ⎪ ⎪(5 + θ3)x1 − 4x 2 ≤ 20 ⎪ ⎪−8x1 + 22x 2 ≤ 121 ⎪ (P1)⎨−4x1 − x 2 ≤ −8 ⎪ ⎪ x1 − 10y1 ≤ 0 ⎪ ⎪ x 2 − 15y1 ≤ 0 ⎪ x ≥ 0, i = 1, 2 ⎪ i ⎪ y ∈ {0, 1}, j = 1, 2 ⎪ j ⎪ 0 ≤ θ ≤ 10, k = 1, 2, 3 ⎩ k

The RIM-mp-MILP problem (RC1) is the partially robust counterpart of (P1). It is formulated as

■

⎧ r(θ ) := min(( −3 + θ1)x1 − 8x 2 + 4y + 2y ) 1 2 ⎪ x ,y ⎪ ⎪ s.t. ⎪ x1 + x 2 ≤ 13 + θ2 ⎪ ⎪15x1 − 4x 2 ≤ 20 ⎪ ⎪−8x1 + 22x 2 ≤ 121 ⎪ (RC1)⎨−4x1 − x 2 ≤ −8 ⎪ ⎪ x1 − 10y1 ≤ 0 ⎪ x − 15y ≤ 0 1 ⎪ 2 ⎪ x ≥ 0, i = 1, 2 ⎪ i ⎪ y ∈ {0, 1}, j = 1, 2 ⎪ j ⎪ 0 ≤ θ ≤ 10, k = 1, 2, 3 ⎩ k

⎧ v := min(7x1 − 8x 2 + 4y + 2y ) 1 2 ⎪ x ,y ⎪ ⎪ s.t. ⎪ x1 + x 2 ≤ 13 ⎪ ⎪15x1 − 4x 2 ≤ 20 ⎪ ⎪−8x1 + 22x 2 ≤ 121 (cvRC1)⎨ ⎪−4x1 − x 2 ≤ −8 ⎪ x − 10y ≤ 0 1 ⎪ 1 ⎪ x − 15y ≤ 0 1 ⎪ 2 ⎪ xi ≥ 0, i = 1, 2 ⎪ ⎪ y ∈ {0, 1}, j = 1, 2 ⎩ j

THE SOLUTION OF RIM-MP-MILP PROBLEMS We outline a decomposition algorithm to solve the partially robust model (RC). The framework is based on the algorithm for the solution of RIM-mp-MILP problems.13 The algorithm iterates between a MINLP master problem and a mp-LP subproblem. The MINLP master problem. The master problem (M) is derived from (RC) by treating the parameter θ as an optimization variable. It is formulated as ⎧ min((c + Hθ )T x + (d + Lθ )T y) ⎪ x ,y,θ ⎪ (M )⎨ s.t. Ax + Ey ≤ b + Fθ ⎪ ⎪ x ∈ n , y ∈ {0, 1} p , θ ∈ Θ ⎩

The master problem contains the bilinear terms θTHx and θTLy in the objective function. Here, we use the global optimization solver BARON21 provided by the software package GAMS22 version 23.3 to solve the MINLP master problem. Let (xopt,yopt,θopt) be the optimal solution of (M), then yopt is input to obtain the following mp-LP subproblem.

Note that we do not require auxiliary variables in the model because of the non-negativity of x. We denote as the conventional robust counterpart of (P) the model whose solutions are immune against all data variation in 8097

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

Table 1. Steps of the Decomposition Algorithm for Problems of Type (RC) Step 0 Step 1 (mp-LP problem)

Step 2 (MINLP problem) Step 3

Initialize a critical region, CR = Θ, solve the MINLP problem (M) to global optimality, identify an initial integer solution, yopt of (M). If (M) is infeasible, go to Step 3. Fix y = yopt, solve (S) in CR to obtain rs(θ), a parametric upper bound of (RC) and the corresponding critical regions CRi, i = 1, ..., l. Identify closed polyhedral convex regions IRi′, I′ = 1, ..., k with ∪i′IRi′ = cl(CR\∪iCRi) where (S) is infeasible. Set rsi′ = ∞ if θ ∈ IRi′. Add rs(θ) and the optimal solution to the envelope of parametric prof iles. Update the set of critical regions of (RC) by incorporating the newly identified regions CRi and IRi′ for all i and i′. For each region CR (a) update the master problem (M) by (i) introducing integer and parametric cuts eq 12, eq 13, (ii) restricting the feasible parameter space to θ ∈ CR and (b) solve (M) and return to Step 1 if a new integer solution, yopt of (M), is found. The iteration terminates in a region where the MINLP master problem is infeasible. The algorithm stops if there are no more critical regions to explore.

The mp-LP Subproblem. By fixing y = yopt, (RC) results in the following RIM-mp-LP problem

when the additive term (d + Lθ)Tyopt in the objective function is neglected, have been identified. The Decomposition Algorithm for RIM-mp-MILP Problems. In the mp-MILP algorithm, the master problem (M) identifies an optimal integer solution of (RC). The subproblem (S) provides the corresponding continuous solution of (RC) at this integer node. Hence, with rs(θ), θ ∈ Θ, an upper bound on the optimal objective value of (RC) has been determined. Further iterations find integer nodes whose optimal profiles yield a tighter upper bound for any subregion of the parameter space. Between every master and subproblem iteration, the MINLP master problem is updated for each currently identified critical region. Integer cuts are introduced into the formulation of (M) in order to exclude previously visited integer solutions. The cuts are as follows,

⎧ r s(θ) := min ((c + Hθ)T x + (d + Lθ)T y opt ) x ⎪ ⎪ s.t. (S)⎨ ⎪ Ax ≤ b′ + Fθ ⎪ n ⎩x ∈  , θ ∈ Θ

with b′ := b − Eyopt. Consider the RIM-mp-LP problem in standard form,

⎧ min c(θ )T x x ⎪ ⎪ (LP)⎨ s.t. Ax = b(θ ) ⎪ ⎪ x ≥ 0, θ ∈ Θ ⎩

∑ yjk j∈J

with all uncertain coefficients depending linearly on θ. Let θ* ∈ Θ be a point such that the problem (LP) has an optimal solution. Then there exists an optimal basis, a nonsingular submatrix of the constraint matrix, for the deterministic problem at θ = θ*. We denote by p′ the index set of the basic variables corresponding to the columns of the optimal basis. For (LP), we define Ap′ to be the matrix whose columns are those of A related to the index p′ and c(θ)p′ to be the vector with entries of c(θ) related to p′. The entries of the optimal solution x(θ)opt of (LP) related to p′ are determined by −1 x(θ )opt p ′ = A p ′ b(θ )

∑ yjk

≤ |J k | − 1,

k = 1, ..., K

k

j∈L

(c + Hθ )T x + (d + Lθ )T y ≤ r s(θ )k − ε ,

(12)

k = 1, ..., K (13)

with ε ≥ 0 sufficiently small. To reduce the need to solve nonconvex optimization problems to global optimality, comparison procedures between objective function values outside the master problems are omitted. Instead, the envelope of parametric prof iles24 for each critical region is retained, storing all solutions that have been identified by the decomposition algorithm. Among all candidate profiles stored in the envelope, the optimal one is always included. The optimal solution of (RC) is then obtained through online function evaluations of the parametric profiles and direct value comparison of the upper bounds for a given realization of parameter values. The steps of the decomposition algorithm for RIM-mp-MILP problems are summarized in Table 1. Theorem 1. If problem (RC) is bounded from below for every θ∈Θ, the decomposition algorithm,Table 1, terminates finitely. Proof. Assume that Step 0 furnishes an integer solution yopt. For the corresponding subproblem (S), let G ⊆ Θ define the union of the LP optimality conditions eq 9 and eq 10 over all optimal basis indices p′ with respect to the standard formulation of type (LP) and the restriction θ ∈ Θ; that is, G:= ∪p′CRp′. The set G is convex which is shown in Theorem 2.2 in the paper by Gal.25 The number of optimal bases of (LP) is bounded by (nm), m ≤ n. Hence, the number of critical regions of (S) is finite. Each critical region CRp′ is a closed polyhedral convex set characterized by eq 11. Thus G is also a closed polyhedral

(8)

(9)

and the dual feasibility condition is described by ATA p−′Tc(θ )p ′ ≤ c(θ )

−

where K denotes the number of identified integer solutions in this critical region, the index sets Jk and Lk are defined by Jk := {j|ykj = 1} and Lk := {j|ykj = 0}, respectively, and |·| denotes the cardinality. Parametric cuts are also introduced to (M). These constraints, which are nonlinear and nonconvex, ensure that only integer nodes that are optimal for a certain realization of the parameters are considered. These cuts are given by

and those corresponding to the nonbasic variables are set to zero. The primal feasibility condition is then characterized by the inequality

A p−′1b(θ ) ≥ 0

k

(10)

The subregion of Θ for which Ap′ remains an optimal basis of (LP) is uniquely identified by the conditions eq 9 and eq 10, which are referred to as LP optimality conditions.10,23 The critical region is given by CR p ′ := {θ ∈ Θ|A p−′1b(θ ) ≥ 0, ATA p−′Tc(θ )p ′ ≤ c(θ )} (11) opt

The solution x(θ) , whose basic variables are defined by eq 8, is an optimal solution of (LP) in CRp′. The parameter value θ* is also contained in this region. The explicit solution of the mp-LP subproblem (S) is found when all optimal bases with respect to its corresponding standard form of type (LP), in which (S) can be transformed 8098

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

convex set. The rest of the region, Θrest = cl(Θ\G), where cl(·) denotes the closure of a set, is either empty or can be represented as union of closed polyhedral convex sets employing the procedure outlined by Bemporad et al.26 or by Dua et al.24 The generated sets may only overlap on their boundaries. The decomposition algorithm proceeds in each identified region, including the regions where (S) is infeasible as (RC) may only be infeasible for the current yopt. If from the solution of an updated master problem a new integer node is identified for a region, this region may again be subdivided into closed polyhedral convex sets as described above. As any MINLP master problem may contain at most 2p integer cuts before rendering infeasible, the decomposition algorithm terminates finitely. ◻ Remark. The solutions stored in the envelope of parametric profiles of (RC) generated by the decomposition algorithm are piecewise affine functions. In the special case that (RC) is an OFC-mp-MILP or RHS-mp-MILP problem, the objective function values related to the solutions from the envelope are also piecewise affine functions. The critical regions obtained with the decomposition algorithm are closed polyhedral convex sets, yet the union of the regions is not necessarily convex. In general, the optimal objective value r(θ) of (RC) is not continuous. Remark. Step 2 of the algorithm is performed on closed polyhedral convex regions. Thus it becomes necessary to represent the region where any mp-LP subproblem is infeasible as union of convex sets whose interiors are nonoverlapping. The characterization of the region of infeasibility for any subproblem influences the envelope of parametric profiles generated by the decomposition algorithm upon termination. For neighboring critical regions, the envelope may contain solutions that are identical such that for the corresponding mp-LP subproblem the same optimal basis is valid in both regions. Example 1 (continued). The partially robust counterpart (RC1) of (P1) is solved with the decomposition algorithm. Step 0: Initial MINLP problem in Θ = {θ ∈3|0≤θk ≤ 10, k = 1,2,.3}. Solve (M) which is equal to problem (RC1) with θ T being an additional optimization variable, yopt 0 = (1,1) . opt Step 1: mp-LP subproblem at y = y0 , θ ∈ Θ. Fix y = (1,1)T and solve the resulting mp-LP problem (S), formulated as

x 2opt =

CR 2 :=

2

( −3 + θ1)x1 − 8x 2 + 4y1 + 2y2 ≤

3

462 8392 θ1 − 149 149

θ ∈ CR1

The resulting MINLP master problem (M1) in CR1 is infeasible as is the master problem related to CR 2 . The algorithm terminates as there are no more critical regions to explore. In each critical region one solution has been identified, which is optimal for (RC1). The optimal objective function value r(θ) of (RC1) is given in Table 2 and the critical regions are depicted in Figure 1. Table 2. Optimal Objective Value and Integer Node of (RC1), Example 1 r(θ)

CR1 CR 2

((462/149)θ1 −(8392/149)) ((55/96)θ1 −(3973/96))

yopt (1,1)T (1,1)T

critical region 0 ≤ θ1 ≤ (65/11), 0 ≤ θ2 ≤ 10, 0 ≤ θ3 ≤ 10 (65/11) ≤ θ1 ≤ 10, 0 ≤ θ2 ≤ 10, 0 ≤ θ3 ≤ 10

Figure 1. Critical regions of (RC1), Example 1.

■

THE SOLUTION OF THE GENERAL MP-MILP PROBLEM (P) The proposed two-stage method is a combined robust optimization and multiparametric programming approach for the approximate solution of (P). As shown in Figure 2, it consists of an approximation stage in which (P) is immunized against

⎛ 462 1975 ⎞T ⎜ ⎟ , ⎝ 149 298 ⎠

{

1

y1 + y2 ≤ 1

The optimal solution x(θ)opt of (S) is given by

CR1 := 0 ≤ θ1 ≤

{ 1165 ≤ θ ≤ 10, 0 ≤ θ ≤ 10, 0 ≤ θ ≤ 10}

The critical regions identified entirely cover the original feasible parameter space, that is, CR1 ∪ CR 2 = Θ. Step 2: MINLP master problems in CR1 and CR 2 . For each critical region CR i , i = 1,2 solve an updated MINLP problem (Mi), which is problem (RC1) with (i) two additional constraints related to the integer and parametric cuts eq 12, eq 13 and (ii) the restriction θ ∈ CR i on the parameters and (iii) θ being treated as an optimization variable. As an example for CR1, this leads to the introduction of the constraints

⎧ r s(θ ) := min(( −3 + θ1)x1 − 8x 2 + 6) ⎪ x ⎪ s.t. ⎪ ⎪ x1 + x 2 ≤ 13 + θ2 ⎪ ⎪15x1 − 4x 2 ≤ 20 ⎪ ⎪−8x1 + 22x 2 ≤ 121 (S )⎨ ⎪−4x1 − x 2 ≤ −8 ⎪ ⎪ x1 ≤ 10 ⎪ x ≤ 15 ⎪ 2 ⎪ xi ≥ 0, i = 1, 2 ⎪ ⎪ 0 ≤ θ ≤ 10, k = 1, 2, 3 ⎩ k

x1opt =

⎛ 55 137 ⎞T ⎜ ⎟ , ⎝ 96 24 ⎠

65 , 0 ≤ θ2 ≤ 10, 0 ≤ θ3 ≤ 10 11

} 8099

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

jeopardizing the applicability of the decomposition algorithm for the solution of (P). The optimal solution of the LHS-mp-LP problem (S) may nevertheless be retrieved, provided all optimal bases and their regions of optimality are identified. In the work by Li and Ierapetritou,10 in which small-size mp-MILP problems with LHS-uncertainty are solved, optimal bases of mp-LP subproblems are obtained through exploration of a sufficiently fine grid of feasible parameter points. Function evaluation determines whether a grid point is already covered by a critical region in which case it can be discarded. Otherwise, if (S) is feasible and bounded at the point in question, a new optimal basis is identified. To the best of our knowledge, with the exception of the single parametric case,14 no algorithm is currently available that attempts to determine the optimal solution of the general mp-MILP problem (P) without employing discretization of the parameter space if A depends on θ. An alternative approach for the solution of (P) is to employ global optimization techniques adapted to the multiparametric framework.27 However, this approach may lead to a large number of critical regions and in an overall increase of the computational requirements. These difficulties are the driving motivation for the proposed two-stage method as an effective and approximate solution strategy for the general mp-MILP problem (P).

Figure 2. Flowchart of the two-stage method.

LHS-uncertainty yielding the partially robust counterpart (RC), and the solution stage in which the decomposition algorithm for the explicit solution of (RC) as outlined in Table 1 is employed. By construction, the optimal objective value r(θ) of (RC) is an upper bound on z(θ) of (P) for every θ ∈ Θ. On the other hand, the partially robust model (RC) is less conservative than the conventional robust counterpart (cvRC) of the general mpMILP problem (P). Hence,

■

COMPUTATIONAL STUDIES In this section, several general mp-MILP problems are presented to illustrate the potential of the two-stage method. Each example problem is solved with the proposed two-stage method. The decomposition algorithm for the solution of RIMmp-MILP problems has been implemented in Matlab with an interface to GAMS and is running on a Linux workstation (Dual 4 Core Intel Xeon processor, 1.6 GHz, 4 GB RAM) with Matlab being single threaded. The mp-LP subproblems are solved with the Parametric Optimization (POP) software.28 Examples 1, 2, and 3 are taken from Li and Ierapetritou.10 The third example was altered to include integer variables in the model. Examples 4 and 5 are variations of process synthesis problems.13,29 Examples 6, 7, and 8 are variations of subproblems of a model predictive control formulation. Examples 9 and 10 are scheduling formulations of batch processes that exhibit price, demand, and processing time uncertainty. Example 11 is a variation of Example 10. The example problems are described in more detail in the appendix. Table 3 provides an overview of the CPU-time requirements of the two-stage method for the example problems in comparison with discretization of the parameter space. Discretization is performed using 100 and 10, respectively, equidistant points for each parameter involved. The deterministic problems are then solved in Matlab with an interface to CPLEX30 version 9. For the two-stage method, the default value for the run time of the master problems is set to 100 s and ε = 1 is chosen for the offset value in the parametric cuts eq 13. A value of 0.01 in Examples 1 to 9 and a value of 0.1 in Examples 10 and 11, respectively, for the relative optimality criterion in the solution of MINLP problems is used. The columns of Table 3 denote the number of constraints with equality constraints given in brackets, the number of continuous variables, the number of binary variables, and the number of parameters. This is followed by a column for the CPU time for the discretization method with respect to 100

r max := sup r(θ ) ≤ v θ∈Θ max

where r and v are both guaranteed upper bounds on the optimal objective value z(θ) of (P). Remark. For the approximate solution of (P) furnished by the two-stage method, not every profile xji(θ), θ ∈ CRi stored in the envelope yields an objective function value rij(θ) that is necessarily less or equal than the conventional robust objective value v everywhere in the critical region, but ri(θ ) = min{rij(θ )} ≤ v j

holds for every θ ∈ CRi. This is illustrated in Example 4, Figure 9. Note that the decomposition algorithm presented in Table 1 is also suitable as a framework for the solution of the general mp-MILP problem (P). If uncertainty in the constraint matrices A and E is present, the constraints of the master problem (M) and those of the subproblem (S) contain bilinear terms. The solution of the LHS-mp-LP subproblem (S) is the bottleneck in the overall solution of (P). From the corresponding LP optimality conditions with respect to the standard formulation of (S) of type (LP), −1 x(θ )opt p ′ = A(θ )p ′ b(θ ) ≥ 0

(14)

and A(θ )T A(θ )−p ′T c(θ )p ′ ≤ c(θ )

(15)

for any optimal basis index p′, it follows that the optimal solution of (S) is a piecewise fractional polynomial function when all coefficients depend linearly on θ. Moreover, the optimal solution is not necessarily continuous. An illustrative example of an LHS-mp-LP problem whose optimal solution is not continuous can be found in the book by Gal.23 The critical regions need not be convex. Likewise, the characterization of the regions in which (S) is infeasible may pose a challenge, 8100

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

Table 3. Computational Requirements of the Two-Stage Method example

m

n

p

q

discretization (CPU h:min:sec)

two-stage method (CPU min:sec)

robust optimization (v)

1 3 4 5 6 7 8 9 10 11

8 9 5 14 26 50 74 254(15) 949(32) 950(32)

2 4 2 4 5 9 13 59 217 217

2 2 2 3 2 4 6 30 96 96

3 3 2 3 2 3 4 2 4 4

7:08:16 (0:00:01) 7:19:30 (0:00:01) 0:00:15 (0:00:01) 8:28:30 (0:00:01) 0:00:08 (0:00:01) 7:38:37 (0:00:01) − (0:00:06) 0:02:09 (0:00:01) − (15:20:21) − (15:12:00)

0:02 0:03 0:03 0:03 0:05 0:23 1:03 0:26 8:10 34:00

−35.7 ∞ 184.6 −27 ∞ ∞ ∞ ∞ −1008.7 −2869.7

Table 4. Number of Iterations, Number of Critical Regions, and Multiplicity of the Envelope of the Two-Stage Method decomposition algorithm

Example 1 2 3 4 5 6 7 8 9 10 11

no. MINLP

no. mpLP

no. CR

mult.

5 2 15 3

1 1 4 2

4 1 9 1

1 1 2 2

constraint is introduced into the model enforcing that at any point the objective value is less than or equal to the optimal objective value v of the conventional robust counterpart, the modified two-stage method remains less conservative than robust optimization in regions with a finite solution. The computational requirements in terms of the number of iterations between master and subproblems of the two-stage method are given in Table 4. The multiplicity of the envelope denotes the maximum of all numbers of solutions stored in the envelope of parametric profiles for any of the critical region. For the small-size problems, Examples 1−4, the optimal solution is also derived using the decomposition algorithm.

two-stage method no. MINLP 3 1 (MILP) 10 3 6 13 73 162 10 4 16

no. mpLP 1 3 2 3 2 10 25 3 1 6

no. CR

mult.

2 1 6 1 3 5 35 69 6 3 10

1 1 2 2 2 1 1 1 2 1 3

■

CONCLUDING REMARKS In this paper, we apply suitable robust optimization and multiparametric programming techniques to derive close to optimal solutions of the mixed integer linear programming problem contaminated with uncertainty. We propose to immunize the general mp-MILP problem against uncertainty in the coefficients of the constraint matrices. We propagate a partially robust counterpart which is of type RIM-mp-MILP, involving varying parameters in the objective function coefficients and right-hand side constraint vector. The optimal partially robust solution serves as solution estimate for the optimal solution of the general mp-MILP problem providing an upper bound on the optimal objective value. It is computationally expensive to solve the general mpMILP problem (P). For the small-size examples we have observed that the two-stage method for the approximate solution of (P) requires less iterations than the decomposition algorithm for the exact solution. Considerable effort may also be needed for the exact solution of mp-LP subproblems of (P) when LHS-uncertainty is present such as discretization of the parameter space, global optimization techniques or enumerating of all bases and subsequent testing for optimality, and for identifying regions where mp-LP subproblems are infeasible. Hence, the decomposition algorithm proves impractical for medium to large-size problems. These issues are overcome by the two-stage method. A second advantage of the two-stage method is the generation of piecewise affine optimal solution estimates and of convex critical regions, which significantly simplifies the characterization of the parameter space. Beneficial of the two-stage method is furthermore the low degree of conservatism of this approach compared to the conventional worst-case oriented robust optimization approach in which the robust model is a deterministic mixed integer linear model. The latter is motivation to apply the two-stage method in pro-active scheduling of short-term batch processes5,17,20 under various types of uncertainty at the optimization stage but that

and 10 grid points, respectively, per parameter. The latter is given in brackets. The second to last column stores the CPU time for the two-stage method. The last column presents the value v of the conventional robust counterpart (cvRC) which is a guaranteed upper bound on the optimal objective value of each example problem. Clearly, an adequate discretization scheme for any problem with a sufficiently large number of parameters is prohibitive. Using 100 grid points per parameter for the test problems with three parameters, the discretization method is outperformed by the two-stage method in terms of computational time. The same scheme is not applicable to problems with four parameters involved as it well exceeds running times of one day. This is indicated by ″−″ in Table 3. The applicability of the two-stage method depends on the success of solving the MINLP master problems to global optimality. In the larger sized Examples 10 and 11, the solver requires a significant amount of time to find the global optimum of the master problems involved. If a globally optimal solution of any master problem is not found within the specified time limit, the current best integer node furnished by the solver is used. In this case, the modified two-stage method identfies a feasible, but not necessarily the optimal partially robust solution of the problem. In other words, an approximate solution is generated for which the corresponding objective function value or the minimum of the objective values if multiple solutions are stored in the envelope of critical regions, respectively, may not be less or equal than the optimal objective value of the conventional robust counterpart for every parameter realization. This is observed in Example 11 around the point of origin and is depicted in Figure 11. If an auxiliary 8101

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

are known to reveal their true values at the time of decision making. Ongoing and future work includes the extension of the two-stage method in embedding different robust models from the open literature31 and in employing suitable relaxations of bilinear terms.32,33 Another application is in robust hybrid control.34,35 As in this work the focus is on the multiparametric case, a comparison of the two-stage method with the algorithms presented by Mitsos and Barton14 for the solution of LHS-p-MILP problems was not performed. For the special case that the partially robust counterpart is a OFC- or a RHS-mp-MILP problem, the performance of the two-stage method employing a branch- and bound-based algorithm1,10 compared to the propagated decomposition algorithm is also a point of general interest and is pursued further. The global solution of mp-MILP problems with LHS-uncertainty is not discussed in depth in this paper. It is a problem of its own right and is subject of forthcoming work.

■

Table 5. Optimal Objective Value and Integer Node of (P1), Example 1 z(θ) CR1 CR2 CR3 CR4

yopt

critical region

(−3 + θ1) (11/2 + 11/15θ2) − (32/15)θ2 − 54 2(−2223 + 231θ1 −242θ3) /(11θ3 + 39) + 6 10θ1 − (1068/11)

(1,1)

(−3973/96) + (55/96)θ1

(1,1)T

T

(1,1)T (1,1)T

0 ≤ θ1 ≤ (65/11), 0 ≤ θ2 ≤ (135/22), 0 ≤ θ3 ≤ ((675 − 78θ2)/(165 + 22θ2)) 0 ≤ θ1 ≤ (65/11), 0 ≤ θ2 ≤ 10, (675 − 78θ2)/(165 + 22θ2) ≤ θ3, (36/55) ≤ θ3 ≤ 10 0 ≤ θ1 ≤ (65/11), (135/22) ≤ θ2 ≤ 10, 0 ≤ θ3 ≤ (36/55) (65/11) ≤ θ1 ≤ 10, 0 ≤ θ2 ≤ 10, 0 ≤ θ3 ≤ 10

APPENDIX

mp-MILP Example Problems

Example 1 (continued). The optimal solution of (P1) obtained with the decomposition algorithm is attained at the integer node yopt = (1,1)T with the corresponding continuous solution given by x1opt =

T ⎛ 11 11 15 4 ⎞ ⎜ + θ2 , + θ2⎟ ⎝2 15 2 15 ⎠

65 135 , 0 ≤ θ2 ≤ , 11 22 675 − 78θ2 ⎫ ⎬ 0 ≤ θ3 ≤ 165 + 22θ2 ⎭

{

CR1 := 0 ≤ θ1 ≤

x 2opt

Figure 3. Critical regions of (P1), Example 1.

v = −35.7

is attained at the point (xopt,yopt)T = (0.5,5.7,1,1)T which is depicted in Figure 5. Example 2. We consider the mp-MILP problem which has only LHS-uncertainty,

T ⎛ 121(θ3 + 5) ⎞ 1 = ⎜462, 80 + ⎟ 11θ3 + 39 ⎝ 2 ⎠

{

CR 2 := 0 ≤ θ1 ≤

65 , 11

675 − 78θ2 ≤ θ3 , 165 + 22θ2

0 ≤ θ2 ≤ 10,

⎧ z(θ ) := min( −2x1 − x 2 + y + y ) 1 2 ⎪ x ,y ⎪ ⎪ s.t. ⎪ x1 + (3 + θ1)x 2 + y ≤ 9 1 ⎪ ⎪ ⎪(2 + θ2)x1 + x 2 − y2 ≤ 8 (P 2)⎨ ⎪ x1 − y1 + y2 ≤ 4 ⎪ ⎪ xi ≥ 0, i = 1, 2 ⎪ ⎪ yj ∈ {0, 1}, j = 1, 2 ⎪ ⎪ 0 ≤ θ ≤ 10, k = 1, 2 ⎩ k

36 ≤ θ3 ≤ 10 55

}

⎛ 201 ⎞⎟T x3opt = ⎜10, ⎝ 22 ⎠

{

CR3 := 0 ≤ θ1 ≤ 0 ≤ θ3 ≤

x4opt =

36 ⎫ ⎬ 55 ⎭

65 , 11

135 ≤ θ2 ≤ 10, 22

⎛ 55 137 ⎞T ⎜ ⎟ , ⎝ 96 24 ⎠

CR 4 :=

{ 1165 ≤ θ ≤ 10, 1

0 ≤ θ2 ≤ 10,

The point y = (0,0)T is the only optimal integer solution of (P2), which is valid in the original feasible parameter set CR1 = {θ ∈ 2|0 ≤ θi ≤ 10, I = 1, 2}. The corresponding optimal noninteger solution and the objective function value are given by

⎫ 0 ≤ θ3 ≤ 10⎬ ⎭

The optimal objective value and the critical regions of (P1) are depicted in Table 5 and Figure 3, respectively. Figures 4 and 5, depicting the optimal objective value z(θ) and optimal partially robust objective value r(θ) of (P1) at θ3 = 0, θ3 = 10, and the nominal value θ3 = 5, clearly show that r(θ) is an upper bound on the optimal objective value z(θ). The optimal objective value v of the conventional robust counterpart (cvRC1) of (P1),

x opt =

1 (15 + 8θ1, 10 + 9θ2)T 5 + 2θ1 + 3θ2 + θ1θ2

and z(θ) = 8102

− 40 − 16θ1 − 9θ2 5 + 2θ1 + 3θ2 + θ1θ2 dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

Figure 4. Optimal objective value z(θ) (surface plot) and optimal partially robust objective value r(θ) (grid plot) of (P1) at θ3 = 0 and θ3 = 10, Example 1.

Figure 6. Optimal objective value z(θ) (surface plot) and optimal partially robust objective value r(θ) (grid plot) of (P2), Example 2.

Figure 5. Optimal objective value z(θ) (surface plot), optimal partially robust objective value r(θ) (red grid plot), and optimal conventional robust objective value v (green grid plot) of (P1) at θ3 = 5, Example 1.

The two-stage method applied to (P3) required a significant lower number of iterations. We depict the corresponding critical regions and the envelope of partially robust objective function values, projected onto the θ1−θ2 space, of problem (P3) in Figure 8. The envelope of partially robust objective values along with the corresponding integer solutions are given in Table 7. Example 4. This is a small process synthesis problem with demand, cost, and conversion rate uncertainty at the time of optimization and the objective is to minimize the costs. The model is given by

respectively. Here, the partially robust and the conventional robust counterpart of (P2) agree. The optimal partially robust solution of (P2) obtained with the two-stage method is attained at (xopt,yopt)T = (0.61,0.64,0,0)T with the corresponding objective value estimate r(θ ) = −1.87,

CR1 = {θ ∈ 2 |0 ≤ θi ≤ 10, i = 1, 2}

In Figure 6, we depict the optimal objective value z(θ) and optimal partially robust objective value r(θ) of (P2). Example 3. This example features LHS-, RHS-, and OFCuncertainty as follows

⎧ z(θ ) := min ((6.4 + 0.25θ )x + 6x x ,y 1 1 2 ⎪ ⎪ + (7.5 + 0.3θ1)y + 5.5y ) 1 2 ⎪ ⎪ s.t. ⎪ ⎪ 0.8x1 + (0.67 + 0.015θ1)x 2 ≥ 10 + θ2 ⎪ x ≤ 40y ⎪ 1 1 (P 4)⎨ ⎪ x 2 ≤ 40y2 ⎪ ⎪ xi ≥ 0, i = 1, 2 ⎪ ⎪ yi ∈ {0, 1}, i = 1, 2 ⎪ 0 ≤ θ ≤ 20 1 ⎪ ⎪ 0 ≤ θ ≤ 10 ⎩ 2

⎧ z(θ ) := min(θx1 + x 2 + y ) 1 ⎪ x ,y ⎪ ⎪ s.t. ⎪−x1 + x 2 + x3 = θ2 + 2y 1 ⎪ ⎪ ⎪ x1 + θ3x 2 + x4 = 1 + θ1y2 (P 3)⎨ ⎪ y2 − y1 ≤ 0 ⎪ ⎪ xi ≥ 0, i = 1, ..., 4 ⎪ ⎪ yj ∈ {0, 1}, j = 1, 2 ⎪ ⎪−5 ≤ θ ≤ 5, k = 1, ..., 3 ⎩ k

Applying the decomposition algorithm, the solutions stored in the envelope of parametric profiles are

The envelope of parametric profiles of (P3) obtained with the decomposition algorithm is presented in Table 6. The critical regions of (P3) are shown in Figure 7.

x1 = (12.5 + 1.25θ2 , 0)T , 8103

y1 = (1, 0)T

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

Table 6. Envelope of Parametric Profiles of (P3), Example 3 z(θ)

y

−∞ θ1 −θ1θ2 0 (−(θ1(θ2θ3 + 1) + (θ2 + 1)))/(θ3−1) 1 (−(θ1(θ2θ3 + 1) + (θ2 + 1)))/(θ3 − 1) θ1(−θ2 − 2) + 1 (−(θ1(θ2θ3 + 1) + (θ2 + 1)))/(θ3 − 1) [(−θ1(θ2θ3 + θ2 + θ3 + 1) −(θ2 + θ3 + 3))/(θ3 − 1)] + 1 1 θ1(−θ2 − 2) + 1 θ1 + 1

CR0 CR1 CR2 CR3 CR4 CR5 CR6 CR7 CR8 CR9

critical region

(0,0)T (0,0)T (0,0)T (0,0)T (1,1)T (0,0)T (1,1)T (0,0)T (1,1)T (1,1)T (1,1)T (1,0)T

−5 ≤ θ1 ≤ 0, −5 ≤ θ2 ≤ 5, θ1θ3 < −1 −5 ≤ θ1 ≤ 0, −1 ≤ θ2 ≤ 5, −5 ≤ θ3, θ1θ3 ≥ −1 0 ≤θ1 ≤ 5, −1 ≤ θ2 ≤ 0, −5≤θ3 ≤ 5 0 ≤ θ1 ≤ 5, 0 ≤ θ2 ≤ 5, −5 ≤θ3 ≤ 5 0 ≤ θ1 ≤ 5, −2 ≤ θ2 ≤ −1, 1 < θ3 ≤ 5 0 ≤ θ1 ≤ 5, −5 ≤ θ2 ≤ −2, −θ1 − 3 ≤ θ2, 1 < θ3 ≤ 5 0 ≤ θ1 ≤ 2, −5 ≤ θ2 ≤ −θ1 −3, 1 < θ3 ≤ 5 0 ≤ θ1 ≤ 5, −2 ≤ θ2 < −1, −5 ≤ θ3 ≤ 1 0 ≤ θ1 ≤ 5, −5 ≤ θ2 ≤ −2, −θ1 − 3 ≤ θ2, −5 ≤ θ3 ≤ 1 −5 ≤ θ1 ≤ 0, −5 ≤ θ2 ≤ −3, −5 ≤ θ3, θ1θ3 ≥ −1

and ⎛ ⎞T 10 + θ2 x 2 = ⎜0, ⎟ , ⎝ 0.67 + 0.025θ1 ⎠

y 2 = (0, 1)T

which are valid in the initial set of feasible parameters. The corresponding objective values, referred to as z1(θ) and z2(θ), respectively, are given in Table 8 and depicted in Figure 9. The parametric profiles obtained with the two-stage method are x1 = (12.5 + 1.25θ2 , 0)T ,

y1 = (1, 0)T

and x 2 = (0, 14.92 + 1.49θ2)T ,

y 2 = (0, 1)T

respectively, which are valid in the initial set of feasible parameters. The envelope of parametric profiles is given in Table 9 and the corresponding partially robust objective function values, represented by r1(θ) and r2(θ), respectively, are depicted in Figure 9.

Figure 7. Critical regions of (P3), Example 3.

Figure 8. (a) Critical regions and (b) objective values from the envelope of parametric profiles of (P3) with two-stage method, Example 3.

Table 7. Envelope of Parametric Profiles of (P3) with Two-Stage Method, Example 3

CR1 CR2 CR3 CR4 CR5 CR6

r(θ)

y

critical region

θ1 −θ1θ2 1 0 1 θ1(−θ2−2) + 1 θ1 + 1

(0,0)T (0,0)T (1,0)T (0,0)T (1,0)T (1,0)T (1,0)T

−5 ≤ θ1 ≤ 0, −1 ≤ θ2 ≤ 5, −5 ≤ θ3 ≤ 5

8104

0 ≤ θ1 ≤ 5, −1 ≤ θ2 ≤ 0, −5 ≤ θ3 ≤ 5 0 ≤ θ1 ≤ 5, 0 ≤ θ2 ≤ 5, −5 ≤ θ3 ≤ 5 0 ≤ θ1 ≤ 5, −2 ≤ θ2 < −1, −5 ≤ θ3 ≤ 5 0 ≤ θ1 ≤ 5, 3 ≤ θ2 ≤ −2, −5 ≤ θ3 ≤ 5 −5 ≤ θ1 ≤ 0, −3 ≤ θ2 < −1, −5 ≤ θ3 ≤ 5 dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

Table 8. Envelope of Parametric Profiles of (P4), Example 4 z(θ) CR1

y

Table 9. Envelope of Parametric Profiles of (P4) with TwoStage Method, Example 4

critical region

3.425θ1 + 8θ2 + 0.3125θ1θ2 + 87.5 (1,0)T 0 ≤ θ1 ≤ 20, 0 ≤ θ2 ≤ 10 (6(10 + θ2)/(0.67 + 0.015θ1) + 5.5 (0,1)T

r(θ) CR1

The region of intersection between z1(θ) and z2(θ) is marked as red contour line and that between r1(θ) and r2(θ) as blue contour line in the graphs of Figure 9. Each intersection splits the set of feasible parameter points into two sub-sets where the lower of the two values equals the optimal objective value of and the optimal partially robust objective value of (P4), respectively. The optimal solution of the conventional robust counterpart of (P4) is attained at the integer node y = (0,1)T and yields the optimal value

y

critical region

3.425θ1 + 8θ2 + 0.3125θ1θ2 + 87.5 (1,0)T 0 ≤ θ1 ≤ 20, 0 ≤ θ2 ≤ 10 8.95θ2 + 95.05 (0,1)T

Table 10. Envelope of Parametric Profiles of (P5) with TwoStage Method, Example 5 continuous solution x = Uθ + v CR1

v = 184.6

critical region

integer solution y

0 ≤ θ1 ≤ 5 2.6 ≤ θ2 ≤ 10 0.75 ≤ θ3 ≤ 0.95

(1,1,0)T

which is also depicted in Figure 9. Example 5. The process synthesis problem involves RHS-, OFC- and LHS-uncertainty due to the presence of cost, demand and conversion rate uncertainty. It is formulated as follows,

(1,0,1)T

⎧ z(θ) := min (− 18(0.82x + θ x ) − 10(0.4x ) 2 3 3 1 ⎪ x ,y ⎪ ⎪ + 2.5x1 + (4 + θ1)x 2 + 5.5x3 + 10y1 + 15y2 ⎪ + 20y ) 3 ⎪ ⎪ s.t. ⎪ ⎪ 2 ≤ 0.82x 2 + θ3x3 ≤ 5 + θ2 ⎪ 0.5x − x − x + x = 0 1 2 3 4 ⎪ ⎪ 0.4x1 ≤ 5 + θ2 ⎪ ⎪ x1 ≤ 16y1 ⎪ ⎪ (P 5)⎨ x 2 ≤ 30y2 ⎪ ⎪ x3 ≤ 30y3 ⎪ ⎪1 ≤ y2 + y3 ⎪ ⎪ x4 ≤ 14 ⎪ x ≥ 0, i = 1, ..., 4 ⎪ i ⎪ y ∈ {0, 1}, i = 1, 2, 3 ⎪i ⎪0 ≤ θ ≤ 5 1 ⎪ ⎪ 0 ≤ θ2 ≤ 10 ⎪ ⎪ 0.75 ≤ θ ≤ 0.95. ⎩ 3

CR2

0 ≤ θ1 ≤ 5 1.56 ≤ θ2 ≤ 2.6 0.75 ≤ θ3 ≤ 0.95

(1,1,0)T

(1,0,1)T

CR3

0 ≤ θ1 ≤ 5 0 ≤ θ2 ≤ 1.56 0.75 ≤ θ3 ≤ 0.95

(1,1,0)T

(1,0,1)T

U

v

0 0 0 0 0 0 0 0

0 1.22 0 1.22 0 0 1.05 1.05

0 0 0 0 0 0 0 0

16 6.09 0 −1.91 16 0 5.26 −2.74

0 0 0 0 0 0 0 0

0 1.22 0 1.22 2.1 0 1.05 0

0 0 0 0 0 0 0 0

16 6.09 0 −1.91 10.53 0 5.26 0

0 0 0 0 0 0 0 0

2.44 1.22 0 0 2.1 0 1.05 0

0 0 0 0 0 0 0 0

12.19 6.09 0 0 10.53 0 5.26 0

Applying the two-stage method, the envelope of partially robust parametric profiles of (P5) is given in Table 10. The optimal objective value v of the conventional robust counterpart of (P5) is given by

Figure 9. (a) Objective values from the envelope of parametric profiles of (P4) with decomposition algorithm and two-stage method, represented by z1(θ),z2(θ), and r1(θ),r2(θ), respectively, and (b) together with the optimal conventional robust objective value v, Example 4. 8105

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

Figure 10. (a) Critical regions of Example 6, and (b) critical regions of Example 9 with two-stage method.

Figure 11. (a) Critical regions, and (b) objective values from the envelope of parametric profiles (surface plot) with modified two-stage method and optimal conventional robust objective value (grid plot), Example 11.

v = −27

Example 11 is a variation of Example 10 in which the objective function has been altered. The parameters θ3 and θ4 affect solely the entries of the constraint matrices A and E. The approximate solution obtained with the modified two-stage method is thus independent of θ3 and θ4. For Example 11, the critical regions and the envelope of objective function values are illustrated in Figure 11. The MAT-files of all example problems, apart from Example 2, and their partially robust counterparts, respectively, are available at http://sites.google.com/site/imperialercmobile/currentresearch/mp-milp-collection. For each of the problems, the envelope of parametric profiles as obtained with the two-stage method is also provided on the website.

and is attained at the integer node y = (1,0,1)T. Examples 6−8. The examples are variations of subproblems related to explicit model predictive control of a system with piecewise affine dynamics which is described in the technical report of Borrelli et al.36 The subproblems arise from a dynamic programming formulation. The parameters represent the states of the system and future control actions. For the example problems, the original quadratic objective functions have been replaced by randomly generated linear objective functions whose coefficients depend linearly on the parameters. The examples are RIM-mpMILP problems and, thus, are in agreement with their partially robust counterparts. The critical regions of Example 6 are depicted in Figure 10a. For these example problems the characterization of the regions where the partially robust model is infeasible, which is a consequence of the requirement of the decomposition algorithm to operate on closed polyhedral convex subsets of the parameter space, is worth mentioning. Upon termination of the two-stage method Example 6 features six regions of infeasibility, Example 7 features 28 regions of infeasibility, and Example 8 features 68 regions of infeasibility. Examples 9−11. Examples 9 and 10 are both short-term batch process scheduling models using a continuous time formulation as described in Ierapetritou and Floudas.37 We have introduced price, demand, and processing time uncertainty. The latter introduces LHS-uncertainty into the models. The partially robust model of Example 9 is modified such that the worst-case with respect to the batch-size related processing time is implemented but not with respect to the fixed processing time. Hence, the partially robust model features uncertain entries in the constraint matrix E which can be equally addressed by the decomposition algorithm, Table 1. The critical regions of Example 9 are depicted in Figure 10b.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support from EPSRC (EP/G059071/1, EP/I014640/1) and from the European Research Council (MOBILE, ERC Advanced Grant, No 226462) is gratefully acknowledged. We would like to thank the anonymous reviewers for their insightful and helpful comments.

■

REFERENCES

(1) Acevedo, J.; Pistikopoulos, E. N. A multiparametric programming approach for linear process engineering problems under uncertainty. Ind. Eng. Chem. Res. 1997, 36, 717−728.

8106

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107

Industrial & Engineering Chemistry Research

Article

(28) ParOS, Parametric Optimization Solutions Ltd (http://www. parostech.co.uk, accessed 2011). (29) Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods of Chemical Process Design; Prentice Hall: NJ, 1997. (30) IBM ILOG CPLEX, (www.ibm.com/software/integration/ optimization/cplex-optimizer, accessed 2011). (31) Li, Z.; Ding, R.; Floudas, C. A. A comparative theoretical and computational study on robust counterpart optimization: I. Robust linear optimization and robust mixed integer linear optimization. Ind. Eng. Chem. Res. 2011, 50, 10567−10603. (32) McCormick, G. P. Computability of global solutions to factorable nonconvex programs: Part I. Convex underestimating problems. Math. Program. 1976, 10, 147−175. (33) Gounaris, C. E.; Misener, R.; Floudas, C. A. Computational comparison of piecewise linear relaxations for pooling problems. Ind. Eng. Chem. Res. 2009, 48, 5742−5766. (34) Goulart, P. J.; Kerrigan, E. C. Input-to-state stability of robust receding horizon control with an expected value cost. Automatica 2008, 44, 1171−1174. (35) Goulart, P. J.; Kerrigan, E. C.; Ralph, D. Efficient robust optimization for robust control with constraints. Math. Program. 2008, 114, 115−147. (36) Borrelli, F.; Baotic, M.; Bemporad, A.; Morari, M. Constrained Optimal Control of Discrete-Time Linear Hybrid Systems; Technical Report; Automatic Control Laboratory, ETH: Zurich, 2003. (37) Ierapetritou, M. G.; Floudas, C. A. Effective continuous-time formulation for short-term scheduling. 1. Multipurpose batch processes. Ind. Eng. Chem. Res. 1998, 37, 4341−4359.

(2) Ryu, J.-H.; Dua, V.; Pistikopoulos, E. N. Proactive scheduling under uncertainty: A Parametric optimization approach. Ind. Eng. Chem. Res. 2007, 46, 8044−8049. (3) Ryu, J.-H.; Pistikopoulos, E. N. A novel approach to scheduling of zero-wait batch processes under processing time variations. Comput. Chem. Eng. 2007, 31, 101−106. (4) Li, Z.; Ierapetritou, M. G. Process scheduling under uncertainty: Review and challenges. Comput. Chem. Eng. 2008, 32, 715−727. (5) Li, Z.; Ierapetritou, M. G. Process scheduling under uncertainty using multiparametric programming. AIChE J. 2007, 53, 3183−3203. (6) Verderame, P. M.; Elia, J. A.; Li, J.; Floudas, C. A. Planning and scheduling under uncertainty: A review across multiple sectors. Ind. Eng. Chem. Res. 2010, 49, 3993−4017. (7) Jenkins, L. Parametric mixed integer programmingan application to solid-waste management. Manage.. Sci. 1982, 28, 1270−1284. (8) Pertsinidis, A. On the parametric optimization of mathematical programs with binary variables and its application in the chemical engineering process synthesis. Ph,D. Thesis, Carnegie Mellon University, Pittsburgh, PA, 1991. (9) Wallace, S. W. Decision making under uncertainty: Is sensitivity analysis of any use? Operat. Res. 2000, 48, 20−25. (10) Li, Z.; Ierapetritou, M. G. A new methodology for the general multiparametric mixed-integer linear programming (MILP) problems. Ind. Eng. Chem. Res. 2007, 46, 5141−5151. (11) Dua, V.; Pistikopoulos, E. N. An algorithm for the solution of multiparametric mixed integer linear programming problems. Ann. Operat. Res. 2000, 99, 123−139. (12) Pertsinidis, A.; Grossmann, I. E.; McRae, G. J. Parametric optimization of MILP programs and a framework for the parametric optimization of MINLPs. Comput. Chem. Eng. 1998, 22, 205−212. (13) Faísca, N. P.; Kosmidis, V. D.; Rustem, B.; Pistikopoulos, E. N. Global optimization of multi-parametric MILP problems. J. Global Opt. 2009, 45, 131−151. (14) Mitsos, A.; Barton, P. I. Parametric mixed-integer 0−1 linear programming: The general case for a single parameter. Eur. J. Operat. Res. 194, 663−686. (15) Dua, V.; Pistikopoulos, E. N. Algorithms for the solution of multiparametric mixed-integer nonlinear optimization problems. Ind. Eng. Chem. Res. 1999, 38, 3976−3987. (16) Ben-Tal, A.; Nemirovski, A. Robust solutions of linear programming problems contaminated with uncertain data. Math. Program. 2000, 88, 411−424. (17) Lin, X.; Janak, S.; Floudas, C. A new robust optimization approach for scheduling under uncertainty: I. Bounded uncertainty. Comput. Chem. Eng. 2004, 28, 1069−1085. (18) Ben-Tal, A.; Ghaoui, L. E.; Nemirovski, A. Robust Optimization; Princeton University Press: Princeton, NJ, 2009. (19) Bertsimas, D.; Sim, M. Robust discrete optimization and network flows. Math. Program. 2003, 98, 49−71. (20) Li, Z.; Ierapetritou, M. G. Robust optimization for process scheduling under uncertainty. Ind. Eng. Chem. Res. 2008, 47, 4148−4157. (21) Sahinidis, N.; Tawarmalani, M. Gams/Baron 9.0.2: Global Optimization of Mixed-Integer Nonlinear Programs, 2009. (22) GAMS, GAMS Development Corporation (http://www.gams.com, accessed 2011). (23) Gal, T. Postoptimal Analyses, Parametric Programming, and Related Topics: Degeneracy, Multicriteria Decision Making, Redundancy; W. de Gruyter: Berlin, 1995. (24) Dua, V.; Bozinis, A.; Pistikopoulos, E. N. A multiparametric programming approach for mixed-integer quadratic engineering problems. Comput. Chem. Eng. 2002, 26, 715−733. (25) Gal, T. Rim multiparametric linear programming. Manage. Sci. 1975, 21, 567−575. (26) Bemporad, A.; Morari, M.; Dua, V.; Pistikopoulos, E. The explicit linear quadratic regulator for constrained systems. (27) Dua, V.; Papalexandri, K. P.; Pistikopoulos, E. N. Global optimization issues in multiparametric continuous and mixed-integer optimization problems. J. Global Opt. 2004, 30, 59−89. 8107

dx.doi.org/10.1021/ie201408p | Ind. Eng. Chem. Res. 2012, 51, 8095−8107