Natural Gas to Liquid Transportation Fuels under Uncertainty Using

thesis under feedstock price, product price, and investment cost uncertainty. ... Chemistry Research. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 1...
1 downloads 0 Views 745KB Size
Subscriber access provided by STEPHEN F AUSTIN STATE UNIV

Process Systems Engineering

Natural Gas to Liquid Transportation Fuels under Uncertainty Using Robust Optimization Logan Ryan Matthews, Yannis A Guzman, Onur Onel, Alexander M Niziolek, and Christodoulos A. Floudas Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b01638 • Publication Date (Web): 18 Jul 2018 Downloaded from http://pubs.acs.org on July 22, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Natural Gas to Liquid Transportation Fuels under Uncertainty Using Robust Optimization Logan R. Matthews,∗,†,‡,¶ Yannis A. Guzman,†,‡,¶ Onur Onel,†,‡,¶ Alexander M. Niziolek,†,‡,¶ and Christodoulos A. Floudas‡,¶ †Department of Chemical and Biological Engineering, Princeton University, Princeton, NJ, USA ‡Artie McFerrin Department of Chemical Engineering, Texas A&M University, College Station, TX, USA ¶Texas A&M Energy Institute, Texas A&M University, College Station, TX, USA E-mail: [email protected]

Abstract The current pricing climate for natural gas and liquid transportation fuels adds a significant amount of uncertainty when designing new natural gas to liquid transportation fuel (GTL) refineries. Robust optimization is a useful tool for optimization under uncertainty, and can be specifically applied to the problem of GTL process synthesis under feedstock price, product price, and investment cost uncertainty. Using historical data to define uncertain price parameters according to an assumed uniform distribution, a process synthesis superstructure with an uncertain objective function was created to maximize the profit of a GTL refinery. Recently developed, tight probabilistic bounds were used a priori and a posteriori in an iterative method to provide solutions with known probabilities of constraint violation for three different uncertainty ∗

To whom correspondence should be addressed.

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

sets. The relative impact of price and investment cost uncertainties are discussed, as is the impact of uncertainty on the overall guaranteed profit of a refinery, the refinery topology, and product distributions. The profitability results with probabilistic guarantees provide useful information for potential GTL refinery investors in the current uncertain energy markets.

1

Introduction

It is often the case that parameters in optimization problems have the potential to realize a wide variety of values. A prime example of this involves price parameters, especially those derived from energy markets whose prices vary moment by moment. When such volatility exists in the energy arena, uncertainty becomes a major factor in optimization and addressing it is a major research priority in order to solve important problems for energy and the environment 1 . Process synthesis is an extremely valuable tool for determining the optimum topologies and product distributions in refineries that produce liquid transportation fuels or chemicals, and includes parameters for key feedstock and product prices when minimizing the refinery cost or maximizing the overall profit. Yet, an unfavorable realization of a price parameter can have an extremely detrimental impact on the profitability of potential refineries in the future. Thus, it is desirable to have optimum solutions to these process synthesis problems that are robust to a variety of parameter realizations. The production of liquid transportation fuels from non-petroleum based feedstocks is an important research area as well. Natural gas is an appealing feedstock to partially offset petroleum consumption in the energy sector due to its increased availability and decreased prices in recent years. The Energy Information Administration predicts a steadily increasing supply of natural gas through 2040, to a point where the United States will become a net exporter of natural gas 2 . Multiple studies have occurred to address the potential of natural gas to liquid transportation fuel (GTL) refineries, specifically in the form of gasoline, diesel, and kerosene so as to supply the current infrastructure 3 . Natural gas is appealing not 2

ACS Paragon Plus Environment

Page 2 of 60

Page 3 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

only as a standalone feedstock, but also in combination with other domestic feedstocks such as biomass and coal 4–6 . A review of liquid transportation fuel production from hybrid and single feedstocks is available from Floudas et al. 7 . Of course, natural gas may also be utilized to produce important chemical products as well, including olefins 8 or aromatics 9 , and potential for new discoveries exists through detailed modeling and process intensification of new technologies as well 10 . Given the established research in deterministic process synthesis and global optimization of natural gas-based processes, a logical advancement in this research is the incorporation uncertainty in parameters such as feedstock and product prices, and investment costs. The idea of incorporating uncertainty into the synthesis of a refinery is not a new problem. For many models, it is possible that a slight deviation from the nominal parameter value could cause an optimal solution to become infeasible. The first applications of process design under uncertainty utilized two-stage nonlinear programming with bounds on the possible uncertain parameter values such that all possible values of the parameters are considered 11,12 . Since then, a variety of stochastic based algorithms have been utilized for these problems 13,14 . Specifically, parametric mixed-integer nonlinear programming (MINLP) 15 , twostage stochastic MINLPs 16,17 , hybrid parametric/stochastic MINLPs 18 , stochastic penalty functions 19 , chance-constrained programming 20 , and even multi-objective stochastic optimization 21 have been utilized to address the inherent uncertainty in process design, synthesis, and optimization. Many studies involving scenario based approaches for incorporating uncertainty have been seen for energy system designs. For instance, Cheali et al. 22 utilized an iterative scenariobased approach for determining optimal biorefinery systems when market uncertainties exists. Rizwan et al. 23 also utilized a scenario-based approach to produce a stochastic MINLP for an algae-based refinery. Energy systems formulated as MILP input-output models in the study of Kasivisvanathan et al. 24 were made robust by accounting for up to K scenarios that the plant may encounter. Scenarios were also generated for a stochastic shale gas to

3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

polymers model by He et al. 25 , who consider material balance and monetary uncertainties. Furthermore, a two-stage stochastic model is utilized by Gao and You 26 to minimize the levelized cost of a shale gas supply chain given uncertainty in the ultimate recovery parameter. Finally, Gargalo et al. 27 utilize Monte Carlo simulations and a two-stage stochastic formulation to account for uncertain parameters in a glycerol-based biorefinery supply chain. Rather than generating uncertain scenarios for inclusion in a process synthesis and global optimization problem, which can quickly increase model complexity, robust optimization techniques may be utilized to handle uncertain price parameters. Robust optimization can be traced back to the worst-case uncertainty sets of Soyster 28 , although modern theory has expanded greatly since the work of Ben-Tal and Nemirovski 29 , El Ghaoui and Lebret 30 , and Bertsimas and Sim 31 . The theory, originally derived for linear programming, has been expanded to include mixed-integer linear programs 32,33 , and in general robust optimization may be applied for any constraint in which uncertain parameters participate linearly. Uncertainty is included via uncertainty sets, which contain the parameter realizations that must be imposed on each constraint. Successfully imposing these uncertainty sets on the uncertain constraints leads to deterministic robust counterparts, and eventually to optimal solutions which remain feasible for all parameter realizations considered in the sets 34,35 . This distinguishes robust optimization from the previously mentioned approaches, as there are no stochastic variables in the model itself and no scenarios must be generated. Beyond this, upper bounds on the probability of constraint violation can be provided when the size of the uncertainty set is varied 36,37 . The powerful developments in probabilistic bounds both a priori and a posteriori by Guzman et al. 38–40 have led to drastic reductions in conservatism in robust solutions, which will be apparent in the GTL refinery. The purpose of this paper is to investigate the impact of feedstock price, product price, and investment cost uncertainty on the profitability, topology, and product distribution of a GTL refinery when robust optimization is applied to a large-scale MINLP GTL superstructure 3 . Process synthesis of a GTL refinery involves rigorous input-output modeling

4

ACS Paragon Plus Environment

Page 4 of 60

Page 5 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

of all potential unit operations, simultaneous heat, power, and water integration, technoeconomic analysis of each feasible solution, and global optimization through a branch-andbound framework to provide a mathematical guarantee of optimality 4,41 . This powerful methodology has been applied to many different feedstock combinations that involve natural gas with both liquid fuels and chemicals as products 5,6,42,43 , and successful incorporation of uncertainty through robust optimization will provide new insights into the feasibility of these refineries in the current economy. The paper is structured as follows. In the next section, background on robust optimization and probabilistic bounds will be provided for important context regarding the GTL mathematical model and optimal solutions. A brief overview of the GTL process superstructure will then be given in Section 3. The transformation of the objective function into its robust counterpart will be discussed in Section 4, along with the defining of uncertain parameters using historical data. Finally, case studies using the box, interval + ellipsoidal, and interval + polyhedral uncertainty sets will be provided in Section 5 to demonstrate the impact of uncertainty on the profitability, topologies, and product distributions of a GTL refinery. Conclusions are provided in Section 6.

2

Robust Optimization and Probabilistic Bounds

Much of the theoretical developments in robust optimization apply to cases where the uncertain parameters of a model can be sequestered into linear or mixed-integer linear inequality constraints. Often in practice, this is easily achievable via simple substitutions or reformulations. The uncertain linear or mixed-integer linear inequality constraint i can thus be generically represented as: X j ∈J / i

aij xj +

X

a ˜ij xj ≤ bi

(1)

j∈Ji

where variables xj can be continuous, binary, or even fixed to accommodate non-coefficient uncertain parameters. The set Ji contains those indices j for which the corresponding pa5

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 60

rameters aij in constraint i are subject to uncertainty. Thus, those aij which are uncertain are denoted with a tilde. Central to robust optimization is the concept that instead of fixing the vector of uncertain parameters of a model to a specific point, multiple possible vectors of parameter realizations will be deterministically imposed upon the model. Sets of parameter realization vectors are called uncertainty sets, and the uncertainty set which is imposed upon constraint i will be denoted as Ξi . A feasible solution to a model with fixed parameters is guaranteed feasible given that all parameters realize their fixed values. Alternately, imposing those realization vectors in Ξi on constraint i then guarantees that constraint i will remain feasible in the case that the true or future realization vector ξi is in Ξi . First, each uncertain parameter a ˜ij of constraint i will be rewritten as a function of random variable ξij :

a ˜ij = aij + ξij a ˆij

(2)

where constant aij is the nominal value of a ˜ij (typically the expected value of a ˜ij ), and a ˆij is a positive constant. Thus, ξij captures random perturbations away from aij which define realizations of a ˜ij . Given bounded uncertainty, that is, when the realizations of a ˜ij can be assumed to fall between lower and upper bounds, the value of a ˆij should be chosen such that ξij ∈ [−1, 1]. The vector ξi can be defined as containing elements ξij , j ∈ Ji , and represents all perturbations away from the nominal parameter realization vector of constraint i. Then, the uncertainty set Ξi can be defined as including realizations of ξi which meet some criteria. Prominent methods in the literature utilize norm-based criteria for including realizations of ξi in Ξi . Using the ∞-norm distance from ξi = 0, uncertainty set Ξi can be defined as such: Ξ∞ i

 =

 ξi : kξi k∞ = max|ξij | ≤ Ψi , j∈Ji

(3)

where Ψi controls the size of Ξ∞ i by restricting the inclusion of realizations of ξi . The ∞ 34 geometry of Ξ∞ i is a hypercube, and Ξi is called a box uncertainty set . Given that ξij

6

ACS Paragon Plus Environment

Page 7 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

are bounded and supported in [−1, 1], selecting Ψi = 1 will include all possible ξi vectors ∞ (and no spurious ones) in Ξ∞ i ; in this case, Ξi is referred to as the worst-case or interval

uncertainty set 28 . With bounded uncertainty, utilizing other norm-based criteria, such as the 1-norm or 2-norm, can include realizations of ξi which have zero probability of occurrence due to the uncertainty set geometries. To prevent this, sets built from other norm-based criteria are intersected with the interval uncertainty set to exclude spurious realizations. Utilizing the 1-norm and intersecting with the interval uncertainty set: )

