Article pubs.acs.org/IECR
Modeling of Hydrogen Networks in a Refinery Using a Stochastic Programming Appraoch Anoop Jagannath and Ali Almansoori* Department of Chemical Engineering, The Petroleum Institute, P.O. Box 2533, Abu Dhabi, United Arab Emirates ABSTRACT: As hydrogen demand in the refineries is increasing, hydrogen management strategies are of great interest for the refiners. Hydrogen management deals with the optimal distribution, routing, and flow allocation of hydrogen gas within a refinery linking the hydrogen producing and consuming units. This is also known as the hydrogen network synthesis or design. Due to uncertainty in the crude feedstock or final product specification in a refinery, there exists variations or uncertainties in the operating conditions (parameters) required for hydrogen network design. Hence there is a need to incorporate these uncertainties while designing hydrogen networks. In this paper, we address the issue of optimal hydrogen network design under uncertainty or uncertain operating parameters. For this, we develop a nonconvex multiscenario Mixed Integer Nonlinear Programming (MINLP) model which is the deterministic equivalent of mathematical two stage stochastic model. The purpose is to develop a network design that is optimal for different operating conditions of a refinery constituted by multiple scenarios of operation. The superstructure for this hydrogen network design under uncertainty includes hydrogen reuse, recycle, and purification options. It also has hydrogen gas stream conditioning equipment, like compressors and valves, to cater to the nonisobaric conditions and heaters and coolers to deal with the nonisothermal conditions. The size of this multiscenario two stage mathematical model builds with increase in the number of scenarios. For this, standard optimization solvers may not be able to guarantee global optimality within reasonable computational time. We propose a specific solution strategy to solve this mathematical model to global optimality by convexifying the nonlinear terms. This strategy is found to exhibit exceptional computational time savings compared to the conventional branch and bound based global optimization approaches. Two case studies are solved to illustrate the applicability of the proposed model and solution strategy.
1. INTRODUCTION Hydrogen is an essential commodity for refinery operations. It is required for producing fuels with cleaner specifications, for breaking down of heavier fuels, and as feed in other hydrocarbon processing operations. It is either imported from the outside or produced in-house in a refinery. This hydrogen is primarily consumed by the processing units, namely, hydrocrackers and hydrotreaters, which have specific demand for hydrogen. The producers of hydrogen are called hydrogen sources, whereas the processing units act as hydrogen consumers. Complex chemical reactions take place in the processing units, and they release purge/off gases which are relatively rich in hydrogen. These purge/off gases may be purified or could be sent to the fuel gas system. The purifiers constitute an integral part of the refinery hydrogen network and help recover hydrogen by purifying the off/purge gases generated from the processing units. The interaction among all the above-mentioned units, namely, hydrogen sources, processing units, purification units, and fuel gas sinks, and developing a network capturing these interactions in an optimal manner constitutes the refinery hydrogen network synthesis problem. This could also be described as the hydrogen distribution problem in a refinery. This is one of the key aspects of hydrogen management strategies adopted in refineries to optimally distribute hydrogen throughout the refinery. Hydrogen network design belongs to the category of application of mass integration principles. Mass integration is a special type of process integration which deals with the analysis of the flow of a particular species within a process. The © XXXX American Chemical Society
process here refers to a holistic system or a group of processes which deal with this species of interest. Process integration aims to study the interactions between these different processes or unit operations dealing with the species of interest to minimize its consumption. For any targeted species in mass integration, two types of units are present. One is the producer (which produces or gives out the species) and other is the consumer (which accepts or consumes the species). The species of interest in the hydrogen network is hydrogen. Producers in a hydrogen network are hydrogen plant or steam methane reformer, catalytic reformer, hydrogen import through pipeline, etc. The consumers are hydrocrackers, hydrotreaters, and fuel gas system. Some interceptor units could be present in the hydrogen network like the purification unit (which upgrades the purity of hydrogen gas) and conditioning equipment like compressor, valve, heater, and cooler (which changes the condition of hydrogen gas like pressure, temperature). These act as both producers and consumers. Many mathematical methods have been developed by several researchers in the past to deal with the issue of hydrogen network synthesis. The pivotal aspect of this is to minimize the pure hydrogen consumption; however, of late even minimizaSpecial Issue: Energy System Modeling and Optimization Conference 2013 Received: March 26, 2014 Revised: June 26, 2014 Accepted: July 14, 2014
A
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
changes in the feedstock or product specifications. The basic motivation for the optimization of hydrogen network design under uncertainty is that the phenomenon of uncertainty in the parameters may not be fully captured in the design phase. This is because during the design phase only the average values of the parameters are used over its entire operating period. The variations or uncertainty in the values of the parameters is only known during the operational phase. In such cases the optimal design predicted during the design phase may be suboptimal or even infeasible for the operational phase. Hence there is a need to incorporate these uncertainties during the design stage of hydrogen networks. The major task that has to be remembered is that the network obtained during the design stage remains feasible and optimal over the entire range of uncertain operating parameters. To cater to the issue of uncertainty, the stochastic programming13−15 technique is used to model the optimization problem. The stochastic programming approach for a process network (hydrogen network as in our case) problem identifies the degree of flexibility with which the process network can be designed such that all realization of uncertain parameters can be balanced and the network can operate optimally. The aim is to find a solution which is feasible, and to some extent optimal, to all possible values of the uncertain parameters by optimizing the objective function. Stochastic programming assumes that even though parameters are uncertain, they are assumed to lie within some given values. To deal with the design of hydrogen networks under uncertainty, in this work we model the optimization problem as a two stage stochastic formulation. In a typical two stage approach, a first stage decision is taken and then random event occurs which could affect the first stage decision taken. Each random event may have some probability of occurrence associated with it. In the second stage, the optimization is carried out in the presence of uncertainty or when the random event occurs. This can be possibly considered as a recourse action taken which optimally hedges or compensates for any possible inconsistency occurring due to first stage decisions. The uncertainty in the second stage decisions is translated into the optimization problem by considering an objective which has the expected value of the second stage decisions. Hence, the objective in a two stage approach involves the first stage decision along with expected second stage or recourse decisions. The deterministic equivalent of a two stage stochastic programming problem16 can be given as follows.
tion of cost of the network seems to be an important aspect. The process integration approaches used in this regard are the pinch analysis and mathematical programming. Since this work deals with mathematical programming, reviewing of literature will focus on the mathematical programming approaches. In the mathematical programming approach, the task will be to construct a mathematical model for all the hydrogen producing (hydrogen sources), consuming (processing units), purifying (purification units), and accepting units (fuel gas sinks) in a refinery and to develop equations showing hydrogen balance for all these units. This mathematical model will depict the hydrogen distribution among these units. Optimization of the above-mentioned model is carried out based on a particular objective, and this gives an optimal network topology for hydrogen distribution among the units. Hallale and Liu1 developed an efficient mathematical programming model for refinery hydrogen distribution for an existing refinery network incorporating pressure constraints, compressors, and purification units in their model. Their objective was to minimize cost and hydrogen consumption in a refinery. Zhang et al.2 integrated hydrogen distribution with refinery planning and a utility system and pointed out the strong interaction among them through a mathematical model. Liu and Zhang3 developed a systematic methodology for appropriate placement of purifier in hydrogen network design. To solve this MINLP problem, they advocated the usage of solution from the relaxation of nonlinear terms as an initial point to the original MINLP model. Fonseca et al.4 developed an adapted LP approach to deal with hydrogen distribution in a refinery. They applied the model on a real case of Porto refinery of GALP Energia network and showed benefits in electricity cost, hydrogen consumption, and gas in fuel system. Khajehpour et al.5 applied engineering insights and developed a reduced superstructure approach to model the hydrogen network distribution in a refinery. They used a genetic algorithm to solve their model. Kumar et al.6 developed and compared LP, NLP, MILP, and MINLP models for the hydrogen network and concluded that MINLP was more suitable and efficient in the design of the hydrogen network. The hydrogen network optimization model along with a strategy for retrofitting an existing network with a purifier was given by Liao et al.7 An improved optimization approach for optimum hydrogen distribution in a refinery involving multiple components along with hydrogen was carried out by Jia and Zhang.8 Their work involved complex structure for hydrogen consumers by having a detailed hydrogen reactor model along with integrated flash calculations for the hydrogen separator. Elkamel et al.9 incorporated a hydrogen network model formulation with petroleum refinery planning and studied combined synergistic interactions with respect to different modes of refinery operation. Jagannath10 improved few modeling aspects of the hydrogen network model of Elkamel et al.9 and solved this model to global optimality. Alhajri et al.11 combined their CO2 emission model with the hydrogen management model from Elkamel et al.9 and presented an overall framework for refinery planning by integrating the refinery planning model with hydrogen management and CO2 emissions. A discrete time MINLP formulation for hydrogen distribution scheduling in a refinery was given by Jiao et al.12 Most of the mathematical optimization models in the literature for hydrogen network design have assumed a fixed set of parameters for all the units in the hydrogen network. These, however, could be uncertain in a real operation due to
min z = c Tx + Eωfω (yω , θω) s.t.
gω(x , yω , θω) = 0
∀ω
hω(x , yω , θω) ≤ 0
∀ω
x ≥ 0, yω ≥ 0
(P1)
In the above problem, x represents the vector of first stage variables, yω and θω are the vector of second stage variables and vector of uncertain parameters, respectively, and ω is the random event or possible outcome. Eω is the mathematical expectation of the random or possible outcome. The first term in the right-hand side of the objective function is the first stage decision, and the other is the second stage decision. gω and hω are the equality and inequality constraints for this problem, respectively. For a typical network design problem, the equality B
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
superstructure synthesis options and alternatives. Second, to the best of our knowledge, there is no work as of now which focuses on the global optimization of the hydrogen network design under uncertainty. The proposed solution strategy helps in solving models to acceptable global solutions within tractable computational time. The layout of this article is as structured as follows. We first describe the problem statement in section 2 and then explain our mathematical model in section 3. In section 4, we describe how the nonlinear terms in our model are convexified and then, with this basis, we present our solution strategy to solve the model to global optimality. In section 5, two case studies are presented to illustrate the use of the mathematical model and the solution strategy to solve the hydrogen network problem under uncertainty to global optimality. Finally, the conclusions are presented in section 6.
constraints are in the form of a mass, component, or energy balance. The inequality constraint could be a capacity constraint or a logical constraint. In problem P1, the uncertain parameters can take finite number of realizations and the uncertainty for each realization can be characterized through different scenarios with corresponding probability pω. Incorporating this into the mathematical model depicted by P1, the multiscenario two stage stochastic model is now given as follows.17 min z = c Tx +
∑ pω fω (yω , θω) ω
s.t.
hω(x , yω , θω) ≤ 0
∀ω
gω(x , yω , θω) = 0
∀ω
x ≥ 0, yω ≥ 0
(P2)
2. PROBLEM STATEMENT In this work, our aim is to develop a mathematical model for the design of a hydrogen network under uncertainty. We extend the mathematical model of the refinery hydrogen network synthesis problem given in Jagannath et al.23 (submitted) to account for uncertainty. In this section, we briefly state the different elements of our model, namely, the units involved in the hydrogen network, uncertain parameters that could exist in the model, and the objective function. The units of the refinery hydrogen network synthesis problem under uncertainty consist of I sources (i = 1, 2, ..., I), M processing units (m = 1, 2, ..., M), N purification units (n = 1, 2, ..., N), and J fuel gas sinks (j = 1, 2, ..., J) all operating under different scenarios of operation S (s = 1, 2, ..., S). The source i (i = 1, 2, ..., I) supplies hydrogen to the entire network, and the stream from source i is characterized by known flow, purity, pressure, and temperature for its hydrogen stream for all s (s = 1, 2, ..., S) scenarios of operation. The processing units m (m = 1, 2, ..., M) are the major consumers of hydrogen in the network. They have known hydrogen flow demand, purity, perpass conversion, and temperature and pressure requirements for s (s = 1, 2, ..., S) scenarios of operation. The stream from the processing unit has known purity, pressure, and temperature for s (s = 1, 2, ..., S) scenarios. A purification unit n (n = 1, 2, ..., N) provides high purity hydrogen stream in the network by purifying the lower purity hydrogen stream. They have an inlet stream and two outlet streams, namely, product (hydrogen rich) and residue stream (lean hydrogen). For the purification unit n (n = 1, 2, ..., N), the hydrogen recovery ratio and product purity are known for s (s = 1, 2, ..., S) scenarios of operation. The temperature and pressure for the inlet feed stream to the purification unit are known for s (s = 1, 2, ..., S) scenarios of operation. The product and residue stream temperature and pressure are also known for s (s = 1, 2, ..., S) scenarios of operation for the same. A fuel gas sink j (j = 1, 2, ..., J) consumes the gas unutilized in the network. They have known limits on flow and purity for all s scenarios of operation (s = 1, 2, ..., S). Also the pressure and temperature of operation of fuel gas sinks are known for s scenarios of operation (s = 1, 2, ..., S). Each source i (i = 1, 2, ..., I) sends its stream to processing units, purification units, and fuel gas sinks. The streams from each of the processing units m (m = 1, 2, ..., M) go to other processing units including itself, purification units, and fuel gas sinks. The product stream from each of the purification unit, n (n = 1, 2, ..., N) goes to other purification units excluding itself and processing units, and the residue stream goes to the fuel gas
The two stage stochastic programming models, similar to the structure shown by P2, have been widely used in the literature17−19 for solving network synthesis problems under uncertainty. The main challenge involved in stochastic network design problems is that the size of the problem is large and increases with scenarios, in comparison to the deterministic problems. Study of flexible or network design under uncertain conditions for hydrogen networks was undertaken by Ahmad et al.,20 Jiao et al.,21 and Lou et al.,22 but all these works were solved only to locally optimal solutions. Although locally optimal solutions seem to be reasonable, there is always a need to look out for the “best feasible” solution to an optimization problem. In this article, we develop a mathematical model for addressing the design of hydrogen networks under uncertainty. This is an extension of an earlier work by Jagannath et al.23 (submitted) which deals with the design of hydrogen networks under steady or constant operating conditions (all parameters being constant). To address the uncertain operating conditions, we segment or divide the operation of refinery into different scenarios and the parameters in the model take different values in each operating scenario. This problem is now formulated as a multiscenario two stage stochastic MINLP for addressing the design of hydrogen network in a refinery under uncertain operational conditions. The objective of this is to minimize the total annualized cost of the entire network. A specific or tailormade solution strategy is implemented to solve this model to global optimality by convexifying the nonlinear terms in the model. Finally two case studies are solved to illustrate the applicability of the proposed model and the solution strategy. The contribution of this article is twofold. First, the mathematical model of Jagannath et al.23 (submitted) has been extended to account for uncertainty. Although the model of Jagannath et al.23 (submitted) has many features in their superstructure such as hydrogen purification, hydrogen reuse and recycle, and presence of compressors, valves, heaters, and coolers to accommodate nonisobaric and nonisothermal conditions, their model is applicable for only a single set of operating conditions. This means that the model of Jagannath et al.23 (submitted) has to be adapted to handle variable or uncertain operating conditions. This involves designing a network which is flexible and resilient in handling such uncertainties. This is the first work, to the best of our knowledge, which has developed a mathematical model for hydrogen network design under uncertainty, incorporating rich C
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
Figure 1. (a) Schematic representation of a transfer line. (b) Schematic representation of connections from the source unit. (c) Schematic representation of connections to/from the processing unit. (d) Schematic representation of connections to/from the purification unit. (e) Schematic representation of connections to the fuel gas sink.
ing units. Pressure swing adsorption and membrane separations are the most common purification units used. Fuel gas sinks could be in the form of flares, gas turbines, boilers, and burners. Gas turbines and boilers consume unutilized gas streams and produce energy for the refinery in the form of steam, heat, and power. In this process, the refinery could have monetary benefit by the surplus energy generated by them. Fuel sinks like flare and burners are mainly used for disposing unused gas stream and usually have operational cost associated with them for disposal. Apart from modeling the flow balances, it is also necessary to model the changes in the stream conditions (pressure and temperature) taking place in each of the transfer lines because of the equipment, namely, valve, compressor, heater, and
sinks. All the above-mentioned connections are represented by a transfer line. Since the network design is nonisothermal and nonisobaric, different gas stream conditioning equipment such as valves, compressors, heaters, and coolers are present on the transfer line. These bring gas streams to the required final conditions of pressure and temperature. The schematic representation of the transfer line is shown in Figure 1a. The connections, in the form of transfer lines, from/to hydrogen sources, processing units, purification units, and fuel gas sinks are shown in Figure 1b−e, respectively. Typical examples of sources would be steam reforming and import of hydrogen from the outside. Catalytic reforming also produces hydrogen as a byproduct stream and can be classified as a source. Hydrocrackers and hydrotreaters are the conventional processD
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
above understanding, we now summarize the hydrogen network synthesis problem under uncertainty as follows. Given: 1. S scenarios (s = 1, 2, ..., S) also known as scenarios of operation of the network 2. I hydrogen sources (i = 1, 2, ..., I) with known profiles for stream flows, temperatures, pressures, and purities for S (s = 1, 2, ..., S) scenarios. 3. M processing units (m = 1, 2, ..., M) with known profiles for inlet flows, purity, reactor per pass conversion, pressure, and temperature for S (s = 1, 2, ..., S) scenarios. The outlet stream is also characterized by known profiles for stream purity, pressure, and temperature for S (s = 1, 2, ..., S) scenarios. 4. At most N purification units (n = 1, 2, ..., N) with known profiles for pressures, temperatures, and flow ranges for its feed for S (s = 1, 2, ..., S) scenarios; known profiles for S (s = 1, 2, ..., S) scenarios for pressures, temperatures, hydrogen recoveries, and hydrogen purities in the hydrogen rich or product streams; temperatures and pressures profile for the residue streams for S (s = 1, 2, ..., S) scenarios. 5. J fuel gas sinks (j = 1, 2, ..., J) with known profiles for temperatures, pressures, flow ranges, and purity ranges for S (s = 1, 2, ..., S) scenarios. 6. Capital expenditure (CAPEX) and operating expenditure (OPEX) for all the conditioning units, transfer lines (pipelines), and purification units. Also given is the cost of power consumption for all the stream conditioning units in all scenarios. 7. Economic data profiles for S (s = 1, 2, ..., S) scenarios for the cost of feed (hydrogen), cost of utilizing (economic returns), and disposing gas streams in fuel gas sinks. Determine: 1. Flows, purities, temperatures, and pressures at all the points in the network for all scenarios of operation. 2. Network topology with existence, duties of all the equipment, conditioning units, and transfer lines. Aim: Minimize the total annualized cost (TAC) of the hydrogen network under uncertainty. Assumptions: 1. For s (s = 1, 2, ..., S) scenarios, there exists probability of occurrence of each scenario denoted by πs. 2. The model is built based on the consideration of total flow in contrast to the component flow. The components in the total flow model are the hydrogen component and the non-hydrogen component. For simplicity and enabling gas property calculations, we assume methane as the non-hydrogen component. 3. All compression processes are single-stage and adiabatic with no restrictions on their inlet and outlet temperatures. 4. All expansions are Joule-Thompson expansions via valves. 5. Zero pressure drops exist across heaters, coolers, and transfer lines. 6. All the flows considered have the unit of MMscf/day, and purity of gas is given by mole fraction. Since the flows are in standard volumetric flows, the model is developed assuming standard conditions of 14.7 psia and 60 °F.
cooler. Since this sequence is generic for all the transfer lines representing all the connections, we need a generic description for each transfer line. To describe a generic transfer line SSpq, we categorize all these units in the network as origin p (p = 1, 2, ..., P) and destination q (q = 1, 2, ..., Q) units. An origin (destination) unit is one in which there exists a product (feed) stream from (to) the unit. Hydrogen sources are the potential origin units whereas fuel gas sinks are potential destination units. Processing and purification units are both potential origin and destination units. The stream conditioning equipment potentially exists along the transfer line. Valves and compressors change the pressure of the stream, and the temperature change of a gas stream is caused by heaters and coolers. The stream conditioning may or may not be required for all the connections in the network. The refinery hydrogen network could be subjected to uncertainty with respect to the hydrogen demand at the inlet of the processing unit and amount of hydrogen exiting from the processing unit. The uncertainty in the former is because of the quality of the crude, whereas the uncertainty in the latter is because of variations in the operating conditions of the processing unit such as product specification, mode of operation of refinery, catalyst activity, etc. To characterize the above-mentioned uncertainties within the framework of our model, we consider two uncertain parameters in our model. For the inlet hydrogen demand at the processing unit, we consider the inlet flow requirement of gas into the processing unit as uncertain. The uncertainty in the amount of hydrogen exiting the processing unit is represented through uncertain reactor per pass conversion. Such uncertainties could affect the optimal operation of the network. To represent a resilient network operation, it is necessary to incorporate the uncertainty design stage itself. The above-mentioned parameters now are considered as uncertain parameters in the model. The network design is formulated as a two stage stochastic program, where in the first stage the network is selected and in the second stage the network operation is carried out. The network operation is divided into many scenarios, and the uncertain parameter takes different values for each operating scenario. The objective function considered for the hydrogen network design or synthesis problem under uncertainty is to minimize the Total Annualized Cost (TAC). The TAC typically consists of two parts, namely, capital/investment expenditure (CAPEX) and operational expenditure (OPEX). As mentioned already, a two stage stochastic formulation is considered in this model which has two stages of decisions, namely, first and second stage; the first stage decisions are made in a way that the expected value of the second stage decision is minimized. In our problem, the first stage involves decisions regarding the design and topology of the network which include existence of equipment such as purification units, transfer lines, and stream conditioning equipment and their design capacities. The second stage decisions cover the aspects regarding the operational performance of the network. This includes flow of gas through transfer lines for each scenario, duty of the entire stream conditioning equipment for each scenario, amount of hydrogen gas used for each scenario, amount of gas sent to fuel gas sinks, etc. From the cost point of view, the first stage decisions represent the capital expenditure (CAPEX) and are independent of the different scenarios of operation, whereas the second stage decisions depict the operational expenditure (OPEX) of the network and are present for all the scenarios. With the E
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article I
3. MODEL FORMULATION 3.1. Material Balance. The two stage stochastic recourse model will be adopted for the design of hydrogen network under uncertainty. The uncertainty in this model will be realized through different scenarios in which uncertain parameters will take different values. This model is based on the MINLP model for the hydrogen network synthesis by Jagannath et al.23 In this section, we first write the material and hydrogen balance for all the units in the hydrogen network. The schematic representation is shown for all the units along with their inlet and outlet connections (Figure 1a−e). The schematic representation for all the units along with their corresponding connections can be visualized as a superstructure for the hydrogen network design. As mentioned previously, the connections are represented by transfer line which will have stream conditioning equipment like cooler, compressor, valve, and heater as a part of its superstructure. For the material balance equations corresponding to different units in the hydrogen network, the units along with their mixers (M) and splitters (S) are lumped together as single entities and equations are written only for these units/entities. In other words, we do not write separate material and component balance equations for splitter and mixers accompanying the units in the network. This is done to decrease the number of variables and constraints in the model. To avoid increasing the length of the article, we only briefly state the material balance equations for all the units in the hydrogen network. For a comprehensive description of the units along with the detailed explanation for the material balance equations, the reader is referred to the works of Jagannath et al.23 (submitted) and Jagannath10 (2012). The description for the indices, parameters and variables are given in the nomenclature section toward the end of this article. Source. The source is the supplier of hydrogen gas to the entire network. It consists of the unit along with a splitter (Figure 1b). The mass balance for source i in scenario s is given in eq 1 with the variable bounds given by eqs 2a− 2c. J
Fis =
M m=1
N
M
∑ Fimszis + ∑ Fnmsypns i=1
∑
+
n=1
Fm ′ msyms
m ′= 1
∀ m = 1, 2, .... , M ;
∀ s = 1, 2, ..., S I j=1
∀ m = 1, 2, .... , M ;
(4)
N
Fmsxms(1 − χms ) = yms (∑ Fmjs +
M
∑ Fmns + ∑ n=1
Fmm ′ s)
m ′= 1
∀ s = 1, 2, ..., S
(5)
⎛⎡ F x (1 − χ ) ⎤ ⎞ ms ms ms ⎥ 0 ≤ Fmjs ≤ min⎜⎜⎢ , F Ujs ⎟⎟ ⎥⎦ yms ⎝⎢⎣ ⎠ ∀ m = 1, 2, .... , M ;
∀ j = 1, 2, ..., J ;
∀ s = 1, 2, ..., S
(6a)
⎛⎡ F x (1 − χ ) ⎤ ⎞ ms ms ms ⎥ 0 ≤ Fmns ≤ min⎜⎜⎢ , FnsU ⎟⎟ ⎥⎦ yms ⎝⎢⎣ ⎠ ∀ m = 1, 2, .... , M ;
∀ n = 1, 2, ..., N ;
∀ s = 1, 2, ..., S
(6b)
⎛⎡ F x (1 − χ ) ⎤ ⎞ ms ms ms ⎥ 0 ≤ Fmm ′ s ≤ min⎜⎜⎢ , Fms⎟⎟ ⎥⎦ yms ⎝⎢⎣ ⎠ ∀ m = 1, 2, .... , M ;
∀ m′ = 1, 2, ..., M ;
∀ s = 1, 2, ..., S
(6c)
Purification Unit. The equations for the purification unit n in scenario s are shown in eqs 7−11. The product (hydrogen rich) stream goes to the purification unit splitter, which in turn directs streams to other purification units and processing units. The residue (hydrogen lean) stream only goes to the fuel gas sinks. I
N
∑ Fijs + ∑ Fims + ∑ Fins j=1
Fmsxms =
Fns =
∀ i = 1, 2, .... , I ;
M
N
M
N
∑ Fins + ∑ Fmns + ∑ Fn′ ns = ∑ Fnms + ∑ Fnn′ s i=1
m=1
n ′= 1
m=1
n ′≠ n
n=1
n ′= 1 n ′≠ n
J
∀ s = 1, 2, ..., S
∑ Fnjs
+
(1)
∀ n = 1, 2, .... , N ;
∀ s = 1, 2, ..., S
j=1
0 ≤ Fijs ≤ min(FisU , F Ujs ) ∀ j = 1, 2, ..., J ;
(7)
∀ i = 1, 2, .... , I ;
∀ s = 1, 2, ..., S
Fnjs = F1njs + F 2njs
(2a)
∀ n = 1, 2, .... , N ;
∀ j = 1, 2, ..., j ; 0 ≤ Fims ≤
min(FisU , Fms)
∀ m = 1, 2, ..., M ;
∀ i = 1, 2, .... , I ; ∀ s = 1, 2, ..., S
I
M
rns(∑ Finszis +
(2b)
∀ s = 1, 2, ..., S
i=1
(8)
N
∑ Fmnsyms + ∑ Fn′ nsypns ) = y m=1
n ′= 1 n ′≠ n
0 ≤ Fins ≤ min(FisU , FnsU ) ∀ n = 1, 2, ..., N ;
M
∀ i = 1, 2, .... , I ; ∀ s = 1, 2, ..., S
N
pns ( ∑ Fnms +
(2c)
m=1
∑ Fn′ ns)
∀ n = 1, 2, .... , N ;
n ′= 1 n ′≠ n
Processing Unit. The flow, hydrogen balance, and bounds for processing unit m in scenario s are shown by eqs 3−6, I
Fms =
N n=1
M
M
∑ Fims + ∑ Fnms + ∑ i=1
∀ s = 1, 2, ..., S
(1 − rns)ypns ( ∑ Fnms +
Fm ′ ms
m=1
m ′= 1
∀ m = 1, 2, .... , M ; ∀ s = 1, 2, ..., S
(9) N
J
∑ Fnn′ s) = rns ∑ F1njs n ′= 1
j=1
n ′≠ n
∀ n = 1, 2, .... , N ;
(3) F
∀ s = 1, 2, ..., S
(10)
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
Figure 2. Generalized superstructure representation for connection between origin and destination units.
0 ≤ Fns ≤ FnsU
∀ n = 1, 2, .... , N ;
I
∀ s = 1, 2, ..., s
FjsxjsL ≤
(11a)
0 ≤ Fnms ≤ min(FnsU , Fms) ∀ m = 1, 2, ..., M ;
0 ≤ Fjs ≤ F Ujs
∀ s = 1, 2, ..., S
∀ n = 1, 2, ..., N ;
∀ s = 1, 2, ..., S
(11d)
∀ j = 1, 2, ..., J ;
∀ s = 1, 2, ..., S
(11e,f)
Fuel Gas Sinks. The material and hydrogen balance for fuel gas sinks for all scenarios is given in eqs 12 and 13 with the bound on fuel gas sinks shown by eq 14. As mentioned in Jagannath et al.,23 the fuel sinks may have inlet gas specifications in terms of density, energy demand, heating value of gas etc., but here we consider only flow and purity requirements. The equations for the fuel gas sink j in scenario s are I
Fjs =
M
N
N
∑ Fijs + ∑ Fmjs + ∑ F1njs + ∑ F 2njs i=1
m=1
∀ j = 1, 2, ..., J ;
n=1
n=1
∀ s = 1, 2, ..., S
∀ s = 1, 2, ..., S
∀ j = 1, 2, ..., J ;
(13)
∀ s = 1, 2, ..., S
3.2. Pressure and Temperature. The description for the indices, parameters, and variables in this section are stated in the nomenclature section toward the end of this article. The units in the hydrogen network operate at different pressures and temperatures. Figure 2 shows the generalized superstructure representation for the stochastic hydrogen network connecting the origin unit p (p = 1, 2, ..., P) with the destination unit q (q = 1, 2, ..., Q) through transfer line SSpq. Since all the units in the network can be categorized as either origin and/or destination units, we adopt a generic way of expressing the mathematical equations of temperature and pressure between all the connections in the hydrogen network. The gas coming out (in) from (to) the origin (destination) unit must have the same temperature and pressure of the unit. When the gas is transported from an origin unit p (p = 1, 2, ..., P) to a destination unit (q = 1, 2, ..., Q), its temperature and pressure must undergo changes. This is done through the presence of the stream conditioning equipment. Since the order of the existence of these equipment along the transfer line is the same, we write a generic equation for a transfer line which is applicable to all the transfer lines in the entire network. Similar to the equations of Jagannath et al.,23 we write down equations necessary for a transfer line connecting origin p (p = 1, 2, ..., P) to destination q (q = 1, 2, ..., Q) in scenario s (s = 1, 2, ..., S). For more details regarding the model equations of the transfer
(11c)
0 ≤ F1njs ≤ min(FnsU , F Ujs ); 0 ≤ F 2njs ≤ min(FnsU , F Ujs ) ∀ n = 1, 2, .... , N ;
n=1
(14)
∀ n = 1, 2, .... , N ;
∀ n′(n′ ≠ n) = 1, 2, ..., N ;
∀ j = 1, 2, ..., J ;
(11b)
N
m=1
∀ j = 1, 2, ..., J ;
∀ s = 1, 2, ..., S
0 ≤ Fnjs ≤ min(FnsU , F Ujs )
i=1
∀ n = 1, 2, .... , N ;
0 ≤ Fnn ′ s ≤ min(FnsU , FnU′ s)
M
∑ Fijszis + ∑ Fmjsyms + ∑ F1njs ≤ FjsxUjs
(12) G
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
⎧C TOUT F − ΔH C if POUT > PIN ps pqs pqs ps qs ⎪ ps L CpsTpqs Fpqs ≤ ⎨ C ⎪CpsTOUTpsFpqs − ΔH pqs if POUTps ≤ PINqs ⎩
line, the reader can refer to the previous works of Jagannath et al.23 and Jagannath et al.24 The equations are shown as follows.
HINpqs
⎧C TOUT F − ΔH C − ΔHV if POUT > PIN ps pqs pqs pqs ps qs ⎪ ps H ⎪ + ΔHpqs ⎪ ⎪ C B = ⎨CpsTOUTpsFpqs − ΔH pqs if POUTps < PINqs + ΔHpqs ⎪ H + Δ H ⎪ pqs ⎪ C H ⎪CpsTOUTpsFpqs − ΔH pqs if POUTps = PINqs + ΔHpqs ⎩
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
∀ p = 1, 2, ..., P ;
HINpqs ≤ CpsT UpqsFpqs ∀ q = 1, 2, ..., Q ; V ΔHpqs
(17)
∀ q = 1, 2, ..., Q ;
∀ s = 1, 2, ..., S
∀ q = 1, 2, ..., Q ;
∀ s = 1, 2, ..., S
p=1
∀ q = 1, 2, ..., Q ;
∀ s = 1, 2, ..., S
∑ HINpqs
∀ q = 1, 2, ..., Q ;
p=1
(20)
To effectively construct the feasible region of an optimization problem, appropriate bounds have to be given to all the variables. This is done as follows. 0 ≤ Fpqs ≤ F Upqs
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
∀ s = 1, 2, ..., S
(21a)
0 ≤ HINpqs ≤ CpsF UpqsT Upqs (19)
∀ q = 1, 2, ..., Q ; B 0 ≤ ΔHpqs ≤
∀ p = 1, 2, ..., P ;
∀ s = 1, 2, ..., S
(21b)
nps ⎤ (CpsTOUTpsF Upqs) ⎡⎛ PINqs ⎞ ⎢⎜ ⎥ ⎟ − 1 ⎢⎣⎜⎝ POUTps ⎟⎠ ⎥⎦ η
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
∀ s = 1, 2, ..., S
(21c)
V 0 ≤ ΔHpqs ≤ μpqs CpsF Upqs(POUTps − PINqs)
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
∀ s = 1, 2, ..., S
(21d)
H L 0 ≤ ΔHpqs ≤ CpsF Upqs(T Upqs − Tpqs )
∀ q = 1, 2, ..., Q ;
∀ p = 1, 2, ..., P ;
∀ s = 1, 2, ..., S
(21e)
C L 0 ≤ ΔH pqs ≤ CpsF Upqs(TOUTps − Tpqs )
⎧C TOUT F − ΔH C + ΔHV if POUT > PIN ps pqs pqs pqs ps qs ⎪ ps H ⎪ + ΔHpqs ⎪ ⎪ C B = ⎨CpsTOUTpsFpqs − ΔH pqs if POUTps < PINqs + ΔHpqs ⎪ H + Δ H ⎪ pqs ⎪ C H ⎪CpsTOUTpsFpqs − ΔH pqs if POUTps = PINqs + ΔHpqs ⎩
∀ p = 1, 2, ..., P ;
(18a)
∀ s = 1, 2, ..., S (18)
When the Joule-Thompson coefficient of the gas is positive, eqs 15-19 hold good for all connections between origin p and destination q. However, this may also tend to go negative under some circumstances. In such conditions, eqs 15a, 16a, and 18a are written in place of eqs 15, 16, and 18. The difference between eq 15 and eq 15a exists in their expression when POUTps > PINqs. For this condition POUTps > PINqs, in eq 15 the coefficient of ΔHvpqs is −1 whereas in eq 15a the coefficient is +1. This is because of the change in the sign of the JouleThompson coefficient of hydrogen gas as seen in eq 18 and eq 18a. In eq 16a, in comparison to eq 16, the term ΔHvpqs is not present in the expression for the condition POUTps > PINqs. This is because the negative coefficient Joule-Thompson coefficient causes an increase in the gas temperature and hence it is excluded in determining the lowest possible temperature attainable in the transfer line connecting origin p (p = 1, 2, ..., P) to destination q (q = 1, 2, ..., Q) in scenario (s = 1, 2, ..., S).
HINpqs
∀ s = 1, 2, ..., S
P
TINqs ∑ CpsFpqs =
⎧ C ⎪ (CpsTOUTpsFpqs − ΔH pqs) if POUT < PIN ps qs ⎪ η ⎪ ⎡ nps ⎤ ⎛ ⎞ PIN ⎨ = qs ⎟⎟ − 1⎥ ⎪ ⎢⎜⎜ ⎥⎦ ⎪ ⎢⎣⎝ POUTps ⎠ ⎪ ⎩ 0 otherwise
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
P
⎧ C F (POUTps − PINqs) if POUTps > PINqs ⎪μ pqs ps pqs =⎨ ⎪ ⎩ 0 otherwise
∀ p = 1, 2, ..., P ;
B ΔHpqs
(16)
∀ p = 1, 2, ..., P ; ∀ s = 1, 2, ...S
(16a)
Equations 15−19 make sure that for all transfer lines connecting origin unit p to destination q, pressure of the stream in all the transfer lines going to the mixer of destination q in all the scenarios is the same as that of the destination unit pressure PINqs. However, the temperatures of the streams in the transfer line going to the inlet mixer of destination q in all scenarios may not be the same as that of the temperature of the destination unit TINqs. To ensure this condition, the following temperature balance constraint is added to the model. Similar to other constraints, this should hold for all the scenarios of operation.
(15)
∀ s = 1, 2, ..., S
∀ s = 1, 2, ..., S
∀ s = 1, 2, ..., S
⎧ C F (POUTps − PINqs) if POUTps > PINqs ⎪− μ pqs ps pqs V ΔHpqs =⎨ ⎪ ⎩ 0 otherwise
⎧C TOUT F − ΔH C − ΔHV if POUT > PIN ps pqs pqs pqs ps qs ⎪ ps L CpsTpqs Fpqs ≤ ⎨ C ⎪ CpsTOUTpsFpqs − ΔH pqs if POUTps ≤ PINqs ⎩ ∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
∀ p = 1, 2, ..., P ; ∀ s = 1, 2, ..., S
∀ q = 1, 2, ..., Q ; (21f)
3.3. Linking Constraints. Similar to other subsections, the description of the parameters and variables in this subsection are stated in the nomenclature section toward the end of this article. The linking constraint connects or links the scenario independent variables with the scenario dependent variables.
(15a) H
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
The design or scenario independent variables are the existence of transfer lines, purification units, and gas stream conditioning equipment like compressors, valves, heaters, and coolers. Their presence does not depend on the different scenarios of operation. All the other variables in the network change according to the different scenarios of network operation considered and hence are called scenario dependent variables. Equation 22 shows the link between the scenario independent variable Fpq and the scenario dependent variable Fpqs. Fpq stands for flow from origin p to destination q whereas Fpqs stands for flow from origin p to destination q in scenario s. This means that the first stage design variable or scenario independent variable Fpq must always be greater than its corresponding second stage scenario dependent variable Fpqs. In other words, this physically means that the capacity of the transfer line connecting origin p and destination q should be large enough such that it can accommodate the flows for every scenario. The existence of the design variable Fpq is denoted by the use of the binary variable XFpq as shown in eq 23, where FUpq and FLpq give the upper and lower limits on the transfer line capacity, respectively. If XFpq = 1, then the pipeline (transfer line) between origin p and destination q is developed or exists; eq 23 then shows the capacity of the pipeline between its bounds FUpq and FLpq. In another case, when XFpq = 0, this forces the flow Fpq to be zero indicating no transfer line exists. FLpq is the value of the minimum flow rate for which decision should be made on transfer line existence. The value of FUpq on the other hand depends on the maximum of bound on the variable Fpqs examined for all scenarios of operation s. The equations are shown as follows. Fpq ≥ Fpqs
∀ p = 1, 2, ..., P ;
B L B B B U B (ΔHpq ) X ΔHpq ≤ ΔHpq ≤ (ΔHpq ) X ΔHpq
∀ p = 1, 2, ..., P ; B U (ΔHpq )
⎤⎞ − 1⎥⎟ ⎥⎦⎟ ⎠
∀ q = 1, 2, ..., Q
(29)
⎛ S (C TOUT F U ) ⎡⎛ PIN ⎞nps ps ps pqs ⎢ qs ⎜⎜ ⎟⎟ = max⎜∑ ⎜ ⎢ η POUT ps ⎠ ⎣⎝ ⎝ s=1 ∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q (30)
V V ΔHpq ≥ ΔHpqs
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
∀ s = 1, 2, ..., S
(31)
V L V V V U V (ΔHpq ) X ΔHpq ≤ ΔHpq ≤ (ΔHpq ) X ΔHpq
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q
(32)
S V U (ΔHpq ) = max(∑ μpqs CpsF Upqs(POUTps − PINqs)) s=1
∀ p = 1, 2, ..., P ; H H ΔHpq ≥ ΔHpqs
∀ q = 1, 2, ..., Q
∀ p = 1, 2, ..., P ;
(33)
∀ q = 1, 2, ..., Q ;
∀ s = 1, 2, ..., S
(34)
H L H H H U H (ΔHpq ) X ΔHpq ≤ ΔHpq ≤ (ΔHpq ) X ΔHpq
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
∀ q = 1, 2, ..., Q
(35)
S
∀ s = 1, 2, ..., S
(22)
H U L (ΔHpq ) = max(∑ Cps(T Upqs − Tpqs )F Upqs) s=1
L Fpq XFpq ≤ Fpq ≤ XFpqF Upq
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q
∀ p = 1, 2, ..., P ; (23)
S
H L = min(∑ CpsΔTpq Fpq)
∀ p = 1, 2, ..., P ;
s=1
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q
s=1
∀ q = 1, 2, ..., Q
C C ΔH pq ≥ ΔH pqs
(24)
Similar to eqs 22−24, the model equation for the scenario independent design variable for the purification unit is shown in eqs 25−27. For other stream conditioning equipment such as compressor, valve, heater, and cooler, the linking constraints are given in eqs 28−30, eqs 31−33, eqs 34−37. and eqs 38−41, respectively. Fn ≥ Fns
(36)
S H L (ΔHpq )
F Upq = max(∑ F Upqs)
∀ q = 1, 2, ..., Q
∀ n = 1, 2, ..., N ;
∀ s = 1, 2, ..., S
(37)
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q ;
∀ s = 1, 2, ..., S
(38)
C L C C C U C (ΔH pq ) X ΔH pq ≤ ΔH pq ≤ (ΔH pq ) X ΔH pq
∀ p = 1, 2, ..., P ;
∀ q = 1, 2, ..., Q
(39)
S C U L (ΔH pq ) = max(∑ Cps(TOUTps − Tpqs )F Upqs)
(25)
s=1
FnLXFn
≤ Fn ≤
XFnFnU
∀ n = 1, 2, ..., N
∀ p = 1, 2, ..., P ;
(26)
∀ q = 1, 2, ..., Q
(40)
S
S
FnU = max(∑ FnsU )
∀ n = 1, 2, ..., N
s=1 B ΔHpq
≥
B ΔHpqs
C L C L (ΔH pq ) = min(∑ CpsΔT pq Fpq)
(27)
∀ q = 1, 2, ..., Q
∀ p = 1, 2, ..., P ;
∀ s = 1, 2, ..., S
∀ p = 1, 2, ..., P ;
s=1
(41)
3.4. Objective Function. The TAC of this multiscenario hydrogen network consists of two segments, namely, CAPEX and OPEX. The CAPEX is multiplied by an annualizing factor
∀ q = 1, 2, ..., Q ; (28) I
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
is incurred for every scenario and is the second stage decision or expected second stage cost and contains second stage variables. The objective is to minimize TAC, and this is given by eq 42.
(AF) and the OPEX is multiplied by on-stream time of refinery operation in a year (OP) and probability of occurrence of each scenario. CAPEX is incurred only at the time of network construction and consists of first stage decision variables. OPEX ⎡ P Q T T T (Fpq)dpq ) + TAC = AF ⎢ ∑ ∑ (apq XFpq + bpq ⎢ p=1 q=1 ⎣
P
Q
P
B
Q
p=1 q=1
p=1 q=1
V⎞ ⎛ P Q P Q ⎡ ΔHV ⎤dpq ⎟ ⎜ V pq C C C C dC V V⎢ pq ⎥ + ∑ ∑ (apqX ΔH pq + bpq(ΔH pq) )+ ∑ ∑ ⎜apqX ΔHpq + bpq ⎟+ ⎢⎣ (|μp |Cp) ⎥⎦ ⎟ p=1 q=1 p=1 q=1 ⎜ ⎝ ⎠
S ⎡ I + OP ∑ πs⎢∑ αisFis + ⎢ s=1 ⎣ i=1 P
+
Q
J
J
I
M
j=1
P
j=1
Q
p=1 q=1
i=1
m=1
Q
V ΔHpqs
dn
P
Q
P
Q
p=1 q=1
p=1 q=1
N
(42)
coefficients are OTpqs for transfer line, OBpqs for compressor, OHpqs for heater, OCpqs for cooler, OVpqs for valve, and on for purification unit.
S
∀ p = 1, 2, ..., P
4. GLOBAL OPTIMIZATION STRATEGY Equations 1−42 give the nonconvex MINLP model for the hydrogen network design under uncertainty. This section discusses the solution strategy required to solve this nonconvex MINLP to global optimality. For this purpose, we present a specialized deterministic global optimization approach. Most of the works on the deterministic global optimization approaches rely on constructing a convex underestimator for the nonconvex function involved in the problem. The optimization problem constructed using the convex underestimator is called as the lower bounding problem (for the case when objective is to be minimized). As the name indicates, the lower bounding problem gives the lower bound or the lowest possible feasible solution to the problem. Since the lower bounding problem has convex underestimators, the solution obtained by this problem is always global. The upper bound for the nonconvex MINLP is obtained by directly solving the MINLP. Any feasible solution to the MINLP forms a valid upper bound on the global optimum. These lower and upper bounds are made to converge within a specified tolerance limit ε. This process is done through a rigorous algorithmic framework where the bounds on the global optimum are constantly updated. Since the bounds are converged within a specified tolerance limit ε, the solution reached through deterministic methods is called as ε-global optimum.25 The commonly used algorithm for deterministic global optimization of MINLP is the spatial branch and bound method.25,26 In this algorithm, the nodes are enumerated and a branch and bound tree is developed. In every node of the branch and bound tree, a lower bounding problem and an upper bounding problem is solved. On the basis of the gap between the lower and upper bounding problem, the nodes are either branched further or are fathomed. This process is continued in an algorithmic framework until no node remains in the branch and bound tree and the algorithm terminates with ε convergence. Other known MINLP solution methods which could be modified to be solved to global optimality are benders decomposition18,27 and specialized outer approximation.28,29 In benders decomposition, the primal problem is made convex by replacing the nonconvex term by its convex underestimator in
s=1 S
Cp = max(∑ Cps)
n=1
⎤ V + ∑ ∑ opqs + ∑ onsFns ⎥ ⎥ (|μps |Cps) p=1 q=1 n=1 ⎦ P
where, μp = max(∑ μps )
N
⎤ ⎥ ∑ (anXFn + bn(Fn) )⎥ ⎥ n=1 ⎦ N
B B T Fpqs + ∑ ∑ opqs ΔHpqs ∑ γjsFjs − (∑ βjsLHVH2(∑ Fijszis+ ∑ Fmjsyms + ∑ F1njs )) + ∑ ∑ opqs
H H C C ΔHpqs + ∑ ∑ opqs ΔH pqs ∑ ∑ opqs p=1 q=1
H
∑ ∑ (apqBX ΔHpqB + bpqB(ΔHpqB)dpq ) + ∑ ∑ (apqHX ΔHpqH + bpqH(ΔHpqH)dpq )
∀ p = 1, 2, ..., P
s=1
The first six terms represent the annualized CAPEX cost for various equipment in the network. Here a, b, and d represent cost coefficients for the fixed investment cost, capital cost coefficient based on capacity of the equipment, and size factor of the equipment, respectively. The remaining terms represent the operating costs of the network which was incurred for all the scenarios, where each scenario has a probability of occurrence. Similar to Jagannath et al.23 (submitted), we have used nonlinear cost correlations to describe the capital cost of equipment at the expense of increasing the complexity of the model. The reason for this is that simple linear correlations are, in general, approximate and tend to overestimate the cost of the equipment. The nonlinear cost functions used, where the purchase cost of the equipment depends on size factor, are of the form, Equipment capital cost = a + b(Q )d
where, a and b are cost coefficients, Q is the size factor for the equipment, and d is the exponent or exponent factor for that equipment. The various cost coefficients are as follows: αis = Cost of source stream i in scenario s. This is positive when hydrogen feed is from hydrogen import or from hydrogen plant. When the feed is from units which produce hydrogen as a byproduct, this term is zero. γjs = Cost of operating sink j in scenario s. This is positive for sinks which burn/incinerate/dispose the gas such as flare or incinerator. This is zero in the case of sinks which produce energy such as boiler, furnace, and turbine. βjs = Economic value by sink j in scenario s. This is negative for sinks which generate energy like boiler, furnace, and turbine and is zero for sinks which burn or dispose the gas reaching fuel sinks such as flare and incinerators. The last six terms are the OPEX cost coefficients for various equipment in the network for different scenarios s. The cost J
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
The parameter γ can be used to attain uniform or nonuniform partitioning of variable domain.
the primal problem and continuing the benders decomposition algorithm. In the specialized outer approximation method, lower and upper bounding problems are solved and the gap between them is checked. The algorithm is terminated when the gap is smaller than a required tolerance criterion. If not then integer cuts are added to the lower bounding problem and the lower bound is updated until the gap is achieved. From the above, we understand that there exist two focal points in any global optimization strategy. First it is necessary to construct an underestimator for the nonconvex term and obtain the lower bounding problem. From this solution, any suitable heuristic has to be applied to obtain the solution of the original problem which serves as the upper bound. Second, the lower and upper bounds have to be updated in a suitable manner within an algorithmic framework. In the following two sub sections, we discuss these two focal points for our multiscenario hydrogen network synthesis problem. We first describe the convex underestimator for the nonlinear term, which gives the lower bounding problem, and then follow it up with the description of the global optimization algorithm. 4.1. Convex Underestimator for Concave Univariate Term. The nonconvex term in our model is the concave univariate term in the objective function which also happens to be nonlinear. This term is usually used to characterize the capital cost of equipment involved in the optimization formulation. The convex underestimator for this term is available in the literature30 and has been widely used in process network synthesis problems31,32 for expressing the capital cost of equipment. For a more clear understanding, the underestimator is represented for a generic concave univariate cost term xα. In the lower bounding problem the concave univariate term xα is replaced by w and is related to its underestimator wu as shown in eq 43. Geometrically eq 43 gives the secant chord of the concave term xα between xL and xU, where xL and xU are the lower and upper bounds of the variable x, respectively, and α is the value of the coefficient (0 < α < 1). ⎛ (xU )α − (x L)α ⎞ w ≥ wu = (x L)α + ⎜ ⎟(x − x L ) ⎝ xU − x L ⎠
⎛ nc − 1 ⎞γ U ⎜ ⎟ (x − x L) ⎝ NC ⎠
anc = x L +
∀ nc = 1, 2, ..., NC , NC + 1
(44) 35
Using the definition of Gounaris et al., we also use binary variable λnc to represent the condition that the binary variable is activated (or equal to one) if it lies within a partition or segment. This is mathematically stated as follows. ⎧ 1 anc ≤ x ≤ anc + 1 λnc = ⎨ ⎩0 otherwise
∀ nc = 1, 2, ...NC
It is known that at any given point of time only one variable is active or the variable x lies within any of the segment or partition. This is ensured by the following condition. NC
∑ λnc = 1
(45)
nc = 1
The convex hull or convex combination equations to model the piecewise partitioning of underestimator equations for the concave univariate term are given by eqs 46−48. dxnc is a continuous variable which is used as a place holder for segmented or partitioned portion of the variable involved in the concave univariate term. For a more information on the different type of formulations used in the piecewise underestimation of the underestimators (big M, convex hull, incremental cost, SOS-I, and SOS-II) and modeling different type of partitions (uniform, nonuniform, linear, and logarithmic), the reader is directed to some of the literature works34−37 in this area. Although refs 34−37 mentions them in context of the bilinear term, the concept can be extended to the concave univariate term in a straightforward manner. NC
x=
∑ dxnc
(46)
nc = 1
(43)
λnc . anc ≤ dxnc ≤ λnc ·anc + 1
As observed from eq 43, the underestimators are linear in nature. Underestimator equations similar to that of eq 43 are written for the concave univariate terms in our formulation T B H C V (Fpq)dpq,(ΔHBpq)dpq,(ΔHHpq)dpq,(ΔHCpq)dpq, (ΔHVpq)dpq and (Fn)dn. 33−35 it can be understood that the lower From recent works, bound (solution of the lower bounding problem) obtained by incorporating linear underestimators can be sometimes be weak and this slows down the convergence in a global optimization algorithm. In an attempt to obtain better (tighter) lower bound, the underestimator equations can be piecewisely partitioned by partitioning the variable involved in the nonconvex term within its domain. This is known to provide better lower bounds in comparison to their linear counterparts. This strong bounds help in accelerating the convergence within the lower and upper bounding problem which eventually helps in reducing the solution time of the algorithm. The variable x in the concave univariate term xα is partitioned within its variable domain [xU, xL] into nc, nc = 1, 2, ..., NC partitions or segments. The variable partitioning is denoted by grid points anc. For NC partitions of the variable domain, NC + 1 [anc,anc+1] grid points exist such that a1 = xL and aNC+1 = xU. The grid points can be generated using eq 44.
(47)
⎛ (anc + 1)α − (anc)α ⎞ ⎟ anc + 1 − anc ⎠ nc = 1 ⎝ NC
NC
w≥
∀ nc = 1, 2, ..., NC
∑ (anc)α ·λnc + ∑ ⎜ nc = 1
(dxnc − (anc ·λnc))
(48)
Incorporating eqs 45−48 into the original mathematical model given by eqs 1−42 results in the formulation for the lower bounding problem obtained by piecewise partitioning of the concave terms. This is tighter and provides better lower bounds in comparison to the lower bounding problem obtained by incorporating eq 43 into the original formulation given by eqs 1−42. Similar to the equations for the linear underestimators, the piecewise convex combination equations (eqs 44−48) for T B H various concave univariate terms ((Fpq)dpq,(ΔHBpq)dpq,(ΔHHpq)dpq, C
V
(ΔHCpq)dpq, (ΔHVpq)dpq and (Fn)dn) in our formulation are also written. Alternatively, we could also use the recently developed technique of multiparametric disaggregation38,39 for piecewise underestimation of concave univariate terms. The formulations based on multiparametric disaggregation or radix based discretization are proven to be superior and provide tighter bounds when compared to the scheme shown in eqs 44−48. K
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
Figure 3. Global optimization strategy.
4.2. Global Optimization Algorithm. The discretization based global optimization strategy is used to solve the mathematical model of hydrogen network design under uncertainty to global optimality. This algorithm is similar to the version used by Castro and Teles40 and Kolodziej et al.39 for bilinear programming problems. The only difference is that we have developed a lower bounding problem by piecewise relaxation of the convex underestimator. Castro and Teles40 and Kolodziej et al. 38 instead used a multiparametric disaggregation method for the lower bounding problem. Our proposed global optimization strategy for solving hydrogen design under uncertainty can be used to solve globally any generic nonconvex MINLP, which has its only nonconvexity within the model in the form of concave univariate terms and the rest of the model being linear. The flowchart for the algorithm is given in Figure 3, and the steps for our algorithm are given as follows. In this subsection to avoid using repeatedly using the term hydrogen network model under uncertainty (eq 1−42), and we instead use the term “original MINLP” or “original full scale MINLP” for clarity purpose.
the variables which are involved in the nonconvex term. Loose or improper bounds can affect the solution quality of the lower bounding problem which results in more solution time for the algorithm to converge. 2. The values for the required tolerance criterion ε, number of partitions or level of discretization NC, maximum computer processing time CPUmax, overall lower bound (OLB), and overall upper bound (OUB) are assigned. Initially, the overall lower bound OLB is set to −∞. For the overall upper bound OUB, the original MINLP problem is solved by means of a local MINLP solver to see if any value for the same can be obtained. In other case, OUB is set to +∞. 3. The lower bounding problem is formulated as outlined in section 4.1. This is solved by adding eqs 44−48 into the original model given by eqs 1−42. This convexified MILP model gives a LB on the global optima. If the value of the LB is greater than that of the OLB, then the solution from this step is set as the new OLB. In the other case, the already existing value prevails as the OLB. 4. To obtain the upper bound UB, from the solution of the lower bounding problem (step 3) the binary design variables of the original model are fixed and the
1. The bounds on all the variables in the model are assigned appropriately, especially for the complicating variables or L
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
Table 1. CAPEX and OPEX Data for the Case Studies equipment
CAPEX parameters aBpq(k$) = 8.4, bBpq(k$/kW) = 3.1, dBpq = 0.6 aVpq(k$) = 0.01, bVpq($ day/MMscf-psia) = 0.0025, dVpq = 1.0 aHpq(k$) = 53, bHpq(k$/MW) = 69, dHpq = 0.8 aCpq(k$) = 53, bCpq(k$/MW)= 69, dCpq = 0.8 aTpq($) = 6000, bTpq($) = 100, dTpq = 0.6 an (k$) = 503.8, bn(k$/MMscf/day)= 347.4, dn = 1.0 OPEX parameters
Compressor Valve Heater Cooler Pipeline (Transfer line) Purification unit stream/operation
ais ($/MMscf) = 2000 ois ($/MMscf) = 5 oTpqs ($/MMscf) = 0.006 oBpqs ($/kWh) = 0.03 oHpqs ($/kWh) = 0.01 oCpqs ($/kWh) = 0.01 oVpqs ($/MMscf-psia) = 0.0001 γjs ($/MMscf) = 75 βjs ($/MMBtu) = 0.25
Hydrogen Purification unit Pipeline (Transfer line) Compression Heating Cooling Valve Expansion Flare/Incineration Surplus revenue
continuous variables are fed as an initial point to solve the original full scale MINLP. Since the binary variables are fixed, the original problem now becomes a nonconvex NLP. This is solved by using any local or global NLP solver. The solution now becomes the upper bound UB. If the value of the UB is lesser than the OUB, then the solution from this step is set as the new OUB or the incumbent value of OUB prevails. Although unlikely, there is a possibility of encountering infeasibility for the upper bounding problem. In such a case, in an attempt to find a solution and continue further iterations, we try to find a suboptimal solution for this. This is done heuristically for each upper bounding problem, either by altering the optimality tolerance or by solving for a specified amount of time until the best solution could be obtained. 5. The gap now between the OUB and OLB is examined and checked for whether it satisfies the necessary tolerance criterion ε as assigned in step 2. If the necessary tolerance criterion ε is reached, then the algorithm terminates and the global solution is given by OUB. The algorithm can also be terminated when the time for the algorithm exceeds its threshold limit given by CPUmax assigned in step 2. In this case, the incumbent OUB becomes the global solution. 6. If the tolerance criterion is not satisfactory enough, the number of partitions is increased by one, the algorithm proceeds to next iteration, and the steps 3−5 are repeated. In every iteration, the OLB and OUB are updated. 7. If the convex relaxation of the original MINLP is found to be infeasible, then the original MINLP is infeasible. Also, in this algorithm it is unlikely that the integer infeasibility in the lower bounding problem can arise because no information on the bounding problems is transmitted among the iterations and that all the bounding problems are solved from scratch for each iteration. Remarks i. In this study, the variables involved in the concave univariate term are partitioned or discretized with the same number of partitions NC chosen for all the variables. But this does not have to be the case and the
number of partitions NC could be different for different variables. However, as the algorithm proceeds to next iteration, the partitions need to be updated. ii. In the algorithm, the numbers of partitions for the variables are updated in steps of one for each iteration. They could, in principle, be updated in any possible manner. Based on heuristics and problem complexity, the number of partitions can be updated such that the algorithm convergence is achieved quickly. For instance, the partitions can be updated in steps of two, three, or any number to achieve fast convergence for the algorithm.
5. CASE STUDIES AND IMPLEMENTATION The proposed model and solution algorithm were implemented in GAMS version 24.1.2.41 GAMS/CPLEX version 12.5.1.0 was used for solving the MILP problem, whereas GAMS/ CONOPT version 3.15 was used for solving the NLP problems. For the global and local optimization of the direct MINLP problem, GAMS/BARON version 12.3.3 and GAMS/ DICOPT were used, respectively. To solve each of the bounding problems, the relative optimality tolerance was adjusted by using the GAMS option optcr. For case study 1 the optcr was set to 0, whereas for case study 2 the optcr was set to 0.06. All the computations were carried out using Dell Optiplex 780 PC having Windows 7 Professional 32-bit operating system with Intel Core 2 Duo processor each containing 3.0 GHz cores. The uncertainty in the parameters is characterized through the use of normal distribution. A simple sampling rule for generating scenarios, as proposed by Li et al.,18 is used in this work. The probability of each scenario is given by the cumulative distribution function of normal distribution. Using the same nomenclature as used by Li et al.,18 for a normal distribution with mean μ, standard deviation σ and b scenarios; then hth scenario sampled between [μ − 3σ, μ + 3σ] for the uncertain parameter e is given as follows (eq 49). eh = −3σ + μ +
⎛ 3σ ⎞ ⎛ 6σ ⎞ ⎜ ⎟ + ⎜ ⎟(h − 1) ⎝ b ⎠ ⎝ b ⎠
(49)
For both the case studies, Table 1 gives the CAPEX and OPEX data for various equipment in the network, Table 2 gives the M
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
Table 2. Specifications on Feeds to Various Destination Units unit
pressure (psia)
purity (%)
temperature (K)
flow (MMscf/day)
conversion (% per pass)
2000 500 200
90.90 87.60 -
320
92.560 75.240 -
36.13 38.01
300 2000 500 600 250 200 300
78.85 80.61 75.14 77.57 ≥20 -
320 320 300
13.030 180.340 52.080 40.080 ≤100 ≤100 ≤200
26.90 32.80 40.90 41.20
Case Study 1 Unit A Unit B FGS Case Study 2 NHT HCU CNHT DHT FL GT PSA a
90a
Recovery of hydrogen in the PSA.
Table 3. Specifications on Product Streams from Various Origin Units unit
a
pressure (psia)
purity (%)
temperature (K)
HP Unit A Unit B
300 1200 350
99 85 76
300 300 300
SRU HI CRU NHT HCU CNHT DHT PSAa PSAb
300 300 300 200 1200 350 400 300 180
93 95 80 75 75 70 73 99 0−50
300 300 300 300 300 300 300 300 280
capacity (MMscf/day) Case Study 1 ≤200 case study 2 45.090 ≤30 30.050 ≤200 ≤200
sp. heat (kJ/MMscf-K)
JT coefficient (K/psia)
adiabatic index
34335 35940 36874
0.000 70 0.000 83 0.001 56
0.2943 0.3007 0.3016
34335 33592 36252 37641 37185 36838 36910 34335 38779
0.000 70 0.000 70 0.004 70 0.009 00 0.003 10 0.007 10 0.005 10 0.000 70 0.007 70
0.2943 0.2943 0.2701 0.2739 0.2967 0.2785 0.2805 0.2943 0.2706
Main hydrogen product. bResidue stream.
Figure 4. Optimal hydrogen network for case study 1.
5.1. Case study 1. For this we consider a hydrogen network optimization synthesis problem which has one hydrogen source (HP) in the form of hydrogen import, two processing units (Unit A and Unit B), and one fuel gas sink (FGS). The fuel gas sink FGS is a boiler which produces surplus electricity, and this reduces the operational cost of the refinery. The inlet hydrogen flow and per pass conversions of Unit A are considered to be uncertain. All the other parameters
specifications on feeds to various destination units, and Table 3 gives the specification of product streams from various origin units. The mathematical model and solution strategy developed in this paper are generic and can be applied to any hydrogen network system. In our study, we chose a small (Case Study 1) and large (Case study 2) case study to demonstrate the applicability of the proposed model and solution strategy. N
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
feasible solution at this stage with minimum computational time, we chose GAMS/DICOPT (MILP − CPLEX, NLP − SNOPT). This is done despite the fact that the computational time of the MINLP is not included as a part of the computational time of the algorithm. After this, the MILP convex piecewise relaxation of the model was solved and the lower bound was obtained as stated in step 3 of the algorithm. This value was 53 863.25 k$/yr. Then the binary first stage design variables were fixed from the above obtained solution and the resulting nonconvex NLP was solved. This solution happens to be 53 890.87 k$/yr. This solution, being better than the solution obtained from DICOPT previously, now becomes the new overall upper bound OUB. Since the overall lower and upper bounds are now within a tolerance of 1% (ε = 0.01), the algorithm terminates and the 53 890.87 k$/yr is the global solution of the problem. The gap between the overall lower and upper bound is now 0.05%, and it can be further reduced to 0.005% by proceeding with the algorithm for more iteration (further partitioning of variable domain), although this can sometimes be computationally intensive. This is similar to the spatial branch and bound where the gap between the lower and upper bounds is reduced further by going down the branch and bound tree set up by branching on the continuous variable. Table 4 shows that the lower bound obtained using piecewise relaxation increases, whereas the gap decreases by increasing the number of partitions NC on the variable domain.
are assumed to have their same values in all the scenarios. The uncertain parameters for the inlet hydrogen flow into Unit A are normally distributed with their mean values shown in Table 2 and their standard deviations being 5 MMscf/day. Similarly even the per-pass conversion of Unit A is normally distributed with mean value given in Table 2 and standard deviation being 0.01. The scenarios are distributed normally following the sampling rule in eq 49 for nine (b = 9) scenarios. All the uncertain parameters are assumed to be independent, and 9 scenarios for each uncertain parameter are expressed. This amounts to having 9 scenarios for 2 independent parameters leading to a total of 81 (92) scenarios assuming all the uncertain parameters to be independent. The minimum (maximum) limit on the temperature attained in a transfer line is 250 K (1000 K). The minimum flow for which a transfer line (pipeline) can exist is set to 0.01 MMscf/day and minimum temperature drop for existence of heater and cooler is 5 K. The minimum duty for existence of compressor in a transfer line is 50 kW and the efficiency of the compression process is assumed to be 75%. The global solution obtained from the proposed algorithm yields a TAC of 53,890.87 k$/yr with first stage annualized capital cost (CAPEX) being 135.18 k$ and second stage expected operating cost (OPEX) being 53755.74 k$/yr. The solution was reached relatively soon with a time of 1.1 CPU seconds and the optimal network for this is given in Figure 4. The problem with 81 scenarios has 4581 continuous variables, 7728 constraints, and 45 discrete variables. We also solved this model with other local solvers to see the quality of the solution obtained. A small perturbation term ϵ (0.0001) was added to the nonlinear term xα for solvers which had convergence issues when the variable x takes the value zero. While GAMS/ DICOPT (CPLEX as MILP solver and CONOPT as NLP solver) and GAMS/SBB returned an infeasible solution within 3.7 and 3.8 CPU seconds, respectively, GAMS/COUENNE and GAMS/BONMIN could not even obtain any solution after 1000 CPU seconds of computation. GAMS/DICOPT (CPLEX as MILP solver and SNOPT as NLP solver) returned a local solution of 54 017.39 k$/yr in 11 CPU seconds. GAMS/ DICOPT (CPLEX as MILP solver and IPOPT as NLP solver), however, provided a solution of 53 890.87 k$/yr in 48 CPU seconds which happened to be the global solution of the problem. This further demonstrates the need to solve these stochastic network synthesis problems globally. GAMS/ BARON was also used to solve the original model to global optimality. BARON also found the same solution of 53 890.87 k$/yr as the global solution but required 5792 CPU seconds. Despite the original model being predominantly linear, we observe that considerable amount of time is required to obtain the global solution to the problem. This also highlights the need for developing specific solution strategies for stochastic network synthesis optimization problems. For the purpose of application of the algorithm, the initial value for number of partitions is set to 1 (NC = 1). The required tolerance criterion is set to 0.01 (ε = 0.01), maximum computational resource usage restricted to 5 h (CPUmax = 5CPU hrs), and OLB is set to −∞. Also uniform partitioning of the variable domain is carried out by setting γ = 1. An overall upper bound OUB found by solving the original model using GAMS/DICOPT (MILP − CPLEX, NLP − SNOPT) yielded a solution of 54 017.39 k$/yr. Although GAMS/DICOPT (MILP − CPLEX, NLP − IPOPT) could also be used in this stage, we preferred not to use the same as it required relatively large computational time. Since our aim was to find a good
Table 4. Progress of Iterations for the Global Optimization Strategy for Case Study 1 number of partitions
overall lower bound
overall upper bound
iteration
NC
OLB
OUB
gap (%)
1 2 3 4 5 6 7
1 2 3 4 5 6 7
53863.25 53878.77 53882.48 53885.70 53888.37 53884.92 53886.78
53890.87 53890.87 53890.87 53890.87 53890.87 53890.87 53890.87
0.051 0.022 0.016 0.010 0.005 -
5.2. Case Study 2. In this case we consider a large system which could exist in a typical refinery. It has three hydrogen sources, four processing units, two fuel gas sinks, and one purification unit. It is necessary to design a hydrogen network for the above-mentioned refinery system at the grassroot level. This case is taken from Alves and Towler,42 however, the units for flows are converted to MMscf/day. Also there was no purification unit considered by Alves and Towler42 for this case study. Hydrogen sources include steam reforming unit (SRU), hydrogen import (HI) from outside, and continuous catalytic reforming unit (CCR). The processing units are the naphtha hydrotreater (NHT), hydrocracking unit (HCU), cracked naphtha hydrotreater (GOHT), and diesel hydrotreater (DHT). The purification unit available for this purpose is a pressure swing adsorption unit (PSA). Two types of fuel gas sinks are present, namely, flare (FL) and gas turbine (GT). Flare (FL) does not have any flow or purity requirement, whereas gas turbine (GT) has a minimum purity requirement of 20% (mole fraction) of hydrogen gas. FL may have disposal charges associated with it whereas GT does not incur any disposal/incineration charges. No surplus revenue is considered O
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
Figure 5. Optimal hydrogen network for case study 2.
lower and upper bounds were within the required tolerance criterion of 1%, the algorithm terminates and 14 451.68 k$/yr is the global solution of the problem. This required 581 CPU seconds. The optimal network has a total cost of 14 451.68 k $/yr with the first stage annualized capital cost (CAPEX) being 693.756 k$ and second stage expected operating cost (OPEX) being 13 757.94 k$/yr. All the solvers namely GAMS/DICOPT (CPLEX as MILP solver and SNOPT as NLP solver), GAMS/ SBB, GAMS/COUENNE, and GAMS/BONMIN could not handle the large problem size and returned infeasible solutions. Even GAMS/DICOPT (CPLEX as MILP solver and IPOPT as NLP solver), which provided the global optimum for case study 1, struggled to provide a feasible solution for this case. GAMS/ BARON could not provide any solution and continued its computation even after 18 000 CPU seconds. The deterministic equivalent of the stochastic hydrogen network design optimization is a large size problem, and its size builds with increasing number of scenarios. The generalized formula for calculating the number of continuous variables for the mathematical model described in eqs 1−42 is given as (5·(| P|·|Q|·|S| − |N|·|S|)) + ((5·(|P|·|Q| − |N|)) + |N|) + (|S|·(|I| + |N| + |J|)) + (|S|·|Q|·|P|) + (|S|·|N|((2·|J|) − 1)), where |P|, |Q|, |I|, | N|, |J|, and |S| are cardinality of the sets containing p origin units, q destination units, i sources, n purification units, j fuel gas sinks, and s scenarios. The symbol “·” is the product operator. The number of discrete or binary variables is given by ((5·(|P|·|Q| − |N|)) + |N|). For these large scale problems,
by burning excess gas in GT. In this case study, we consider only the cost of hydrogen associated with HI. Hydrogen produced from SRU is not considered for cost since it is an inhouse hydrogen producer in the refinery, albeit it may involve feed (natural gas) and operational expenses. Hydrogen from CCR is a byproduct and is not of very high purity, so its cost is considered zero. All the other parameters are chosen similar to that used in case study 1. Uncertainty is considered in the inlet flow requirement and per pass conversion of HCU. Uncertainty is also seen in the inlet flow requirement of CNHT. Five scenarios for 3 uncertain parameters are considered, and this accounts for a total of 125 (53) scenarios assuming all the uncertain parameters to be independent. The uncertain parameters are normally distributed with their mean values shown in Table 2. Standard deviation of 5 MMscf/day was considered for inlet flow requirements into HCU and CNHT. For the per pass conversion in HCU, the standard deviation considered was 0.01. The model for this case study has 43 026 continuous variables, 64 375 equations, and 276 discrete variables. The number of partitions NC, the required tolerance criterion ε, maximum computational resource usage CPUmax, and overall lower bound OLB were all set to their values as in case study 1. The initial value of the overall upper bound was set to +∞. The proposed algorithm yielded a lower bound of 14 403.89 k$/yr and an upper bound of 14 451.68 k$/yr in the first iteration (NC = 1) with a gap of 0.3%. Since the overall P
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
information regarding the uncertainty is known, the original deterministic models with their corresponding probabilities are solved for all the scenarios. The solution to this is called expected value using a perfect prediction (EVPP), and this value is 14 308.82 k$/yr. The difference between this and the solution of the stochastic formulation (also called expected value under uncertainty (EVUU)) is called the expected value of perfect information (EVPI) which is 142.86 k$/yr. The relatively low value of EVPI shows that the stochastic model minimizes the effect of uncertainty of in operating aspects of HCU and CNHT by optimally hedging against different scenarios of operation and thereby produces a minimum TAC.
relying on solvers to directly solve this type of problems may not be a suitable option. The local solvers in GAMS may provide a solution for a small scale problem but may not provide a solution or return infeasible solutions for large scale problems. Global solvers may require a lot of time to solve even small scale problems. Hence from the case studies, we infer that a specific solution strategy similar to the one proposed in this paper can not only handle problems of varying sizes but also ensure global optimum up to the required tolerance criterion. From the case studies, we see that a strategy which applies NC = 1 (or no discretization) was enough to give a lower bound which in turn provided a gap which was well within the tolerance criterion. This is an important conclusion because increasing the partitions can be sometimes computationally expensive given the substantial increase in the binary and continuous variables, albeit it increases the relaxation quality by providing tighter lower bounds. Hence when the nonlinearity is only in the form of concave univariate term for any general network synthesis problem, following a strategy which applies NC = 1 (or no discretization) can usually provide strong lower bounds on the global optimum. The formulation corresponding to NC = 1 can be obtained by incorporating linear constraints in eq 43 into the original MINLP model given by eqs 1−42. The optimal network for case study 2 is shown in Figure 5. Apart from providing the overall network topology, the stochastic programming approach also provides optimal decisions which need to be practiced for all individual scenarios. We shall now illustrate the benefits of optimization under uncertainty. From the solution of our stochastic model, we observe that the maximum capacity or the gas flow handled by the purification unit is 10.829 MMscf/day. The deterministic version of the model solved using the values in Table 2 and Table 3 and this constitutes the nominal or base case for this case study. Solving the problem for nominal or base case shows that the gas flow handled by purification unit is 6.707 MMscf/ day. Due to any uncertainty in hydrogen demand at the processing unit, the gas flow into the PSA exceeds 4.122 MMscf/day; then the network operation will be infeasible or nonoptimal. This difference is called the flexibility or flexible network design. The variation in the operating conditions inlet flow requirement at HCU and CNHT and per pass conversion at HCU can result in fluctuations of gas flow handled by the PSA. This is compensated or hedged by designing the purification unit PSA by an additional capacity of 4.122 MMscf/day. Such flexibility in design cannot be accommodated in single scenario optimization approaches thereby revealing the importance of incorporating uncertainty in network optimization problems. This explanation of flexibility in capacity of the units can be straightforwardly extended to all the other units such as transfer line, compressor, valve, heater, and cooler. Similarly, the nominal case predicts HI to be 14.48 MMscf/day, whereas the stochastic optimization approach provided the HI varying from 4.744 MMscf/day to 25.226 MMscf/day depending on the different scenarios of operation. In order to hedge against this variation of hydrogen requirement, the refinery has to make itself equipped with additional hydrogen import of 10.746 MMscf/day to handle the uncertainties in the operation of HCU and CNHT. This information is very crucial and serves as an important aspect for the capacity and utility planning for a refinery. It may also be interesting to compare the compensatory or hedged decisions with the case when the perfect information for the uncertainty is known a priori. Since the perfect
6. CONCLUSION This paper provides a mathematical two stage stochastic programming approach for the design of hydrogen networks under uncertainty. This model was solved to global optimality using a solution strategy which is based on the piecewise convexification of the concave univariate term, which was the only nonconvex term in the model. The model and the solution strategy were implemented in the optimization package GAMS. Two case studies were solved to show the applicability of the proposed model and solution strategy. The model and solution strategy when applied on the two case study problems resulted in a resilient and flexible network design which was capable of handling uncertainties in operating conditions. It was shown through the two case studies that generally these stochastic programming models are of large size and that relying on solvers available in optimization packages may guarantee obtaining a solution or the solution may not be achieved within tractable computational time. Hence there is a need for developing solution strategies which help in solving these large problems within reasonable computational time. For a presented case study which comprised of 276 discrete variables and 43 026 continuous variables, the proposed solution strategy was able to solve the model to accepted global optimality within 10 min of the computational time. This shows the ability of our proposed model and solution strategy in handling industrial scale problems of larger size and complexity. Also the solution strategy, for the case studies, required solution times which were much faster and superior compared to that of one of the state-of-art global solver GAMS/BARON. From our study, we infer that tremendous potential exists for future research in this area. Although our model is generic enough to handle uncertainties in other parameters, we preferred to restrict ourselves to only a few parameters (like inlet flow requirement and per pass conversion for the processing unit) to prevent increasing the size of the problem. It would be interesting to include uncertainty in many other parameters, so that this can give a realistic network design. The challenge in this, however, would be the tremendous increase in the problem size. Next would be to use advanced stochastic programming techniques such as scenario reduction to handle the increasing size of the problem and sample average approximation to characterize the uncertainty more accurately instead of a simple sampling rule used in this study. Also, the proposed model developed by us is risk neutral. It would be interesting to include the component of risk in the model and develop a robust mathematical model to cater to both the best and the worst possible scenarios. Q
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
■
Article
dBpq
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Phone: +971 2 607 5583. Fax: +971 2 607 5200.
dCpq
Notes
dHpq
The authors declare no competing financial interest.
■
NOTATION
dTpq
Nomenclature
CAPEX LP MILP MINLP NLP OPEX TAC
Capital expenditure Linear program Mixed integer linear program Mixed integer nonlinear program Nonlinear program Operating expenditure Total annualized cost
dVpq eh FUn , FLn
Indices
i m n nc j p q s h
FUpq, FLpq
Source Processing unit Purification unit Number of partitions on the variable domain Fuel gas sink Origin unit Destination unit Scenario Scenario under consideration
FUis FUjs Fms FUns
Parameters
an anc aBpq aCpq aHpq aTpq aVpq AF b bn bBpq bCpq bHpq bTpq bVpq Cps dn
FUpqs
Fixed capital expenditure of purification unit n Grid point for NC partitions of variable domain Fixed capital expenditure of compressor present between origin p and destination q Fixed capital expenditure of cooler present between origin p and destination q Fixed capital expenditure of heater present between origin p and destination q Fixed capital expenditure of a transfer line present between origin p and destination q Fixed capital expenditure of valve present between origin p and destination q Annualization factor Number of scenarios Variable capital expenditure of purification unit n Variable capital expenditure of compressor present between origin p and destination q Variable capital expenditure of cooler present between origin p and destination q Variable capital expenditure of heater present between origin p and destination q Variable capital expenditure of a transfer line present between origin p and destination q Variable capital expenditure of valve present between origin p and destination q Constant specific heat capacity of stream from origin p in scenario s Exponent factor for variable capital expenditure of purification unit n
h LHVH2 ncps ons OBpqs OCpqs OHpqs OTpqs OVpqs OP PINqs POUTps rns TLpqs, OUpqs
R
Exponent factor for variable capital expenditure of compressor present between origin p and destination q Exponent factor for variable capital expenditure of cooler present between origin p and destination q Exponent factor for variable capital expenditure of heater present between origin p and destination q Exponent factor for variable capital expenditure of a transfer line present between origin p and destination q Exponent factor for variable capital expenditure of valve present between origin p and destination q hth scenario sampled for uncertain parameter e Upper and lower limits on the purification unit n capacity Upper and lower limits on the transfer line capacity connecting origin p and destination q Known upper bound on flow from hydrogen source i in scenario s Known upper bound on flow into fuel gas sink j in scenario s Known required feed flow into processing unit m in scenario s Known upper bound on flow into purification unit n in scenario s Upper bound on flow from origin p to destination q in scenario s Scenario under consideration Lower heating value of hydrogen Adiabatic index of gas stream from origin p in scenario s Operating expenditure of purification unit n in scenario s Operating expenditure of compressor present between origin p and destination q in scenario s Operating expenditure of cooler present between origin p and destination q in scenario s Operating expenditure of heater present between origin p and destination q in scenario s Operating expenditure of a transfer line present between origin p and destination q in scenario s Operating expenditure of valve present between origin p and destination q in scenario s Operating period of a refinery Known inlet pressure of destination unit q in scenario s Known exit pressure of origin unit p in scenario s Known recovery of hydrogen in purification unit n for scenario s Lower and upper bounds on temperature attained in the transfer line connecting origin p and destination q in scenario s dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research TINqs TOUTps wu x L, x U xLjs, xUjs xms yms ypns zis α αis βjs γ γjs ηs μ μps πs σ χms (ΔHBpq)U,
(ΔHBpq)L
(ΔHHpq)U, (ΔHHpq)L (ΔHCpq)U, (ΔHCpq)L (ΔHVpq)U, (ΔHVpq)L ΔTCpq ΔTHpq
Article
Fis Fijs Fims Fins Fjs Fm′ms
Known inlet temperature of destination unit q in scenario s Known exit temperature of origin unit p in scenario s Symbol to show the underestimator expression for concave univariate term Lower and upper bounds on a general continuous variable x Lower and upper bounds on purity of gas stream into fuel gas sink j in scenario s Known purity (mole fraction) of hydrogen for inlet stream into processing unit m in scenario s Known purity (mole fraction) of hydrogen for exit stream from processing unit m in scenario s Known purity (mole fraction) of hydrogen in product or hydrogen rich stream from purification unit n in scenario s Known purity (mole fraction) of hydrogen from source i in scenario s Generic value for exponent factor (0 < α < 1) Cost associated with source i in scenario s Cost coefficient for revenue generated from surplus output by fuel gas sink j in scenario s Parameter to show uniform or nonuniform partitioning of variable domain Cost coefficient associated with disposal of gas streams in fuel gas sink j in scenario s Efficiency of the compression process in scenario s Mean Joule-Thompson expansion coefficient of stream from origin p in scenario s Known probability of occurrence for scenario s Standard deviation Known per pass conversion of hydrogen in the processing unit m in scenario s Upper and lower bound on the maximum duty of compressor in transfer line SSpq connecting origin p and destination q Upper and lower bound on the maximum duty of heater in transfer line SS pq connecting origin p and destination q Upper and lower bound on the maximum duty of cooler in transfer line SS pq connecting origin p and destination q Upper and lower bound on the maximum duty of valve in transfer line SS p q connecting origin p and destination q Minimum temperature drop for the existence of cooler in transfer line SSpq connecting origin p and destination q Minimum temperature drop for the existence of heater in transfer line SSpq connecting origin p and destination q
Fmm′s Fnms Fmjs Fmns Fn Fns Fn′ns Fnn′s Fnjs F1njs F2njs Fpq Fpqs HINpqs
TAC x w ΔTBpq ΔHBpqs ΔHHpq ΔHHpqs ΔHCpq ΔHCpqs ΔHVpq ΔHVpqs
Binary Variables
Non-Negative Variables
dxnc
Flow from source i in scenario s Flow from source i to fuel gas sink j in scenario s Flow from source i to processing unit m in scenario s Flow from source i to purification unit n in scenario s Feed into fuel gas sink j in scenario s Flow from other processing unit m′ to processing unit m in scenario s Flow from processing unit m to other processing unit m′ in scenario s Flow from purification unit n to processing unit m in scenario s Flow from processing unit m to fuel gas sink j in scenario s Flow from processing unit m to purification unit n in scenario s Maximum capacity of the purification unit n Feed into purification unit n in scenario s Flow from other purification unit n′ to purification unit n in scenario s Flow from purification unit n to other purification unit n′ in scenario s Flow of the residue stream from purification unit n to fuel gas sink j in scenario s Hydrogen flow in the residue stream from purification unit n to fuel gas sink j in scenario s Non-hydrogen flow in the residue stream from purification unit n to fuel gas sink j in scenario s Maximum flow from origin p to destination q Flow from origin unit p to destination q in scenario s Variable to represent the product of flow (Fpqs), specific heat (Cps), and temperature (Tpqs) of a gas stream in transfer line SSpq connecting origin p and destination q in scenario s Total annualized cost Continuous variable involved in concave univariate term Auxiliary variable to represent convex underestimator of concave univariate term Maximum duty of the compressor in transfer line SSpq connecting origin p and destination q Product of Fpqs·Cps and temperature change (ΔTpqs) during compression in transfer line SSpq from origin p to destination q in scenario s Maximum duty of the heater in transfer line SSpq connecting origin p and destination q Product of Fpqs·Cps and temperature change (ΔTpqs) during heating in transfer line SSpq from origin p to destination q in scenario s Maximum duty of the cooler in transfer line SSpq connecting origin p and destination q Product of Fpqs·Cps and temperature change (ΔTpqs) during cooling in transfer line SSpq from origin p to destination q in scenario s Maximum duty of the valve in transfer line SSpq connecting origin p and destination q Product of Fpqs·Cps and temperature change (ΔTpqs) during expansion in transfer line SSpq from origin p to destination q in scenario s
XFn
Continuous variable used as a place holder for partitioned portion of the variable involved in concave univariate term
XFpq S
Binary variable to denote existence of purification unit n Binary variable to denote existence of flow from origin p to destination q dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research XHBpq XΔHHpq XΔHCpq XΔHVpq λnc
■
Article
(19) Al-Redhwan, S. A.; Crittenden, B. D.; Lababidi, H. M. S. Wastewater minimization under uncertain operational conditions. Comput. Chem. Eng. 2005, 29 (5), 1009−1021. (20) Ahmad, M. I.; Zhang, N.; Jobson, M. Modelling and optimization for design of hydrogen networks for multi-period operation. J. Clean. Prod. 2010, 18 (9), 889−899. (21) Jiao, Y.; Su, H.; Hou, W.; Li, P. Design and Optimization of Flexible Hydrogen Systems in Refineries. Ind. Eng. Chem. Res. 2013, 52 (11), 4113−4131. (22) Lou, J.; Liao, Z.; Jiang, B.; Wang, J.; Yang, Y. Robust optimization of hydrogen network. Int. J. Hydrogen Energy 2014, 39 (3), 1210−1219. (23) Jagannath, A.; Elkamel, A.; Karimi, I. A. Improved synthesis of hydrogen networks for refineries. Ind. Eng. Chem. Res. 2014, DOI: 10.1021/ie5005042. (24) Jagannath, A.; Hasan, M. M. F.; Al-Fadhli, F. M.; Karimi, I. A.; Allen, D. T. Minimize Flaring through Integration with Fuel Gas Networks. Ind. Eng. Chem. Res. 2012, 51 (39), 12630−12641. (25) Floudas, C. A. Nonlinear and Mixed Integer Optimization: Fundamentals and Applications; Oxford University Press: New York, NY, 1995. (26) Misener, R.; Floudas, C. A. Global optimization of large-scale generalized pooling problems: Quadratically constrained MINLP models. Ind. Eng. Chem. Res. 2010, 49 (11), 5424−5438. (27) Armagan, E. Decomposition Algorithms for Global Solution of Deterministic and Stochastic Pooling Problems in Natural Gas Value Chains. M.Sci. Thesis, Massachusetts Institute of Technology: Cambridge, MA, 2009. (28) Karuppiah, R.; Furman, K. C.; Grossmann, I. E. Global optimization for scheduling refinery crude oil operations. Comput. Chem. Eng. 2008, 32 (11), 2745−2766. (29) Bergamini, M. L.; Grossmann, I. E.; Scenna, N.; Aguirre, P. An improved piecewise outer-approximation algorithm for the global optimization of MINLP models involving concave and bilinear terms. Comput. Chem. Eng. 2008, 32 (3), 477−493. (30) Zamora, J.; Grossmann, I. E. A Branch and Contract Algorithm for Problems with Concave Univariate, Bilinear and Linear Fractional Terms. J. Global Optim. 1999, 14 (3), 217−249. (31) Karuppiah, R.; Grossmann, I. E. Global optimization for the synthesis of integrated water systems in chemical processes. Comput. Chem. Eng. 2006, 30 (4), 650−673. (32) Saif, Y.; Elkamel, A.; Pritzker, M. Global optimization of reverse osmosis network for wastewater treatment and minimization. Ind. Eng. Chem. Res. 2008, 47 (9), 3060−3070. (33) Bergamini, M. L.; Aguirre, P.; Grossmann, I. E. Logic-based outer approximation for globally optimal synthesis of process networks. Comput. Chem. Eng. 2005, 29 (9), 1914−1933. (34) Wicaksono, D. S.; Karimi, I. A. Piecewise MILP under- and overestimators for global optimization of bilinear programs. AIChE J. 2008, 54 (4), 991−1008. (35) Gounaris, C. E.; Misener, R.; Floudas, C. A. Computational comparison of piecewise-linear relaxations for pooling problems. Ind. Eng. Chem. Res. 2009, 48 (12), 5742−5766. (36) Hasan, M. M. F.; Karimi, I. A. Piecewise linear relaxation of bilinear programs using bivariate partitioning. AIChE J. 2010, 56 (7), 1880−1893. (37) Misener, R.; Thompson, J. P.; Floudas, C. A. Apogee: Global optimization of standard, generalized, and extended pooling problems via linear and logarithmic partitioning schemes. Comput. Chem. Eng. 2011, 35 (5), 876−892. (38) Kolodziej, S.; Castro, P.; Grossmann, I. E. Global optimization of bilinear programs with a multiparametric disaggregation technique. J. Global Optim. 2013, 57 (4), 1039−1063. (39) Teles, J. P.; Castro, P. M.; Matos, H. A. Global optimization of water networks design using multiparametric disaggregation. Comput. Chem. Eng. 2012, 40 (0), 132−147. (40) Castro, P. M.; Teles, J. P. Comparison of global optimization algorithms for the design of water-using networks. Comput. Chem. Eng. 2013, 52 (0), 249−261.
Binary variable to denote existence of compressor in transfer line SSpq connecting origin p and destination q Binary variable to denote existence of heater in transfer line SSpq connecting origin p and destination q Binary variable to denote existence of cooler in transfer line SSpq connecting origin p and destination q Binary variable to denote existence of valve in transfer line SSpq connecting origin p and destination q Binary variable to represent the condition that continuous variable x lies within two successive grid points anc, anc+1
REFERENCES
(1) Hallale, N.; Liu, F. Refinery hydrogen management for clean fuels production. Adv. Environ. Res. 2001, 6 (1), 81−98. (2) Zhang, J.; Zhu, X. X.; Towler, G. P. A simultaneous optimization strategy for overall integration in refinery planning. Ind. Eng. Chem. Res. 2001, 40 (12), 2640−2653. (3) Liu, F.; Zhang, N. Strategy of purifier selection and integration in hydrogen networks. Chem. Eng. Res. Des. 2004, 82 (10), 1315−1330. (4) Fonseca, A.; Sa, V.; Bento, H.; Tavares, M. L. C.; Pinto, G.; Gomes, L. A. C. N. Hydrogen distribution network optimization: a refinery case study. J. Cleaner Prod. 2008, 16 (16), 1755−1763. (5) Khajehpour, M.; Farhadi, F.; Pishvaie, M. R. Reduced superstructure solution of MINLP problem in refinery hydrogen management. Int. J. Hydrogen Energy 2009, 34 (22), 9233−9238. (6) Kumar, A.; Gautami, G.; Khanam, S. Hydrogen distribution in the refinery using mathematical modeling. Energy 2010, 35 (9), 3763− 3772. (7) Liao, Z.; Wang, J.; Yang, Y.; Rong, G. Integrating purifiers in refinery hydrogen networks: a retrofit case study. J. Cleaner Prod. 2010, 18 (3), 233−241. (8) Jia, N.; Zhang, N. Multi-component optimization for refinery hydrogen networks. Energy 2011, 36 (8), 4663−4670. (9) Elkamel, A.; Alhajri, I.; Almansoori, A.; Saif, Y. Integration of hydrogen management in refinery planning with rigorous process models and product quality specifications. Int. J. Proc. Syst. Eng. 2011, 1 (3/4), 302−330. (10) Jagannath, A. Modeling and Optimization of Gas Networks in a refinery. M.Eng. Thesis, National University of Singapore, Singapore, 2012. (11) Alhajri, I.; Saif, Y.; Elkamel, A.; Almansoori, A. Overall integration of the management of H2 and CO2 within refinery planning using rigorous process models. Chem. Eng. Commun. 2013, 200 (1), 139−161. (12) Jiao, Y.; Su, H.; Hou, W.; Liao, Z. A multiperiod optimization model for hydrogen system scheduling in refinery. Ind. Eng. Chem. Res. 2012, 51 (17), 6085−6098. (13) Sahinidis, N. V. Optimization under uncertainty: state-of-the-art and opportunities. Comput. Chem. Eng. 2004, 28 (6−7), 971−983. (14) Almansoori, A.; Shah, N. Design and operation of a future hydrogen supply chain: Multi-period model. Int. J. Hydrogen Energy. 2009, 34 (19), 7883−7897. (15) Betancourt-Torcat, A.; Almansoori, A.; Elkamel, A.; RicardezSandoval, L. Stochastic Modeling of the Oil Sands Operations under Greenhouse Gas Emission Restrictions and Water Management. Energy Fuels 2013, 27 (9), 5559−5578. (16) Birge, J. R.; Louveaux, F. V. Introduction to Stochastic Programming; Springer New York: New York, NY, 1997. (17) Karuppiah, R.; Grossmann, I. E. Global optimization of multiscenario mixed integer nonlinear programming models arising in the synthesis of integrated water networks under uncertainty. Comput. Chem. Eng. 2008, 32 (1−2), 145−160. (18) Li, X.; Armagan, E.; Tomasgard, A.; Barton, P. I. Stochastic pooling problem for natural gas production network design and operation under uncertainty. AIChE J. 2011, 57 (8), 2120−2135. T
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Industrial & Engineering Chemistry Research
Article
(41) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R. GAMS: A Users Guide; 2005. (42) Alves, J. J.; Towler, G. P. Analysis of refinery hydrogen distribution systems. Ind. Eng. Chem. Res. 2002, 41 (23), 5759−5769.
■
NOTE ADDED AFTER ASAP PUBLICATION After this paper was published ASAP July 23, 2014, corrections were made to Figure 1a, eq 11b, eq 42, and several places in the text. The corrected version was reposted August 20, 2014.
U
dx.doi.org/10.1021/ie5011004 | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
