- Email: [email protected]
Adaptive robust optimization with minimax regret criterion: Multiobjective optimization framework and computational algorithm for planning and scheduling under uncertainty
Accepted Manuscript Title: Data-driven adaptive robust optimization with minimax regret criterion: Multiobjective optimization framework and computational algorithm for planning and scheduling under uncertainty Authors: Chao Ning, Fengqi You PII: DOI: Reference:
S0098-1354(17)30345-9 https://doi.org/10.1016/j.compchemeng.2017.09.026 CACE 5906
To appear in:
Computers and Chemical Engineering
Received date: Revised date: Accepted date:
14-7-2017 27-9-2017 28-9-2017
Please cite this article as: Ning, Chao., & You, Fengqi., Data-driven adaptive robust optimization with minimax regret criterion: Multiobjective optimization framework and computational algorithm for planning and scheduling under uncertainty.Computers and Chemical Engineering https://doi.org/10.1016/j.compchemeng.2017.09.026 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Data-Driven Adaptive Robust Optimization with Minimax Regret Criterion: Multiobjective Optimization Framework and Computational Algorithm for Planning and Scheduling under Uncertainty
Chao Ning, Fengqi You*
Robert Frederick Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, New York 14853, USA *
Corresponding author. Phone: (607) 255-1162; Fax: (607) 255-9166; E-mail:
[email protected] September 28, 2017 Submitted to Computers & Chemical Engineering
Highlights
ARO framework that incorporates conventional robustness and minimax regret criteria Solution algorithms for bi-criterion multilevel mixed-integer programming problem Innovative applications on planning and scheduling under uncertainty Leveraging big data analytics for data-driven decision-making under uncertainty
Abstract Regret is defined as the deviation of objective value from the perfect information solution, and serves as an important evaluation metric for decision-making under uncertainty. This paper proposes a novel framework that effectively incorporates the minimax regret criterion into twostage adaptive robust optimization (ARO). In addition to the conventional robustness criterion, this ARO framework also simultaneously optimizes the worst-case regret to push the 1
performance of the resulting solution towards the utopia one under perfect information. By using a data-driven uncertainty set, we formulate a multiobjective ARO problem that generates a set of Pareto-optimal solutions to reveal the systematic trade-offs between the conventional robustness and minimax regret criteria. The resulting multi-level mixed-integer programming problem cannot be solved directly by any off-the-shelf optimization solvers, so we further propose tailored column-and-constraint generation algorithms to address the computational challenge. Two applications on process network planning and batch process scheduling are presented to demonstrate the applicability of the proposed framework and the efficiency of the proposed solution algorithms.
Key words: Adaptive robust optimization, minimax regret, big data, multiobjective optimization, planning and scheduling
Nomenclature DDBCARO Framework The main notations used in the DDBCARO modeling framework are listed below. Parameters b
vector of cost coefficients corresponding to second-stage decisions
c
vector of cost coefficients corresponding to first-stage decisions
u
vector of uncertainties
μk
vector of parameters in normal Wishart distribution
ηk
mean vector in the Dirichlet process mixture model of Gaussian distributions
Cw
worst-case cost
li
index of mixture components assigned to the observation oi
m
total number of basic uncertainty sets
M
truncation level in the Variational inference
M0
constant with a very large value
Rw
worst-case regret 2
RHw
optimal highest worst-case regret
RLw
lowest worst-case regret
U
uncertainty set
Ui
basic uncertainty set
oi
the observation generated from the Dirichlet process mixture model
α
concentration parameter in the Dirichlet process
λk
parameter in normal Wishart distribution
θk
random variable drawn from the base distribution
πk
parameter of weight in the Dirichlet process
ωi
parameter in normal Wishart distribution
βk
proportion of stick being broken from the remaining stick
Λi
scaling factor in the data-driven uncertainty set
i
a hyper-parameter for Beta distribution in the Variational inference
vi
a hyper-parameter for Beta distribution in the Variational inference
γi
weight of mixture component
γ*
threshold value for the weight of mixture component
Θ
support of the base distribution
F0
base distribution in the Dirichlet process
Φi
uncertainty budget in data-driven uncertainty set
A
coefficient matrix corresponding to the first-stage decisions
Hk
precision matrix in the Dirichlet process mixture model of Gaussian distributions
M
coefficient matrix corresponding to the uncertainties
T
technology matrix in the DDBCARO framework
W
recourse matrix in the DDBCARO framework
Ψk
parameter in the normal Wishart distribution Variables
x
vector of first-stage decisions 3
y
vector of second-stage decisions
Application to planning of process networks The sets, parameters, and variables used in the application part are summarized below. Sets/indices I
set of processes indexed by i
J
set of chemicals indexed by j
T
set of time periods indexed by t Parameters
cei
maximum number of expansions for process i over the planning horizon
dujt
demand of chemical j in time period t
sujt
supply of chemical j in time period t
c1it
variable investment cost for process i in time period t
c2it
fixed investment cost for process i in time period t
c3it
unit operating cost for process i in time period t
c4it
purchase price of chemical j in time period t
cbt
maximum allowable investment in time period t
qeitL
lower bound for capacity expansion of process i in time period t
qeitU
upper bound for capacity expansion of process i in time period t
vjt
sale price of chemical j in time period t
κij
mass balance coefficient for chemical j in process i Binary variables
Yit
binary variable that indicates whether process i is chosen for expansion in time period t Continuous variables
Pjt
purchase amount of chemical j in time period t
Qit
total capacity of process i in time period t
QEit
capacity expansion of process i in time period t
Sjt
sale amount of chemical j in time period t
Wit
operation level of process i in time period t 4
Application to batch process scheduling The sets, parameters, and variables used in the data-driven robust scheduling of batch processes are summarized below. Sets/indices I
set of tasks indexed by i
J
set of equipment units indexed by j
N
set of time points indexed by n
S
set of states indexed by s Parameters
bimin
lower bound of batch size of task i
bimax
upper bound of batch size of task i
cs
storage capacity for state s
ds
demand for state s
fci
fixed cost of task i
ht
time horizon
mis
mass balance coefficient for consumption or production of state s in task i
ps
price of state s
sints
initial amount of state s
sisi
binary parameter that indicates whether task i uses state s as input
snps
binary parameter that indicates whether state s is a product
sosi
binary parameter that indicates whether task i produces state s as output
vdi
variable processing time of task i
vci
variable cost of task i Binary variables
Wsin
binary variable that indicates whether task i starts at time point n
Wfin
binary variable that indicates whether task i finishes at or before time point n
5
Continuous variables Bsin
batch size of task i that starts at time point n
Bpin
batch size of task i being processed at time point n
Bfin
batch size of task i that finishes at or before time point n
BIisn
amount of state s used as input of task i at time point n
BOisn
amount of state s produced from task i at time point n
Din
processing time of task i that starts at time point n
SAs
amount of state s sold at time point n
STsn
inventory level of state s at time point n
Tn
time corresponding to time point n
Tfin
finish time of task i that starts at time point n
Tsin
start time of task i that starts at time point n
1. Introduction In recent years, decision making under uncertainty has received tremendous attention from both academia and industry [1-4]. As an inherent characteristic of process systems, uncertainty could render the deterministic solution suboptimal or even infeasible [5-7]. To this end, extensive research efforts have been made towards the development of systematic methodologies for handling uncertainty, such as stochastic programming approach and robust optimization [812]. Robust optimization has attracted considerable interest from the process systems engineering community due to its good balance between solution quality and computational tractability [13-17]. By introducing recourse actions, adaptive robust optimization (ARO) captures the sequential nature of decision-making processes, and is typically less conservative than static robust optimization [18, 19]. Therefore, it has a fruitful of applications, including process scheduling [20, 21], energy systems [22-24], and supply chain optimization [25]. Although conventional ARO has a number of attractive features, it does not account for an important evaluation metric, known as regret, in decision-making theory [26]. In practice, the solution of an optimization problem under the presence of uncertainty is usually evaluated ex post by comparing its performance with that of the perfect information solution [27]. The objective value of the perfect information solution is the best possible value that could be 6
achieved by assuming that the decision maker had known the exact uncertainty realization a priori. While the perfect information solution is a utopia in the sequential decision making under uncertainty, this solution serves as an ideal benchmark for the quality of the optimized solution. The minimax regret criterion holds the promise to significantly improve the aforementioned evaluation performance, since it aims for the least deviation from the ideal benchmark [28, 29]. Due to its practical relevance, the minimax regret criterion has demonstrated a broad array of successful applications, such as dynamic pricing [30], energy supply systems [31-33], and sustainability [34]. Additionally, a decision maker may hold different degrees of preference towards the regret metric. For instance, one might be willing to sacrifice the performance in minimax regret criterion for a better performance in the conventional robustness criterion or vice versa. Unfortunately, ARO models with a single criterion fail to meet this research need, because they are not capable to reveal the trade-offs between the conventional robustness and minimax regret criteria. Based on all the concerns above, it is imperative to develop an ARO framework that effectively accounts for the conventional robustness and minimax regret criteria. To fill the knowledge gap, we are facing a number of research challenges. The first challenge is how to formulate an ARO model that simultaneously optimizes the conventional robustness and minimax regret criteria, and also reveals the systematic trade-offs between them. Besides, the multi-level optimization structure, coupled with the bi-objective feature, makes the resulting optimization problem computationally challenging. Therefore, the second research challenge is to develop an efficient computational algorithm for the global optimization of the resulting ARO problem, which cannot be solved directly by any off-the-shelf optimization solvers. The third challenge lies in the innovative applications of the general framework to address various planning and scheduling problems. This paper proposes a general bi-criterion ARO (BCARO) framework that effectively accounts for the conventional robustness and minimax regret criteria. By minimizing the worstcase regret over all uncertainty realizations within an uncertainty set, the proposed framework pushes the performance of the resulting solution towards the utopia one. In this way, the evaluation performance against the ideal benchmark is significantly improved. Since the worstcase uncertainty realization is not uniquely defined for stochastic multiobjective optimization, different concepts of robust efficiency have been proposed [35]. In this work, we employ the concept of point-based robust efficiency due to its clear interpretation and wide applications [367
38]. The proposed BCARO framework generates a set of Pareto-optimal solutions to reveal the systematic trade-offs between the conventional robustness and minimax regret criteria. This framework is general enough to apply to various types of uncertainty sets, such as the polyhedral and data-driven uncertainty sets [11, 39]. A data-driven ARO model of this framework is developed based on the Dirichlet process mixture model [39]. The resulting optimization problem is formulated as a bi-criterion multi-level mixed-integer programming problem, which cannot be solved directly by any off-the-shelf optimization solvers. Therefore, we further propose tailored column-and-constraint generation (C&CG) algorithms for its efficient solution. To illustrate the applicability of the proposed framework and efficiency of the proposed solution algorithms, we present two applications on process network planning and batch process scheduling. The major novelties of this paper are summarized as follows.
A general BCARO framework that effectively accounts for the conventional robustness and minimax regret criteria, and reveals the trade-offs between them;
Tailored decomposition-based C&CG algorithms that efficiently solve the resulting bicriterion multi-level optimization problems;
Innovative applications of the general framework on strategic planning of process networks and multipurpose batch process scheduling under uncertainty.
The rest of this paper is organized as follows. Section 2 provides relevant background materials on conventional ARO and the minimax regret criterion. In Section 3, we present the general BCARO modeling framework and solution algorithms. Section 4 provides a numerical example. Two applications on process network planning and batch process scheduling are illustrated in Section 5 and Section 6, respectively. Conclusions are drawn in Section 7.
2. Background In this section, we briefly review relevant background on conventional ARO and the minimax regret criterion.
2.1. Conventional adaptive robust optimization The general form of a conventional two-stage adaptive robust mixed-integer programming model is given as follows [19]. 8
min max min cT x bT y x
uU y ( x ,u )
s.t. Ax d, x Rn1 Z n2
(1)
x, u y Rn3 : Wy h Tx Mu
where x is the first-stage decision made “here-and-now” before the uncertainty u is realized, while the second-stage decision or recourse decision y is postponed in a “wait-and-see” manner after the uncertainties are revealed. x includes both continuous and integer variables, while y only includes continuous variables. c and b are the vectors of the cost coefficients. U is an uncertainty set that characterizes the region in which uncertainty realizations reside. The worst-case cost associated with a certain first-stage decision x over all the uncertainty realizations in the uncertainty set U is given by [32], C w x max min cT x bT y
(2)
uU y ( x ,u )
With the definition of worst-case cost in (2), the conventional ARO problem in (1) can be interpreted as an optimization problem that aims to find a first-stage decision x with the lowest worst-case cost. Note that the conventional robustness criterion in (1) is also known as the minimax cost criterion. After observing the uncertainty realization, the decision maker takes the “wait-and-see” decision y based on the fully adaptive decision rule, which is determined by the following optimization problem. min bT y y
s.t. y Rn3
(3)
Wy h Tx* Mu
where u is the observed uncertainty realization, and x* is the optimal first-stage decision in the conventional ARO problem. Note that x* in (3) should be treated as a parameter since it has been made at the initial stage of the decision-making process.
2.2. Adaptive robust optimization under the minimax regret criterion As mentioned in Section 1, the perfect information solution could serve as an ideal benchmark for the quality of an optimized solution [27]. The optimization problem under perfect information of uncertainty realization is presented in (4).
9
u min cT xˆ bT yˆ xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2
(4)
Wyˆ Txˆ h Mu, yˆ Rn3
where Π(u) is the best objective value that could be achieved by assuming that the decision maker had known the uncertainty realization a priori. Note that, in optimization problem (4), all the decisions are made in a “wait-and-see” mode with full information of uncertainty realization. Given a first-stage decision x, the regret under a specific uncertainty realization u is defined by, R x, u min cT x bT y u
(5)
y ( x ,u )
where R(x, u) represents the regret. It characterizes the deviation, in terms of total costs, from the perfect information solution under the uncertainty realization u. Under the minimax regret criterion, the decision maker aims to minimize the worst-case regret [40], which is shown below: R w x max R x, u uU
max min cT x bT y u y ( x,u ) uU
(6)
max min cT x bT y u uU y ( x ,u )
where Rw(x) denotes the worst-case regret over all uncertainty realizations within the uncertainty set U. Based on the above definitions, a general two-stage adaptive robust mixed-integer programming model with the minimax regret objective is then given by, min max min cT x bT y u x uU y ( x ,u ) s.t. Ax d, x Rn1 Z n2
x, u y Rn3 : Wy h Tx Mu u min c xˆ b yˆ T
(7)
T
xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2 Wyˆ Txˆ h Mu, yˆ Rn3
The optimization problem in (7) focuses on the relative robustness with respect to the perfect information solution, while the optimization problem in (1) aims for the abosulte robustness. Both optimization problems have multi-level optimization structures. However, the 10
ARO problem under the minimax regret criterion has one more level due to the innermost perfect information optimization problem. After uncertainty is realized, the “wait-and-see” decision y is made based on the same optimization problem in (3). This is because u is a constant at the second stage of the decision-making process. Based on the above conventional ARO and minimax regret criterion, the proposed BCARO framework is developed in the next section.
3. Bi-Criterion Adaptive Robust Optimization Framework In this section, we propose a novel framework that incorporates the minimax regret criterion into conventional ARO through the lens of multiobjective optimization. Notably, the proposed BCARO framework is general enough to apply to various uncertainty sets of conventional ARO. For the ease of exposition, the BCARO framework is developed with a data-driven uncertainty set [39].
3.1. Bi-criterion adaptive robust optimization model To model uncertainty, we adopt a data-driven uncertainty set using both l1 and l∞ norms [39], as presented in (8) and (9). U
U i U1
U2
(8)
Um
i: i *
Ui u u μi i Ψ1/2 z i i z,
1, z 1 i
(9)
where U is the data-driven uncertainty set, Ui is the basic uncertainty set corresponding to the i th component, Λi is a scaling factor, γi is the weight of the i th mixture component,
i
i
i
i 1
v i
j 1
vj j
vj
and M 1 i 1 i . The value of γi indicates the probability of the M 1
component i, and γ* is a threshold value. m is the total number of mixture components having probabilities larger than γ*. i
i 1
i i 1 dim u
[41], and τi, νi, μi, λi, ωi, ψi are the
variational inference results. Details of the data-driven uncertainty set are presented in the Appendix.
11
We further propose a general ARO framework that accounts for the conventional robustness and minimax regret criteria simultaneously with the data-driven uncertainty set. The proposed data-driven bi-criterion adaptive robust mixed-integer programming model, denoted as (DDBCARO), is presented as follows. min Worst-case cost C w x x
min Worst-case regret R w x x
s.t. Ax d, x Rn1 Z n2 C w x max max min cT x bT y i1, , m uU i y ( x ,u )
(DDBCARO)
R w x max max min cT x bT y u i1, , m uU i y ( x ,u )
x, u y Rn3 : Wy h Tx Mu u min cT xˆ bT yˆ xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2 Wyˆ Txˆ h Mu, yˆ Rn3
where Ui is the data-driven basic uncertainty set defined in (9). (DDBCARO) optimizes the objectives of the worst-case cost Cw(x) and the worst-case regret Rw(x) simultaneously. By explicitly incorporating the minimax regret criterion, it pushes the performance of the resulting solution towards the utopia one under perfect information of uncertainty realization. One way to handle multiobjective optimization is the weighted sum method. Essentially, the method attaches weights to different objective functions, and eventually optimizes over a weighed sum of all the objectives. Following the weighted sum approach [42], we further propose a unified optimization model (DDBCARO-λ) shown as follows.
12
min 1 C w x R w x x s.t. Ax d, x Rn1 Z n2 C w x max max min cT x bT y i1, , m uU i y ( x ,u ) R w x max max min cT x bT y u i1, , m uU i y ( x ,u )
(DDBCARO-λ)
x, u y Rn3 : Wy h Tx Mu u min cT xˆ bT yˆ xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2 Wyˆ Txˆ h Mu, yˆ Rn3
where λ is a parameter ranging from 0 to 1. (DDBCARO-λ) could provide a wide spectrum of solutions by changing the value of λ from 0 to 1. Specifically, it reduces to the conventional ARO model when λ equals to 0, and it becomes the ARO with minimax regret criterion when λ equals to 1. Thus, the proposed modeling framework (DDBCARO-λ) unifies the ARO models in Section 2.1 and Section 2.2. If the value of λ is chosen between 0 and 1, the proposed model makes a trade-off between the conventional robustness and minimax regret criteria. Another popular method for multiobjective optimization is the ε-constraint method that holds one of the objective functions as a constraint. Following the ε-constraint method [42], the objective on the worst-case regret is transformed into an extra constraint that specifies the upper bound of Rw(x). The resulting optimization model (DDBCARO-ε) is presented by, min C w x x
s.t. R w x Ax d, x Rn1 Z n2 C w x max max min cT x bT y i1, , m uU i y ( x ,u )
(DDBCARO-ε)
R w x max max min cT x bT y u i1, , m uU i y ( x ,u )
x, u y Rn3 : Wy h Tx Mu u min c xˆ b yˆ T
T
xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2 Wyˆ Txˆ h Mu, yˆ Rn3
13
Similarly, (DDBCARO-ε) can reveal the trade-offs between the conventional robustness and minimax regret criteria by following three steps of the ɛ-constraint method [43]. The first step is to minimize the worst-case regret subject to the same constraints in (DDBCARO), i.e. without the ɛ-constraint, to obtain the lowest worst-case regret RLw . The second one is to minimize the worst-case cost subject to the same constraints in (DDBCARO) to obtain its optimal solution x*, which in turn yields the optimal highest worst-case regret RHw R w x* . The last step is to divide the interval between RLw and RHw into small intervals, and solve problem (DDBCARO-ɛ) by fixing ɛ to the end points of each small interval. The resulting optimization problems (DDBCARO-λ) and (DDBCARO-ε) are both multilevel mixed-integer programming problems. They cannot be solved directly by any off-the-shelf optimization solvers due to the multi-level optimization structures. To this end, computational algorithms are further developed to efficiently solve the optimization problems.
3.2. Computational algorithm In this section, we address the computational challenge by developing two efficient solution algorithms for global optimization of problems (DDBCARO-λ) and (DDBCARO-ε) in Section 3.2.1 and 3.2.2, respectively. To facilitate the solution procedure, we first reformulate the worst-case regret to reduce the number of levels. R w x max max min cT x bT y u i1, , m uU i y ( x ,u ) max max min cT x bT y min cT xˆ bT yˆ xˆ , yˆ u i1, , m uU i y ( x ,u ) cT x max max min bT y cT xˆ bT yˆ i1, , m xˆ , yˆ u ,uU i y ( x ,u )
(10)
where u xˆ Rn Z n , y Rn : Axˆ d, Wyˆ Txˆ h Mu . 1
2
3
3.2.1. Solution algorithm for (DDBCARO-λ) Based on the above reformulation of the worst-case regret, (DDBCARO-λ) can be transformed into (11)-(15) using the epigraph reformulation. min 1 1 2
(11)
x ,1 ,2
14
s.t. Ax d, x Rn1 Z n2
(12)
1 cT x max max min bT y i1, , m uU y ( x ,u )
(13)
i
2 cT x max
min bT y cT xˆ bT yˆ
max
i1, , m xˆ , yˆ u ,uUi y ( x ,u )
x, u y Rn3 : Wy h Tx Mu
(14) (15)
where η1 and η2 are the auxiliary variables. Note the reformulated problem above is still a multilevel optimization problem due to the max-max-min structures. This multi-level optimization problem can be reformulated to an equivalent single-level optimization problem by enumerating all the extreme points, and introducing their corresponding recourse variables [19]. However, this single-level optimization problem could be a large-scale optimization problem due to the potentially large number of extreme points. Therefore, we partially enumerate the extreme points, and construct the master problem (MP-λ) shown below. min 1 1 2
x ,1 ,2 , y1l , y l2
s.t. Ax d
1 cT x bT y1l , l 1, , r (MP-λ)
Tx Wy1l h Mu1l , l 1,
,r
Tx Wy l2 h Mul2 , l 1,
,r
2 cT x bT y l2 cT xˆ l bT yˆ l , l 1, , r x R Z , y , y R , l 1, n1
n2
l 1
n3
l 2
,r
where r is the current iteration number. Since the master problem (MP-λ) contains only a subset of the extreme points, it provides a relaxation of the original problem. To enumerate the important uncertainty realizations for both (13) and (14) on-the-fly, two groups of sub-problems are solved in each iteration. The sub-problem corresponding to (13), denoted as (SUP1,i), is presented by Q1i x max min bT y1 u1U i
y1
s.t. Wy1 h Tx Mu1 y1 Rn3
where i is the index of basic uncertainty sets. The optimal u1 is the identified uncertainty realization that should be added to the (MP-λ).
15
(16)
Note that (16) is a max-min optimization problem, which can be transformed to a singlelevel optimization problem by employing strong duality or KKT conditions [19, 44]. The resulting problem after taking the dual of inner problem is shown in (17). T Q1i x max h Tx Mμi φ i MΨ1/2 i i tj z j φt z ,φ t j T s.t. W φ b
(17)
φ0 z
1, z 1 i
where φ is the vector of dual variables associated with the constraints Wy h Tx Mu , φt represents the t th component in vector φ, and zj is the j th component in vector z. Since there are bilinear terms zjφt in the objective function of (17), we need to further reformulate the sub-problem using a linearization technique. To facilitate the derivation, zj can be divided into two parts z j z j z j , where z j and z j are two non-negative variables. Thus, we have u μi i Ψ1/2 i i z z
(18)
The uncertainty budget parameter Φi is typically set as an integer value due to its physical meaning, as reported in the existing literature [44, 45]. For the integer uncertainty budget, the optimal z j and z j take the values of 0 or 1, and these two variables can be restricted to binary
variables. By using Glover’s linearization for the bilinear terms G tj z j φt and G tj z j φt [46], which are the products of a binary variable z j (and z j ) and a continuous variable φt, we can
reformulate (17) into the following optimization problem in (19).
16
Q1i x
h Tx Mμ T φ Tr MΨ1/2 T G G i i i i s.t. WT φ b max
φ , z , z ,G ,G
φ0 0 G φ eT 0 G φ eT G e M 0 z
T
G e M 0 z
T
e M e z
G φ eT e M 0 e z
T
G φ eT
T
0
eT z z i z z e,
z , z 0,1
K
(19)
where e and e both represent column vectors whose elements are all ones, and they have different dimensions, i.e., dim e dim u K and dim e dim φ . The t th row and j th column elements of matrices G and G are denoted as G tj and G tj , respectively. The notation Tr is the trace of a matrix, e.g. for matrix Q qij R nn , Tr Q in1 qii . Note that the vector or matrix inequalities should be interpreted elementwise. The sub-problem corresponding to (14), denoted as (SUP2,i), is presented by Q2i x
max
xˆ , yˆ u 2 ,u 2 U i
min bT y 2 cT xˆ bT yˆ y2
s.t. Wy 2 h Tx Mu 2 y 2 Rn3
Similarly, we can reformulate optimization problem (20) as follows
17
(20)
h Tx Mμ T φ Tr MΨ1/2 T G G cT xˆ bT yˆ Q2i x max i i i i φ,z ,z G ,G , xˆ , yˆ
s.t. WT φ b φ0 0 G φ eT 0 G φ eT G e M 0 z
T
G e M 0 z
T
e M e z
G φ eT e M 0 e z
T
G φ eT
T
0
eT z z i z z e,
z , z 0,1
K
Axˆ d Wyˆ Txˆ h Mμ i i MΨ1/2 i i z z
(21)
xˆ Rn1 Z n2 , yˆ Rn3
The proposed solution algorithm iteratively solves a master problem and two groups of subproblems until the relative optimality gap reduces to a predefined tolerance. The pseudocode of the algorithm is given in Figure 1, where ζ represents the predefined tolerance, and m represents the total number of basic uncertainty sets. As the optimality and feasibility cuts are added to the master problem at each iteration, the lower bound provided by (MP-λ) becomes tighter as the iteration moves on. The computational algorithm terminates in finite iterations because of the finite number of extreme points.
18
3.2.2. Solution algorithm for (DDBCARO-ε) We employ the techniques developed in the previous subsection to address problem (DDBCARO-ε), which can be first reformulated into (22) using the epigraph reformulation. min x ,
s.t. Ax d, x Rn1 Z n2
cT x max max min bT y i1, , m uU y ( x ,u ) i
cT x max
(22)
min b y c xˆ b yˆ T
max
i1, , m xˆ , yˆ u ,uU i
T
T
y ( x ,u )
x, u y Rn3 : Wy h Tx Mu
By enumerating a subset of the extreme points, we arrive at its corresponding master problem (MP-ε), which is shown below. min
x , , y1l , y l2
s.t. Ax d
cT x bT y1l , l 1, , r (MP-ε)
Tx Wy1l h Mu1l , l 1,
,r
Tx Wy l2 h Mul2 , l 1,
,r
cT x bT y l2 cT xˆ l bT yˆ l , l 1, , r x R Z , y , y R , l 1, n1
n2
l 1
n3
l 2
,r
Note that the sub-problems for (DDBCARO-ε) are the same as in the previous subsection. The pseudocode of the proposed solution algorithm is directly provided in Figure 2. To handle the ɛ-constraint, the solution algorithm in Figure 2 incorporates additional steps to check feasibility. From Steps 13-15, we can see that it updates the upper bound only when the ɛconstraint is feasible.
19
4. Numerical Example In this section, a small numerical example of the proposed BCARO framework, as shown in (23), is presented to illustrate the conventional robustness and minimax regret criteria and the trade-off between them. We consider the following numerical example: min max min x
uU y ( x ,u )
3x1 5 x2 12 x3 15y1 30 y2 22 y3
min max min 3 x1 5 x2 12 x3 15y1 30 y2 22 y3 u x uU y ( x ,u ) s.t. x1 x2 x3 200 xi 0, i 1, 2,3
(x, u) y Rn3 : yi ui xi , i 1, 2,3
(23)
u min 3xˆ1 5 xˆ2 12 xˆ3 15yˆ 1 30 yˆ 2 22 yˆ3 xˆ , yˆ
s.t. xˆ1 xˆ2 xˆ3 200 xˆi yˆi ui , i 1, 2,3 xˆi , yˆi 0, i 1, 2,3
where x1, x2 and x3 are the first-stage or “here-and-now” decision variables, u1, u2 and u3 are the uncertain parameters, and y1, y2 and y3 are the second-stage or “wait-and-see” decision variables. Note that the “here-and-now” decisions are made prior to the resolution of uncertainties, whereas “wait-and-see” decisions could be made after the uncertainty realizations. U is the data-driven uncertainty set. Π(u) is the objective value of the perfect information solution. Problem (23) is an instance of (BCARO). By using the weighted sum method, it can be transformed into optimization problem (BCARO-λ), where λ=0 corresponds to the conventional robustness criterion, and λ=1 corresponds to the minimax regret criterion. A trade-off can be achieved by choosing a value of λ between 0 and 1. The scatter plot of uncertainty data and the data-driven uncertainty set are shown in Figure 3, where the red dots represent uncertainty realizations of u1, u2 and u3, and the union of the
20
green polytopes is the data-driven uncertainty set. A total of 500 uncertainty data points is used for uncertainty modeling.
The optimization models and the proposed tailored C&CG algorithm are coded in GAMS 24.7.3. The computational experiments are implemented on a computer with an Intel (R) Core (TM) i7-6700 CPU @ 3.40 GHz and 32 GB RAM. The relative optimality tolerance for the tailored C&CG algorithm is set to be 0.01%. The numerical results are summarized in Table 1. We can observe that the conventional robustness criterion generates the lowest worst-case cost, and has the highest worst-case regret. Since the minimax regret criterion only focuses on the regret metric, it reduces the worst-case regret by 48.7% compared with the conventional robustness criterion. However, it suffers from the highest worst-case cost of 1,833.2. A trade-off between the conventional robustness and minimax regret criteria can be made when choosing λ to be 0.4 in (BCARO-λ). In the two-stage ARO setting, the cost is calculated as a sum of the first-stage cost and the second-stage cost. The costs can be different for different uncertainty realizations. The first-stage cost is 3x1+5x2+12x3, which does not depend on the uncertainty realization. Hence, the first-stage cost can be calculated based on the values of the first-stage decisions provided in Table 1. The second-stage cost is 15y1+30y2+22y3, which depends on the uncertainty realizations. To obtain the values of y1, y2 and y3, the following optimization problem in (24) is solved for each uncertainty realization. min 15y1 30 y2 22 y3 y
s.t. yi ui xi , i 1, 2,3
(24)
yi 0, i 1, 2,3
For each uncertainty realization, the cost can then be calculated as 3x1+5x2+12x3+15y1+30y2+22y3. We use an uncertainty realization in Figure 3 to illustrate the concept of regret. For this uncertainty realization, u1=55.2, u2=41.8, u3=74.4. The perfect information solution corresponding to this uncertainty realization is xˆ1 55.2, xˆ2 41.8, xˆ3 74.4, yˆ1 0, yˆ2 0, yˆ3 0 , and the corresponding cost is 1,267.4. The (BCARO-λ) solution at λ=1.0 is x1=58.1, x2=48.9, x3=49.8, y1=0, y2=0, y3=24.6, and the corresponding cost is 1,557.6. The regret is the difference between these two costs, which is 290.2. 21
To further demonstrate the trade-off between the conventional robustness and minimax regret criteria, we present the Pareto-optimal curve in Figure 4. The X-axis represents the worstcase regret, while the Y-axis denotes the worst-case cost. The Pareto-optimal curve is plotted following three steps of the ɛ-constraint method [43]. First, we solve problem (BCARO-λ) at λ=1 to obtain the lowest worst-case regret RLw . Next, we solve problem (BCARO-λ) at λ=0 to obtain the optimal solution x*, which in turn yields the optimal highest worst-case regret RHw R w x* . The last step is dividing the interval between RLw and RHw into small intervals, and solving problem (BCARO-ɛ) by fixing ɛ to the end points of each small interval. In Figure 4, the yellow dot represents a wise trade-off solution. Compared with the conventional robustness criterion, it greatly reduces the worst-case regret by 21.1% at the expense of a tiny increase (3.9%) in worstcase cost. The result illustrates the advantage of the proposed framework in revealing a systematic trade-off between the conventional robustness and minimax regret criteria.
5. Application to Strategic Planning of Chemical Process Networks under Supply and Demand Uncertainties To demonstrate the application of the proposed modeling framework and computational algorithms, we consider a case study on a multi-period process network planning under uncertainty. Large-scale chemical complexes often consist of interconnected processes and various chemicals [47-49]. These interconnections allow for the synergies between different processes in chemical production [49, 50]. The chemicals in a process network include feedstocks, intermediates and products. The process network planning problem determines the purchase levels of feedstocks, sales of products, capacity expansions, and production profiles of processes at each time period over the entire planning horizon [51]. In this application, the supply of feedstocks and demand of products are uncertain.
5.1. Bi-criterion adaptive robust planning of process networks model The bi-criterion ARO model for process network planning under uncertainty is formulated as follows. The two objectives are to maximize the worst-case net present value (NPV) and to 22
minimize the worst-case regret given in (25) and (26). The definition of worst-case NPV and worst-case regret are presented in (27)-(28). The constraints can be classified into capacity expansion constraints (29)-(30), budget constraints (31)-(32), production level constraint (33), mass balance constraint (34), supply and demand constraints (35)-(36), non-negativity constraints (37)-(38), and integrity constraint (39). A list of indices/sets, parameters and variables is given in the Nomenclature section, where all the parameters are denoted in lowercase symbols, and all the variables are denoted in upper-case symbols. The objective functions are defined by (25) and (26), as shown below max NPV w Qit , QEit , Yit
(25)
min R w Qit , QEit , Yit
(26)
Qit ,QEit ,Yit
Qit ,QEit ,Yit
where Qit is a decision on total capacity of process i in time period t, QEit is a decision variable for capacity expansion level of process i in time period t, Yit is a binary variable to determine on whether process i is expanded in time period t, NPVw represents the worst-case NPV and Rw denotes the worst-case regret. The worst-case NPV and the worst-case regret are presented in (27) and (28), respectively. NPV w mindem max jt S jt c1it QEit c2it Yit c3it Wit c 4 jt Pjt (27) du jt U , Pjt ,Qit , iI tT iI tT iI tT jJ tT jJ tT su U sup S ,W jt
jt
it
R w maxdem min du jt , su jt jt S jt c1it QEit c2it Yit c3it Wit c 4 jt Pjt du jt U , Pjt ,Qit , iI tT iI tT iI tT jJ tT jJ tT su U sup S ,W jt
jt
it
(28) where Sjt is the amount of chemical j sold to the market in time period t, Wit represents the operating level of process i during time period t, and Pjt is the amount of chemical j purchased in time period t. Parameters c1it, c2it, c3it and c4jt are coefficients for variable investment cost, fixed investment cost, operating cost and purchase cost, respectively. νjt is the sale price of chemical j in time period t. dujt is the demand of chemical j in time period t and sujt is the market availability of chemical j in time period t. Π(dujt, sujt) is the NPV of the perfection information solution. Constraint (29) enforces the upper and lower bounds of capacity expansions. Specifically, if Yit=1, i.e. process i is chosen to be expanded in time period t, the expanded capacity should be within the range [qeitL, qeitU]. 23
qeitL Yit QEit qeitU Yit ,
i, t
(29)
Constraint (30) specifies the updates of total capacity of each process i. Qit Qit 1 QEit ,
i, t
(30)
Constraints (31) and (32) specify the maximum number of capacity expansions for each process and investment budget for each time period, respectively.
Y
it
cei ,
i
(31)
t
c1
it
QEit c2it Yit cbt ,
t
(32)
i
where cei is the maximum expansion times for process i, and cbt is the investment budget for time period t. Constraint (33) enforces that the production level of a process cannot exceed its total capacity.
Wit Qit ,
i, t
(33)
Constraint (34) specifies the mass balance for chemicals. To be more specific, the purchase amount plus the production amount is equal to the total amount sold to the market and consumed by processes. Pjt ij Wit S jt 0,
j , t
(34)
i
where κij is the mass balance coefficient for chemical j in process i. Constraint (35) enforces that the purchase amount of a feedstock should not be higher than its available amount provided by the supplier. Constraint (36) specifies that the sale amount of a product cannot exceed its market demand. Pjt su jt ,
j, t
(35)
S jt du jt ,
s, j, t
(36)
Constraints (37)-(38) indicate the non-negativity of continuous decisions. QEit 0, Qit 0 i, t
(37)
Pjt , S jt , Wit 0,
(38)
j, t
Constraint (39) ensures that Yit is a binary variable. Yit 0,1 ,
i, t
(39)
24
The data-driven uncertainty sets for demand and supply uncertainties is shown in (40)-(41). U dem k : kdem *
du du μ
dem k
kdem Ψ dem dem z, z k k
1, z 1 dem k
(40)
dem where μ dem , Ψ dem are the inference results learned from uncertain product demand data k , k k
using the Variational inference algorithm. U sup k : ksup *
su su μ
sup k
sup ksup Ψsup z k k z,
1, z 1 sup k
(41)
sup sup where μsup are the inference results learned from uncertain supply data employing k , k , Ψk
the Variational inference algorithm. The resulting optimization problem on strategic planning of process network under uncertainty is a bi-objective multi-level mixed-integer programming problem. Due to the multilevel optimization structure, it cannot be solved directly by any off-the-shelf optimization solvers. We employ the proposed computational algorithms to efficiently solve the resulting large-scale optimization problems.
5.2. Case study overview In this subsection, a specific case study is presented to demonstrate the applicability of the proposed approach. The process network in this case study consists of 28 chemicals, 38 processes, 10 suppliers and 16 markets [48], as shown in Figure 5. Chemicals A-J represent feedstocks, which can be purchased from suppliers or manufactured by certain processes. Chemicals K-Z are products, which are sold to the markets. The intermediates are denoted as AA and AB. This complex process network has such flexibility that many manufacturing options are available. For example, Chemical E can be manufactured by Process 22, Process 25 or purchased from a supplier. In this case study, we consider 5 time periods over the planning horizon, and the duration of each time period is 2 years. It is assumed that processes 12, 13, 16 and 38 have initial capacities of 40, 25, 300 and 200 kt/y, respectively, at the beginning of time period 1. These processes are not allowed to be expanded until time period 2. The other chemical processes can be installed at the beginning of time period 1. Due to market fluctuations, availabilities of raw materials and demands of final products are subject to uncertainties in the process network planning problem [52]. In this case study, a total of 26 parameters, including supply of all feedstocks A-J and demands for all products K-Z are 25
uncertain. For the supply uncertainty, a total of 2,000 uncertainty data are used. Each data point has 10 dimensions for a combination of all feedstocks. For the demand uncertainty, 3,000 uncertainty data points are used, and each data point is for a combination of all the 16 products.
5.3. Results and discussions We solve the multiobjective multi-level optimization problem of process network planning using the computational algorithms proposed in Section 3.2. The optimization methods are modeled in GAMS 24.7.3 [53], and are implemented on a computer with an Intel (R) Core (TM) i7-6700 CPU @ 3.40 GHz and 32 GB RAM. The master problem and sub-problems are solved using CPLEX 12 with an optimality tolerance of 0. The relative optimality tolerance for the tailored C&CG algorithm is set to be 0.01%. The uncertainty budgets for supply and demand are set to be 3 and 5, respectively. A sensitivity analysis on the values of the uncertainty budgets is performed at the end of this subsection. We first solve optimization problem (BCARO-λ) using the computational algorithm in Section 3.2.1. To show how the two objective values change with parameter λ, we calculate the worst-case NPV and the worst-case regret under different choices of parameter λ and present them in Figure 6. The X axis represents parameter λ ranging from 0 to 1, and the Y axis represents the worst-case NPV or worst-case regret. As can be seen from Figure 6, when parameter λ increases, both the worst-case NPV and regret decreases. The proposed BCARO framework provides a wide spectrum of solutions, including the ARO solutions under the conventional robustness criterion (λ=0) and the minimax regret criterion (λ=1). Based on the decision maker’ preference to conventional robustness and minimax regret criteria, a wise-choice solution can be obtained by selecting the value of λ between 0 and 1. We select three instances of (BCARO-λ) corresponding to λ=0, λ=0.2 and λ=1, and summarize the computational results and model statistics in Table 2. Since the optimality cuts are added to the mater problem at each iteration, its problem size increases as the computational algorithm iterates. The master problem at the first iteration involves 190 integer variables, 1,323 26
continuous variables, and 1,793 constraints. The size of the master problem at the last iteration of the computational algorithm is listed in Table 2. Compared with the ARO problem under the conventional robustness criterion, the ARO problem under the minimax regret criterion costs a much longer computational time (409 seconds v.s. 180 seconds). Notably, all the three instances can be solved to the specific optimality tolerance within only 7 minutes.
To further investigate the trade-off between these two robustness criteria, we employ the εconstraint based computational algorithm in Section 3.2.2, and demonstrate the Pareto-optimal curve with cost breakdown in Figure 7. The Pareto-optimal curve is plotted following the same procedure as described in Section 4. On the Pareto-optimal curve in Figure 7, three points are selected for further discussions. Point A represents the ARO solution under the conventional robustness criterion. It has the highest worst-case NPV of $544.0 MM, and the largest worst-case regret of $705.9 MM. In contrast to Point A, Point C represents the ARO solution under the minimax regret criterion, having the smallest worst-case regret of $90.1 MM, and the lowest worst-case NPV of $464.3 MM. Comparing Point A with Point B, we can see that the worst-case regret is significantly reduced by 52.3%, while the worst-case NPV decreases by merely 3.2%. This suggests that Point B could be a wise-choice solution based on the trade-off between the worst-case NPV and regret. The details on cost breakdowns, including fixed investment cost, variable investment cost, operating cost, and purchase cost, are shown as pie charts in Figure 7. From the pie charts, we can observe that almost half of the total cost comes from purchasing feedstocks for Points A, B and C. Moreover, the percentage of the investment costs of Point C is higher than those of the other two points due to larger process capacities and more frequent capacity expansions. The optimal design and planning decisions at time period 4 determined by the solutions of Points A, B and C are shown in Figure 8, Figure 9 and Figure 10, respectively. In these three figures, the optimal total capacities are displayed under operating processes.
To illustrate the optimal capacity expansion activities, we present the capacity profiles during the entire planning horizon for Points A, B and C in Figure 11, Figure 12 and Figure 13, respectively. From Figure 11, we can see that a total of 14 processes are selected in the optimal process network for Point A. By comparing Figure 11 and Figure 13, we can conclude that the 27
optimal expansion frequency of Point C is higher than that of Point A. Additionally, the optimal total capacities of some processes for Points B and C are different. For example, the optimal total capacity of Process 32 is 245 kt/y determined by Point B, while the corresponding capacity is 55 kt/y larger determined by Point C.
To demonstrate the efficiency of the computational algorithms in Sections 3.2.1 and 3.2.2, we select two instances and display the upper and lower bounds in each iteration of the algorithms in Figure 14 and Figure 15, respectively. In these two figures, the green dots stand for the upper bounds, and the yellow circles represent the lower bounds. The X-axis represents the iteration number, and the Y-axis represents the objective function values. The first instance corresponds to optimization problem (BCARO-λ) with λ=0. It can be seen from Figure 14 that the relative optimality gap decreases significantly from 89.46% to 3.90% during the first three iterations. The proposed computational algorithm takes only 7 iterations to reduce the relative optimality gap to 0.70%, and requires 9 iterations to reach a relative optimality gap below the predefined tolerance of 0.01%. The result demonstrates the efficiency of the proposed computational algorithm in Section 3.2.1.
The second instance corresponds to Point B in Figure 7. The performance of the computational algorithm is shown in Figure 15. The lower bound is not updated in the first iteration because the first-stage decision does not satisfy the ɛ-constraint. The lower bound is updated after the feasibility cuts are added to the master problem in the second iteration. The proposed computational algorithm takes only 4 iterations to converge within a reasonable solution time. The result demonstrates the efficiency of the proposed computational algorithm in Section 3.2.2.
iteration of the tailored C&CG algorithm in the instance corresponding to Point B.
To further investigate how the values of uncertainty budgets influence the trade-off between the conventional robustness and minimax regret criteria, we present the Pareto-optimal curves under different uncertainty budgets in Figure 16. For all these Pareto-optimal curves, we can see 28
a similar trend that worst-case NPV decreases as the worst-case regret reduces. For a fixed worst-case regret, increasing either the supply or demand uncertainty budget leads to a lower worst-case NPV. By comparing the red, purple and green curves, we can see that the Paretooptimal curve lies close to the upper left corner as we reduce the uncertainty budget for demand. From the blue, purple and orange curves, we can conclude that the Pareto-optimal curve is less sensitive to the change in the supply uncertainty budget. Therefore, the impact of the demand uncertainty budget is more significant than that of the supply uncertainty budget.
6. Application to Batch Process Scheduling under Demand Uncertainty This section presents another application of the proposed approach to multipurpose batch process scheduling under demand uncertainty. Multipurpose batch process scheduling is an operational problem solved on a daily or weekly basis [54]. In the scheduling problem, the following information is given: the structure of the facility, recipe data, equipment units and their sizes, storage capacities of intermediates, product prices and processing time [55, 56]. The objective is to generate the highest profit by determining the optimal batch sizes, the assignment of tasks to equipment units, the sequencing and timing of tasks. The solution of a batch scheduling problem is often represented by Gantt charts. In this application, product demands are subject to uncertainty.
6.1. Bi-criterion adaptive robust scheduling of batch process model The multi-criterion ARO model for batch process scheduling is formulated as follows. Following the conventional ARO model of batch process scheduling [21], the task assignment decision is made “here-and-now”, while other decisions are made in a “wait-and-see” manner after knowing the demand uncertainty realization. The two objectives are to maximize the worstcase profit and to minimize the worst-case regret given in (42)-(43). The definition of worst-case profit and worst-case regret are presented in (44)-(45). The constraints can be classified into assignment constraints (46)-(51), time constraints (52)-(62), batch size constraints (63)-(79), mass balance and storage constraints (80)-(84), and non-negativity and integrity constraints (85) -(86). A list of indices/sets, parameters and variables is given in the Nomenclature section. The objective functions are defined by (42)-(43), as shown below max Profit w Wsin ,Wfin
Wsin ,Wfin
29
(42)
min R w Wsin ,Wfin
(43)
Wsin ,Wfin
The worst-case profit and regret are presented in (44) and (45), respectively.
Profit w min
max
dU Bfin , Bpin , Bsin , BIisn , BOisn , STin , SAs ,Tn , Din ,Tfin ,Tsin
R w max dU
ps STs|n| vci Bsin fci Wsin i n i n sP
d ps STs|n| vci Bsin fci Wsin Bfin , Bpin , Bsin , BIisn , i n i n sP BO , ST , SA ,T , min
isn
in
s
(44)
(45)
n
Din ,Tfin ,Tsin
Constraints (46)-(47) specify that at most one task can start or finish at each equipment unit at time point n. Constraint (48) specifies that a task must finish in the future once it starts. Constraint (49) enforces that at most one task can be processed at each unit at any time point n.
Ws
in
1
j, n
(46)
Wf
in
1
j, n
(47)
iI j
iI j
Ws
in
n
i
(48)
n
Ws iI j n ' n
Wfin ,
in '
Wfin ' 1,
j , n
(49)
Constraint (50) enforces that no tasks could finish at the first time point. Similarly, constraint (51) specifies that no task could begin at the last time point.
Wfi1 0
i
(50)
Wsi|n| 0
i
(51)
Constraints (52)-(61) are for the timing of tasks. Constraint (52) describes the processing time of each task. The relationship between starting time, finish time and processing time of a task is given in constraints (53)-(56). Constraint (57) specifies that task i starting at time point n1 should finish at or before time point n if Wfin equals to 1. Constraint (58) describes the relationship between starting time and the time corresponding to time point n. fdi Wsin +vdi Bsin Din Tsin Din ht 1 Wsin Tfin ,
30
i, n i, n
(52) (53)
Tfin Tsin Din ht 1 Wsin , Tfin Tfi n1 Wsin ht ,
i, n
Din Tfin Tfi n1 ,
i, n
Tfi n1 Tn ht 1 Wfin ,
Tsin Tn ,
i, n
(54) (55) (56)
i, n 2
i, n
(57) (58)
Constraints (59)-(62) defines the initial condition and upper bound of related variables.
T1 0 Tn T n1 ,
(59) n 2
(60)
T|n| ht
Tfin ht ,
(61)
i, n
(62)
Constraints (63)-(69) specify the lower and upper bounds of batch sizes, and linking variables associated with batches.
Bsin biminWsin ,
i, n
(63)
Bsin bimaxWsin ,
i, n
(64)
Bfin biminWfin ,
i, n
(65)
Bfin bimaxWfin ,
i, n
(66)
Bpin bimin Wsin ' Wf in ' , n ' n n ' n
i, n
(67)
Bpin bimax Wsin ' Wf in ' , n ' n n ' n
i, n
(68)
i, n 2
(69)
Bsi n1 Bpi n1 Bpin Bfin ,
Constraints (70)-(73) describe the input and output of batch operations, as well as their upper bounds.
sisi BIisn sisi mis Bsin 0, i, s, n
(70)
sosi BOisn sosi mis Bfin 0, i, s, n
(71)
sisi BIisn sisi bimax mis Wsin , i, s, n
(72)
31
sosi BOisn sosi bimax mis Wfin , i, s, n
(73)
Constraints (74)-(79) are the initial and final conditions for the variables associated with batch sizes.
Bpin 0, i, n 1
(74)
Bfin 0, i, n 1
(75)
BOisn 0, i, s, n 1
(76)
Bsin 0, i, n n
(77)
Bpin 0, i, n n
(78)
BIisn 0, i, s, n n
(79)
Constraint (80) describes the initial inventory of different materials.
STsn sisi BIisn sints , s, n 1
(80)
i
Constraint (81) is for the mass balance. STsn STs n1
B
iTOs
O isn
I Bisn ,
s, n
(81)
iTI s
Constraint (82) specifies the storage capacity. If a chemical is sold out as a final product, constraint (83) states that its sale amount is equal to its inventory level at the last time point. Constraint (84) enforces that sale amount of each final product should meet its demand. Constraints (85) and (86) defines nonnegative continuous variables and binary variables.
STin cs , s, n
(82)
snps STin snps SAs 0, s, n n
(83)
snps SAs snps ds , s
(84)
Bfin , Bpin , Bsin , BIisn , BOisn , STin , SAs , Tn , Din , Tfin , Tsin 0
(85)
Wsin ,Wfin 0,1
(86)
Uncertainty set for product demands, denoted as U, can be directly defined as a polytope or constructed from uncertainty data. The resulting optimization problem on batch process scheduling is a bi-objective multi-level mixed-integer programming problem. It explicitly
32
accounts for the conventional robustness and minimax regret criteria to benefit batch process operations.
6.2. Case study 1 In this subsection, we consider a specific case study on the short-term scheduling of multipurpose batch process [57]. The state-task network of this batch process is shown in Figure 17. As can be seen from Figure 17, this multipurpose batch process involves one heating task, three reaction tasks, and one separation task. This network involves three raw materials (A-C), four intermediates (Hot A, Impure E, and I1-I2) and two products (P1-P2). There are four equipment units: one heater, two reactors, and one separator. The mass balance coefficients are also listed in Figure 17. The scheduling horizon is 8 hours, and there are 6 time points.
The demands for Product 1 and Product 2 are subject to uncertainty. In this case study, the uncertainty set of product demands is directly given in Figure 18. This uncertainty set is a polytope with four extreme points (Points A-D). The coordinates of these extreme points are given as follows: Point A (40, 30), Point B (52, 28), Point C (24, 42) and Point D (12, 44).
By using the computational algorithm in Section 3.2.1, we solve the corresponding problem by changing λ from 0 to 1 with an increment of 0.05. The Gantt charts associated with λ=0 and λ=1 are demonstrated in Figure 19 (a) and (b), respectively. By comparing the two Gantt charts, we can observe that the assignment decisions determined by the two approaches are different. Specifically, Reaction 3 starts at a unique time point in the scheduling solution determined by the conventional robustness criterion, while Reaction 3 shares the same starting time point with Reaction 1 in the scheduling solution determined by the minimax regret criterion. The values of worst-case profit and regret are listed in Figure 19. The scheduling solution determined by the minimax regret criterion has the lowest worst-case regret of $13.6, and only decreases the worstcase profit by 0.3% compared with solution under the conventional robustness criterion. Due to the discrete nature of first-stage decisions and problem setup, there is no trade-off solution, when λ takes other values between 0 and 1.
33
6.3. Case study 2 In this subsection, another case study is considered to demonstrate that the solutions with the conventional robustness criterion and the minimax regret criterion could have the same worstcase profit and worst-case regret. This case study is originated from an industrial multipurpose batch process in The Dow Chemical Company [58-61]. Figure 20 displays the state-task network of the batch process in this case study. This multipurpose batch process has one preparation task, three reaction tasks, two packing tasks and two drumming tasks. The chemicals include four raw materials (M1-M4), six intermediates (I1-I6) and four products (P1-P4). The equipment units include one mixer, two reactors, one finishing system and one drumming line. The scheduling horizon is one week, i.e. 168 hours, and there are 11 time points.
Due to market fluctuations, product demands for P1-P4 are subject to uncertainties in this case study. A total of 200 uncertainty data are used for constructing the uncertainty set, and each data point has 4 dimensions for a combination of all final products. We employ the computational algorithm in Section 3.2.1 to solve the resulting scheduling problems with different values of parameter λ. The Gantt charts of the ARO solutions under the conventional robustness criterion and the minimax regret criterion are shown in Figure 21 (a) and (b), respectively. From the Gantt chart, we can observe that the assignment decisions are different, but the worst-case profit and regret remain unchanged for the two solutions. This is because of a different problem setup from Case study 1.
7. Conclusions In this paper, a novel ARO framework was proposed to account for the minimax regret criterion. Apart from the conventional robustness criterion, the worst-case regret was introduced as another objective to minimize the deviation from the perfect information solution. A set of Paretooptimal solutions was generated to reveal the systematic trade-offs between the conventional robustness and minimax regret criteria. The proposed optimization framework was general enough to apply to various kinds of uncertainty sets, such as the polyhedral and data-driven uncertainty sets. A data-driven version of this framework, named DDBCARO, leveraged the advantages of the minimax regret criterion, and at the same time exploited the value of big data via the Dirichlet process mixture model. Since the resulting optimization problem has multi-level 34
optimization structures, we further proposed the tailored C&CG algorithms to address the computational challenge. We presented a numerical example and applications on strategic planning of process networks and batch process scheduling. The results showed that the proposed framework could make a systematic trade-off between the conventional robustness criterion and minimax regret criterion.
Acknowledgments The authors acknowledge financial support from the National Science Foundation (NSF) CAREER Award (CBET-1643244).
Appendix A. Details of the Data-Driven Uncertainty Set In data-driven robust optimization paradigm, one paramount ingredient is the uncertainty set, which is learned directly from uncertainty data. To model uncertainty data with high fidelity, we adopt the Dirichlet process mixture model in this work. The Dirichlet process mixture model is a powerful nonparametric Bayesian statistical model, which has an infinite-dimensional parameter space. The resulting data-driven uncertainty set could adjust its complexity based on the data structure and complexity, thereby truthfully capturing the nature of uncertainty [39]. The Dirichlet process is defined as follows [62]. A random distribution F is said to be distributed as a Dirichlet process, denoted as F ~ DP (α, F0), if for any fixed number of partitions (Θ1, ..., Θr) of Θ, (F(Θ1),…, F(Θr)) obeys a Dirichlet distribution, i.e.,
F (1 ),
, F (r ) ~ Dir F0 (1 ), , F0 (r )
(A1)
where F0 is a base distribution on Θ, α is the concentration parameter that indicates the users’ belief in F0. Dir is the notation for the Dirichlet distribution. Besides, a random draw from the Dirichlet process has the following representation based on a stick-breaking construction [63].
F k 1 kk
(A2)
where the weight k k j 1 1 j , and βk is the proportion of stick being broken from the k 1
remaining stick. βk follows a Beta distribution, denoted as βk ~ Beta (1, α). The random variable
35
drawn from the base distribution F0 is denoted as θk. k is a delta function which is equal to 1 at θk, and 0 otherwise. The Dirichlet process is a distribution on distributions. A random draw from a Dirichlet process, i.e. DP (α, F0), is a distribution F. We can still have a random draw from the distribution F. The Dirichlet process mixture model adds one more level to this hierarchy. It utilizes θk as the parameters of the data distribution. As shown in eq. (A2), the Dirichlet process mixture model can be interpreted as a mixture model that has a potentially infinite number of mixture components. The general Dirichlet process mixture model is summarized as follows [64].
k ~ Beta(1, ) k F0 ~ F0 li 1 , 2 ,
~ Mult ( )
(A3)
oi li ~ p (oi li )
where Mult denotes a multinomial distribution, and li is the index of mixture components to which the observation oi is assigned. Variational inference algorithm is employed for the Dirichlet process mixture model [64, 65]. In Bayesian statistics, the model parameters are considered as random variables, and the inference is to obtain the posterior distribution of these latent random variables. The mixtures of Gaussian distributions are adopted in this work. θk = (ηk, Hk) includes mean vector and precision matrix. It is assumed that θk ~ NW (μk, λk, ωk, ψk), where NW represents the normal Wishart distribution. In the variational inference algorithm, the posterior distribution is approximated by a factorized variational distribution q (l, β, α, η, H), as given below q l, β, , η, H
N
M 1
M
i 1
k 1
k 1
q li q k q q ηk , Hk
(A4)
where M is the truncation level. Suppose Θ= (l, β, α, η, H) is set of all latent variables and parameters of the Dirichlet process mixture model, and U= (u1, …, uN) are the observations of uncertain parameters. The Dirichlet process mixture model, specifies the distribution p(U, Θ). From a Bayesian perspective, the goal of inference is to find the posterior distribution of Θ given the observation data U, that is p(Θ|U). In the variational inference, the log marginal probability is decomposed into two parts. log p U ELBO q KL q p
36
(A5)
where the first part is called the evidence lower bound (ELBO), and the second part is KullbackLeibler divergence or KL divergence; q is the variational distribution to approximate the posterior distribution p(Θ|U). They are defined as follows [65, 66]. p U, Θ ELBO q q Θ log dΘ q Θ
(A6)
p Θ U KL q p q Θ log dΘ q Θ
(A7)
Given the KL(q||p) is a distance between the distributions q and p, the goal of the variational algorithm is to minimize KL(q||p) to achieve the best approximation. Since log p(U)is a constant with respect to q, minimizing the KL divergence levels to maximizing the ELBO. To obtain a tractable inference, the variational inference utilize a factorized variational distribution q(l, β, α, η, H). q Gamma 1* , 2*
(A8)
where Gamma represents the Gamma distribution, and the hyper-parameters are updated by
1* 1 M 1
(A9)
2* 2 i 1 ln 1 i M 1
(A10)
The notation <.> represents the expectation of corresponding random variable. All closed forms of <.> are listed at the end of Appendix A. The posterior of βk is as follows. q k Beta k , vk
(A11)
where Beta represents the Beta distribution, and the hyper-parameters are updated by k* 1 n1 q ln k
(A12)
vk* n1 i k 1 q ln i
(A13)
N
N
M
The posterior of q(ηk, Hk)is the Normal-Wishart distribution, and the related parameters are calculated in (A15)-(A19).
q ηk , Hk N ηk μ*k , k*Ψ*k1 W H k Ψ*k1 , k* 1
Nk n1 q ln k N
37
(A14) (A15)
k* k Nk
(A16)
k* k Nk
(A17)
μ*k 1 Ψ Ψ k n 1 q ln k u n Nk N
* k
+
k N k 1 k N k N k
N n '1
μ 1
k
* k
k
N n 1
q ln k u n
1 n '1 q ln ' k un ' un N k N
1 q ln ' k u n ' μ k N k
q ln k
N n '1
(A18) n '1 q ln ' k un '
T
N
q ln ' k u n ' μ k
nk*
(A19)
T
(A20)
* i 1 ni M
where ρnk* is defined as follows. k 1 1 K 1 T nk* exp ln k i 1 ln 1 i ln H k ln 2 un ηk H k un ηk 2 2 2
(A21)
where K= dim (u) is the dimension of uncertainty data. The expressions involved < > are listed below.
un ηk
T
1* * 2
(A22)
ln k k k vk
(A23)
ln 1 k vk k vk
(A24)
* 1 i K * ln H k i 1 k K ln 2 ln Ψ k 2
(A25)
H k u n ηk
K
k
k u n μ*k Ψ*k1 u n μ*k T
where ψ is the digamma function. The variational inference algorithm is summarized below.
38
(A26)
References [1] [2]
[3] [4] [5]
[6]
[7]
[8]
[9] [10] [11] [12] [13] [14] [15]
[16]
[17] [18]
N. V. Sahinidis, "Optimization under uncertainty: State-of-the-art and opportunities," Computers & Chemical Engineering, vol. 28, pp. 971-983, 2004. I. E. Grossmann, R. M. Apap, B. A. Calfa, P. García-Herreros, and Q. Zhang, "Recent advances in mathematical programming techniques for the optimization of process systems under uncertainty," Computers & Chemical Engineering, vol. 91, pp. 3-14, 2016. E. N. Pistikopoulos, "Uncertainty in process design and operations," Computers & Chemical Engineering, vol. 19, pp. 553-563, 1995. J. M. Mulvey, R. J. Vanderbei, and S. A. Zenios, "Robust Optimization of Large-Scale Systems," Operations Research, vol. 43, pp. 264-281, 1995. Y. A. Guzman, L. R. Matthews, and C. A. Floudas, "New a priori and a posteriori probabilistic bounds for robust counterpart optimization: I. Unknown probability distributions," Computers & Chemical Engineering, vol. 84, pp. 568-598, 2016. K. Tong, F. You, and G. Rong, "Robust design and operations of hydrocarbon biofuel supply chain integrating with existing petroleum refineries considering unit cost objective," Computers & Chemical Engineering, vol. 68, pp. 128-139, 2014. A. C. S. Amaro and A. P. F. D. Barbosa-Póvoa, "The effect of uncertainty on the optimal closed-loop supply chain planning under different partnerships structure," Computers & Chemical Engineering, vol. 33, pp. 2144-2158, 2009. Y. Yuan, Z. Li, and B. Huang, "Robust optimization under correlated uncertainty: Formulations and computational study," Computers & Chemical Engineering, vol. 85, pp. 58-71, 2016. J. R. Birge and F. Louveaux, Introduction to stochastic programming: Springer Science & Business Media, 2011. A. Ben-Tal and A. Nemirovski, "Robust optimization – methodology and applications," Mathematical Programming, vol. 92, pp. 453-480, 2002. D. Bertsimas, D. B. Brown, and C. Caramanis, "Theory and applications of robust optimization," SIAM Review, vol. 53, pp. 464-501, 2011. Z. Li and Z. Li, "Optimal robust optimization approximation for chance constrained optimization problem," Computers & Chemical Engineering, vol. 74, pp. 89-99, 2015. K. Shang, Z. Feng, L. Ke, and F. T. S. Chan, "Comprehensive Pareto Efficiency in robust counterpart optimization," Computers & Chemical Engineering, vol. 94, pp. 75-91, 2016. N. H. Lappas and C. E. Gounaris, "Multi-stage adjustable robust optimization for process scheduling under uncertainty," AIChE Journal, vol. 62, pp. 1646-1667, 2016. D. E. Majewski, M. Wirtz, M. Lampe, and A. Bardow, "Robust multi-objective optimization for sustainable design of distributed energy supply systems," Computers & Chemical Engineering, vol. 102, pp. 26-39, 2017. C. Ning and F. You, "A data-driven multistage adaptive robust optimization framework for planning and scheduling under uncertainty," AIChE Journal, vol. 63, pp. 4343-4369, 2017. C. Shang, X. Huang, and F. You, "Data-Driven Robust Optimization Based on Kernel Learning," Computers & Chemical Engineering, 2017. A. Ben-Tal, A. Goryashko, E. Guslitzer, and A. Nemirovski, "Adjustable robust solutions of uncertain linear programs," Mathematical Programming, vol. 99, pp. 351-376, 2004.
39
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26] [27]
[28] [29] [30] [31]
[32]
[33] [34]
B. Zeng and L. Zhao, "Solving two-stage robust optimization problems using a columnand-constraint generation method," Operations Research Letters, vol. 41, pp. 457-461, 2013. Q. Zhang, M. F. Morari, I. E. Grossmann, A. Sundaramoorthy, and J. M. Pinto, "An adjustable robust optimization approach to scheduling of continuous industrial processes providing interruptible load," Computers & Chemical Engineering, vol. 86, pp. 106-119, 2016. H. Shi and F. You, "A computational framework and solution algorithms for two-stage adaptive robust scheduling of batch manufacturing processes under uncertainty," AIChE Journal, vol. 62, pp. 687-703, 2016. J. Gong and F. You, "Optimal processing network design under uncertainty for producing fuels and value-added bioproducts from microalgae: Two-stage adaptive robust mixed integer fractional programming model and computationally efficient solution algorithm," AIChE Journal, vol. 63, pp. 582-600, 2017. J. Gong, D. J. Garcia, and F. You, "Unraveling Optimal Biomass Processing Routes from Bioconversion Product and Process Networks under Uncertainty: An Adaptive Robust Optimization Approach," ACS Sustainable Chemistry & Engineering, vol. 4, pp. 31603173, 2016. D. E. Majewski, M. Lampe, P. Voll, and A. Bardow, "TRusT: A Two-stage Robustness Trade-off approach for the design of decentralized energy supply systems," Energy, vol. 118, pp. 590-599, 2017. D. Yue and F. You, "Optimal supply chain design and operations under multi-scale uncertainties: Nested stochastic robust optimization modeling framework and solution algorithm," AIChE Journal, vol. 62, pp. 3041-3055, 2016. D. E. Bell, "Regret in Decision Making under Uncertainty," Operations Research, vol. 30, pp. 961-981, 1982. T. Assavapokee, M. J. Realff, and J. C. Ammons, "Min-Max Regret Robust Optimization Approach on Interval Data Uncertainty," Journal of Optimization Theory and Applications, vol. 137, pp. 297-316, 2008. P. Kouvelis and G. Yu, Robust discrete optimization and its applications vol. 14: Springer Science & Business Media, 2013. I. R. Dunning, "Advances in robust and adaptive optimization: algorithms, software, and insights," Massachusetts Institute of Technology, 2016. Z.-L. Chen, N. G. Hall, and H. Kellerer, "Dynamic Pricing to Minimize Maximum Regret," Production and Operations Management, vol. 26, pp. 47-63, 2017. R. Yokoyama, K. Fujiwara, M. Ohkura, and T. Wakui, "A revised method for robust optimal design of energy supply systems based on minimax regret criterion," Energy Conversion and Management, vol. 84, pp. 196-208, 2014. B. Chen, J. Wang, L. Wang, Y. He, and Z. Wang, "Robust Optimization for Transmission Expansion Planning: Minimax Cost vs. Minimax Regret," IEEE Transactions on Power Systems, vol. 29, pp. 3069-3077, 2014. R. Jiang, J. Wang, M. Zhang, and Y. Guan, "Two-Stage Minimax Regret Robust Unit Commitment," IEEE Transactions on Power Systems, vol. 28, pp. 2271-2282, 2013. R. Loulou and A. Kanudia, "Minimax regret strategies for greenhouse gas abatement: methodology and application," Operations Research Letters, vol. 25, pp. 219-230, 1999.
40
[35] [36]
[37] [38] [39]
[40]
[41] [42] [43] [44] [45]
[46] [47]
[48]
[49]
[50]
[51] [52] [53]
J. Ide and A. Schöbel, "Robustness for uncertain multi-objective optimization: a survey and analysis of different concepts," OR Spectrum, vol. 38, pp. 235-271, 2016. L. Wang, Q. Li, R. Ding, M. Sun, and G. Wang, "Integrated scheduling of energy supply and demand in microgrids under uncertainty: A robust multi-objective optimization approach," Energy, vol. 130, pp. 1-14, 2017. J. Fliege and R. Werner, "Robust multiobjective optimization & applications in portfolio optimization," European Journal of Operational Research, vol. 234, pp. 422-433, 2014. D. Kuroiwa and G. M. Lee, "On robust multiobjective optimization," Vietnam J. Math, vol. 40, pp. 305-317, 2012. C. Ning and F. You, "Data-Driven Adaptive Nested Robust Optimization: General Modeling Framework and Efficient Computational Algorithm for Decision Making Under Uncertainty," AIChE Journal, vol. 63, pp. 3790-3817, 2017. H. Aissi, C. Bazgan, and D. Vanderpooten, "Min–max and min–max regret versions of combinatorial optimization problems: A survey," European Journal of Operational Research, vol. 197, pp. 427-438, 2009. T. Campbell and J. P. How, "Bayesian nonparametric set construction for robust optimization," in 2015 American Control Conference (ACC), 2015, pp. 4216-4221. G. P. Rangaiah and A. Bonilla-Petriciolet, Multi-objective optimization in chemical engineering: developments and applications: John Wiley & Sons, 2013. F. You and I. E. Grossmann, "Design of responsive supply chains under demand uncertainty," Computers & Chemical Engineering, vol. 32, pp. 3090-3111, 2008. A. Thiele, T. Terry, and M. Epelman, "Robust linear optimization with recourse," Rapport technique, pp. 4-37, 2009. D. Bertsimas, E. Litvinov, X. A. Sun, J. Zhao, and T. Zheng, "Adaptive Robust Optimization for the Security Constrained Unit Commitment Problem," IEEE Transactions on Power Systems, vol. 28, pp. 52-63, 2013. F. Glover, "Improved linear integer programming formulations of nonlinear integer problems," Management Science, vol. 22, pp. 455-460, 1975. N. V. Sahinidis, I. E. Grossmann, R. E. Fornari, and M. Chathrathi, "Optimization model for long range planning in the chemical industry," Computers & Chemical Engineering, vol. 13, pp. 1049-1063, 1989. F. You and I. E. Grossmann, "Stochastic Inventory Management for Tactical Process Planning Under Uncertainties: MINLP Models and Algorithms," AIChE Journal, vol. 57, pp. 1250-1277, 2011. D. Yue and F. You, "Planning and Scheduling of Flexible Process Networks Under Uncertainty with Stochastic Inventory: MINLP Models and Algorithm," AIChE Journal, vol. 59, pp. 1511-1532, 2013. A. A. Levis and L. G. Papageorgiou, "A hierarchical solution approach for multi-site capacity planning under uncertainty in the pharmaceutical industry," Computers & Chemical Engineering, vol. 28, pp. 707-725, 2004. S. Ahmed and N. V. Sahinidis, "Robust process planning under uncertainty," Industrial & Engineering Chemistry Research, vol. 37, pp. 1883-1892, 1998. M. L. Liu and N. V. Sahinidis, "Optimization in Process Planning under Uncertainty," Industrial & Engineering Chemistry Research, vol. 35, pp. 4154-4165, 1996. E. Rosenthal, GAMS-A user’s guide. Washington, DC: GAMS Development Corporation, 2008. 41
[54]
[55]
[56]
[57]
[58]
[59]
[60]
[61]
[62] [63] [64] [65] [66]
S. Moniz, A. P. Barbosa-Póvoa, and J. P. de Sousa, "Simultaneous regular and nonregular production scheduling of multipurpose batch plants: A real chemical– pharmaceutical case study," Computers & Chemical Engineering, vol. 67, pp. 83-102, 2014. Y. Chu and F. You, "Model-based integration of control and operations: Overview, challenges, advances, and opportunities," Computers & Chemical Engineering, vol. 83, pp. 2-20, 2015. E. Muñoz, E. Capón-García, M. Moreno-Benito, A. Espuña, and L. Puigjaner, "Scheduling and control decision-making under an integrated information environment," Computers & Chemical Engineering, vol. 35, pp. 774-786, 2011. E. Kondili, C. C. Pantelides, and R. W. H. Sargent, "A general algorithm for short-term scheduling of batch operations—I. MILP formulation," Computers & Chemical Engineering, vol. 17, pp. 211-227, 1993. Y. Chu, J. M. Wassick, and F. You, "Efficient scheduling method of complex batch processes with general network structure via agent-based modeling," AIChE Journal, vol. 59, pp. 2884-2906, 2013. J. M. Wassick, A. Agarwal, N. Akiya, J. Ferrio, S. Bury, and F. You, "Addressing the operational challenges in the development, manufacture, and supply of advanced materials and performance products," Computers & Chemical Engineering, vol. 47, pp. 157-169, 2012. Y. Chu, F. You, J. M. Wassick, and A. Agarwal, "Integrated planning and scheduling under production uncertainties: Bi-level model formulation and hybrid solution method," Computers & Chemical Engineering, vol. 72, pp. 255-272, 2015. Y. Chu, F. You, and J. M. Wassick, "Hybrid method integrating agent-based modeling and heuristic tree search for scheduling of complex batch processes," Computers & Chemical Engineering, vol. 60, pp. 277-296, 2014. W. Y. Teh, "Dirichlet Process," in Encyclopedia of Machine Learning, C. Sammut and G. I. Webb, Eds., ed Boston, MA: Springer US, 2010, pp. 280-287. J. Sethuraman, "A constructive definition of Dirichlet priors," Statistica Sinica, vol. 4, pp. 639-650, 1994. D. M. Blei and M. I. Jordan, "Variational inference for Dirichlet process mixtures," pp. 121-143, 2006. K. Kurihara, M. Welling, and N. A. Vlassis, "Accelerated variational Dirichlet process mixtures," in Advances in neural information processing systems, 2006, pp. 761-768. X. Wei and C. Li, "The infinite Student's t-mixture for robust modeling," Signal Processing, vol. 92, pp. 224-234, 2012.
42
Algorithm. Tailored C&CG algorithm for (DDBCARO-λ) 1: 2:
Set LB , UB , k 0 , and UB LB do while UB
;
Solve (MP-λ) to obtain x*k 1 ,1k 1* ,2k 1* ;
3:
Update LB 1 1k 1* 2k 1* ;
4: 5: 6:
for i 1 to m do Solve (SUP1,i) to obtain u1,k i 1 and Q1i x*k 1 ;
Solve (SUP2,i) to obtain uk2,i1 , xˆ ik 1 , yˆ ik 1 and Q2i x*k 1 ;
7: 8:
end
i* arg max Q1i x*k 1 ;
9:
i
u1k 1 u1,k i*1 and Q1 x*k 1 Q1i x*k 1 ; *
10: 11:
j * arg max Q2j x*k 1 ; j
, xˆ k 1 xˆ kj*1 , yˆ k 1 yˆ kj*1 and Q2 x*k 1 Q2j x*k 1 ;
12:
u
13:
Update UB min UB, cT x*k 1 1 Q1 x*k 1 Q2 x*k 1 ;
k 1 2
u
k 1 2, j*
*
k 1 1
Create second-stage variables y 14:
1 c x b y , Tx Wy T
T
k 1 1
2 c x b y T
T
k 1 2
c xˆ T
k 1 1
k 1
,y
k 1 2
, and add cuts
k 1 1
h Mu
b yˆ k 1 , Tx Wy k2 1 h Mu k2 1 T
to (MP-λ);
k k 1 ;
15: 16:
end
17:
return UB ;
Figure 1. The pseudocode of the proposed computational algorithm for (DDBCARO-λ).
43
Algorithm. Tailored C&CG algorithm for (DDBCARO-ε) 1: 2:
Set LB , UB , k 0 , and UB LB do while UB
;
3:
Solve (MP-ε) to obtain x*k 1 , k 1* ;
4: 5: 6:
Update LB ; for i 1 to m do Solve (SUP1,i) to obtain u1,k i 1 and Q1i x*k 1 ; k 1*
Solve (SUP2,i) to obtain uk2,i1 , xˆ ik 1 , yˆ ik 1 and Q2i x*k 1 ;
7: 8:
end
i* arg max Q1i x*k 1 ;
9:
i
u1k 1 u1,k i*1 and Q1 x*k 1 Q1i x*k 1 ; *
10: 11:
j * arg max Q2j x*k 1 ; j
12:
u
k 1 2
u
, xˆ k 1 xˆ kj*1 , yˆ k 1 yˆ kj*1 and Q2 x*k 1 Q2j x*k 1 ; *
if cT x*k 1 Q2 x*k 1 then
13:
Update UB min UB, cT x*k 1 Q1 x*k 1 ;
14: 15: 16:
k 1 2, j*
end Create second-stage variables y1k 1 , y k2 1 , and add cuts
cT x bT y1k 1 , Tx Wy1k 1 h Mu1k 1
cT x bT y k2 1 cT xˆ k 1 bT yˆ k 1 , Tx Wy k2 1 h Mu k2 1
to (MP-ε);
k k 1 ;
17: 18:
end
19:
return UB ;
Figure 2. The pseudocode of the proposed computational algorithm for (DDBCARO-ε).
44
Figure 3. The scatter plot of uncertainty data and the data-driven uncertainty set.
Figure 4. The Pareto-optimal curve in the numerical example.
45
Figure 5. The chemical process network in the case study.
46
Figure 6. The worst-case NPV and worst-case regret under different values of parameter λ.
Figure 7. Pareto-optimal curve of the case study with breakdown of the total cost.
47
Figure 8. The optimal design and planning decisions for time period 4 determined by the solution of Point A. The optimal production capacity of each process is displayed under each operating process. 48
Figure 9. The optimal design and planning decisions for time period 4 determined by the solution of Point B. The optimal production capacity of each process is displayed under each operating process. 49
Figure 10. The optimal design and planning decisions for time period 4 determined by the solution of Point C. The optimal production capacity of each process is displayed under each operating process. 50
Figure 11. Optimal capacity expansion decisions over the entire planning horizon for Point A in the case study.
Figure 12. Optimal capacity expansion decisions over the entire planning horizon for Point B in the case study.
51
Figure 13. Optimal capacity expansion decisions over the entire planning horizon for Point C in the case study.
Figure 14. Upper and lower bounds in each iteration of the tailored C&CG algorithm for optimization problem (BCARO-λ) with λ=0.
52
Figure 15. Upper and lower bounds in each
Figure 16. Pareto-optimal curves under different uncertainty budgets in the case study.
53
Figure 17. State-task network of the multipurpose batch process in case study 1.
Figure 18. Uncertainty set of product demands in case study 1.
54
Figure 19. Gantt charts determined by (a) the conventional robustness criterion, and (b) the minimax regret criterion in case study 1.
55
Figure 20. The state-task network of the multipurpose batch process in case study 2.
56
Figure 21. Gantt charts determined by (a) the conventional robustness criterion, and (b) the minimax regret criterion in case study 2.
57
Algorithm. Variational inference algorithm N Input uncertainty data set U un n1 ; 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13:
Set ζe and q(l, β, α, η, H); ELBOt q ELBOt 1 q e do while ELBOt 1 q Update q by the expression in (A8); for k=1 to M do Update q k by the expression in (A11); Update q ηk , H k by the expression in (A14); end for n=1 to N do Update q ln by the expression in (A20); end end N
M 1
M
i 1
k 1
k 1
Output q l, β, , η, H q li q k q q ηk , H k
Figure A1. The variational inference algorithm for the Dirichlet process mixture model.
58
Table 1. Results of the BCARO-λ approach in the numerical example. BCARO-λ (λ=0) *
BCARO-λ (λ=0.4) *
BCARO-λ (λ=1) *
Worst-case cost Worst-case regret Iterations Total CPU (s)
1,555.7 878.6 5 7
1,615.5 694.2 5 8
1,833.2 451.0 5 8
First-stage decision x
x2 56.3
x2 48.1
x2 48.9
x3 81.3
x3 70.2
x3 49.8
x1 62.4
x1 58.9
x1 58.1
* λ is the weight of the minimax regret criterion, and (1-λ) is the weight of the conventional robustness criterion in model (BCARO-λ).
Table 2. Comparisons of model and solution statistics for (BCARO-λ).
Bin. Var.
BCARO-λ (λ=0)* Master problem 190
Sub-problem
BCARO-λ (λ=0.2)* Master problem
Sub-problem
BCARO-λ (λ=1)* Master problem
Sub-problem
Cont. Var.
8,843
302 17,141
190
302
190
302
5,083
17,141
8,843
Constraints
11,569
48,749
17,141
6,681
48,749
11,569
Worst-case NPV ($MM)
544.0
505.6
464.3
48,749
Worst-case regret ($MM)
705.9
207.7
90.1
Total CPU (s)
180
214
409
Iterations
9
5
9
* λ is the weight of the conventional robustness criterion, and (1-λ) is the weight of the conventional robustness criterion in model (BCARO-λ).
59
S0098-1354(17)30345-9 https://doi.org/10.1016/j.compchemeng.2017.09.026 CACE 5906
To appear in:
Computers and Chemical Engineering
Received date: Revised date: Accepted date:
14-7-2017 27-9-2017 28-9-2017
Please cite this article as: Ning, Chao., & You, Fengqi., Data-driven adaptive robust optimization with minimax regret criterion: Multiobjective optimization framework and computational algorithm for planning and scheduling under uncertainty.Computers and Chemical Engineering https://doi.org/10.1016/j.compchemeng.2017.09.026 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Data-Driven Adaptive Robust Optimization with Minimax Regret Criterion: Multiobjective Optimization Framework and Computational Algorithm for Planning and Scheduling under Uncertainty
Chao Ning, Fengqi You*
Robert Frederick Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, New York 14853, USA *
Corresponding author. Phone: (607) 255-1162; Fax: (607) 255-9166; E-mail:
[email protected] September 28, 2017 Submitted to Computers & Chemical Engineering
Highlights
ARO framework that incorporates conventional robustness and minimax regret criteria Solution algorithms for bi-criterion multilevel mixed-integer programming problem Innovative applications on planning and scheduling under uncertainty Leveraging big data analytics for data-driven decision-making under uncertainty
Abstract Regret is defined as the deviation of objective value from the perfect information solution, and serves as an important evaluation metric for decision-making under uncertainty. This paper proposes a novel framework that effectively incorporates the minimax regret criterion into twostage adaptive robust optimization (ARO). In addition to the conventional robustness criterion, this ARO framework also simultaneously optimizes the worst-case regret to push the 1
performance of the resulting solution towards the utopia one under perfect information. By using a data-driven uncertainty set, we formulate a multiobjective ARO problem that generates a set of Pareto-optimal solutions to reveal the systematic trade-offs between the conventional robustness and minimax regret criteria. The resulting multi-level mixed-integer programming problem cannot be solved directly by any off-the-shelf optimization solvers, so we further propose tailored column-and-constraint generation algorithms to address the computational challenge. Two applications on process network planning and batch process scheduling are presented to demonstrate the applicability of the proposed framework and the efficiency of the proposed solution algorithms.
Key words: Adaptive robust optimization, minimax regret, big data, multiobjective optimization, planning and scheduling
Nomenclature DDBCARO Framework The main notations used in the DDBCARO modeling framework are listed below. Parameters b
vector of cost coefficients corresponding to second-stage decisions
c
vector of cost coefficients corresponding to first-stage decisions
u
vector of uncertainties
μk
vector of parameters in normal Wishart distribution
ηk
mean vector in the Dirichlet process mixture model of Gaussian distributions
Cw
worst-case cost
li
index of mixture components assigned to the observation oi
m
total number of basic uncertainty sets
M
truncation level in the Variational inference
M0
constant with a very large value
Rw
worst-case regret 2
RHw
optimal highest worst-case regret
RLw
lowest worst-case regret
U
uncertainty set
Ui
basic uncertainty set
oi
the observation generated from the Dirichlet process mixture model
α
concentration parameter in the Dirichlet process
λk
parameter in normal Wishart distribution
θk
random variable drawn from the base distribution
πk
parameter of weight in the Dirichlet process
ωi
parameter in normal Wishart distribution
βk
proportion of stick being broken from the remaining stick
Λi
scaling factor in the data-driven uncertainty set
i
a hyper-parameter for Beta distribution in the Variational inference
vi
a hyper-parameter for Beta distribution in the Variational inference
γi
weight of mixture component
γ*
threshold value for the weight of mixture component
Θ
support of the base distribution
F0
base distribution in the Dirichlet process
Φi
uncertainty budget in data-driven uncertainty set
A
coefficient matrix corresponding to the first-stage decisions
Hk
precision matrix in the Dirichlet process mixture model of Gaussian distributions
M
coefficient matrix corresponding to the uncertainties
T
technology matrix in the DDBCARO framework
W
recourse matrix in the DDBCARO framework
Ψk
parameter in the normal Wishart distribution Variables
x
vector of first-stage decisions 3
y
vector of second-stage decisions
Application to planning of process networks The sets, parameters, and variables used in the application part are summarized below. Sets/indices I
set of processes indexed by i
J
set of chemicals indexed by j
T
set of time periods indexed by t Parameters
cei
maximum number of expansions for process i over the planning horizon
dujt
demand of chemical j in time period t
sujt
supply of chemical j in time period t
c1it
variable investment cost for process i in time period t
c2it
fixed investment cost for process i in time period t
c3it
unit operating cost for process i in time period t
c4it
purchase price of chemical j in time period t
cbt
maximum allowable investment in time period t
qeitL
lower bound for capacity expansion of process i in time period t
qeitU
upper bound for capacity expansion of process i in time period t
vjt
sale price of chemical j in time period t
κij
mass balance coefficient for chemical j in process i Binary variables
Yit
binary variable that indicates whether process i is chosen for expansion in time period t Continuous variables
Pjt
purchase amount of chemical j in time period t
Qit
total capacity of process i in time period t
QEit
capacity expansion of process i in time period t
Sjt
sale amount of chemical j in time period t
Wit
operation level of process i in time period t 4
Application to batch process scheduling The sets, parameters, and variables used in the data-driven robust scheduling of batch processes are summarized below. Sets/indices I
set of tasks indexed by i
J
set of equipment units indexed by j
N
set of time points indexed by n
S
set of states indexed by s Parameters
bimin
lower bound of batch size of task i
bimax
upper bound of batch size of task i
cs
storage capacity for state s
ds
demand for state s
fci
fixed cost of task i
ht
time horizon
mis
mass balance coefficient for consumption or production of state s in task i
ps
price of state s
sints
initial amount of state s
sisi
binary parameter that indicates whether task i uses state s as input
snps
binary parameter that indicates whether state s is a product
sosi
binary parameter that indicates whether task i produces state s as output
vdi
variable processing time of task i
vci
variable cost of task i Binary variables
Wsin
binary variable that indicates whether task i starts at time point n
Wfin
binary variable that indicates whether task i finishes at or before time point n
5
Continuous variables Bsin
batch size of task i that starts at time point n
Bpin
batch size of task i being processed at time point n
Bfin
batch size of task i that finishes at or before time point n
BIisn
amount of state s used as input of task i at time point n
BOisn
amount of state s produced from task i at time point n
Din
processing time of task i that starts at time point n
SAs
amount of state s sold at time point n
STsn
inventory level of state s at time point n
Tn
time corresponding to time point n
Tfin
finish time of task i that starts at time point n
Tsin
start time of task i that starts at time point n
1. Introduction In recent years, decision making under uncertainty has received tremendous attention from both academia and industry [1-4]. As an inherent characteristic of process systems, uncertainty could render the deterministic solution suboptimal or even infeasible [5-7]. To this end, extensive research efforts have been made towards the development of systematic methodologies for handling uncertainty, such as stochastic programming approach and robust optimization [812]. Robust optimization has attracted considerable interest from the process systems engineering community due to its good balance between solution quality and computational tractability [13-17]. By introducing recourse actions, adaptive robust optimization (ARO) captures the sequential nature of decision-making processes, and is typically less conservative than static robust optimization [18, 19]. Therefore, it has a fruitful of applications, including process scheduling [20, 21], energy systems [22-24], and supply chain optimization [25]. Although conventional ARO has a number of attractive features, it does not account for an important evaluation metric, known as regret, in decision-making theory [26]. In practice, the solution of an optimization problem under the presence of uncertainty is usually evaluated ex post by comparing its performance with that of the perfect information solution [27]. The objective value of the perfect information solution is the best possible value that could be 6
achieved by assuming that the decision maker had known the exact uncertainty realization a priori. While the perfect information solution is a utopia in the sequential decision making under uncertainty, this solution serves as an ideal benchmark for the quality of the optimized solution. The minimax regret criterion holds the promise to significantly improve the aforementioned evaluation performance, since it aims for the least deviation from the ideal benchmark [28, 29]. Due to its practical relevance, the minimax regret criterion has demonstrated a broad array of successful applications, such as dynamic pricing [30], energy supply systems [31-33], and sustainability [34]. Additionally, a decision maker may hold different degrees of preference towards the regret metric. For instance, one might be willing to sacrifice the performance in minimax regret criterion for a better performance in the conventional robustness criterion or vice versa. Unfortunately, ARO models with a single criterion fail to meet this research need, because they are not capable to reveal the trade-offs between the conventional robustness and minimax regret criteria. Based on all the concerns above, it is imperative to develop an ARO framework that effectively accounts for the conventional robustness and minimax regret criteria. To fill the knowledge gap, we are facing a number of research challenges. The first challenge is how to formulate an ARO model that simultaneously optimizes the conventional robustness and minimax regret criteria, and also reveals the systematic trade-offs between them. Besides, the multi-level optimization structure, coupled with the bi-objective feature, makes the resulting optimization problem computationally challenging. Therefore, the second research challenge is to develop an efficient computational algorithm for the global optimization of the resulting ARO problem, which cannot be solved directly by any off-the-shelf optimization solvers. The third challenge lies in the innovative applications of the general framework to address various planning and scheduling problems. This paper proposes a general bi-criterion ARO (BCARO) framework that effectively accounts for the conventional robustness and minimax regret criteria. By minimizing the worstcase regret over all uncertainty realizations within an uncertainty set, the proposed framework pushes the performance of the resulting solution towards the utopia one. In this way, the evaluation performance against the ideal benchmark is significantly improved. Since the worstcase uncertainty realization is not uniquely defined for stochastic multiobjective optimization, different concepts of robust efficiency have been proposed [35]. In this work, we employ the concept of point-based robust efficiency due to its clear interpretation and wide applications [367
38]. The proposed BCARO framework generates a set of Pareto-optimal solutions to reveal the systematic trade-offs between the conventional robustness and minimax regret criteria. This framework is general enough to apply to various types of uncertainty sets, such as the polyhedral and data-driven uncertainty sets [11, 39]. A data-driven ARO model of this framework is developed based on the Dirichlet process mixture model [39]. The resulting optimization problem is formulated as a bi-criterion multi-level mixed-integer programming problem, which cannot be solved directly by any off-the-shelf optimization solvers. Therefore, we further propose tailored column-and-constraint generation (C&CG) algorithms for its efficient solution. To illustrate the applicability of the proposed framework and efficiency of the proposed solution algorithms, we present two applications on process network planning and batch process scheduling. The major novelties of this paper are summarized as follows.
A general BCARO framework that effectively accounts for the conventional robustness and minimax regret criteria, and reveals the trade-offs between them;
Tailored decomposition-based C&CG algorithms that efficiently solve the resulting bicriterion multi-level optimization problems;
Innovative applications of the general framework on strategic planning of process networks and multipurpose batch process scheduling under uncertainty.
The rest of this paper is organized as follows. Section 2 provides relevant background materials on conventional ARO and the minimax regret criterion. In Section 3, we present the general BCARO modeling framework and solution algorithms. Section 4 provides a numerical example. Two applications on process network planning and batch process scheduling are illustrated in Section 5 and Section 6, respectively. Conclusions are drawn in Section 7.
2. Background In this section, we briefly review relevant background on conventional ARO and the minimax regret criterion.
2.1. Conventional adaptive robust optimization The general form of a conventional two-stage adaptive robust mixed-integer programming model is given as follows [19]. 8
min max min cT x bT y x
uU y ( x ,u )
s.t. Ax d, x Rn1 Z n2
(1)
x, u y Rn3 : Wy h Tx Mu
where x is the first-stage decision made “here-and-now” before the uncertainty u is realized, while the second-stage decision or recourse decision y is postponed in a “wait-and-see” manner after the uncertainties are revealed. x includes both continuous and integer variables, while y only includes continuous variables. c and b are the vectors of the cost coefficients. U is an uncertainty set that characterizes the region in which uncertainty realizations reside. The worst-case cost associated with a certain first-stage decision x over all the uncertainty realizations in the uncertainty set U is given by [32], C w x max min cT x bT y
(2)
uU y ( x ,u )
With the definition of worst-case cost in (2), the conventional ARO problem in (1) can be interpreted as an optimization problem that aims to find a first-stage decision x with the lowest worst-case cost. Note that the conventional robustness criterion in (1) is also known as the minimax cost criterion. After observing the uncertainty realization, the decision maker takes the “wait-and-see” decision y based on the fully adaptive decision rule, which is determined by the following optimization problem. min bT y y
s.t. y Rn3
(3)
Wy h Tx* Mu
where u is the observed uncertainty realization, and x* is the optimal first-stage decision in the conventional ARO problem. Note that x* in (3) should be treated as a parameter since it has been made at the initial stage of the decision-making process.
2.2. Adaptive robust optimization under the minimax regret criterion As mentioned in Section 1, the perfect information solution could serve as an ideal benchmark for the quality of an optimized solution [27]. The optimization problem under perfect information of uncertainty realization is presented in (4).
9
u min cT xˆ bT yˆ xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2
(4)
Wyˆ Txˆ h Mu, yˆ Rn3
where Π(u) is the best objective value that could be achieved by assuming that the decision maker had known the uncertainty realization a priori. Note that, in optimization problem (4), all the decisions are made in a “wait-and-see” mode with full information of uncertainty realization. Given a first-stage decision x, the regret under a specific uncertainty realization u is defined by, R x, u min cT x bT y u
(5)
y ( x ,u )
where R(x, u) represents the regret. It characterizes the deviation, in terms of total costs, from the perfect information solution under the uncertainty realization u. Under the minimax regret criterion, the decision maker aims to minimize the worst-case regret [40], which is shown below: R w x max R x, u uU
max min cT x bT y u y ( x,u ) uU
(6)
max min cT x bT y u uU y ( x ,u )
where Rw(x) denotes the worst-case regret over all uncertainty realizations within the uncertainty set U. Based on the above definitions, a general two-stage adaptive robust mixed-integer programming model with the minimax regret objective is then given by, min max min cT x bT y u x uU y ( x ,u ) s.t. Ax d, x Rn1 Z n2
x, u y Rn3 : Wy h Tx Mu u min c xˆ b yˆ T
(7)
T
xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2 Wyˆ Txˆ h Mu, yˆ Rn3
The optimization problem in (7) focuses on the relative robustness with respect to the perfect information solution, while the optimization problem in (1) aims for the abosulte robustness. Both optimization problems have multi-level optimization structures. However, the 10
ARO problem under the minimax regret criterion has one more level due to the innermost perfect information optimization problem. After uncertainty is realized, the “wait-and-see” decision y is made based on the same optimization problem in (3). This is because u is a constant at the second stage of the decision-making process. Based on the above conventional ARO and minimax regret criterion, the proposed BCARO framework is developed in the next section.
3. Bi-Criterion Adaptive Robust Optimization Framework In this section, we propose a novel framework that incorporates the minimax regret criterion into conventional ARO through the lens of multiobjective optimization. Notably, the proposed BCARO framework is general enough to apply to various uncertainty sets of conventional ARO. For the ease of exposition, the BCARO framework is developed with a data-driven uncertainty set [39].
3.1. Bi-criterion adaptive robust optimization model To model uncertainty, we adopt a data-driven uncertainty set using both l1 and l∞ norms [39], as presented in (8) and (9). U
U i U1
U2
(8)
Um
i: i *
Ui u u μi i Ψ1/2 z i i z,
1, z 1 i
(9)
where U is the data-driven uncertainty set, Ui is the basic uncertainty set corresponding to the i th component, Λi is a scaling factor, γi is the weight of the i th mixture component,
i
i
i
i 1
v i
j 1
vj j
vj
and M 1 i 1 i . The value of γi indicates the probability of the M 1
component i, and γ* is a threshold value. m is the total number of mixture components having probabilities larger than γ*. i
i 1
i i 1 dim u
[41], and τi, νi, μi, λi, ωi, ψi are the
variational inference results. Details of the data-driven uncertainty set are presented in the Appendix.
11
We further propose a general ARO framework that accounts for the conventional robustness and minimax regret criteria simultaneously with the data-driven uncertainty set. The proposed data-driven bi-criterion adaptive robust mixed-integer programming model, denoted as (DDBCARO), is presented as follows. min Worst-case cost C w x x
min Worst-case regret R w x x
s.t. Ax d, x Rn1 Z n2 C w x max max min cT x bT y i1, , m uU i y ( x ,u )
(DDBCARO)
R w x max max min cT x bT y u i1, , m uU i y ( x ,u )
x, u y Rn3 : Wy h Tx Mu u min cT xˆ bT yˆ xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2 Wyˆ Txˆ h Mu, yˆ Rn3
where Ui is the data-driven basic uncertainty set defined in (9). (DDBCARO) optimizes the objectives of the worst-case cost Cw(x) and the worst-case regret Rw(x) simultaneously. By explicitly incorporating the minimax regret criterion, it pushes the performance of the resulting solution towards the utopia one under perfect information of uncertainty realization. One way to handle multiobjective optimization is the weighted sum method. Essentially, the method attaches weights to different objective functions, and eventually optimizes over a weighed sum of all the objectives. Following the weighted sum approach [42], we further propose a unified optimization model (DDBCARO-λ) shown as follows.
12
min 1 C w x R w x x s.t. Ax d, x Rn1 Z n2 C w x max max min cT x bT y i1, , m uU i y ( x ,u ) R w x max max min cT x bT y u i1, , m uU i y ( x ,u )
(DDBCARO-λ)
x, u y Rn3 : Wy h Tx Mu u min cT xˆ bT yˆ xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2 Wyˆ Txˆ h Mu, yˆ Rn3
where λ is a parameter ranging from 0 to 1. (DDBCARO-λ) could provide a wide spectrum of solutions by changing the value of λ from 0 to 1. Specifically, it reduces to the conventional ARO model when λ equals to 0, and it becomes the ARO with minimax regret criterion when λ equals to 1. Thus, the proposed modeling framework (DDBCARO-λ) unifies the ARO models in Section 2.1 and Section 2.2. If the value of λ is chosen between 0 and 1, the proposed model makes a trade-off between the conventional robustness and minimax regret criteria. Another popular method for multiobjective optimization is the ε-constraint method that holds one of the objective functions as a constraint. Following the ε-constraint method [42], the objective on the worst-case regret is transformed into an extra constraint that specifies the upper bound of Rw(x). The resulting optimization model (DDBCARO-ε) is presented by, min C w x x
s.t. R w x Ax d, x Rn1 Z n2 C w x max max min cT x bT y i1, , m uU i y ( x ,u )
(DDBCARO-ε)
R w x max max min cT x bT y u i1, , m uU i y ( x ,u )
x, u y Rn3 : Wy h Tx Mu u min c xˆ b yˆ T
T
xˆ , yˆ
s.t. Axˆ d, xˆ Rn1 Z n2 Wyˆ Txˆ h Mu, yˆ Rn3
13
Similarly, (DDBCARO-ε) can reveal the trade-offs between the conventional robustness and minimax regret criteria by following three steps of the ɛ-constraint method [43]. The first step is to minimize the worst-case regret subject to the same constraints in (DDBCARO), i.e. without the ɛ-constraint, to obtain the lowest worst-case regret RLw . The second one is to minimize the worst-case cost subject to the same constraints in (DDBCARO) to obtain its optimal solution x*, which in turn yields the optimal highest worst-case regret RHw R w x* . The last step is to divide the interval between RLw and RHw into small intervals, and solve problem (DDBCARO-ɛ) by fixing ɛ to the end points of each small interval. The resulting optimization problems (DDBCARO-λ) and (DDBCARO-ε) are both multilevel mixed-integer programming problems. They cannot be solved directly by any off-the-shelf optimization solvers due to the multi-level optimization structures. To this end, computational algorithms are further developed to efficiently solve the optimization problems.
3.2. Computational algorithm In this section, we address the computational challenge by developing two efficient solution algorithms for global optimization of problems (DDBCARO-λ) and (DDBCARO-ε) in Section 3.2.1 and 3.2.2, respectively. To facilitate the solution procedure, we first reformulate the worst-case regret to reduce the number of levels. R w x max max min cT x bT y u i1, , m uU i y ( x ,u ) max max min cT x bT y min cT xˆ bT yˆ xˆ , yˆ u i1, , m uU i y ( x ,u ) cT x max max min bT y cT xˆ bT yˆ i1, , m xˆ , yˆ u ,uU i y ( x ,u )
(10)
where u xˆ Rn Z n , y Rn : Axˆ d, Wyˆ Txˆ h Mu . 1
2
3
3.2.1. Solution algorithm for (DDBCARO-λ) Based on the above reformulation of the worst-case regret, (DDBCARO-λ) can be transformed into (11)-(15) using the epigraph reformulation. min 1 1 2
(11)
x ,1 ,2
14
s.t. Ax d, x Rn1 Z n2
(12)
1 cT x max max min bT y i1, , m uU y ( x ,u )
(13)
i
2 cT x max
min bT y cT xˆ bT yˆ
max
i1, , m xˆ , yˆ u ,uUi y ( x ,u )
x, u y Rn3 : Wy h Tx Mu
(14) (15)
where η1 and η2 are the auxiliary variables. Note the reformulated problem above is still a multilevel optimization problem due to the max-max-min structures. This multi-level optimization problem can be reformulated to an equivalent single-level optimization problem by enumerating all the extreme points, and introducing their corresponding recourse variables [19]. However, this single-level optimization problem could be a large-scale optimization problem due to the potentially large number of extreme points. Therefore, we partially enumerate the extreme points, and construct the master problem (MP-λ) shown below. min 1 1 2
x ,1 ,2 , y1l , y l2
s.t. Ax d
1 cT x bT y1l , l 1, , r (MP-λ)
Tx Wy1l h Mu1l , l 1,
,r
Tx Wy l2 h Mul2 , l 1,
,r
2 cT x bT y l2 cT xˆ l bT yˆ l , l 1, , r x R Z , y , y R , l 1, n1
n2
l 1
n3
l 2
,r
where r is the current iteration number. Since the master problem (MP-λ) contains only a subset of the extreme points, it provides a relaxation of the original problem. To enumerate the important uncertainty realizations for both (13) and (14) on-the-fly, two groups of sub-problems are solved in each iteration. The sub-problem corresponding to (13), denoted as (SUP1,i), is presented by Q1i x max min bT y1 u1U i
y1
s.t. Wy1 h Tx Mu1 y1 Rn3
where i is the index of basic uncertainty sets. The optimal u1 is the identified uncertainty realization that should be added to the (MP-λ).
15
(16)
Note that (16) is a max-min optimization problem, which can be transformed to a singlelevel optimization problem by employing strong duality or KKT conditions [19, 44]. The resulting problem after taking the dual of inner problem is shown in (17). T Q1i x max h Tx Mμi φ i MΨ1/2 i i tj z j φt z ,φ t j T s.t. W φ b
(17)
φ0 z
1, z 1 i
where φ is the vector of dual variables associated with the constraints Wy h Tx Mu , φt represents the t th component in vector φ, and zj is the j th component in vector z. Since there are bilinear terms zjφt in the objective function of (17), we need to further reformulate the sub-problem using a linearization technique. To facilitate the derivation, zj can be divided into two parts z j z j z j , where z j and z j are two non-negative variables. Thus, we have u μi i Ψ1/2 i i z z
(18)
The uncertainty budget parameter Φi is typically set as an integer value due to its physical meaning, as reported in the existing literature [44, 45]. For the integer uncertainty budget, the optimal z j and z j take the values of 0 or 1, and these two variables can be restricted to binary
variables. By using Glover’s linearization for the bilinear terms G tj z j φt and G tj z j φt [46], which are the products of a binary variable z j (and z j ) and a continuous variable φt, we can
reformulate (17) into the following optimization problem in (19).
16
Q1i x
h Tx Mμ T φ Tr MΨ1/2 T G G i i i i s.t. WT φ b max
φ , z , z ,G ,G
φ0 0 G φ eT 0 G φ eT G e M 0 z
T
G e M 0 z
T
e M e z
G φ eT e M 0 e z
T
G φ eT
T
0
eT z z i z z e,
z , z 0,1
K
(19)
where e and e both represent column vectors whose elements are all ones, and they have different dimensions, i.e., dim e dim u K and dim e dim φ . The t th row and j th column elements of matrices G and G are denoted as G tj and G tj , respectively. The notation Tr is the trace of a matrix, e.g. for matrix Q qij R nn , Tr Q in1 qii . Note that the vector or matrix inequalities should be interpreted elementwise. The sub-problem corresponding to (14), denoted as (SUP2,i), is presented by Q2i x
max
xˆ , yˆ u 2 ,u 2 U i
min bT y 2 cT xˆ bT yˆ y2
s.t. Wy 2 h Tx Mu 2 y 2 Rn3
Similarly, we can reformulate optimization problem (20) as follows
17
(20)
h Tx Mμ T φ Tr MΨ1/2 T G G cT xˆ bT yˆ Q2i x max i i i i φ,z ,z G ,G , xˆ , yˆ
s.t. WT φ b φ0 0 G φ eT 0 G φ eT G e M 0 z
T
G e M 0 z
T
e M e z
G φ eT e M 0 e z
T
G φ eT
T
0
eT z z i z z e,
z , z 0,1
K
Axˆ d Wyˆ Txˆ h Mμ i i MΨ1/2 i i z z
(21)
xˆ Rn1 Z n2 , yˆ Rn3
The proposed solution algorithm iteratively solves a master problem and two groups of subproblems until the relative optimality gap reduces to a predefined tolerance. The pseudocode of the algorithm is given in Figure 1, where ζ represents the predefined tolerance, and m represents the total number of basic uncertainty sets. As the optimality and feasibility cuts are added to the master problem at each iteration, the lower bound provided by (MP-λ) becomes tighter as the iteration moves on. The computational algorithm terminates in finite iterations because of the finite number of extreme points.
18
3.2.2. Solution algorithm for (DDBCARO-ε) We employ the techniques developed in the previous subsection to address problem (DDBCARO-ε), which can be first reformulated into (22) using the epigraph reformulation. min x ,
s.t. Ax d, x Rn1 Z n2
cT x max max min bT y i1, , m uU y ( x ,u ) i
cT x max
(22)
min b y c xˆ b yˆ T
max
i1, , m xˆ , yˆ u ,uU i
T
T
y ( x ,u )
x, u y Rn3 : Wy h Tx Mu
By enumerating a subset of the extreme points, we arrive at its corresponding master problem (MP-ε), which is shown below. min
x , , y1l , y l2
s.t. Ax d
cT x bT y1l , l 1, , r (MP-ε)
Tx Wy1l h Mu1l , l 1,
,r
Tx Wy l2 h Mul2 , l 1,
,r
cT x bT y l2 cT xˆ l bT yˆ l , l 1, , r x R Z , y , y R , l 1, n1
n2
l 1
n3
l 2
,r
Note that the sub-problems for (DDBCARO-ε) are the same as in the previous subsection. The pseudocode of the proposed solution algorithm is directly provided in Figure 2. To handle the ɛ-constraint, the solution algorithm in Figure 2 incorporates additional steps to check feasibility. From Steps 13-15, we can see that it updates the upper bound only when the ɛconstraint is feasible.
19
4. Numerical Example In this section, a small numerical example of the proposed BCARO framework, as shown in (23), is presented to illustrate the conventional robustness and minimax regret criteria and the trade-off between them. We consider the following numerical example: min max min x
uU y ( x ,u )
3x1 5 x2 12 x3 15y1 30 y2 22 y3
min max min 3 x1 5 x2 12 x3 15y1 30 y2 22 y3 u x uU y ( x ,u ) s.t. x1 x2 x3 200 xi 0, i 1, 2,3
(x, u) y Rn3 : yi ui xi , i 1, 2,3
(23)
u min 3xˆ1 5 xˆ2 12 xˆ3 15yˆ 1 30 yˆ 2 22 yˆ3 xˆ , yˆ
s.t. xˆ1 xˆ2 xˆ3 200 xˆi yˆi ui , i 1, 2,3 xˆi , yˆi 0, i 1, 2,3
where x1, x2 and x3 are the first-stage or “here-and-now” decision variables, u1, u2 and u3 are the uncertain parameters, and y1, y2 and y3 are the second-stage or “wait-and-see” decision variables. Note that the “here-and-now” decisions are made prior to the resolution of uncertainties, whereas “wait-and-see” decisions could be made after the uncertainty realizations. U is the data-driven uncertainty set. Π(u) is the objective value of the perfect information solution. Problem (23) is an instance of (BCARO). By using the weighted sum method, it can be transformed into optimization problem (BCARO-λ), where λ=0 corresponds to the conventional robustness criterion, and λ=1 corresponds to the minimax regret criterion. A trade-off can be achieved by choosing a value of λ between 0 and 1. The scatter plot of uncertainty data and the data-driven uncertainty set are shown in Figure 3, where the red dots represent uncertainty realizations of u1, u2 and u3, and the union of the
20
green polytopes is the data-driven uncertainty set. A total of 500 uncertainty data points is used for uncertainty modeling.
The optimization models and the proposed tailored C&CG algorithm are coded in GAMS 24.7.3. The computational experiments are implemented on a computer with an Intel (R) Core (TM) i7-6700 CPU @ 3.40 GHz and 32 GB RAM. The relative optimality tolerance for the tailored C&CG algorithm is set to be 0.01%. The numerical results are summarized in Table 1. We can observe that the conventional robustness criterion generates the lowest worst-case cost, and has the highest worst-case regret. Since the minimax regret criterion only focuses on the regret metric, it reduces the worst-case regret by 48.7% compared with the conventional robustness criterion. However, it suffers from the highest worst-case cost of 1,833.2. A trade-off between the conventional robustness and minimax regret criteria can be made when choosing λ to be 0.4 in (BCARO-λ). In the two-stage ARO setting, the cost is calculated as a sum of the first-stage cost and the second-stage cost. The costs can be different for different uncertainty realizations. The first-stage cost is 3x1+5x2+12x3, which does not depend on the uncertainty realization. Hence, the first-stage cost can be calculated based on the values of the first-stage decisions provided in Table 1. The second-stage cost is 15y1+30y2+22y3, which depends on the uncertainty realizations. To obtain the values of y1, y2 and y3, the following optimization problem in (24) is solved for each uncertainty realization. min 15y1 30 y2 22 y3 y
s.t. yi ui xi , i 1, 2,3
(24)
yi 0, i 1, 2,3
For each uncertainty realization, the cost can then be calculated as 3x1+5x2+12x3+15y1+30y2+22y3. We use an uncertainty realization in Figure 3 to illustrate the concept of regret. For this uncertainty realization, u1=55.2, u2=41.8, u3=74.4. The perfect information solution corresponding to this uncertainty realization is xˆ1 55.2, xˆ2 41.8, xˆ3 74.4, yˆ1 0, yˆ2 0, yˆ3 0 , and the corresponding cost is 1,267.4. The (BCARO-λ) solution at λ=1.0 is x1=58.1, x2=48.9, x3=49.8, y1=0, y2=0, y3=24.6, and the corresponding cost is 1,557.6. The regret is the difference between these two costs, which is 290.2. 21
To further demonstrate the trade-off between the conventional robustness and minimax regret criteria, we present the Pareto-optimal curve in Figure 4. The X-axis represents the worstcase regret, while the Y-axis denotes the worst-case cost. The Pareto-optimal curve is plotted following three steps of the ɛ-constraint method [43]. First, we solve problem (BCARO-λ) at λ=1 to obtain the lowest worst-case regret RLw . Next, we solve problem (BCARO-λ) at λ=0 to obtain the optimal solution x*, which in turn yields the optimal highest worst-case regret RHw R w x* . The last step is dividing the interval between RLw and RHw into small intervals, and solving problem (BCARO-ɛ) by fixing ɛ to the end points of each small interval. In Figure 4, the yellow dot represents a wise trade-off solution. Compared with the conventional robustness criterion, it greatly reduces the worst-case regret by 21.1% at the expense of a tiny increase (3.9%) in worstcase cost. The result illustrates the advantage of the proposed framework in revealing a systematic trade-off between the conventional robustness and minimax regret criteria.
5. Application to Strategic Planning of Chemical Process Networks under Supply and Demand Uncertainties To demonstrate the application of the proposed modeling framework and computational algorithms, we consider a case study on a multi-period process network planning under uncertainty. Large-scale chemical complexes often consist of interconnected processes and various chemicals [47-49]. These interconnections allow for the synergies between different processes in chemical production [49, 50]. The chemicals in a process network include feedstocks, intermediates and products. The process network planning problem determines the purchase levels of feedstocks, sales of products, capacity expansions, and production profiles of processes at each time period over the entire planning horizon [51]. In this application, the supply of feedstocks and demand of products are uncertain.
5.1. Bi-criterion adaptive robust planning of process networks model The bi-criterion ARO model for process network planning under uncertainty is formulated as follows. The two objectives are to maximize the worst-case net present value (NPV) and to 22
minimize the worst-case regret given in (25) and (26). The definition of worst-case NPV and worst-case regret are presented in (27)-(28). The constraints can be classified into capacity expansion constraints (29)-(30), budget constraints (31)-(32), production level constraint (33), mass balance constraint (34), supply and demand constraints (35)-(36), non-negativity constraints (37)-(38), and integrity constraint (39). A list of indices/sets, parameters and variables is given in the Nomenclature section, where all the parameters are denoted in lowercase symbols, and all the variables are denoted in upper-case symbols. The objective functions are defined by (25) and (26), as shown below max NPV w Qit , QEit , Yit
(25)
min R w Qit , QEit , Yit
(26)
Qit ,QEit ,Yit
Qit ,QEit ,Yit
where Qit is a decision on total capacity of process i in time period t, QEit is a decision variable for capacity expansion level of process i in time period t, Yit is a binary variable to determine on whether process i is expanded in time period t, NPVw represents the worst-case NPV and Rw denotes the worst-case regret. The worst-case NPV and the worst-case regret are presented in (27) and (28), respectively. NPV w mindem max jt S jt c1it QEit c2it Yit c3it Wit c 4 jt Pjt (27) du jt U , Pjt ,Qit , iI tT iI tT iI tT jJ tT jJ tT su U sup S ,W jt
jt
it
R w maxdem min du jt , su jt jt S jt c1it QEit c2it Yit c3it Wit c 4 jt Pjt du jt U , Pjt ,Qit , iI tT iI tT iI tT jJ tT jJ tT su U sup S ,W jt
jt
it
(28) where Sjt is the amount of chemical j sold to the market in time period t, Wit represents the operating level of process i during time period t, and Pjt is the amount of chemical j purchased in time period t. Parameters c1it, c2it, c3it and c4jt are coefficients for variable investment cost, fixed investment cost, operating cost and purchase cost, respectively. νjt is the sale price of chemical j in time period t. dujt is the demand of chemical j in time period t and sujt is the market availability of chemical j in time period t. Π(dujt, sujt) is the NPV of the perfection information solution. Constraint (29) enforces the upper and lower bounds of capacity expansions. Specifically, if Yit=1, i.e. process i is chosen to be expanded in time period t, the expanded capacity should be within the range [qeitL, qeitU]. 23
qeitL Yit QEit qeitU Yit ,
i, t
(29)
Constraint (30) specifies the updates of total capacity of each process i. Qit Qit 1 QEit ,
i, t
(30)
Constraints (31) and (32) specify the maximum number of capacity expansions for each process and investment budget for each time period, respectively.
Y
it
cei ,
i
(31)
t
c1
it
QEit c2it Yit cbt ,
t
(32)
i
where cei is the maximum expansion times for process i, and cbt is the investment budget for time period t. Constraint (33) enforces that the production level of a process cannot exceed its total capacity.
Wit Qit ,
i, t
(33)
Constraint (34) specifies the mass balance for chemicals. To be more specific, the purchase amount plus the production amount is equal to the total amount sold to the market and consumed by processes. Pjt ij Wit S jt 0,
j , t
(34)
i
where κij is the mass balance coefficient for chemical j in process i. Constraint (35) enforces that the purchase amount of a feedstock should not be higher than its available amount provided by the supplier. Constraint (36) specifies that the sale amount of a product cannot exceed its market demand. Pjt su jt ,
j, t
(35)
S jt du jt ,
s, j, t
(36)
Constraints (37)-(38) indicate the non-negativity of continuous decisions. QEit 0, Qit 0 i, t
(37)
Pjt , S jt , Wit 0,
(38)
j, t
Constraint (39) ensures that Yit is a binary variable. Yit 0,1 ,
i, t
(39)
24
The data-driven uncertainty sets for demand and supply uncertainties is shown in (40)-(41). U dem k : kdem *
du du μ
dem k
kdem Ψ dem dem z, z k k
1, z 1 dem k
(40)
dem where μ dem , Ψ dem are the inference results learned from uncertain product demand data k , k k
using the Variational inference algorithm. U sup k : ksup *
su su μ
sup k
sup ksup Ψsup z k k z,
1, z 1 sup k
(41)
sup sup where μsup are the inference results learned from uncertain supply data employing k , k , Ψk
the Variational inference algorithm. The resulting optimization problem on strategic planning of process network under uncertainty is a bi-objective multi-level mixed-integer programming problem. Due to the multilevel optimization structure, it cannot be solved directly by any off-the-shelf optimization solvers. We employ the proposed computational algorithms to efficiently solve the resulting large-scale optimization problems.
5.2. Case study overview In this subsection, a specific case study is presented to demonstrate the applicability of the proposed approach. The process network in this case study consists of 28 chemicals, 38 processes, 10 suppliers and 16 markets [48], as shown in Figure 5. Chemicals A-J represent feedstocks, which can be purchased from suppliers or manufactured by certain processes. Chemicals K-Z are products, which are sold to the markets. The intermediates are denoted as AA and AB. This complex process network has such flexibility that many manufacturing options are available. For example, Chemical E can be manufactured by Process 22, Process 25 or purchased from a supplier. In this case study, we consider 5 time periods over the planning horizon, and the duration of each time period is 2 years. It is assumed that processes 12, 13, 16 and 38 have initial capacities of 40, 25, 300 and 200 kt/y, respectively, at the beginning of time period 1. These processes are not allowed to be expanded until time period 2. The other chemical processes can be installed at the beginning of time period 1. Due to market fluctuations, availabilities of raw materials and demands of final products are subject to uncertainties in the process network planning problem [52]. In this case study, a total of 26 parameters, including supply of all feedstocks A-J and demands for all products K-Z are 25
uncertain. For the supply uncertainty, a total of 2,000 uncertainty data are used. Each data point has 10 dimensions for a combination of all feedstocks. For the demand uncertainty, 3,000 uncertainty data points are used, and each data point is for a combination of all the 16 products.
5.3. Results and discussions We solve the multiobjective multi-level optimization problem of process network planning using the computational algorithms proposed in Section 3.2. The optimization methods are modeled in GAMS 24.7.3 [53], and are implemented on a computer with an Intel (R) Core (TM) i7-6700 CPU @ 3.40 GHz and 32 GB RAM. The master problem and sub-problems are solved using CPLEX 12 with an optimality tolerance of 0. The relative optimality tolerance for the tailored C&CG algorithm is set to be 0.01%. The uncertainty budgets for supply and demand are set to be 3 and 5, respectively. A sensitivity analysis on the values of the uncertainty budgets is performed at the end of this subsection. We first solve optimization problem (BCARO-λ) using the computational algorithm in Section 3.2.1. To show how the two objective values change with parameter λ, we calculate the worst-case NPV and the worst-case regret under different choices of parameter λ and present them in Figure 6. The X axis represents parameter λ ranging from 0 to 1, and the Y axis represents the worst-case NPV or worst-case regret. As can be seen from Figure 6, when parameter λ increases, both the worst-case NPV and regret decreases. The proposed BCARO framework provides a wide spectrum of solutions, including the ARO solutions under the conventional robustness criterion (λ=0) and the minimax regret criterion (λ=1). Based on the decision maker’ preference to conventional robustness and minimax regret criteria, a wise-choice solution can be obtained by selecting the value of λ between 0 and 1. We select three instances of (BCARO-λ) corresponding to λ=0, λ=0.2 and λ=1, and summarize the computational results and model statistics in Table 2. Since the optimality cuts are added to the mater problem at each iteration, its problem size increases as the computational algorithm iterates. The master problem at the first iteration involves 190 integer variables, 1,323 26
continuous variables, and 1,793 constraints. The size of the master problem at the last iteration of the computational algorithm is listed in Table 2. Compared with the ARO problem under the conventional robustness criterion, the ARO problem under the minimax regret criterion costs a much longer computational time (409 seconds v.s. 180 seconds). Notably, all the three instances can be solved to the specific optimality tolerance within only 7 minutes.
To further investigate the trade-off between these two robustness criteria, we employ the εconstraint based computational algorithm in Section 3.2.2, and demonstrate the Pareto-optimal curve with cost breakdown in Figure 7. The Pareto-optimal curve is plotted following the same procedure as described in Section 4. On the Pareto-optimal curve in Figure 7, three points are selected for further discussions. Point A represents the ARO solution under the conventional robustness criterion. It has the highest worst-case NPV of $544.0 MM, and the largest worst-case regret of $705.9 MM. In contrast to Point A, Point C represents the ARO solution under the minimax regret criterion, having the smallest worst-case regret of $90.1 MM, and the lowest worst-case NPV of $464.3 MM. Comparing Point A with Point B, we can see that the worst-case regret is significantly reduced by 52.3%, while the worst-case NPV decreases by merely 3.2%. This suggests that Point B could be a wise-choice solution based on the trade-off between the worst-case NPV and regret. The details on cost breakdowns, including fixed investment cost, variable investment cost, operating cost, and purchase cost, are shown as pie charts in Figure 7. From the pie charts, we can observe that almost half of the total cost comes from purchasing feedstocks for Points A, B and C. Moreover, the percentage of the investment costs of Point C is higher than those of the other two points due to larger process capacities and more frequent capacity expansions. The optimal design and planning decisions at time period 4 determined by the solutions of Points A, B and C are shown in Figure 8, Figure 9 and Figure 10, respectively. In these three figures, the optimal total capacities are displayed under operating processes.
To illustrate the optimal capacity expansion activities, we present the capacity profiles during the entire planning horizon for Points A, B and C in Figure 11, Figure 12 and Figure 13, respectively. From Figure 11, we can see that a total of 14 processes are selected in the optimal process network for Point A. By comparing Figure 11 and Figure 13, we can conclude that the 27
optimal expansion frequency of Point C is higher than that of Point A. Additionally, the optimal total capacities of some processes for Points B and C are different. For example, the optimal total capacity of Process 32 is 245 kt/y determined by Point B, while the corresponding capacity is 55 kt/y larger determined by Point C.
To demonstrate the efficiency of the computational algorithms in Sections 3.2.1 and 3.2.2, we select two instances and display the upper and lower bounds in each iteration of the algorithms in Figure 14 and Figure 15, respectively. In these two figures, the green dots stand for the upper bounds, and the yellow circles represent the lower bounds. The X-axis represents the iteration number, and the Y-axis represents the objective function values. The first instance corresponds to optimization problem (BCARO-λ) with λ=0. It can be seen from Figure 14 that the relative optimality gap decreases significantly from 89.46% to 3.90% during the first three iterations. The proposed computational algorithm takes only 7 iterations to reduce the relative optimality gap to 0.70%, and requires 9 iterations to reach a relative optimality gap below the predefined tolerance of 0.01%. The result demonstrates the efficiency of the proposed computational algorithm in Section 3.2.1.
The second instance corresponds to Point B in Figure 7. The performance of the computational algorithm is shown in Figure 15. The lower bound is not updated in the first iteration because the first-stage decision does not satisfy the ɛ-constraint. The lower bound is updated after the feasibility cuts are added to the master problem in the second iteration. The proposed computational algorithm takes only 4 iterations to converge within a reasonable solution time. The result demonstrates the efficiency of the proposed computational algorithm in Section 3.2.2.
iteration of the tailored C&CG algorithm in the instance corresponding to Point B.
To further investigate how the values of uncertainty budgets influence the trade-off between the conventional robustness and minimax regret criteria, we present the Pareto-optimal curves under different uncertainty budgets in Figure 16. For all these Pareto-optimal curves, we can see 28
a similar trend that worst-case NPV decreases as the worst-case regret reduces. For a fixed worst-case regret, increasing either the supply or demand uncertainty budget leads to a lower worst-case NPV. By comparing the red, purple and green curves, we can see that the Paretooptimal curve lies close to the upper left corner as we reduce the uncertainty budget for demand. From the blue, purple and orange curves, we can conclude that the Pareto-optimal curve is less sensitive to the change in the supply uncertainty budget. Therefore, the impact of the demand uncertainty budget is more significant than that of the supply uncertainty budget.
6. Application to Batch Process Scheduling under Demand Uncertainty This section presents another application of the proposed approach to multipurpose batch process scheduling under demand uncertainty. Multipurpose batch process scheduling is an operational problem solved on a daily or weekly basis [54]. In the scheduling problem, the following information is given: the structure of the facility, recipe data, equipment units and their sizes, storage capacities of intermediates, product prices and processing time [55, 56]. The objective is to generate the highest profit by determining the optimal batch sizes, the assignment of tasks to equipment units, the sequencing and timing of tasks. The solution of a batch scheduling problem is often represented by Gantt charts. In this application, product demands are subject to uncertainty.
6.1. Bi-criterion adaptive robust scheduling of batch process model The multi-criterion ARO model for batch process scheduling is formulated as follows. Following the conventional ARO model of batch process scheduling [21], the task assignment decision is made “here-and-now”, while other decisions are made in a “wait-and-see” manner after knowing the demand uncertainty realization. The two objectives are to maximize the worstcase profit and to minimize the worst-case regret given in (42)-(43). The definition of worst-case profit and worst-case regret are presented in (44)-(45). The constraints can be classified into assignment constraints (46)-(51), time constraints (52)-(62), batch size constraints (63)-(79), mass balance and storage constraints (80)-(84), and non-negativity and integrity constraints (85) -(86). A list of indices/sets, parameters and variables is given in the Nomenclature section. The objective functions are defined by (42)-(43), as shown below max Profit w Wsin ,Wfin
Wsin ,Wfin
29
(42)
min R w Wsin ,Wfin
(43)
Wsin ,Wfin
The worst-case profit and regret are presented in (44) and (45), respectively.
Profit w min
max
dU Bfin , Bpin , Bsin , BIisn , BOisn , STin , SAs ,Tn , Din ,Tfin ,Tsin
R w max dU
ps STs|n| vci Bsin fci Wsin i n i n sP
d ps STs|n| vci Bsin fci Wsin Bfin , Bpin , Bsin , BIisn , i n i n sP BO , ST , SA ,T , min
isn
in
s
(44)
(45)
n
Din ,Tfin ,Tsin
Constraints (46)-(47) specify that at most one task can start or finish at each equipment unit at time point n. Constraint (48) specifies that a task must finish in the future once it starts. Constraint (49) enforces that at most one task can be processed at each unit at any time point n.
Ws
in
1
j, n
(46)
Wf
in
1
j, n
(47)
iI j
iI j
Ws
in
n
i
(48)
n
Ws iI j n ' n
Wfin ,
in '
Wfin ' 1,
j , n
(49)
Constraint (50) enforces that no tasks could finish at the first time point. Similarly, constraint (51) specifies that no task could begin at the last time point.
Wfi1 0
i
(50)
Wsi|n| 0
i
(51)
Constraints (52)-(61) are for the timing of tasks. Constraint (52) describes the processing time of each task. The relationship between starting time, finish time and processing time of a task is given in constraints (53)-(56). Constraint (57) specifies that task i starting at time point n1 should finish at or before time point n if Wfin equals to 1. Constraint (58) describes the relationship between starting time and the time corresponding to time point n. fdi Wsin +vdi Bsin Din Tsin Din ht 1 Wsin Tfin ,
30
i, n i, n
(52) (53)
Tfin Tsin Din ht 1 Wsin , Tfin Tfi n1 Wsin ht ,
i, n
Din Tfin Tfi n1 ,
i, n
Tfi n1 Tn ht 1 Wfin ,
Tsin Tn ,
i, n
(54) (55) (56)
i, n 2
i, n
(57) (58)
Constraints (59)-(62) defines the initial condition and upper bound of related variables.
T1 0 Tn T n1 ,
(59) n 2
(60)
T|n| ht
Tfin ht ,
(61)
i, n
(62)
Constraints (63)-(69) specify the lower and upper bounds of batch sizes, and linking variables associated with batches.
Bsin biminWsin ,
i, n
(63)
Bsin bimaxWsin ,
i, n
(64)
Bfin biminWfin ,
i, n
(65)
Bfin bimaxWfin ,
i, n
(66)
Bpin bimin Wsin ' Wf in ' , n ' n n ' n
i, n
(67)
Bpin bimax Wsin ' Wf in ' , n ' n n ' n
i, n
(68)
i, n 2
(69)
Bsi n1 Bpi n1 Bpin Bfin ,
Constraints (70)-(73) describe the input and output of batch operations, as well as their upper bounds.
sisi BIisn sisi mis Bsin 0, i, s, n
(70)
sosi BOisn sosi mis Bfin 0, i, s, n
(71)
sisi BIisn sisi bimax mis Wsin , i, s, n
(72)
31
sosi BOisn sosi bimax mis Wfin , i, s, n
(73)
Constraints (74)-(79) are the initial and final conditions for the variables associated with batch sizes.
Bpin 0, i, n 1
(74)
Bfin 0, i, n 1
(75)
BOisn 0, i, s, n 1
(76)
Bsin 0, i, n n
(77)
Bpin 0, i, n n
(78)
BIisn 0, i, s, n n
(79)
Constraint (80) describes the initial inventory of different materials.
STsn sisi BIisn sints , s, n 1
(80)
i
Constraint (81) is for the mass balance. STsn STs n1
B
iTOs
O isn
I Bisn ,
s, n
(81)
iTI s
Constraint (82) specifies the storage capacity. If a chemical is sold out as a final product, constraint (83) states that its sale amount is equal to its inventory level at the last time point. Constraint (84) enforces that sale amount of each final product should meet its demand. Constraints (85) and (86) defines nonnegative continuous variables and binary variables.
STin cs , s, n
(82)
snps STin snps SAs 0, s, n n
(83)
snps SAs snps ds , s
(84)
Bfin , Bpin , Bsin , BIisn , BOisn , STin , SAs , Tn , Din , Tfin , Tsin 0
(85)
Wsin ,Wfin 0,1
(86)
Uncertainty set for product demands, denoted as U, can be directly defined as a polytope or constructed from uncertainty data. The resulting optimization problem on batch process scheduling is a bi-objective multi-level mixed-integer programming problem. It explicitly
32
accounts for the conventional robustness and minimax regret criteria to benefit batch process operations.
6.2. Case study 1 In this subsection, we consider a specific case study on the short-term scheduling of multipurpose batch process [57]. The state-task network of this batch process is shown in Figure 17. As can be seen from Figure 17, this multipurpose batch process involves one heating task, three reaction tasks, and one separation task. This network involves three raw materials (A-C), four intermediates (Hot A, Impure E, and I1-I2) and two products (P1-P2). There are four equipment units: one heater, two reactors, and one separator. The mass balance coefficients are also listed in Figure 17. The scheduling horizon is 8 hours, and there are 6 time points.
The demands for Product 1 and Product 2 are subject to uncertainty. In this case study, the uncertainty set of product demands is directly given in Figure 18. This uncertainty set is a polytope with four extreme points (Points A-D). The coordinates of these extreme points are given as follows: Point A (40, 30), Point B (52, 28), Point C (24, 42) and Point D (12, 44).
By using the computational algorithm in Section 3.2.1, we solve the corresponding problem by changing λ from 0 to 1 with an increment of 0.05. The Gantt charts associated with λ=0 and λ=1 are demonstrated in Figure 19 (a) and (b), respectively. By comparing the two Gantt charts, we can observe that the assignment decisions determined by the two approaches are different. Specifically, Reaction 3 starts at a unique time point in the scheduling solution determined by the conventional robustness criterion, while Reaction 3 shares the same starting time point with Reaction 1 in the scheduling solution determined by the minimax regret criterion. The values of worst-case profit and regret are listed in Figure 19. The scheduling solution determined by the minimax regret criterion has the lowest worst-case regret of $13.6, and only decreases the worstcase profit by 0.3% compared with solution under the conventional robustness criterion. Due to the discrete nature of first-stage decisions and problem setup, there is no trade-off solution, when λ takes other values between 0 and 1.
33
6.3. Case study 2 In this subsection, another case study is considered to demonstrate that the solutions with the conventional robustness criterion and the minimax regret criterion could have the same worstcase profit and worst-case regret. This case study is originated from an industrial multipurpose batch process in The Dow Chemical Company [58-61]. Figure 20 displays the state-task network of the batch process in this case study. This multipurpose batch process has one preparation task, three reaction tasks, two packing tasks and two drumming tasks. The chemicals include four raw materials (M1-M4), six intermediates (I1-I6) and four products (P1-P4). The equipment units include one mixer, two reactors, one finishing system and one drumming line. The scheduling horizon is one week, i.e. 168 hours, and there are 11 time points.
Due to market fluctuations, product demands for P1-P4 are subject to uncertainties in this case study. A total of 200 uncertainty data are used for constructing the uncertainty set, and each data point has 4 dimensions for a combination of all final products. We employ the computational algorithm in Section 3.2.1 to solve the resulting scheduling problems with different values of parameter λ. The Gantt charts of the ARO solutions under the conventional robustness criterion and the minimax regret criterion are shown in Figure 21 (a) and (b), respectively. From the Gantt chart, we can observe that the assignment decisions are different, but the worst-case profit and regret remain unchanged for the two solutions. This is because of a different problem setup from Case study 1.
7. Conclusions In this paper, a novel ARO framework was proposed to account for the minimax regret criterion. Apart from the conventional robustness criterion, the worst-case regret was introduced as another objective to minimize the deviation from the perfect information solution. A set of Paretooptimal solutions was generated to reveal the systematic trade-offs between the conventional robustness and minimax regret criteria. The proposed optimization framework was general enough to apply to various kinds of uncertainty sets, such as the polyhedral and data-driven uncertainty sets. A data-driven version of this framework, named DDBCARO, leveraged the advantages of the minimax regret criterion, and at the same time exploited the value of big data via the Dirichlet process mixture model. Since the resulting optimization problem has multi-level 34
optimization structures, we further proposed the tailored C&CG algorithms to address the computational challenge. We presented a numerical example and applications on strategic planning of process networks and batch process scheduling. The results showed that the proposed framework could make a systematic trade-off between the conventional robustness criterion and minimax regret criterion.
Acknowledgments The authors acknowledge financial support from the National Science Foundation (NSF) CAREER Award (CBET-1643244).
Appendix A. Details of the Data-Driven Uncertainty Set In data-driven robust optimization paradigm, one paramount ingredient is the uncertainty set, which is learned directly from uncertainty data. To model uncertainty data with high fidelity, we adopt the Dirichlet process mixture model in this work. The Dirichlet process mixture model is a powerful nonparametric Bayesian statistical model, which has an infinite-dimensional parameter space. The resulting data-driven uncertainty set could adjust its complexity based on the data structure and complexity, thereby truthfully capturing the nature of uncertainty [39]. The Dirichlet process is defined as follows [62]. A random distribution F is said to be distributed as a Dirichlet process, denoted as F ~ DP (α, F0), if for any fixed number of partitions (Θ1, ..., Θr) of Θ, (F(Θ1),…, F(Θr)) obeys a Dirichlet distribution, i.e.,
F (1 ),
, F (r ) ~ Dir F0 (1 ), , F0 (r )
(A1)
where F0 is a base distribution on Θ, α is the concentration parameter that indicates the users’ belief in F0. Dir is the notation for the Dirichlet distribution. Besides, a random draw from the Dirichlet process has the following representation based on a stick-breaking construction [63].
F k 1 kk
(A2)
where the weight k k j 1 1 j , and βk is the proportion of stick being broken from the k 1
remaining stick. βk follows a Beta distribution, denoted as βk ~ Beta (1, α). The random variable
35
drawn from the base distribution F0 is denoted as θk. k is a delta function which is equal to 1 at θk, and 0 otherwise. The Dirichlet process is a distribution on distributions. A random draw from a Dirichlet process, i.e. DP (α, F0), is a distribution F. We can still have a random draw from the distribution F. The Dirichlet process mixture model adds one more level to this hierarchy. It utilizes θk as the parameters of the data distribution. As shown in eq. (A2), the Dirichlet process mixture model can be interpreted as a mixture model that has a potentially infinite number of mixture components. The general Dirichlet process mixture model is summarized as follows [64].
k ~ Beta(1, ) k F0 ~ F0 li 1 , 2 ,
~ Mult ( )
(A3)
oi li ~ p (oi li )
where Mult denotes a multinomial distribution, and li is the index of mixture components to which the observation oi is assigned. Variational inference algorithm is employed for the Dirichlet process mixture model [64, 65]. In Bayesian statistics, the model parameters are considered as random variables, and the inference is to obtain the posterior distribution of these latent random variables. The mixtures of Gaussian distributions are adopted in this work. θk = (ηk, Hk) includes mean vector and precision matrix. It is assumed that θk ~ NW (μk, λk, ωk, ψk), where NW represents the normal Wishart distribution. In the variational inference algorithm, the posterior distribution is approximated by a factorized variational distribution q (l, β, α, η, H), as given below q l, β, , η, H
N
M 1
M
i 1
k 1
k 1
q li q k q q ηk , Hk
(A4)
where M is the truncation level. Suppose Θ= (l, β, α, η, H) is set of all latent variables and parameters of the Dirichlet process mixture model, and U= (u1, …, uN) are the observations of uncertain parameters. The Dirichlet process mixture model, specifies the distribution p(U, Θ). From a Bayesian perspective, the goal of inference is to find the posterior distribution of Θ given the observation data U, that is p(Θ|U). In the variational inference, the log marginal probability is decomposed into two parts. log p U ELBO q KL q p
36
(A5)
where the first part is called the evidence lower bound (ELBO), and the second part is KullbackLeibler divergence or KL divergence; q is the variational distribution to approximate the posterior distribution p(Θ|U). They are defined as follows [65, 66]. p U, Θ ELBO q q Θ log dΘ q Θ
(A6)
p Θ U KL q p q Θ log dΘ q Θ
(A7)
Given the KL(q||p) is a distance between the distributions q and p, the goal of the variational algorithm is to minimize KL(q||p) to achieve the best approximation. Since log p(U)is a constant with respect to q, minimizing the KL divergence levels to maximizing the ELBO. To obtain a tractable inference, the variational inference utilize a factorized variational distribution q(l, β, α, η, H). q Gamma 1* , 2*
(A8)
where Gamma represents the Gamma distribution, and the hyper-parameters are updated by
1* 1 M 1
(A9)
2* 2 i 1 ln 1 i M 1
(A10)
The notation <.> represents the expectation of corresponding random variable. All closed forms of <.> are listed at the end of Appendix A. The posterior of βk is as follows. q k Beta k , vk
(A11)
where Beta represents the Beta distribution, and the hyper-parameters are updated by k* 1 n1 q ln k
(A12)
vk* n1 i k 1 q ln i
(A13)
N
N
M
The posterior of q(ηk, Hk)is the Normal-Wishart distribution, and the related parameters are calculated in (A15)-(A19).
q ηk , Hk N ηk μ*k , k*Ψ*k1 W H k Ψ*k1 , k* 1
Nk n1 q ln k N
37
(A14) (A15)
k* k Nk
(A16)
k* k Nk
(A17)
μ*k 1 Ψ Ψ k n 1 q ln k u n Nk N
* k
+
k N k 1 k N k N k
N n '1
μ 1
k
* k
k
N n 1
q ln k u n
1 n '1 q ln ' k un ' un N k N
1 q ln ' k u n ' μ k N k
q ln k
N n '1
(A18) n '1 q ln ' k un '
T
N
q ln ' k u n ' μ k
nk*
(A19)
T
(A20)
* i 1 ni M
where ρnk* is defined as follows. k 1 1 K 1 T nk* exp ln k i 1 ln 1 i ln H k ln 2 un ηk H k un ηk 2 2 2
(A21)
where K= dim (u) is the dimension of uncertainty data. The expressions involved < > are listed below.
un ηk
T
1* * 2
(A22)
ln k k k vk
(A23)
ln 1 k vk k vk
(A24)
* 1 i K * ln H k i 1 k K ln 2 ln Ψ k 2
(A25)
H k u n ηk
K
k
k u n μ*k Ψ*k1 u n μ*k T
where ψ is the digamma function. The variational inference algorithm is summarized below.
38
(A26)
References [1] [2]
[3] [4] [5]
[6]
[7]
[8]
[9] [10] [11] [12] [13] [14] [15]
[16]
[17] [18]
N. V. Sahinidis, "Optimization under uncertainty: State-of-the-art and opportunities," Computers & Chemical Engineering, vol. 28, pp. 971-983, 2004. I. E. Grossmann, R. M. Apap, B. A. Calfa, P. García-Herreros, and Q. Zhang, "Recent advances in mathematical programming techniques for the optimization of process systems under uncertainty," Computers & Chemical Engineering, vol. 91, pp. 3-14, 2016. E. N. Pistikopoulos, "Uncertainty in process design and operations," Computers & Chemical Engineering, vol. 19, pp. 553-563, 1995. J. M. Mulvey, R. J. Vanderbei, and S. A. Zenios, "Robust Optimization of Large-Scale Systems," Operations Research, vol. 43, pp. 264-281, 1995. Y. A. Guzman, L. R. Matthews, and C. A. Floudas, "New a priori and a posteriori probabilistic bounds for robust counterpart optimization: I. Unknown probability distributions," Computers & Chemical Engineering, vol. 84, pp. 568-598, 2016. K. Tong, F. You, and G. Rong, "Robust design and operations of hydrocarbon biofuel supply chain integrating with existing petroleum refineries considering unit cost objective," Computers & Chemical Engineering, vol. 68, pp. 128-139, 2014. A. C. S. Amaro and A. P. F. D. Barbosa-Póvoa, "The effect of uncertainty on the optimal closed-loop supply chain planning under different partnerships structure," Computers & Chemical Engineering, vol. 33, pp. 2144-2158, 2009. Y. Yuan, Z. Li, and B. Huang, "Robust optimization under correlated uncertainty: Formulations and computational study," Computers & Chemical Engineering, vol. 85, pp. 58-71, 2016. J. R. Birge and F. Louveaux, Introduction to stochastic programming: Springer Science & Business Media, 2011. A. Ben-Tal and A. Nemirovski, "Robust optimization – methodology and applications," Mathematical Programming, vol. 92, pp. 453-480, 2002. D. Bertsimas, D. B. Brown, and C. Caramanis, "Theory and applications of robust optimization," SIAM Review, vol. 53, pp. 464-501, 2011. Z. Li and Z. Li, "Optimal robust optimization approximation for chance constrained optimization problem," Computers & Chemical Engineering, vol. 74, pp. 89-99, 2015. K. Shang, Z. Feng, L. Ke, and F. T. S. Chan, "Comprehensive Pareto Efficiency in robust counterpart optimization," Computers & Chemical Engineering, vol. 94, pp. 75-91, 2016. N. H. Lappas and C. E. Gounaris, "Multi-stage adjustable robust optimization for process scheduling under uncertainty," AIChE Journal, vol. 62, pp. 1646-1667, 2016. D. E. Majewski, M. Wirtz, M. Lampe, and A. Bardow, "Robust multi-objective optimization for sustainable design of distributed energy supply systems," Computers & Chemical Engineering, vol. 102, pp. 26-39, 2017. C. Ning and F. You, "A data-driven multistage adaptive robust optimization framework for planning and scheduling under uncertainty," AIChE Journal, vol. 63, pp. 4343-4369, 2017. C. Shang, X. Huang, and F. You, "Data-Driven Robust Optimization Based on Kernel Learning," Computers & Chemical Engineering, 2017. A. Ben-Tal, A. Goryashko, E. Guslitzer, and A. Nemirovski, "Adjustable robust solutions of uncertain linear programs," Mathematical Programming, vol. 99, pp. 351-376, 2004.
39
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26] [27]
[28] [29] [30] [31]
[32]
[33] [34]
B. Zeng and L. Zhao, "Solving two-stage robust optimization problems using a columnand-constraint generation method," Operations Research Letters, vol. 41, pp. 457-461, 2013. Q. Zhang, M. F. Morari, I. E. Grossmann, A. Sundaramoorthy, and J. M. Pinto, "An adjustable robust optimization approach to scheduling of continuous industrial processes providing interruptible load," Computers & Chemical Engineering, vol. 86, pp. 106-119, 2016. H. Shi and F. You, "A computational framework and solution algorithms for two-stage adaptive robust scheduling of batch manufacturing processes under uncertainty," AIChE Journal, vol. 62, pp. 687-703, 2016. J. Gong and F. You, "Optimal processing network design under uncertainty for producing fuels and value-added bioproducts from microalgae: Two-stage adaptive robust mixed integer fractional programming model and computationally efficient solution algorithm," AIChE Journal, vol. 63, pp. 582-600, 2017. J. Gong, D. J. Garcia, and F. You, "Unraveling Optimal Biomass Processing Routes from Bioconversion Product and Process Networks under Uncertainty: An Adaptive Robust Optimization Approach," ACS Sustainable Chemistry & Engineering, vol. 4, pp. 31603173, 2016. D. E. Majewski, M. Lampe, P. Voll, and A. Bardow, "TRusT: A Two-stage Robustness Trade-off approach for the design of decentralized energy supply systems," Energy, vol. 118, pp. 590-599, 2017. D. Yue and F. You, "Optimal supply chain design and operations under multi-scale uncertainties: Nested stochastic robust optimization modeling framework and solution algorithm," AIChE Journal, vol. 62, pp. 3041-3055, 2016. D. E. Bell, "Regret in Decision Making under Uncertainty," Operations Research, vol. 30, pp. 961-981, 1982. T. Assavapokee, M. J. Realff, and J. C. Ammons, "Min-Max Regret Robust Optimization Approach on Interval Data Uncertainty," Journal of Optimization Theory and Applications, vol. 137, pp. 297-316, 2008. P. Kouvelis and G. Yu, Robust discrete optimization and its applications vol. 14: Springer Science & Business Media, 2013. I. R. Dunning, "Advances in robust and adaptive optimization: algorithms, software, and insights," Massachusetts Institute of Technology, 2016. Z.-L. Chen, N. G. Hall, and H. Kellerer, "Dynamic Pricing to Minimize Maximum Regret," Production and Operations Management, vol. 26, pp. 47-63, 2017. R. Yokoyama, K. Fujiwara, M. Ohkura, and T. Wakui, "A revised method for robust optimal design of energy supply systems based on minimax regret criterion," Energy Conversion and Management, vol. 84, pp. 196-208, 2014. B. Chen, J. Wang, L. Wang, Y. He, and Z. Wang, "Robust Optimization for Transmission Expansion Planning: Minimax Cost vs. Minimax Regret," IEEE Transactions on Power Systems, vol. 29, pp. 3069-3077, 2014. R. Jiang, J. Wang, M. Zhang, and Y. Guan, "Two-Stage Minimax Regret Robust Unit Commitment," IEEE Transactions on Power Systems, vol. 28, pp. 2271-2282, 2013. R. Loulou and A. Kanudia, "Minimax regret strategies for greenhouse gas abatement: methodology and application," Operations Research Letters, vol. 25, pp. 219-230, 1999.
40
[35] [36]
[37] [38] [39]
[40]
[41] [42] [43] [44] [45]
[46] [47]
[48]
[49]
[50]
[51] [52] [53]
J. Ide and A. Schöbel, "Robustness for uncertain multi-objective optimization: a survey and analysis of different concepts," OR Spectrum, vol. 38, pp. 235-271, 2016. L. Wang, Q. Li, R. Ding, M. Sun, and G. Wang, "Integrated scheduling of energy supply and demand in microgrids under uncertainty: A robust multi-objective optimization approach," Energy, vol. 130, pp. 1-14, 2017. J. Fliege and R. Werner, "Robust multiobjective optimization & applications in portfolio optimization," European Journal of Operational Research, vol. 234, pp. 422-433, 2014. D. Kuroiwa and G. M. Lee, "On robust multiobjective optimization," Vietnam J. Math, vol. 40, pp. 305-317, 2012. C. Ning and F. You, "Data-Driven Adaptive Nested Robust Optimization: General Modeling Framework and Efficient Computational Algorithm for Decision Making Under Uncertainty," AIChE Journal, vol. 63, pp. 3790-3817, 2017. H. Aissi, C. Bazgan, and D. Vanderpooten, "Min–max and min–max regret versions of combinatorial optimization problems: A survey," European Journal of Operational Research, vol. 197, pp. 427-438, 2009. T. Campbell and J. P. How, "Bayesian nonparametric set construction for robust optimization," in 2015 American Control Conference (ACC), 2015, pp. 4216-4221. G. P. Rangaiah and A. Bonilla-Petriciolet, Multi-objective optimization in chemical engineering: developments and applications: John Wiley & Sons, 2013. F. You and I. E. Grossmann, "Design of responsive supply chains under demand uncertainty," Computers & Chemical Engineering, vol. 32, pp. 3090-3111, 2008. A. Thiele, T. Terry, and M. Epelman, "Robust linear optimization with recourse," Rapport technique, pp. 4-37, 2009. D. Bertsimas, E. Litvinov, X. A. Sun, J. Zhao, and T. Zheng, "Adaptive Robust Optimization for the Security Constrained Unit Commitment Problem," IEEE Transactions on Power Systems, vol. 28, pp. 52-63, 2013. F. Glover, "Improved linear integer programming formulations of nonlinear integer problems," Management Science, vol. 22, pp. 455-460, 1975. N. V. Sahinidis, I. E. Grossmann, R. E. Fornari, and M. Chathrathi, "Optimization model for long range planning in the chemical industry," Computers & Chemical Engineering, vol. 13, pp. 1049-1063, 1989. F. You and I. E. Grossmann, "Stochastic Inventory Management for Tactical Process Planning Under Uncertainties: MINLP Models and Algorithms," AIChE Journal, vol. 57, pp. 1250-1277, 2011. D. Yue and F. You, "Planning and Scheduling of Flexible Process Networks Under Uncertainty with Stochastic Inventory: MINLP Models and Algorithm," AIChE Journal, vol. 59, pp. 1511-1532, 2013. A. A. Levis and L. G. Papageorgiou, "A hierarchical solution approach for multi-site capacity planning under uncertainty in the pharmaceutical industry," Computers & Chemical Engineering, vol. 28, pp. 707-725, 2004. S. Ahmed and N. V. Sahinidis, "Robust process planning under uncertainty," Industrial & Engineering Chemistry Research, vol. 37, pp. 1883-1892, 1998. M. L. Liu and N. V. Sahinidis, "Optimization in Process Planning under Uncertainty," Industrial & Engineering Chemistry Research, vol. 35, pp. 4154-4165, 1996. E. Rosenthal, GAMS-A user’s guide. Washington, DC: GAMS Development Corporation, 2008. 41
[54]
[55]
[56]
[57]
[58]
[59]
[60]
[61]
[62] [63] [64] [65] [66]
S. Moniz, A. P. Barbosa-Póvoa, and J. P. de Sousa, "Simultaneous regular and nonregular production scheduling of multipurpose batch plants: A real chemical– pharmaceutical case study," Computers & Chemical Engineering, vol. 67, pp. 83-102, 2014. Y. Chu and F. You, "Model-based integration of control and operations: Overview, challenges, advances, and opportunities," Computers & Chemical Engineering, vol. 83, pp. 2-20, 2015. E. Muñoz, E. Capón-García, M. Moreno-Benito, A. Espuña, and L. Puigjaner, "Scheduling and control decision-making under an integrated information environment," Computers & Chemical Engineering, vol. 35, pp. 774-786, 2011. E. Kondili, C. C. Pantelides, and R. W. H. Sargent, "A general algorithm for short-term scheduling of batch operations—I. MILP formulation," Computers & Chemical Engineering, vol. 17, pp. 211-227, 1993. Y. Chu, J. M. Wassick, and F. You, "Efficient scheduling method of complex batch processes with general network structure via agent-based modeling," AIChE Journal, vol. 59, pp. 2884-2906, 2013. J. M. Wassick, A. Agarwal, N. Akiya, J. Ferrio, S. Bury, and F. You, "Addressing the operational challenges in the development, manufacture, and supply of advanced materials and performance products," Computers & Chemical Engineering, vol. 47, pp. 157-169, 2012. Y. Chu, F. You, J. M. Wassick, and A. Agarwal, "Integrated planning and scheduling under production uncertainties: Bi-level model formulation and hybrid solution method," Computers & Chemical Engineering, vol. 72, pp. 255-272, 2015. Y. Chu, F. You, and J. M. Wassick, "Hybrid method integrating agent-based modeling and heuristic tree search for scheduling of complex batch processes," Computers & Chemical Engineering, vol. 60, pp. 277-296, 2014. W. Y. Teh, "Dirichlet Process," in Encyclopedia of Machine Learning, C. Sammut and G. I. Webb, Eds., ed Boston, MA: Springer US, 2010, pp. 280-287. J. Sethuraman, "A constructive definition of Dirichlet priors," Statistica Sinica, vol. 4, pp. 639-650, 1994. D. M. Blei and M. I. Jordan, "Variational inference for Dirichlet process mixtures," pp. 121-143, 2006. K. Kurihara, M. Welling, and N. A. Vlassis, "Accelerated variational Dirichlet process mixtures," in Advances in neural information processing systems, 2006, pp. 761-768. X. Wei and C. Li, "The infinite Student's t-mixture for robust modeling," Signal Processing, vol. 92, pp. 224-234, 2012.
42
Algorithm. Tailored C&CG algorithm for (DDBCARO-λ) 1: 2:
Set LB , UB , k 0 , and UB LB do while UB
;
Solve (MP-λ) to obtain x*k 1 ,1k 1* ,2k 1* ;
3:
Update LB 1 1k 1* 2k 1* ;
4: 5: 6:
for i 1 to m do Solve (SUP1,i) to obtain u1,k i 1 and Q1i x*k 1 ;
Solve (SUP2,i) to obtain uk2,i1 , xˆ ik 1 , yˆ ik 1 and Q2i x*k 1 ;
7: 8:
end
i* arg max Q1i x*k 1 ;
9:
i
u1k 1 u1,k i*1 and Q1 x*k 1 Q1i x*k 1 ; *
10: 11:
j * arg max Q2j x*k 1 ; j
, xˆ k 1 xˆ kj*1 , yˆ k 1 yˆ kj*1 and Q2 x*k 1 Q2j x*k 1 ;
12:
u
13:
Update UB min UB, cT x*k 1 1 Q1 x*k 1 Q2 x*k 1 ;
k 1 2
u
k 1 2, j*
*
k 1 1
Create second-stage variables y 14:
1 c x b y , Tx Wy T
T
k 1 1
2 c x b y T
T
k 1 2
c xˆ T
k 1 1
k 1
,y
k 1 2
, and add cuts
k 1 1
h Mu
b yˆ k 1 , Tx Wy k2 1 h Mu k2 1 T
to (MP-λ);
k k 1 ;
15: 16:
end
17:
return UB ;
Figure 1. The pseudocode of the proposed computational algorithm for (DDBCARO-λ).
43
Algorithm. Tailored C&CG algorithm for (DDBCARO-ε) 1: 2:
Set LB , UB , k 0 , and UB LB do while UB
;
3:
Solve (MP-ε) to obtain x*k 1 , k 1* ;
4: 5: 6:
Update LB ; for i 1 to m do Solve (SUP1,i) to obtain u1,k i 1 and Q1i x*k 1 ; k 1*
Solve (SUP2,i) to obtain uk2,i1 , xˆ ik 1 , yˆ ik 1 and Q2i x*k 1 ;
7: 8:
end
i* arg max Q1i x*k 1 ;
9:
i
u1k 1 u1,k i*1 and Q1 x*k 1 Q1i x*k 1 ; *
10: 11:
j * arg max Q2j x*k 1 ; j
12:
u
k 1 2
u
, xˆ k 1 xˆ kj*1 , yˆ k 1 yˆ kj*1 and Q2 x*k 1 Q2j x*k 1 ; *
if cT x*k 1 Q2 x*k 1 then
13:
Update UB min UB, cT x*k 1 Q1 x*k 1 ;
14: 15: 16:
k 1 2, j*
end Create second-stage variables y1k 1 , y k2 1 , and add cuts
cT x bT y1k 1 , Tx Wy1k 1 h Mu1k 1
cT x bT y k2 1 cT xˆ k 1 bT yˆ k 1 , Tx Wy k2 1 h Mu k2 1
to (MP-ε);
k k 1 ;
17: 18:
end
19:
return UB ;
Figure 2. The pseudocode of the proposed computational algorithm for (DDBCARO-ε).
44
Figure 3. The scatter plot of uncertainty data and the data-driven uncertainty set.
Figure 4. The Pareto-optimal curve in the numerical example.
45
Figure 5. The chemical process network in the case study.
46
Figure 6. The worst-case NPV and worst-case regret under different values of parameter λ.
Figure 7. Pareto-optimal curve of the case study with breakdown of the total cost.
47
Figure 8. The optimal design and planning decisions for time period 4 determined by the solution of Point A. The optimal production capacity of each process is displayed under each operating process. 48
Figure 9. The optimal design and planning decisions for time period 4 determined by the solution of Point B. The optimal production capacity of each process is displayed under each operating process. 49
Figure 10. The optimal design and planning decisions for time period 4 determined by the solution of Point C. The optimal production capacity of each process is displayed under each operating process. 50
Figure 11. Optimal capacity expansion decisions over the entire planning horizon for Point A in the case study.
Figure 12. Optimal capacity expansion decisions over the entire planning horizon for Point B in the case study.
51
Figure 13. Optimal capacity expansion decisions over the entire planning horizon for Point C in the case study.
Figure 14. Upper and lower bounds in each iteration of the tailored C&CG algorithm for optimization problem (BCARO-λ) with λ=0.
52
Figure 15. Upper and lower bounds in each
Figure 16. Pareto-optimal curves under different uncertainty budgets in the case study.
53
Figure 17. State-task network of the multipurpose batch process in case study 1.
Figure 18. Uncertainty set of product demands in case study 1.
54
Figure 19. Gantt charts determined by (a) the conventional robustness criterion, and (b) the minimax regret criterion in case study 1.
55
Figure 20. The state-task network of the multipurpose batch process in case study 2.
56
Figure 21. Gantt charts determined by (a) the conventional robustness criterion, and (b) the minimax regret criterion in case study 2.
57
Algorithm. Variational inference algorithm N Input uncertainty data set U un n1 ; 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13:
Set ζe and q(l, β, α, η, H); ELBOt q ELBOt 1 q e do while ELBOt 1 q Update q by the expression in (A8); for k=1 to M do Update q k by the expression in (A11); Update q ηk , H k by the expression in (A14); end for n=1 to N do Update q ln by the expression in (A20); end end N
M 1
M
i 1
k 1
k 1
Output q l, β, , η, H q li q k q q ηk , H k
Figure A1. The variational inference algorithm for the Dirichlet process mixture model.
58
Table 1. Results of the BCARO-λ approach in the numerical example. BCARO-λ (λ=0) *
BCARO-λ (λ=0.4) *
BCARO-λ (λ=1) *
Worst-case cost Worst-case regret Iterations Total CPU (s)
1,555.7 878.6 5 7
1,615.5 694.2 5 8
1,833.2 451.0 5 8
First-stage decision x
x2 56.3
x2 48.1
x2 48.9
x3 81.3
x3 70.2
x3 49.8
x1 62.4
x1 58.9
x1 58.1
* λ is the weight of the minimax regret criterion, and (1-λ) is the weight of the conventional robustness criterion in model (BCARO-λ).
Table 2. Comparisons of model and solution statistics for (BCARO-λ).
Bin. Var.
BCARO-λ (λ=0)* Master problem 190
Sub-problem
BCARO-λ (λ=0.2)* Master problem
Sub-problem
BCARO-λ (λ=1)* Master problem
Sub-problem
Cont. Var.
8,843
302 17,141
190
302
190
302
5,083
17,141
8,843
Constraints
11,569
48,749
17,141
6,681
48,749
11,569
Worst-case NPV ($MM)
544.0
505.6
464.3
48,749
Worst-case regret ($MM)
705.9
207.7
90.1
Total CPU (s)
180
214
409
Iterations
9
5
9
* λ is the weight of the conventional robustness criterion, and (1-λ) is the weight of the conventional robustness criterion in model (BCARO-λ).
59