( Ξ1∩∞ = i

ξi : kξi k1 =

X

|ξij | ≤ Γi , kξi k∞ ≤ 1

(4)

j∈Ji

where ξij are bounded and supported in [−1, 1]. Size parameter Γi controls the inclusion of is thus called the realizations of ξi via the 1-norm, which is polyhedral in geometry; Ξ1∩∞ i interval + polyhedral uncertainty set 31 . Selecting Γi = |Ji | will include all possible ξi vectors (and no spurious ones) in Ξ1∩∞ . Utilizing the 2-norm and intersecting with the interval i uncertainty set: Ξ2∩∞ = i

 

ξi : kξi k2 =



sX

ξij2 ≤ Ωi , kξi k∞

j∈Ji

  ≤1 

(5)

where ξij are bounded and supported in [−1, 1]. Size parameter Ωi controls the inclusion is thus called of realizations of ξi via the 2-norm, which is ellipsoidal in geometry; Ξ2∩∞ i p the interval + ellipsoidal uncertainty set 44 . Selecting Ωi = |Ji | will include all possible ξi vectors (and no spurious ones) in Ξ2∩∞ . i These three uncertainty sets will be utilized for investigating price and investment cost uncertainty, as bounded uncertain parameters will be assumed. For completeness, it is noted that uncertainty sets may also be created for constraints with only unbounded parameters using ellipsoidal and polyhedral geometries. In this case, the sets would be defined using each geometry’s respective norm, without intersecting with the interval set. Furthermore, new uncertainty sets were recently developed that allow both bounded and unbounded random 7

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 60

variables to exist in the same constraint 35 . To impose a given uncertainty set Ξi onto constraint i, we first rewrite constraint (1) in terms of ξi : X

aij xj +

j

X

ξij a ˆij xj ≤ bi ,

(6)

j∈Ji

and then impose Ξi via the maximization operator: ( X

aij xj + max ξi ∈Ξi

j

) X

ξij a ˆij xj

≤ bi .

(7)

j∈Ji

The inner maximization can be solved in its context of an inequality constraint. This yields a deterministic replacement of constraint (1), called a robust counterpart, which imposes Ξi onto constraint i and can substituted into the original model. For box uncertainty sets, the robust counterpart of constraint (1) is 28,34 : X

aij xj + Ψi

j

X

a ˆij |xj | ≤ bi .

(8)

j∈Ji

For interval + polyhedral uncertainty sets, the robust counterpart of constraint (1) is 31 : X

aij xj +

j

X

pij + Γi zi ≤ bi

j∈Ji

zi + pij ≥ a ˆij |xj | ∀j ∈ Ji

(9)

zi ≥ 0 pij ≥ 0

∀j ∈ Ji .

For interval + ellipsoidal uncertainty sets, the robust counterpart of constraint (1) is 44 : X j

aij xj +

X

a ˆij |xj − zij | + Ωi

j∈Ji

sX

a ˆ2ij zij2 ≤ bi .

(10)

j∈Ji

With bounded uncertainty, setting the values of size parameters Ψi , Γi , or Ωi , which will

8

ACS Paragon Plus Environment

Page 9 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

be generically referred to as ∆i , to their maximum values such that all possible realizations of ξi are included in Ξi will result in zero probability that constraint i will be rendered infeasible by the true or future realization of ξi : (

) X

Pr

aij xj +

j

X

(max)

ξij a ˆij xj > bi : ξij ∈ [−1, 1], ∆i = ∆i

= 0.

(11)

j∈Ji

However, imposing every possible realization onto constraint i will result in very conservative solutions of its model. Less conservative solutions can be obtained by reducing the value of ∆i ; the consequence of this action is that the probability of constraint violation of constraint i becomes nonzero: )

( Pr

X

aij xj +

j

X

ξij a ˆij xj > bi : ξij ∈ [−1, 1], ∆i
0.

(12)

j∈Ji

Say there is an acceptable probability of constraint violation of constraint i, denoted as i . Prior to solving the uncertain model via its robust formulation, size parameter ∆i should be set a priori so that i is guaranteed, regardless of the optimal solution x∗ . The probabilistic expressions B(∆i ) which relate probability of constraint violation to ∆i are called a priori probabilistic bounds; all provide upper bounds on the probability because they must characterize all x∗ : ( Pr

) X j

aij xj +

X

ξij a ˆij xj > bi

≤ B(∆i ).

(13)

j∈Ji

Then the smallest parameter ∆i should be selected such that B(∆i ) ≤ i . A tighter a priori bound B(∆i ) will allow the selection of a smaller value of ∆i which guarantees i and thus yields a less conservative robust solution. An alternative approach is to characterize a posteriori an optimal robust solution x∗ . The probabilistic expressions B(x∗ ) which relate probability of constraint violation to x∗ are called a posteriori probabilistic bounds, although

9

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 60

there exist some which are not bounds but provide an exact probability 40 : ( Pr

) X j

aij x∗j +

X

ξij a ˆij x∗j > bi

≤ B(x∗ ).

(14)

j∈Ji

A tighter a posteriori bound B(x∗ ) will assign a lower probability of constraint violation to constraint i with a given x∗ .

3

Gas to Liquids Process Synthesis Superstructure

The main focus of this paper is the incorporation of uncertainty into a GTL process synthesis superstructure. Thus, it is important to understand the basic processes which are considered through input-output modeling and inclusion in the MINLP model. The following sections will highlight key components in the superstructure which are used for a) natural gas conversion to synthesis gas (syngas), b) syngas cleaning, c) hydrocarbon production through Fischer-Tropsch (FT) or methanol synthesis, d) hydrocarbon upgrading to liquid fuels, e) hydrogen and oxygen production, and f) wastewater treatment. For brevity and emphasis on robust optimization, some details and figures regarding unit operations are removed from these sections, but can be found in previous work by Baliban et al. 3 . A general process overview by section is shown in Figure 1. The full mathematical MINLP model and associated nomenclature may be found in the Supporting Information.

3.1

Natural Gas Conversion

The natural gas in the refinery is assumed to arrive at pipeline conditions. Thus, the assumed molar composition is based on an average of pipeline samples taken by the NETL; it contains 93.1% methane, 3.2% ethane, 0.7% propane, 0.4% butane, 1.0% carbon dioxide, and 1.6% nitrogen 45 . Arriving at 31 bar and 25◦ C, the natural gas is first fed through a zinc oxide sulfur guard bed to protect catalysts that exist downstream 45 .

10

ACS Paragon Plus Environment

Page 11 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Auto Thermal Reforming

Natural Gas Synthesis Gas Production

Hydrogen

Methanol Synthesis

Main Products Utilities

Acid Gas Sulfur

Fischer-Tropsch Synthesis

Methanol to Gasoline

Methanol to Olefins Hydrocarbon Upgrading

Natural Gas

Synthesis Gas Treatment

Oxygen

Hydrocarbon Production

Steam Reforming

Mobil Olefins to Gasoline and Distillate

ZSM-5 Upgrading

Standard Upgrading

LPG-Gasoline Separation

Gasoline

Diesel

Wastewater Treatment

Kerosene

Heat, Power, and Water Integration

LPG Hydrogen and Oxygen Production

Figure 1: Overview of the GTL superstructure by general section, highlighting the main products produced from natural gas.

11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 60

Natural gas can be converted in two primary ways: autothermal reforming (ATR) and steam reforming (SMR). While other conversion methods exist, like direct conversion to methanol or direct conversion to olefins, these routes have been shown to currently lack economic competitiveness, and were not considered viable routes in this model 3 . The ATR and SMR units both create synthesis gas for use in FT or methanol synthesis. The SMR operates at 30 bar and 800◦ C. Reactions 15 and 16 constrain the SMR at equilibrium. H2 O + CO ↔ H2 + CO2

(15)

CH4 + H2 O ↔ CO + 3H2

(16)

The resulting reformed gas is sent to the syngas cleaning section. The SMR operates with only a steam source required, unlike the ATR which needs a pure oxygen stream. Combustion of natural gas or recycled gas provides the heat necessary to drive the endothermic reaction in the SMR, and thus ambient air is also sent to the reactor. Alternatively, the ATR produces synthesis gas through partial combustion of oxygen and steam reforming in the same reactor, and thus a high purity source of oxygen is required to feed oxygen and steam simultaneously. This occurs either through a cryogenic air separation unit that provides a 99.5 wt% oxygen stream, or through electrolysis of water to gain a 100 wt% source of oxygen. The oxygen and natural gas are both fed at 300◦ C, while the steam is sent at 550◦ C. The reactor itself operates at 30 bar and 900◦ C. Just as with the SMR, the reformed gases from the ATR are sent to the syngas cleaning section. While the majority of the input natural gas will be converted through either the ATR or SMR, two other options exist in the process superstructure. A fuel combustor is available to provide heat for processes, or a gas turbine can be used to produce electricity. A summary of the possible natural gas conversion routes is shown in Figure 2.

12

ACS Paragon Plus Environment

Page 13 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Steam

Reformed Gas

Steam Reformer Recycle Natural Gas

Input Natural Gas

Natural Gas for Utility Generation

Sulfur Guard

Steam

Auto Thermal Reformer

Recycle Light Gases

Reformed Gas

Oxygen

Figure 2: Natural gas conversion options. Natural gas is combined with recycled gases and can be converted to synthesis gas either through autothermal reforming or steam reforming. Electricity can also be created through a gas turbine, or heat is provided for the processes through a fuel combustor.

13

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.2

Syngas Cleaning

The reformed gas leaving the ATR or SMR may require a compositional shift for use as synthesis gas, and thus it can be routed to a reverse water-gas shift reactor. This reactor always operates at 26 bar, but the temperature can be varied to 400, 500, or 600◦ C depending on the syngas composition required downstream. CO2 in the plant may also be recycled and consumed in the reverse water-gas shift reactor. After the composition is adjusted, the synthesis gas is cooled to 35◦ C for water knockout; the vapor effluent continues either to a Rectisol unit for CO2 removal or directly to the upgrading section if the CO2 level is acceptable. The Rectisol unit, operating at 25◦ C and 20.1 bar, provides a stream of CO2 which can be compressed to 31 bar for recycle to reformers or the reverse water-gas shift units, or compressed to 150 bar for sequestration. If recycle or sequestration is not necessary, CO2 can also be vented to the atmosphere. At this point, the syngas is prepared for upgrading into hydrocarbons.

3.3

Hydrocarbon Production

Fischer-Tropsch Production of Hydrocarbons Six different Fischer-Tropsch reactors exist in the hydrocarbon production section in order to convert the synthesis gas into raw FT hydrocarbons for upgrading. These six reactors have varying catalysts and temperatures that produce varying levels of hydrocarbons and wax with or without the presence of a watergas shift reaction. An iron-based catalyst exists in four of the FT reactors. Two of these reactors facilitate the reverse water-gas shift reaction, one with minimum wax (320◦ C) and one with nominal amounts of wax (240◦ C) 46,47 . The other two facilitate a forward water-gas shift reactor at 260◦ C, but with different levels of wax production 48,49 . Finally, no water-gas shift reaction occurs in the cobalt-based nominal (240◦ C) and minimal (320◦ C) wax FischerTropsch reactors. After hydrocarbon production is complete, the raw FT hydrocarbons are sent downstream to one of two FT-upgrading areas in the refinery.

14

ACS Paragon Plus Environment

Page 14 of 60

Page 15 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Methanol Synthesis The methanol synthesis route begins with the conversion of synthesis gas in a methanol synthesis reactor at 300◦ C and 50 bar. The conversion of carbon monoxide is approximately 30-40%, and thus the effluent is flashed, with methanol continuing downstream and light gases recycled back to the methanol synthesis reactor. Five percent of the light gas is purged and recycled. The flashed methanol is degassed and is ready to be further converted to liquid fuels through one of two methanol converting processes in the hydrocarbon upgrading section. Because the methanol synthesis process will be heavily involved in the robust results, a flowsheet of the methanol synthesis and upgrading processes is shown in Figure 3. Clean Syngas

Methanol Synthesis

Output Gasoline Offgas

Offgas

Methanol Flash

Methanol to Olefins

Methanol Degasser

Offgas

Wastewater

Olefin Fractionation

Distillate to Distillate Hydrotreater

Olefins to Gasoline and Distillate

Hydrocarbon Fractionation

Raw MTG Product to LPG/Gasoline Separation

Methanol to Gasoline

Kerosene to Distillate Hydrotreater

Raw MOGD Product to LPG/Gasoline Separation

Wastewater

Figure 3: Overview of the methanol synthesis and upgrading portions of the process synthesis superstructure. Reprinted with permission from Matthews et al. 50 (Ind. Eng. Chem. Res., 2016, 55 (12), pp 32033225). Copyright 2015 American Chemical Society.

15

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.4

Hydrocarbon Upgrading

FT Hydrocarbon Upgrading to Fuels If Fischer-Tropsch reactors are utilized in the process, the effluent from these reactors can be upgraded in a variety of ways. There is only one option for the wax upgrading, as it is sent to a wax hydrocracker which produces diesel, kerosene, and C5 -C6 cuts. On the other hand, the hydrocarbon stream can be upgraded in one of two ways. First, a ZSM-5 catalytic unit exists to upgrade these hydrocarbons directly into gasoline and distillate which are fractionated and separated from the sour water produced in the reaction, which is sent to the wastewater treatment facility 48,49 . The resulting gasoline is sent to the LPG-Gasoline separation section while the distillate is hydrotreated to the final diesel product. On the other hand, standard upgrading can be utilized based on a Bechtel design; first, water and oxygenates need to be removed from the raw FT product via a sequence of separating units. 51,52 The water-lean hydrocarbons can then be fractionated and sent through a series of units to produce the final products. Naphtha, kerosene, and distillate from the fractionation unit are hydrotreated, while the light C3 -C5 cut is alkylated. Light gases and offgas produced in these processes are sent to a saturated gas plant to either produce LPG or be recycled back in the refinery.

Methanol Upgrading to Fuels Once methanol has been produced and purified via a methanol synthesis reactor, the methanol can be upgraded through one of two reaction pathways. The first utilizes a methanol-to-gasoline (MTG) reactor which creates gasoline and byproduct liquefied petroleum gas (LPG). The reaction occurs using a ZSM-5 catalyst and occurs at 400◦ C

53,54

; the effluent is separated in the LPG-gasoline separation section.

Alternatively, methanol can first be converted to olefins via a methanol-to-olefins reaction that occurs at 482◦ C with a SAPO-34 catalyst 54,55 . Then, the olefins are converted to gasoline, diesel, and kerosene over a ZSM-5 catalyst in an olefins-to-gasoline/distillate reactor at 300◦ C

56,57

. This sequence, abbreviated in short as the MOGD reactors, can operate at

different levels to produce essentially pure gasoline, or gasoline and distillate at a ratio

16

ACS Paragon Plus Environment

Page 16 of 60

Page 17 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

of 0.12:1. The ability to produce a variety of products at differing ratios may be beneficial when considering product price uncertainty. The modeling of the olefins-to-gasoline/distillate reactor assumes the subsequent fractionation and final hydrotreating that would be required to produce the fuels, which are sent to the gasoline, diesel, and kerosene blenders respectively.

LPG-Gasoline Separation The raw products from the methanol to gasoline reactor and ZSM-5 FT hydrocarbon upgrading unit produce mixtures of gasoline and LPG that must be separated in this section. Briefly, the raw product proceeds through a condensation knockout column to concentrate the hydrocarbons and remove light gases before proceeding to a deethanizer column 53 . The light gas from this deethanizer unit proceeds to an absorber column and is removed, while the C3 + hydrocarbons continue forward to a stabilizer column. The C3 /C4 hydrocarbons from this stabilizer can be alkylized in a unit with input butane to produce a LPG/alkylate mixture that is separated. The LPG is sold as a byproduct and the iso-octane product is blended with other gasoline product from the stabilizer column 53 . Some of the effluent from the stabilizer unit is recycled as lean oil for the absorber column; the light gases recovered in this unit are recycled into the refinery, while any higher hydrocarbons recovered are sent back to the deethanizer as reflux. The products exiting this section include the gasoline product, byproduct LPG, and light gases recycled into the refinery.

3.5

Other Supporting Refinery Sections

The GTL refinery contains many supporting sections which will be touched on briefly. First, wastewater treatment is critical and is addressed onsite with a biological digestor and reverse osmosis unit available as options in the superstructure. The wastewater treatment facility also provides steam for process units and cooling water. Hydrogen and oxygen production is also critical for the refinery. Hydrogen can be produced either through a pressure swing adsorption unit or electrolysis of water. If an electrolyzer is used, oxygen is also produced; if not, a cryogenic air separation unit is required in the facility as well for oxygen production.

17

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 60

Finally, a considerable amount of light gas is generated in the refinery. An interior gas loop can recycle light gas to units that will convert the carbon to valuable products; these units include the autothermal and steam reformers. Alternatively, the light gas could be sent to the pressure swing adsorption unit in the internal recycle loop. An external recycle loop will provide a purge of the light gases and can either consume them as fuel in a fuel combustor for heat in the refinery, or create electricity by using the light gases in a gas turbine. Together, the refinery sections avoid unnecessary loss of value from the useful materials, and ensure that the refinery can operate properly.

3.6

Refinery Economics and Objective Function

The total direct costs for each unit in the GTL refinery are estimated from literature sources and take the form of Equation (17),  T DC = (1 + BOP ) · C0 ·

S Sf

sf ,

(17)

where C0 is the base cost, S is the actual capacity, Sf is the base capacity, and the investment scales with scaling factor sf . The values of these parameters for each unit may be found in the Supporting Information. The balance of plant (BOP) percentage is assumed to be 20%. The investment costs are normalized to $2015 using the Chemical Engineering Plant Cost Index. The indirect costs are assumed to be 32% of the TDC, and when summed with the TDC gives the total overnight cost (TOC). These costs are levelized according to Kreutz et al. 58 , with capital charges (CC) calculated for the refinery based on a levelized capital charge rate (LCCR) of 14.38%/year, and the interest during construction factor (IDCF) of 1.0716.

CC = LCCR · IDCF · T OC

(18)

To gain a consistent basis across the refinery, all costs and revenue are normalized based on 18

ACS Paragon Plus Environment

Page 19 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

the amount of product produced. For the investment costs, this takes the form LCCR · IDCF · T OCu . CAP · P rod · LHVP rod

V CostIN = u

(19)

It is assumed that the yearly operating capacity (CAP ) is 330 days/year, with the production capacity P rod converted to an energy basis using the lower heating value LHVP rod . The operating and maintenance (OM) is assumed to be 5% of the TOC; knowing this, the overall levelized unit cost (CostUu ) is CostUu

=

V CostIN u

+

CostOM u

 =

LCCR · IDCF OM + CAP 365

 ·

T OCu . P rod · LHVP rod

(20)

Substituting the expression for T OCu into Equation (20) shows the nonlinear nature of the unit cost functions due to capacity variables Su : CostUu

 =

LCCR · IDCF OM + CAP 365



1.32 · (1 + BOP ) · C0 · · P rod · LHVP rod



Su Sf,u

sf .

(21)

This unit cost is one of many costs that appears in the objective function of minimizing the negative profit (i.e. maximizing the profit) of the GTL refinery. All costs are levelized by the energy content in GJ of the products produced, so that the objective function takes the form: min

X f ∈F

Costf −

X

Costp + CostCO2 +

p∈P

X

CostUu

(22)

u∈UIN V

Here, the set F represents the feedstocks of natural gas, butane, input electricity, and water, while the set P represents the products gasoline, diesel, kerosene, propane/LPG, and output electricity. The Cost variables are defined to always be positive. Thus, Costf represents the overall purchase costs for feedstock f , and Costp represents the revenue for selling product p. The cost of carbon sequestration (CostSeq ) is also included; this may be necessary to ensure that the overall greenhouse gas emissions remain at or below the levels of equivalent petroleum-based processes. Finally, the investment and operating and maintenance costs are 19

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

included through CostUu , which is defined over the set of units with associated investment cost parameters, UIN V . The minimization of the objective function requires global optimization of a large-scale mixed-integer nonlinear problem. A rigorous branch-and-bound algorithm is utilized to accomplish this minimization. Initially, 300 starting points are created at the root node for the branch-and-bound algorithm using solution pool functionality in CPLEX 59 , solving lower bound mixed-integer linear relaxations of the original MINLP. These initial points are used to branch into child nodes from the root node; at each child node, ten MILP solutions points are generated, and then binary variables are fixed in order solve the resulting NLP using CONOPT 3 60 . The NLP solutions serve as upper bounds, and if the NLP solution at a node is less than the existing upper bound, the upper bound is replaced. Nodes are fathomed if the lower bound solution is within a tolerance  of the current upper bound, and branching occurs to the nodes with the lowest current lower bound. The algorithm is allowed to run until all nodes are visited, or until a time period of 100 hours is complete. At this point, an optimum solution is reported as the upper bound, and the optimality gap is known based on the lowest value of the MILP at the remaining nodes. This complete framework for the GTL refinery can provide nominal, globally optimal solutions to the MINLP process synthesis problem. When the global optimization of the refinery is complete, the topology is used to calculate the minimum investment costs for the heat exchanger network. Using the information from the optimum plant topology, the minimum number of heat exchanger matches are calculated 61 , and the final minimum annualized investment costs are then calculated. 62 Altogether, this provides a globally optimum GTL refinery at the desired capacity with a known optimality gap for nominal parameter values. The novel inclusion of uncertainty via robust optimization takes place in the next section.

20

ACS Paragon Plus Environment

Page 20 of 60

Page 21 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

4

GTL Robust Counterpart

4.1

Robust Objective Function Formulation

In a GTL refinery with product and feedstock price uncertainty as well as investment cost uncertainty, the uncertain parameters appear in the objective function. The nominal objective function (22) can be rewritten according to (23), where ζ represents the negative profit which is minimized. min ζ s.t. ζ =

X

Costf −

f ∈F

X

X

Costp + CostCO2 +

p∈P

CostUu

(23)

u∈UIN V

Furthermore, the new constraint in (23) is converted into an inequality for use with robust optimization,

−ζ +

X f ∈F

Costf −

X

X

Costp + CostCO2 +

p∈P

CostUu ≤ 0

(24)

u∈UIN V

Note that at an optimal solution, this is equivalent to (23) as ζ will be minimized against the bound provided by the other terms on the left-hand side of (24). The price uncertainty impacts the variables Costf and Costp for the feedstocks and products. These variables represent the costs due to the flowrates of the input and output streams from the refinery, as reflected in Equations (25) and (26) below.

Costf = Cf · Ff

∀f ∈F

(25)

Costp = Cp · Fp

∀p∈P

(26)

The total feedstock and product species flowrates Ff and Fp , normalized based on the total energy of product produced, are multiplied by the respective feedstock and product prices Cf and Cp . The costs Cf and Cp are uncertain in this superstructure, and can be denoted for 21

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 60

general component j as C˜j . If Equations (25) and (26) were substituted into the objective function, it would cause the now uncertain constraint (24) to take the form of −ζ +

X f ∈F

C˜f · Ff −

X

C˜p · Fp + CostCO2 +

p∈P

X

CostUu ≤ 0

(27)

u∈UIN V

Then, the uncertain prices can be defined as

C˜f = Cf + ξf Cˆf

∀f ∈F

(28)

C˜p = Cp + ξp Cˆp

∀p∈P

(29)

and

according to Section 2. It is appealing, however, to keep the form of (24) where Equations (25) and (26) are not directly substituted into the objective function, but the total feedstock and product costs are included to remain consistent with previous presentations of the superstructure in literature 3 . This can be done with robust optimization through a transformation of the uncertain parameters. The fractional perturbation of the actual prices from their nominal values may be calculated as a ˆj = Cˆj /Cj , and this fraction represents the maximum fractional change in the contribution of variables Costj to the objective function. That is, a new uncertain parameter a ˜j may be defined representing the uncertain coefficient of the feedstock and product pricing costs in the objective function:

a ˜f = 1 + ξf a ˆf

∀f ∈F

(30)

a ˜p = 1 + ξp a ˆp

∀ p ∈ P.

(31)

and

These definitions allow a more intuitive form of the robust objective function to be defined

22

ACS Paragon Plus Environment

Page 23 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

in Equation (32).

−ζ +

X

a ˜f · Costf −

f ∈F

X

X

a ˜p · Costp + CostCO2 +

p∈P

CostUu ≤ 0

(32)

u∈UIN V

Equations (25) and (26) remain in the model with the cost parameters Cj at their nominal values, as defined in the following section. Thus, the variables Costf and Costp still represent the nominal cost and revenue from the feedstocks and products, and the uncertainty is only included due to the parameter a ˜. This formulation is mathematically equivalent to (27). There is considerable potential for uncertainty in the investment costs for a refinery of this scale as well. Thus, uncertain parameters will be defined to represent this uncertainty. The expression for the cost of a unit, Equation (21), is nonlinear, which initially might imply complications for use with robust optimization. However, the assumption of this paper is that the uncertainty will exist in the base costs, C0,u . Since these parameters participate linearly in the expression, the theory of robust optimization will still hold with some simple reformulations. Designating an uncertain base cost C˜0,u = C0,u + ξu Cˆ0,u , it is easy to see that the overall cost of a unit, Costu , is proportional to the base cost, just as Costj is to C˜j . Thus, any fractional change in C˜0,u would lead to a fractional change in the cost of the unit, and thus we can consider the coefficient of the unit costs to be uncertain in the objective function. That is, we can define the coefficients of the unit costs,

a ˜u = 1 + ξu a ˆu ,

(33)

where a ˆu represents the maximum fractional perturbation, i.e.

ˆ0,u C . C0,u

Finally, we can utilize these definitions to gain a final version of the objective function with uncertain parameters specified.

−ζ +

X f ∈F

a ˜f · Costf −

X

a ˜p · Costp + CostCO2 +

p∈P

X u∈UIN V

23

ACS Paragon Plus Environment

a ˜u CostUu ≤ 0

(34)

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 60

This constraint still contains uncertain parameters which are defined using random variables, and thus (34) is not fully deterministic at this point. The robust counterparts outlined in Section 2 will be applied to create tractable forms of the objective function. The box, interval + ellipsoidal, and interval + polyhedral counterparts will be utilized in the following case studies. In each of the formulations, the set J will represent the total set of eight uncertain price parameters participating in the constraint. The price of water is not considered uncertain in these trials, and thus J contains all products and feedstocks except water. A total of 45 units exist in the set UIN V , so there are 45 uncertain investment cost parameters as well. The box counterpart is presented first in (35), with parameter Ψ determining the size of the box: −ζ +

X

Costf −

f ∈F

X

X

Costp + CostCO2 +

p∈P

 X  CostU a ˆj Costj + u +Ψ

u∈UIN V

j∈J

 X

a ˆu Costu  ≤ 0 (35)

u∈UIN V

It should be noted that due to definition of the uncertain parameters a ˜, the nominal parameter values are always one and thus only a ˆ appears in the robust counterparts. This straightforward box counterpart remains linear and is a one-for-one exchange with the original constraint. That is, the original constraint (24) is represented robustly by replacing it with new constraint, (35). The interval + ellipsoidal counterpart requires new variables in the model, and takes the form of constraints (36) and (37). The size of the uncertainty set is determined using Ω. −ζ+

X

f ∈F

Costf −

X

p∈P

Costp +CostCO2 +

X

u∈UIN V

X

CostU u+

a ˆj uj +

j∈J

X

sX a ˆ2j zj2 + a ˆj uj +Ω

u∈UIN V

j∈J

X

a ˆ2u zu2 ≤ 0

u∈UIN V

(36)

−uj ≤ Costj − zj ≤ uj

∀j∈J (37)

−uu ≤ Costu − zu ≤ uu

∀ u ∈ UIN V

Variables uj and zj are auxiliary variables required by the counterpart and represent the impact on variables due to uncertainty by the interval (uj ) and ellipsoidal (zj ) portions of the counterpart. The complexity of the counterpart is immediately evident in the nonlinearities 24

ACS Paragon Plus Environment

Page 25 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

due to the ellipsoidal contribution; as MILP relaxations are required for lower bound solves in the branch-and-bound algorithm, these nonlinearities appear to pose a critical concern. However, the constraint is a second-order conic constraint, and with careful formulation the lower-bound problem can be solved in CPLEX as an MIQCP. 59 Thus, the branch-and-bound algorithm is only slightly influenced. Finally, the interval + polyhedral robust counterpart is included for comparison through reformulation in (38) and (39).

−ζ +

X f ∈F

Costf −

X p∈P

Costp + CostCO2 +

X

CostUu +

X

u∈UIN V

z + wj ≥ a ˆj Costj

wj +

j∈J

X

wu + Γz ≤ 0 (38)

u∈UIN V

∀j∈J (39)

z + wu ≥ a ˆu Costu

∀ u ∈ UIN V

Like the box counterpart, the robust formulation remains linear. However, it does require auxiliary variables z and w, which respectively define the polyhedral and interval contributions to the constraint.

4.2

Price Parameter Definition

The variability in natural gas and liquid transportation fuel price from a five-year period between June 2011 to May 2016 can be seen in Figures 4 and 5. This volatility in the prices can be considered using robust optimization, but the parameters must be identified with a nominal and perturbation value to meet the requirements of a robust formulation. While there are many approaches to this assignment of parameters, a simple solution was utilized for this study. The recent advances by Guzman et al. 38 ,39 ,40 provide drastically improved bounds both a priori and a posteriori, especially when the parameters are subject to uniform, normal, or discrete distributions. Thus, a simple assumption will be made that the prices will vary uniformly between upper and lower bounds defined by historical data. The definition of the uncertain parameters a ˜j from Equations (30) and (31) involves 25

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

7 Natural Gas Price ($/TSCF)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 60

6.5 6 5.5 5 4.5 4 3.5 3 2.5 12/11

12/12

12/13 Date

12/14

12/15

Figure 4: Pricing data from June 2011 to May 2016 for natural gas 63 . historical data over the five year period from June 2011 to May 2016. The lowest and highest prices are assumed to be the bounds on the uniform distribution, and the nominal value is defined as the midpoint. Using the definitions from Section 4.1, if Cjmin and Cjmax represent the lowest and highest prices over the five year period for feedstock or product j, then the nominal price Cj =

Cjmax +Cjmin 2

ˆj and Cˆj = Cjmax − Cj . Thus, the fractional perturbation a

used in the condensed robust counterpart is defined as,

a ˆj =

Cˆj . Cj

(40)

Defining a ˆj in this manner, and defining random variable ξj to be bounded with a uniform distribution in [−1, 1], ensures that all price values for the uncertain parameters in the five year period are included in the worst-case uncertainty set. The values for a ˆj can be seen in Table 1. A uniform assumption with these parameter definitions is typically a conservative definition. As seen in Figure 5, the prices for gasoline, diesel, and kerosene have tended

26

ACS Paragon Plus Environment

Page 27 of 60

3.5

Gasoline Diesel Kerosene Propane Butane

3 Price ($/Gal)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

2.5 2 1.5 1 0.5 12/11

12/12

12/13 Date

12/14

12/15

Figure 5: Pricing data from June 2011 to May 2016 for gasoline, diesel, kerosene, LPG, and butane 64,65 . to be higher than or around the midpoint value, yet a uniform approximation gives equal probability that the price will fall at any point between the highest and lowest values. Thus, using probabilistic bounds with assumed uniform distributions will provide a conservative approximation of the final uncertain profits, yet will greatly improve the results relative to using probabilistic bounds without a known distribution function. The current probabilistic bounds assume each random variable is independent from one another, which is certainly an assumption in the case of these prices, as some correlation likely exists between products such as gasoline, diesel, and kerosene. In the very worst case, this assumption of independence may allow some price scenarios to occur that are very unlikely in the energy market. For instance, on a $/MMBtu scale, the worst case price for natural gas as a feedstock (it’s highest price) is higher than the worst case price for propane, a product (thus it’s lowest price). While it is unlikely that this scenario will ever actually occur, the robust optimization implementation can still provide what the optimal solution

27

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 60

Table 1: Characterization of uncertain feedstock and product price parameters for use in a robust counterpart formulation. Cost Parameter Natural Gas Butane Input Electricity Gasoline Diesel Kerosene LPG Output Electricity

Unit

Lowest Price

Highest Price

Nominal Value

a ˆ

$/TSCF $/gal ¢/kWh $/gal $/gal $/gal $/gal ¢/kWh

2.91 0.50 6.38 1.05 1.02 1.02 0.40 6.38

6.20 2.05 7.62 3.20 3.31 3.30 1.65 7.62

4.75 1.28 7.00 2.13 2.16 2.16 1.03 7.00

0.387 0.606 0.089 0.508 0.530 0.528 0.609 0.089

would be for this case, and the use of probabilistic bounds will help to remove the influence of these unlikely price realizations. There is never a scenario where natural gas is more valuable than the main products, gasoline, diesel, and kerosene. The use of the ellipsoidal or polyhedral shaped uncertainty sets will also help to alleviate some of the extreme cases for prices that may be very unlikely in an actual energy market. In future work, recent advances in the incorporation of correlated uncertainty may be considered 66–69 . For instance, the work of Ning and You 68 utilizes principle component analysis and kernel smoothing methods to construct “data-driven” uncertainty sets that take shapes specific to uncertain parameter data rather than using traditional uncertainty set geometries. This is another effective means for addressing conservatism in robust optimization, albeit with more complicated uncertainty sets and robust counterparts. These methods are not focused on in this study due to the desire to use the recent advances in tight probabilistic bounds 38–40 , which efficiently provide high quality solutions via simple uncertainty sets and basic assumptions on parameters. In this study, assuming independence among each uncertain parameter allows conservative approximations on the upper bound of the probability of constraint violation to be gained. It should also be noted from Table 1 that input and output electricity are given distinct uncertain parameters, providing a potential uncertain scenario where electricity must be purchased and sold at different prices. At the worst case, this would mean that electricity is purchased for more than it is sold for. This is a conservative assumption. However, given

28

ACS Paragon Plus Environment

Page 29 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

that electricity is never both purchased and sold simultaneously in a solution, the use of a posteriori bounds will effectively only include one uncertain price for electricity. Finally, assumptions must be made regarding the level of uncertainty in the investment cost parameters. As all investment costs are calculated using parameters from the literature, a reasonable level of accuracy may be assumed. However, in order to gain a full understanding of the impact of investment cost uncertainty relative to price uncertainty, a serious deviation in investment cost is allowed from their nominal values. It is assumed that the investment costs could vary by up to 50% from their nominal value; that is, a ˆu = 0.50, ∀u ∈ UIN V . This level of deviation is at the same order of magnitude as the a ˆj values for many of the price parameters. Multiple case studies will be presented in Section 5. The first will consider the impact of both price uncertainty and investment cost uncertainty, while the second will consider price uncertainty alone. For the case when investment cost uncertainty is not considered, the a ˆu values will all be set to zero, effectively removing their influence from the robust counterparts. Before discussing the case studies, however, a procedure for gaining very high quality robust solutions via probabilistic bounds will be discussed.

4.3

Improved Solutions with Probabilistic Bounds

The goal of these investigations is to acquire very high quality robust solutions to the GTL process synthesis problem, such that topologies may be gained at desired levels of risk as measured by the probability of constraint violation. In this case, the probability of constraint violation corresponds to the probability that the actual realization of the feedstock prices, product prices, or investment costs will lead to a lower profit than reported by the optimization model. As explained in Section 2, probabilities of constraint violation can be defined either a priori and a posteriori. Bound GMF6 from Guzman et al. 39 is used to set the values of Ψ and Γ a priori for the box and interval + polyhedral sets, while bound GMF7’ 39 is used to define Ω for the interval + ellipsoidal uncertainty set. Bound GMF12 from Guz29

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

man et al. 40 , specific to the uniform bound, will be used to characterize the uncertainty a posteriori. Ideally, there will be a low level of risk associated with the optimum topologies for the GTL refinery. To demonstrate how the level of risk impacts the refinery topologies and profits, solutions will be gained with target bounds on the probabilities of constraint violation ranging from 5% to 30% in 5% increments. It is very easy to achieve these increments using a priori bounds, which can be utilized to define the size of the uncertainty set based on these target values. It is known, however, that a posteriori bounds are significantly tighter than their analogous a priori bounds 38,40 ; given that an exact expression for the probability of constraint violation, GMF12, exists for the case when all parameters are subject to a uniform distribution 40 , better solutions and higher profits will be gained by utilizing a posteriori bounds. At the same time, it is much more difficult to gain an a posteriori bound at a specific probability, since these bounds characterize solutions after an optimization problem is solved. Li and Floudas 37 addressed this issue through an iterative method, which takes advantage of the generally monotonic relationship between ∆ and post . Given this relationship, a bisection method or similar algorithm may be applied to iteratively solve the model, starting from a ∆ value gained via a priori bounds, and adjusting ∆ based on post values until it is within the required tolerance of the desired probability of constraint violation, desired . This iterative method is very effective for finding high quality solutions whose probabilities of constraint violation are bounded at specific values a posteriori. However, it may involve many solves of the optimization problem, and thus requires that the model may be solved very quickly. In the case of a GTL refinery, the global optimization of the process synthesis superstructure requires 100 hours of wall time. Furthermore, it is desirable to have solutions at six different a posteriori probabilities of constraint violation from 5% to 30% for comparison purposes. Performing this iterative algorithm six times, with a problem so complex that it requires 100 hours to find globally optimum solutions, is clearly not computationally

30

ACS Paragon Plus Environment

Page 30 of 60

Page 31 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

feasible. Thus, a new iterative approach is required for this case when the model runtime is extensive. This approach takes advantage of the ability to conduct multiple solves of the model in parallel, gaining a wide spread of data points connecting ∆ to post . First, the values of ∆ at probabilities of constraint violation ranging from 5% to 95% will be calculated using a priori probabilistic bounds, either GMF6 and GMF7’ depending on the type of uncertainty set 39 . By using 5% increments, this corresponds to 19 ∆ values and thus 19 different model instances, which may be solved in parallel, gaining distinct solutions to characterize using a posteriori bounds. Combined with the nominal and worst case solutions, this first iteration will thus lead to 21 points relating ∆ and post . The goal is to gain high quality solutions within a certain tolerance or percent error of some desired probabilities of constraint violation, i.e. 5%, 10%, 15%, 20%, 25%, and 30%. For these cases studies, the solutions will be within one percent error: desired  − post < 0.01. desired

(41)

It is quite possible that after this first iteration, some of these 19 initial solutions may be within this percent error of a desired probability between 5% and 30%. However, it is more likely that this does not happen, and a second iteration will be required. In this second iteration, the existing data points from the first iteration are utilized to predict the appropriate values of ∆ at each desired . This will be done via linear interpolation. If ∆1 and ∆2 correspond to post and post values such that post ≤ desired ≤ post 1 2 1 2 , then the value of ∆ in the second iteration may be calculated as

(2)



desired − post 1 = ∆1 + (∆2 − ∆1 ) · post . 2 − post 1

(42)

These ∆(2) values will be calculated for each of the six desired probabilities for which a 1% percent error was not met in the first iteration. 31

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Once these new values have been calculated, the models are solved again, and the method will terminate if the tolerance is met at each of the six desired probabilities; it is often the case that this occurs in the second iteration for the GTL synthesis problem. However, if any of the desired probabilities are not met, then more linear interpolation will occur utilizing the new data gained in the second iteration relating ∆ to post . Once again, the model would be solved at the new approximations ∆(3) , and the process would repeat until a solution is gained with an post value within the desired tolerance. It is clear that utilizing solves in parallel with linear interpolation for ∆ values will drastically decrease the amount of time required to gain high quality solutions compared to previous iterative methods. It is unlikely that the bisection method would converge in just few solves; if the bisection method required even 10 solves, there is a drastic reduction in both the computational resources and the time required to gain high quality solutions. Thus, this method will be implemented for the GTL process synthesis problem under various types of uncertainty. The values of ∆ for the interval + polyhedral uncertainty set and the interval + ellipsoidal uncertainty set are dependent on the number of uncertain parameters that exist in the constraint. Thus, the initial values will be different depending on whether investment costs are considered uncertain in the model. To demonstrate this, Table S2 is included in the Supporting Information and contains the initial ∆ values for the starting point of the algorithm, including the two endpoints when pri = 0 (worst-case) and pri = 1 (nominal case). In total, there are 21 models per uncertainty set which will be solved in the first iteration. Finally, once a solution has been found, it is possible to characterize the solution based on the change in the objective function value due to the increased robustness. In the case of the GTL refinery, the change in profit due to robustness is quantified using the Price of

32

ACS Paragon Plus Environment

Page 32 of 60

Page 33 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Robustness (PoR) 31 and Normalized Price of Robustness (NPoR) 38 defined in (43) and (44). ζN − ζR . PoR = ζN

(43)

ζN − ζR . NPoR = ζN − ζW C

(44)

In these expressions, ζN and ζW C represent the nominal and worst case objective function values, respectively, while ζR represents the robust value at a specific solution. With both the methodology of gaining solutions and the quantification of solution quality defined, it is now possible to discuss the computational studies in the following section.

5

Computational Studies

The computational studies in this section demonstrate the impact of uncertainty on the process synthesis of GTL refineries that produce gasoline, diesel, and kerosene at a capacity equivalent to 50 thousand barrels per day (kBPD) of gasoline by energy content. The relative impact of price uncertainty compared to investment cost uncertainty will be seen through cases with and without investment cost uncertainty. It will become obvious through these cases studies that the volatility in prices will have a larger impact than the investment cost uncertainty. In both cases, optimal profits will be provided for a priori probabilities of constraint violation ranging from 5% to 95% in step sizes of 5% for each of the three main uncertainty sets. This will provide a clear demonstration of the impact of uncertainty on the optimal solutions. For the case when only price uncertainty is considered, the iterative algorithm will be implemented to gain specific solutions at low probabilities of constraint violation. The studies were performed on the Ada Supercomputing cluster at Texas A&M University. The GTL superstructure with a box uncertainty set is an MINLP with 14,891 variables and 15,116 constraints. The interval + ellipsoidal robust counterpart adds 107 variables and

33

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

107 constraints to the formulation, and the interval + polyhedral robust counterpart adds 54 variables and 53 constraints to the formulation. In each of the implementations, 28 of the variables are binary. There also exists 569 bilinear terms, one trilinear term, and six quadrilinear terms in each formulation. Thus, the branch-and-bound global optimization algorithm described in Section 3.6 is given 100 hours to execute in order to provide solutions with upper bounds, lower bounds, and an optimality gap. The results discussed in the following sections will focus on the upper bound solutions from each of these runs, which as discussed in Section 3.6, represent the best solution found for the MINLP when binary variables are fixed. The lower bound solutions are not as informative as they represent relaxations, and will not be discussed in detail; for reference, in the runs with price uncertainty alone, the maximum gap between an upper bound solution and its lower bound is $3.05/GJ and the average gap is $2.21/GJ. Furthermore, the following sections will focus on the profits of the refineries. Since the objective function was written as a minimization of the negative profit, the optimum profits are simply calculated by multiplying this optimal objective function value by −1. To begin the analysis of the computational results, the nominal and worst case solutions will first be discussed. These points provide context for the various optimal topologies and the broader discussion of the use of probabilistic bounds in process synthesis.

5.1

Nominal and Worst Case Analysis

The difference between the nominal and worst case scenarios is striking, and heavily emphasizes the need for tight probabilistic bounds when using robust optimization. The nominal profit from the GTL refinery is $4.93/GJ, normalized by the total energy of the products produced. This is a reasonable profit that could spur investment in a GTL facility. In this nominal case, the natural gas is upgraded to reformed gas through an autothermal reformer. A methanol intermediate is produced from the reformed gas via a methanol synthesis reactor, and the primary product gasoline is produced along with byproduct LPG via a methanol 34

ACS Paragon Plus Environment

Page 34 of 60

Page 35 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

to gasoline reactor. An autothermal reformer is actually the optimal natural gas upgrading process regardless of the level of robustness, but the products produced and reactors required may change based on the probability of constraint violation, as seen in future sections. Under the assumptions for the uncertain parameters, the nominal profit corresponds to an a posteriori probability of constraint violation of 50%, meaning with this topology the profit would be at least $4.93/GJ approximately 50% of the time. On the other hand, the worst case scenario presents a very dire picture for the refinery. When only feedstock and product price are considered, the worst case value is -$5.51/GJ, a $10.44/GJ decrease in profit, meaning the refinery would operate at a loss if all of the prices reach their worst case values simultaneously. This worst case value is found using the box uncertainty set with Ψ = 1. In this case, gasoline and byproduct LPG are still produced via the methanol to gasoline reactor. If the investment costs are allowed to vary by up to 50% of their nominal value, then the worst case drops even further to -$-7.56/GJ. The same optimal topology occurs in this case as well. These worst case results provide a disappointing outlook for the potential GTL refinery. Of course, they are also vastly too conservative. The likelihood of each uncertain parameter reaching its worst case value simultaneously is exceedingly low. Yet, the change in the objective function values between the nominal and worst case values still clearly demonstrate the need to take uncertainty into account. Fortunately, the following sections will demonstrate how probabilistic bounds will provide solutions which are robust to most uncertain realizations without overly compromising the profit of the GTL refineries.

5.2

Robustness and Profit

The relationship between the overall profit of the refinery and the probability of constraint violation may be seen in Figures 6 through 8 for each of the three main uncertainty sets. In these figures, only feedstock and product price uncertainty are considered, and two data series are included to compare the profit when an a priori probabilistic bound is utilized compared 35

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

to an a posteriori probabilistic bound. These data series contain a horizontal correspondence between points, as each profit value is from a solution which was characterized a priori and a posteriori. It is immediately obvious that regardless of whether a box, interval + ellipsoidal, or interval + polyhedral uncertainty set is utilized, the a posteriori probability GMF12 is significantly tighter than the a priori bound, as expected. This is due to the fact that GMF12 is an exact expression for the probability of constraint violation for a specific solution when all parameters are subject to uniform distributions, whereas only bounds are available a priori that must consider all possible solutions to the optimization problem. The closest the a priori and a posteriori bounds come is in Figure 7 for the interval + ellipsoidal set, for which the very tight a priori bound GMF7’ may be utilized.

4 Overall Profit ($/GJ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 60

2 0 −2 −4 a priori GMF6 a posteriori GMF12

−6 0

0.2 0.4 0.6 0.8 Probability of Constraint Violation

1

Figure 6: Comparison of the optimal profits using a priori and a posteriori probabilities of constraint violation for the box uncertainty set. The improvement in profit which is gained due to the information from probabilistic bounds is immediately evident when comparing to the worst case profit of the previous section. Rather than a worst case profit of -$5.51/GJ, positive profits may occur in a GTL

36

ACS Paragon Plus Environment

Page 37 of 60

Overall Profit ($/GJ)

4 2 0 −2 −4 a priori GMF7’ a posteriori GMF12

−6 0

0.2 0.4 0.6 0.8 Probability of Constraint Violation

1

Figure 7: Comparison of the optimal profits using a priori and a posteriori probabilities of constraint violation for the interval + ellipsoidal uncertainty set.

4 Overall Profit ($/GJ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

2 0 −2 −4 a priori GMF6 a posteriori GMF12

−6 0

0.2 0.4 0.6 0.8 Probability of Constraint Violation

1

Figure 8: Comparison of the optimal profits using a priori and a posteriori probabilities of constraint violation for the interval + polyhedral uncertainty set.

37

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

refinery if small levels of risk are considered. For instance, when the interval + ellipsoidal uncertainty set is utilized, the “breakeven” probability of constraint violation is approximately 11.6%. This is determined by the approximate probability of constraint violation at which the profit of a solution is $0.00/GJ. Because the uncertainty only appears in the objective function, this probabilistic information implies that under the assumptions of these case studies, a topology exists for a GTL refinery which would remain profitable 88% of the time. Similar results may be gained from the box or interval + polyhedral uncertainty sets. The “breakeven” probability of constraint violation would be approximately 15.7% for box uncertainty set and 12.5% for the interval + polyhedral uncertainty set. This indicates that the interval + ellipsoidal uncertainty set tends to give better results for these GTL refineries, especially at low probabilities of constraint violation. The relative performances of each uncertainty set, and the reasons for it, will be discussed in more detail in Section 5.5. Before this, however, the focus will shift towards the inclusion of investment cost uncertainty into these process synthesis problems. Very similar figures to Figures 6 through 8 could be constructed for the case with investment cost uncertainty, demonstrating how a posteriori bounds are superior to a priori bounds in this context as well. Given this is true, it is more informative to simply focus on the relative profits with and without investment cost uncertainty when only a posteriori bounds are utilized. This will be discussed in the following section.

5.3

The Impact of Investment Cost Uncertainty

In the discussion of worst case solutions, it was seen that the inclusion of investment cost uncertainty lowers the overall profit by $2.05/GJ when each uncertain investment cost parameter is allowed to vary by 50% of its nominal value and realizes its worst case value. It follows then that the investment cost uncertainty will likely make a large impact at very low probabilities of constraint violation, but the question remains as to whether it would 38

ACS Paragon Plus Environment

Page 38 of 60

Page 39 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

substantially decrease the profits at reasonably tolerable levels of risk. This question is answered in Figure 9. In this figure the overall profit is plotted against the a posteriori probability of constraint violation when investment uncertainty is and is not included and the box uncertainty set is utilized. All data points were gained by characterizing solutions whose ∆ values were set before each run using a priori bounds as shown in Table S2 of the Supporting Information. The red squares in Figure 9 correspond to the case when investment cost uncertainty is included, and as is expected, these red squares appear much lower on the left side of the figure when the probability of constraint violation is very close to zero. But, under the assumption of uniform distributions for each of the parameters, the difference between the two data series disappears almost immediately once some risk is allowed. Similar figures may be generated with the interval + ellipsoidal and interval + polyhedral uncertainty sets, but are omitted here for brevity. These results clearly indicate that the inclusion of investment cost uncertainty does not significantly impact the results gained from considering price uncertainty alone, if an investor allows even a small amount of uncertainty in their investment via a posteriori probabilistic bounds. It is immediately evident that the most important uncertain parameters are those that deal with the price of natural gas, gasoline, diesel, and kerosene; the results imply that uncertainty in these price parameters have the largest potential to harm the objective function value. Intuitively, this may be seen from the cost and revenue contributions at the nominal solutions. The sum of the Costu variables for all units u ∈ UIN V for the nominal solution is $4.14/GJ. In contrast, the total revenue from gasoline is $13.76/GJ and the total cost of natural gas is $7.15/GJ. Especially considering that each unit has its own uncertain investment cost parameter, it is clear that the price of natural gas and gasoline will impact the overall profit substantially more than one, or even all, of the investment cost parameters. While the presence of investment cost uncertainty is significant in the worst case scenario, it becomes negligible as the sizes of the uncertainty sets shrink and the price uncertainty

39

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

6 4 2 Profit ($/GJ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 60

0 −2 −4 −6

With Investment Uncertainty No Investment Uncertainty

−8 0

0.1 0.2 0.3 0.4 A Posteriori Probability of Constraint Violation

0.5

Figure 9: Comparison of the profit vs a posteriori probability of constraint violation for the box uncertainty set with and without investment cost uncertainty. becomes dominant. Thus, the rest of the results discussed in the following sections will focus on applying the iterative method to the case when only uncertain feedstock and product prices are imposed on the model.

5.4

Improving Results via the Iterative Method

The iterative method described in Section 4.3 is applied to the case when only feedstock and product prices are uncertain, in order to find high quality solutions within a given tolerance of desired a posteriori probabilities of constraint violation. This was implemented for each of the uncertainty sets, with the results provided in Tables 2 and 3 for the interval + ellipsoidal set and interval + polyhedral set, respectively; results for the box set may be found in Tables S3 of the Supporting Information. For the box set, one of the desired a posteriori probabilities (post = 0.30) was actually gained by chance from the first iteration of runs, as were results for post = 0.20 and post = 0.25 with the interval + polyhedral uncertainty set.

40

ACS Paragon Plus Environment

Page 41 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

The remaining probabilities for each uncertainty set were gained within the desired tolerance after two iterations, except for one case with the interval + ellipsoidal set for desired = 0.20 which required four iterations. It is clear then, that this approach will acquire high quality solutions very quickly. The tables each contain the percent error, PoR, and NPoR for each of the iterations. It is clear that the iterative method will greatly improve all of these metrics compared to the use of a priori bounds alone, and these metrics provide valuable insights into the effect of robustness on profit. Consider the case when at most a 15% probability of constraint violation is desired. The interval + polyhedral uncertainty set provides a profit of $0.34/GJ; from the PoR it is known that this is a 93% reduction in profit from the nominal case. This profit is still closer to the nominal profit than the worst case, however, as the NPoR is approximately 44%. That is, using an interval + polyhedral uncertainty set and a 15% probability of constraint violation will move the profit approximately 44% of the way from the nominal to the worst case value. The performance of the interval + polyhedral uncertainty set is better than the box uncertainty set at a 15% probability of constraint violation, but worse than the interval + ellipsoidal set. Interestingly, the NPoR is below 50% when post ≥ 0.15 for all of the sets considered, indicating that the profit will be closer to the nominal than the worst case if only 15% risk is allowed. For the interval + ellipsoidal set, NPoR drops below 50% even for the case when post = 0.10. The tail representing that 10% of risk clearly makes a tremendous difference on the profit, but is unlikely enough that it can reasonably be neglected when making the final investment decisions. At the highest levels of risk considered, 30%, the profits are $2.17/GJ, $2.13/GJ, and $1.88/GJ for the box, interval + ellipsoidal, and interval + polyhedral uncertainty sets respectively. This corresponds to a 56% reduction in profit, or about 26% of the distance from the nominal case to the worst case, but it also indicates that there exists a topology such that for 70% of the possible uncertainty realizations, the profit will be at least $2.17/GJ

41

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

or higher. This provides valuable information for a potential investor in the GTL refinery regarding what profits could occur, and with what likelihood. As will be seen in the next section, even more information may be gathered when reviewing the specific refinery topology that provided this profit, and how it compares to the topologies at the nominal case and other probabilities of constraint violation.

5.5

Robustness and Product Distribution

As it has become obvious through Tables 2, 3, and S3 of the Supporting Information, the choice of uncertainty set leads to very different results regarding the profit of the robust solutions at various probabilities of constraint violation. This is made even more clear through Figure 10, which shows the profit vs probability of constraint violation for each of the sets in the same figure. It is important then to understand the reason why the robust profits differ, and how this should impact the final decision an investor might make regarding a GTL refinery. The primary difference in the robust profits is due to the different product distributions chosen by the various uncertainty sets. Figure 11 provides the product distribution at each of the six a posteriori probabilities gained through the iterative method with the box uncertainty set. It is immediately clear that regardless of the probability chosen, the product distribution will be the same; only gasoline will be produced as the primary product, with some byproduct LPG coming from the methanol to gasoline reactor. The box uncertainty set, scaled by the ∞-norm, does not incentivize diversification of products, and the nominal and worst-case topology is chosen for any probability. The interval + polyhedral uncertainty set does incentivize diversification of products as it is scaled by the 1-norm, as seen in Figure 13. For the interval + polyhedral set, all three major products are produced with byproduct LPG. This would require a completely different plant topology than the nominal solution. A methanol to olefins reactor and an olefins to gasoline/distillate reactor would be included, allowing the production of diesel 42

ACS Paragon Plus Environment

Page 42 of 60

Page 43 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Table 2: Results from the iterative method for the interval + ellipsoidal robust counterpart. The final a posteriori probabilities of constraint violation in the final iteration are within one percent error of the desired probabilities. desired Ω

Profit ($/GJ)

Error post GM F 12 (%)

PoR (%)

NPoR (%)

169.58 152.29 141.17 132.76 124.76 118.44 112.22 106.37 101.33 95.98 91.34 85.26 80.22 74.18 66.59 58.79 50.32 42.01 30.01

80.07 71.91 66.66 62.69 58.91 55.93 52.99 50.23 47.85 45.32 43.13 40.26 37.88 35.03 31.44 27.76 23.76 19.84 14.17

126.17 105.11 90.38 79.73 70.86 56.68

59.58 49.63 42.68 37.65 33.46 26.76

78.85

37.23

79.00

37.30

First Iteration 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95

1.3589 1.2026 1.0975 1.0147 0.9445 0.8823 0.8255 0.7726 0.7223 0.6739 0.6266 0.5798 0.5330 0.4854 0.4364 0.3846 0.3285 0.2647 0.1848

-3.43 -2.58 -2.03 -1.61 -1.22 -0.91 -0.60 -0.31 -0.07 0.20 0.43 0.73 0.97 1.27 1.65 2.03 2.45 2.86 3.45

0.0061 0.0158 0.0269 0.0389 0.0518 0.0655 0.0801 0.0957 0.1121 0.1295 0.1484 0.1688 0.1911 0.2376 0.2641 0.2921 0.3224 0.3568 0.4001

87.80 84.20 82.07 80.55 79.28 78.17 77.11 76.08 75.09 74.10 73.02 71.87 70.60 66.06 64.79 63.49 62.07 60.36 57.88

Second Iteration 0.05 0.10 0.15 0.20 0.25 0.30

0.9543 0.7594 0.6229 0.5239 0.4625 0.3700

-1.29 -0.25 0.47 1.00 1.44 2.13

0.0498 0.0999 0.1499 0.1957 0.2500 0.3000

0.40 0.10 0.07 2.15 0.00 0.00

Third Iteration 0.20

0.5200

1.04

0.1976

1.20

Fourth Iteration 0.20

0.5179

1.03

0.1987

0.65

and kerosene, and the methanol to gasoline reactor would remain to efficiently produce gasoline and LPG. At low probabilities of constraint violation, producing all of the products is clearly superior to producing just gasoline; the interval + polyhedral set has superior

43

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 60

Table 3: Results from the iterative method for the interval + polyhedral robust counterpart. The final a posteriori probabilities of constraint violation in the final iteration are within one percent error of the desired probabilities. PoR (%)

NPoR (%)

212.30 211.20 205.51 202.11 189.89 179.08 169.39 163.06 153.91 145.10 137.57 129.98 122.44 114.85 106.85 98.92 91.23 82.25 70.86 58.59 1.72

100.24 99.73 97.04 95.43 89.67 84.56 79.98 77.00 72.68 68.51 64.96 61.38 57.82 54.23 50.45 46.71 43.08 38.84 33.46 27.67 0.81

1.7561 -1.44 0.0498 0.40 129.13 1.3881 -0.41 0.0998 0.20 108.22 1.1326 0.34 0.1496 0.27 93.03 Found from Γ = 0.9291 in first iteration. Found from Γ = 0.7486 in first iteration. 0.5804 1.88 0.2993 0.23 61.89

60.97 51.10 43.93

desired Γ

Profit ($/GJ)

post GM F 12

Error (%)

First Iteration 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00

8.0000 3.8435 3.4014 3.1042 2.8700 2.6714 2.4955 2.3349 2.1851 2.0429 1.9059 1.7722 1.6400 1.5076 1.3731 1.2342 1.0879 0.9291 0.7486 0.5227 0.0000

-5.53 -5.48 -5.20 -5.03 -4.43 -3.90 -3.42 -3.11 -2.66 -2.22 -1.85 -1.48 -1.11 -0.73 -0.34 0.05 0.43 0.87 1.44 2.04 4.84

0.0000 0.0000 0.0000 0.0004 0.0021 0.0046 0.0078 0.0133 0.0196 0.0275 0.0359 0.0482 0.0630 0.0809 0.1024 0.1282 0.1596 0.1984 0.2481 0.3178 0.5000

0.00 100.00 100.00 99.73 98.95 98.16 97.40 96.20 95.10 93.89 92.82 91.24 89.50 87.55 85.37 82.91 80.05 76.66 72.43 66.55 50.00

Second Iteration 0.05 0.10 0.15 0.20 0.25 0.30

29.22

profits for post ≤ 0.20. Yet, for probabilities of constraint violation of 25% or higher, the box set and nominal topology lead to higher profits. Despite the diversification, the interval + polyhedral uncertainty set does not vary the ratios of products at the various probabilities considered in the iterative method. It would appear that an uncertainty set which provides solutions with different product ratios at different probabilities of constraint violation would

44

ACS Paragon Plus Environment

Page 45 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

perform best. This is exactly what is seen with the interval + ellipsoidal uncertainty set, as shown in Figure 12. At a 5% probability of constraint violation, all of the products are produced, again utilizing both a methanol to gasoline route, and a methanol to olefins and olefins to gasoline/distillate process. Compared to the interval + polyhedral set, however, more gasoline and LPG are produced. As the probability of constraint violation is increased, more gasoline is produced, progressing back towards the nominal distribution of only gasoline and LPG. By a 25% probability of constraint violation, the same refinery topology is chosen with both the interval + ellipsoidal uncertainty set and the box uncertainty set. This change in ratios is what one might expect as more uncertainty is added. Diversification is very useful, but the ratio should not be the same for all levels of risk, and the interval + ellipsoidal uncertainty set confirms this intuition. It is for this reason that the interval + ellipsoidal uncertainty set generally provides the best results across all probabilities of constraint violation. Through the use of the 2norm for scaling the set, it provides the appropriate level of diversity in the products at low probabilities of constraint violation, but also is able to identify the point when the nominal topology becomes the most advantageous. If one was to use only one uncertainty set for future investigations, the interval + ellipsoidal uncertainty set is recommended. However, it is very useful to have the information from both the box and the interval + ellipsoidal sets at low probabilities of constraint violation, due to the fact that the nominal product distribution appears for both sets at relatively low probabilities of constraint violation. Since these product distributions are identical and higher performing than the interval + polyhedral set at a 25% probability of constraint violation, it is known that the nominal solution is the optimal solution for at least 75% of the possible uncertainty realizations. Furthermore, it can be stated that if the nominal topology is built, with a methanol to gasoline reactor, then the profit will be at least $1.44/GJ 75% of the time, and at least $4.93/GJ 50% of the time. Finally, because the nominal topology is the only topology that is gained

45

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

from the box uncertainty set, more guarantees can be made on the profit based on each of the results of the iterative method. That is, it can be stated that the profit will be above $0.65/GJ 80% of the time and that a refinery with this topology will not be profitable in 15% of the possible uncertainty realizations. All of this information can be very useful to an investor seeking to understand the risk associated with building a GTL refinery.

4 Overall Profit ($/GJ)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 60

2 0 −2 −4

Box Interval + Ellipsoidal Interval + Polyhedral

−6 0

0.1 0.2 0.3 0.4 a posteriori Probability of Constraint Violation

0.5

Figure 10: A comparison of the profits vs a posteriori probability of constraint violation for each of the uncertainty sets when only price uncertainty is considered.

5.6

Performance at Nominal Parameter Values

A final metric for understanding the robust solutions involves the performance of the topologies when the parameters are at their nominal values. The robust solutions report a profit which will be guaranteed no matter what parameters realizations in the uncertainty set actually occur. This allows the ability to make guarantees about the profit being as high or higher than the robust profit, but it also emphasizes the “worst-case” possibility of that uncertainty set. If the parameters actually realize their nominal values, effectively their ex46

ACS Paragon Plus Environment

Page 47 of 60

Product Produced (kBPD)

60

LPG Kerosene Diesel Gasoline

50 40 30 20 10 0 0.05 0.10 0.15 0.20 0.25 0.30 a posteriori Probability of Constraint Violation

Figure 11: Product distribution at low probabilities of constraint violation using the box robust counterpart.

60 Product Produced (kBPD)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

LPG Kerosene Diesel Gasoline

50 40 30 20 10 0 0.05 0.10 0.15 0.20 0.25 0.30 a posteriori Probability of Constraint Violation

Figure 12: Product distribution at low probabilities of constraint violation using the interval + ellipsoidal robust counterpart.

47

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research

60 Product Produced (kBPD)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 48 of 60

LPG Kerosene Diesel Gasoline

50 40 30 20 10 0 0.05 0.10 0.15 0.20 0.25 0.30 a posteriori Probability of Constraint Violation

Figure 13: Product distribution at low probabilities of constraint violation using the interval + polyhedral robust counterpart. pected values given a uniform distribution, then the profits will be much higher than those reported by the robust solution. Table 4 reports the robust profits and profits at the nominal parameter values for each of the results gained through the iterative method. In all cases, as would be expected, the profits calculated with the parameters at their nominal values are significantly higher than those reported by the robust solutions. Given that the box set always results in the nominal topology, it is not surprising that its “nominal equivalent” column contains profits that are effectively at the nominal solution of approximately $4.93/GJ. This is true for the interval + ellipsoidal set at post ≥ 0.25 as well. The most interesting results involve the impact of diversification when the parameters are at their nominal values. As seen in the previous section, the interval + polyhedral uncertainty set provides the same product distribution for each of the iterative method results, and thus the column of profits at the nominal parameter values are all approximately the same, at around $3.52/GJ.

48

ACS Paragon Plus Environment

Page 49 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

This is a very significant drop in profit, representing over 27% of the nominal value of $4.93/GJ. The diversification of the interval + ellipsoidal set is not as drastic, meaning its profits at the nominal parameter values are not as poor as the interval + polyhedral uncertainty set. Yet, in all cases of diversification there is at least a $1/GJ drop in profit. Given the relatively low amount of uncertain parameter realizations which would cause diversification to be beneficial, it is perhaps difficult to justify choosing a robust topology that produces all of the products unless the investor is extremely risk averse. All of this analysis points to the value of utilizing probabilistic bounds for GTL refineries under uncertainty. With probabilistic bounds, it can be seen that the nominal topology and product distribution would be the robust optimal solution for at least 75% of the uncertain parameter realizations. Furthermore, this topology performs significantly better than the other robust optimal topologies if parameter values are realized at or near the nominal values. Thus, with the assumptions made in this study, one can reasonably invest in a GTL refinery which uses a methanol intermediate to make gasoline and LPG, knowing that this topology is both a reasonably robust and profitable decision. Table 4: Comparison of the profit guaranteed by the robust solutions to the nominally equivalent profits provided by the topologies when all parameters are at their nominal values. The nominal solution from running the box uncertainty set at Ψ = 0 reports a profit of approximately $4.93/GJ post Robust $/GJ 0.05 0.10 0.15 0.20 0.25 0.30

-2.08 -0.95 -0.11 0.65 1.43 2.17

Box Nominal Equivalent $/GJ 4.85 4.85 4.83 4.85 4.93 4.96

Interval + Ellipsoidal Robust Nominal Equivalent $/GJ $/GJ -1.29 -0.25 0.47 1.03 1.44 2.13

3.75 3.80 3.82 3.88 4.93 4.93

49

ACS Paragon Plus Environment

Interval + Polyhedral Robust Nominal Equivalent $/GJ $/GJ -1.44 -0.41 0.34 0.87 1.44 1.88

3.53 3.52 3.54 3.50 3.55 3.52

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

6

Conclusion

Robust optimization provides a convenient and efficient means for deterministically including uncertainty in optimization problems. The development of tight a priori and a posteriori probabilistic bounds to characterize robust solutions provides meaningful answers to important problems, in this case centering on the production of liquid transportation fuels via natural gas under feedstock price, product price, and investment cost uncertainty. The analysis conducted has demonstrated the differences that robust counterpart formulations have on the overall profits and product distributions of a GTL refinery. An iterative algorithm to gain high quality a posteriori solutions was implemented successfully on a large-scale MINLP. Robust solutions with guaranteed profits were found for probabilities ranging from 5-30%, with the nominal price of robustness dropping below 50% at probabilities of constraint violation of 15% and higher. The impact of uncertainty on optimal topologies and product distributions were discussed in detail. As a result, it is seen that the nominal topology is also considered robust for probabilities of constraint violation of 25% and higher, as the nominal product distribution is found at this point using the interval + ellipsoidal set. Thus, given the assumptions of a uniform distribution for uncertain parameters, a potential investor can build a GTL refinery estimating that the profit will remain above $0.65/GJ 80% the time in a refinery that will also give a profit of $4.93/GJ if the parameters remain at their nominal values. The work conducted here sets a precedent for the widespread inclusion of uncertainty via robust optimization for many large-scale process systems applications.

7

Acknowledgments

Logan R. Matthews, Yannis A. Guzman, Onur Onel, and Alexander M. Niziolek are all grateful for the leadership and support of Professor Christodoulos A. Floudas. Professor Floudas passed away on August 14th, 2016. He provided comments on an early draft of this work before he passed away. The authors gratefully acknowledge financial support from 50

ACS Paragon Plus Environment

Page 50 of 60

Page 51 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

the National Science Foundation (NSF CBET-1548540). This research was conducted with Government support under and awarded by DoD, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a, which LRM gratefully acknowledges.

8

Supporting Information

Supporting Information is available online and contains the following items: a table containing investment cost parameters for each unit, a table providing the initial ∆ values for each case study as calculated via a priori bounds, a table with the complete results for the box uncertainty set, and the full MINLP model for the GTL superstructure. This information is available free of charge via the Internet at http://pubs.acs.org/.

51

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Floudas, C. A.; Niziolek, A. M.; Onel, O.; Matthews, L. R. Multi-scale systems engineering for energy and the environment: Challenges and opportunities. AIChE Journal 2016, 62, 602–623. (2) U.S. Energy Information Administration, Annual Energy Outlook 2015 with projections to 2040 ; 2015. (3) Baliban, R. C.; Elia, J. A.; Floudas, C. A. Novel Natural Gas to Liquids Processes: Process Synthesis and Global Optimization Strategies. AIChE Journal 2013, 59, 505– 531. (4) Baliban, R. C.; Elia, J. A.; Misener, R.; Floudas, C. A. Global optimization of a MINLP process synthesis model for thermochemical based conversion of hybrid coal, biomass, and natural gas to liquid fuels. Computers & Chemical Engineering 2012, 42, 64 – 86. (5) Baliban, R. C.; Elia, J. A.; Floudas, C. A. Biomass and Natural Gas to Liquid Transportation Fuels: Process Synthesis, Global Optimization, and Topology Analysis. Industrial & Engineering Chemistry Research 2013, 52, 3381–3406. (6) Baliban, R. C.; Elia, J. A.; Weekman, V.; Floudas, C. A. Process synthesis of hybrid coal, biomass, and natural gas to liquids via Fischer-Tropsch synthesis, ZSM5 catalytic conversion, methanol synthesis, methanol-to-gasoline, and methanol-toolefins/distillate technologies. Computers & Chemical Engineering 2012, 47, 29 – 56. (7) Floudas, C. A.; Elia, J. A.; Baliban, R. C. Hybrid and single feedstock energy processes for liquid transportation fuels: A critical review. Computers & Chemical Engineering 2012, 41, 24 – 51. (8) Onel, O.; Niziolek, A. M.; Floudas, C. A. Optimal Production of Light Olefins from Nat-

52

ACS Paragon Plus Environment

Page 52 of 60

Page 53 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

ural Gas via the Methanol Intermediate. Industrial & Engineering Chemistry Research 2016, 55, 3043–3063. (9) Niziolek, A. M.; Onel, O.; Floudas, C. A. Production of benzene, toluene, and xylenes from natural gas via methanol: Process synthesis and global optimization. AIChE Journal 2016, 62, 1531–1556. (10) Onel, O.; Niziolek, A. M.; Butcher, H.; Wilhite, B. A.; Floudas, C. A. Multi-scale approaches for gas-to-liquids process intensification: CFD modeling, process synthesis, and global optimization. Computers & Chemical Engineering 2017, 105, 276 – 296. (11) Grossmann, I. E.; Sargent, R. W. H. Optimum design of chemical plants with uncertain parameters. AIChE Journal 1978, 24, 1021–1028. (12) Halemane, K. P.; Grossmann, I. E. Optimal process design under uncertainty. AIChE Journal 1983, 29, 425–433. (13) Pistikopoulos, E.; Ierapetritou, M. Novel approach for optimal process design under uncertainty. Computers & Chemical Engineering 1995, 19, 1089 – 1110. (14) Sahinidis, N. V. Optimization under uncertainty: State-of-the-art and opportunities. Computers & Chemical Engineering 2004, 28, 971 – 983. (15) Acevedo, J.; Pistikopoulos, E. N. A Parametric MINLP Algorithm for Process Synthesis Problems under Uncertainty. Industrial & Engineering Chemistry Research 1996, 35, 147–158. (16) Acevedo, J.; Pistikopoulos, E. N. Stochastic optimization based algorithms for process synthesis under uncertainty. Computers & Chemical Engineering 1998, 22, 647 – 671. (17) Steimel, J.; Harrmann, M.; Schembecker, G.; Engell, S. A framework for the modeling and optimization of process superstructures under uncertainty. Chemical Engineering Science 2014, 115, 225 – 237. 53

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(18) Hen´e, T. S.; Dua, V.; Pistikopoulos, E. N. A Hybrid Parametric/Stochastic Programming Approach for Mixed-Integer Nonlinear Problems under Uncertainty. Industrial & Engineering Chemistry Research 2002, 41, 67–77. (19) Chaudhuri, P. D.; Diwekar, U. M. Process synthesis under uncertainty: A penalty function approach. AIChE Journal 1996, 42, 742–752. (20) Wendt, M.; Li, P.; Wozny, G. Nonlinear Chance-Constrained Process Optimization under Uncertainty. Industrial & Engineering Chemistry Research 2002, 41, 3621–3629. (21) Geraili, A.; Romagnoli, J. A. A multiobjective optimization framework for design of integrated biorefineries under uncertainty. AIChE Journal 2015, 61, 3208–3222. (22) Cheali, P.; Quaglia, A.; Gernaey, K. V.; Sin, G. Effect of Market Price Uncertainties on the Design of Optimal Biorefinery Systems: A Systematic Approach. Industrial & Engineering Chemistry Research 2014, 53, 6021–6032. (23) Rizwan, M.; Zaman, M.; Lee, J. H.; Gani, R. Optimal processing pathway selection for microalgae-based biorefinery under uncertainty. Computers & Chemical Engineering 2015, 82, 362 – 373. (24) Kasivisvanathan, H.; Ubando, A. T.; Ng, D. K. S.; Tan, R. R. Robust Optimization for Process Synthesis and Design of Multifunctional Energy Systems with Uncertainties. Industrial & Engineering Chemistry Research 2014, 53, 3196–3209. (25) He, C.; Pan, M.; Zhang, B.; Chen, Q.; You, F.; Ren, J. Monetizing shale gas to polymers under mixed uncertainty: Stochastic modeling and likelihood analysis. AIChE Journal 2018, 64, 2017–2036. (26) Gao, J.; You, F. Deciphering and handling uncertainty in shale gas supply chain design and optimization: Novel modeling framework and computationally efficient solution algorithm. AIChE Journal 2015, 61, 3739–3755. 54

ACS Paragon Plus Environment

Page 54 of 60

Page 55 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(27) Gargalo, C. L.; Carvalho, A.; Gernaey, K. V.; Sin, G. Optimal Design and Planning of Glycerol-Based Biorefinery Supply Chains under Uncertainty. Industrial & Engineering Chemistry Research 2017, 56, 11870–11893. (28) Soyster, A. L. Convex Programming with Set-Inclusive Constraints and Applications to Inexact Linear Programming. Operations Research 1972, 21, pp. 1154–1157. (29) Ben-Tal, A.; Nemirovski, A. Robust solutions of uncertain linear programs. Operations Research Letters 1999, 25, 1 – 13. (30) El Ghaoui, L.; Lebret, H. Robust Solutions to Least-Squares Problems with Uncertain Data. SIAM Journal on Matrix Analysis and Applications 1997, 18, 1035–1064. (31) Bertsimas, D.; Sim, M. The Price of Robustness. Operations Research 2004, 52, 35–53. (32) Lin, X.; Janak, S. L.; Floudas, C. A. A new robust optimization approach for scheduling under uncertainty: I. Bounded uncertainty. Computers & Chemical Engineering 2004, 28, 1069 – 1085. (33) Janak, S. L.; Lin, X.; Floudas, C. A. A new robust optimization approach for scheduling under uncertainty: II. Uncertainty with known probability distribution. Computers & Chemical Engineering 2007, 31, 171 – 195. (34) Li, Z.; Ding, R.; Floudas, C. A. A Comparative Theoretical and Computational Study on Robust Counterpart Optimization: I. Robust Linear Optimization and Robust Mixed Integer Linear Optimization. Industrial & Engineering Chemistry Research 2011, 50, 10567–10603. (35) Matthews, L. R.; Guzman, Y. A.; Floudas, C. A. Generalized robust counterparts for constraints with bounded and unbounded uncertain parameters. Computers & Chemical Engineering 2017, In Press.

55

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(36) Li, Z.; Tang, Q.; Floudas, C. A. A Comparative Theoretical and Computational Study on Robust Counterpart Optimization: II. Probabilistic Guarantees on Constraint Satisfaction. Industrial & Engineering Chemistry Research 2012, 51, 6769–6788. (37) Li, Z.; Floudas, C. A. A Comparative Theoretical and Computational Study on Robust Counterpart Optimization: III. Improving the Quality of Robust Solutions. Industrial & Engineering Chemistry Research 2014, 53, 13112–13124. (38) Guzman, Y. A.; Matthews, L. R.; Floudas, C. A. New a priori and a posteriori probabilistic bounds for robust counterpart optimization: I. Unknown probability distributions. Computers & Chemical Engineering 2016, 84, 568 – 598. (39) Guzman, Y. A.; Matthews, L. R.; Floudas, C. A. New a priori and a posteriori probabilistic bounds for robust counterpart optimization: II. A priori bounds for known symmetric and asymmetric probability distributions. Computers & Chemical Engineering 2017, 101, 279 – 311. (40) Guzman, Y. A.; Matthews, L. R.; Floudas, C. A. New a priori and a posteriori probabilistic bounds for robust counterpart optimization: III. Exact and near-exact a posteriori expressions for known probability distributions. Computers & Chemical Engineering 2017, 103, 116 – 143. (41) Baliban, R. C.; Elia, J. A.; Floudas, C. A. Simultaneous process synthesis, heat, power, and water integration of thermochemical hybrid biomass, coal, and natural gas facilities. Computers & Chemical Engineering 2012, 37, 297 – 327. (42) Niziolek, A. M.; Onel, O.; Elia, J. A.; Baliban, R. C.; Xiao, X.; Floudas, C. A. Coal and Biomass to Liquid Transportation Fuels: Process Synthesis and Global Optimization Strategies. Industrial & Engineering Chemistry Research 2014, 53, 17002–17025. (43) Niziolek, A. M.; Onel, O.; Elia, J. A.; Baliban, R. C.; Floudas, C. A. Coproduction of

56

ACS Paragon Plus Environment

Page 56 of 60

Page 57 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

liquid transportation fuels and C6 C8 aromatics from biomass and natural gas. AIChE Journal 2015, 61, 831–856. (44) Ben-Tal, A.; Nemirovski, A. Robust solutions of linear programming problems contaminated with uncertain data. Mathematical programming 2000, 88, 411–424. (45) Chou, V.; Kuehn, N. Assessment of Hydrogen Production with CO2 Capture, Volume 1: Baseline State of the Art Plants; 2010. (46) de Klerk, A. Fischer-Tropsch Refining; Wiley, 2011. (47) Steynberg, A.; Dry, M. Fischer-Tropsch Technology; Studies in Surface Science and Catalysis; Elsevier Science, 2004. (48) Mobil Research and Development Corporation, Slurry Fischer-Tropsch/Mobil Two Stage Process of Converting Syngas to High Octane Gasoline; 1983. (49) Research, M.; Corporation, D. Two-Stage Process For Conversion of Synthesis Gas to High Quality Transportation Fuels; 1985. (50) Matthews, L. R.; Niziolek, A. M.; Onel, O.; Pinnaduwage, N.; Floudas, C. A. Biomass to Liquid Transportation Fuels via Biological and Thermochemical Conversion: Process Synthesis and Global Optimization Strategies. Industrial & Engineering Chemistry Research 2016, 55, 3203–3225. (51) Bechtel, Baseline Design/Economics for Advanced Fischer-Tropsch Technology; 1992. (52) Bechtel, Fuels Process Flowsheet Simulation Model of a Battelle Biomass-Based Gasification, Fischer-Tropsch Liquefaction, and Combined-Cycle Power Plant; 1998. (53) National Reneweable Energy Laboratory, Gasoline from wood via integrated gasification, synthesis, and methanol-to-gasoline technologies; 2011.

57

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(54) Tabak, S.; Yurchak, S. Conversion of methanol over ZSM-5 to fuels and chemicals. Catalysis Today 1990, 6, 307 – 327. (55) Mokrani, T.; Scurrell, M. Gas Conversion to Liquid Fuels and Chemicals: The Methanol Route-Catalysis and Processes Development. Catalysis Reviews 2009, 51, 1–145. (56) Tabak, S.; Krambeck, F. Shaping process makes fuels. Hydrocarbon Process 1985, 64 . (57) Tabak, S. A.; Krambeck, F. J.; Garwood, W. E. Conversion of propylene and butylene over ZSM-5 catalyst. AIChE Journal 1986, 32, 1526–1531. (58) Kreutz, T. G.; Larson, E. D.; Liu, G.; Williams, R. Fischer-Tropsch fuels from coal and biomass. Proceedings of the 25th Annual International Pittsburgh Coal Conference. 2008. (59) IBM, CPLEX User’s Manual Version 12 Release 6. 2013. (60) ARKI Consulting & Development A/S, CONOPT 3. 2015. (61) Baliban, R. C.; Elia, J. A.; Floudas, C. A. Optimization framework for the simultaneous process synthesis, heat and power integration of a thermochemical hybrid biomass, coal, and natural gas facility. Computers & Chemical Engineering 2011, 35, 1647 – 1690. (62) Elia, J. A.; Baliban, R. C.; Floudas, C. A. Toward Novel Hybrid Biomass, Coal, and Natural Gas Processes for Satisfying Current Transportation Fuel Demands, 2: Simultaneous Heat and Power Integration. Industrial & Engineering Chemistry Research 2010, 49, 7371–7388. (63) U.S. Energy Information Administration, United States Natural Gas Industrial Price. 2018. (64) U.S. Energy Information Administration, Refiner Petroleum Product Prices by Sales Type. 2018. 58

ACS Paragon Plus Environment

Page 58 of 60

Page 59 of 60 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(65) Wall Street Journal, Cash Prices - Market Data Center. 2018. (66) Yuan, Y.; Li, Z.; Huang, B. Robust optimization under correlated uncertainty: Formulations and computational study. Computers & Chemical Engineering 2016, 85, 58 – 71. (67) Shang, C.; Huang, X.; You, F. Data-driven robust optimization based on kernel learning. Computers & Chemical Engineering 2017, 106, 464 – 479. (68) Ning, C.; You, F. Data-driven decision making under uncertainty integrating robust optimization with principal component analysis and kernel smoothing methods. Computers & Chemical Engineering 2018, 112, 190 – 210. (69) Ning, C.; You, F. Data-driven stochastic robust optimization: General computational framework and algorithm leveraging machine learning for optimization under uncertainty in the big data era. Computers & Chemical Engineering 2018, 111, 115 – 133.

59

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table of Contents Graphic

Figure 14: Table of Contents Graphic

60

ACS Paragon Plus Environment

Page 60 of 60