Financial Risk Management in Heat Exchanger Networks Considering

May 16, 2018 - This article is part of the 2017 European Symposium on ... To this end, three different risk metrics are employed: variability index, d...
0 downloads 0 Views 3MB Size
Subscriber access provided by Kaohsiung Medical University

Process Systems Engineering

Financial risk management in heat exchanger networks considering multiple utility sources with uncertain costs L. V. Pavão, Carlos Pozo-Fernandez, Laureano Jiménez, Mauro A. S. S. Ravagnani, and Caliane Bastos Borba Costa Ind. Eng. Chem. Res., Just Accepted Manuscript • Publication Date (Web): 16 May 2018 Downloaded from http://pubs.acs.org on May 16, 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 55 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

Financial risk management in heat exchanger networks considering multiple utility sources with uncertain costs Leandro V. Pavãoa, Carlos Pozob, Laureano Jiménezc, Caliane B. B. Costaa, Mauro A. S. S. Ravagnania1 a

Department of Chemical Engineering, State University of Maringá,

Av. Colombo, 5790, Bloco D90, CEP 87020900, Maringá, PR, Brazil b

Centre for Process Systems Engineering (CPSE), Imperial College London SW7 2AZ, United Kingdom c

Departament d’Enginyeria Química, Universitat Rovira i Virgili Av. Països Catalans, 26, 43007, Tarragona, Spain

ABSTRACT In a plant, different utilities may be produced from different sources whose costs vary singularly. Therefore, evaluating the financial risk associated to heat exchanger networks (HEN) with multiple utilities is essential. This work aims to address the HEN synthesis under uncertainty by presenting several optimization models for managing financial risk. To this end, three different risk metrics are employed: Variability Index (VI), Downside risk (DR) and Risk Area Ratio (RAR). The resulting models are based on an enhanced version of the stage-wise superstructure (SWS), which incorporates the potential to place utilities at all HEN stages. The study considers fixed utility sources as well as the possibility to vary those yearly. Furthermore, we propose a novel “pseudointersection” calculation scheme for the RAR metric, which proves efficient in assessing how risks can be reduced at minimal opportunity loss. The proposed models are able to efficiently identify solutions appealing to different investor profiles.

1

Corresponding author. Tel: +55 (44) 3011-4774, Fax: +55 (44) 3011-4793

E-mail address: [email protected]

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

Keywords: Uncertainty; Financial risk management; Multi-objective optimization; Heat Exchanger Networks; Meta-heuristics 1

Introduction

Energy integration is a recurrent subject in the Process Engineering literature. Industrial processes might be highly energy-demanding in order to maintain appropriate operating conditions. For that reason, efficient energy recovery systems able to reduce the use of auxiliary energy exchange sources (i.e., heating and cooling via utilities) can yield substantial economic savings. Such potential has drawn much attention to energy integration strategies such as the optimal synthesis of heat exchanger networks (HEN). Essentially, the HEN synthesis problem consists of finding a configuration of heat exchangers integrating process streams with optimal trade-off between operating cost (i.e., auxiliary heating/cooling requirements) and investment cost (i.e., number and areas of heat exchanger units). Several works have tackled the problem with mathematical programming formulations minimizing total annual costs (TAC), with noteworthy contributions from Floudas et al.1 and of Yee and Grossmann.2 The former proposed a hyperstructure that gave rise to a mixed integer nonlinear programming (MINLP) problem comprising heat exchangers with all possible stream matches as well as several piping options, such as splits, by-passes, crossflows and sequential heat exchangers in single stream split branches. The latter presented the stage-wise superstructure (SWS) concept, from which an MINLP model was also derived. This model is able to lead to cost-effective HEN configurations despite being structurally simpler than the hyperstructure1 (e.g., one unit per stream split branch, utilities placed at streams’ ends and crossflows not possible). Other authors took advantage of insights derived from the hyperstructure and the SWS models and combined them with new concepts to develop novel formulations. For instance, the SWS-based model of PonceOrtega et al.3 comprised utilities placement in an additional stream split branch in all stages. Also based on the SWS, Huo et al.4 presented a model that had splits in some of the stages, while others were composed by serial units only. Huang and Karimi5 presented a hybrid superstructure,

ACS Paragon Plus Environment

Page 2 of 55

Page 3 of 55 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

combining the stages concept of Yee and Grossmann2 with structural insights based on the hyperstructure of Floudas et al.1 Kim et al.6 presented a superstructure containing stages and substages, enabling the use of serial heat exchangers in stream split branches. Pavão et al.7 introduced an extended superstructure model based on the SWS, allowing multiple heaters or coolers in intermediate and ending stages of the superstructure. Those units could be placed one per stream split branch, and a branch was present for each available utility. Given their nonlinearities, non-convexities and combinatorial nature, solving HEN synthesis MINLP models is an onerous task, which requires efficient optimization techniques for achieving near-optimal solutions. While deterministic solution techniques, such as reduced-gradient algorithms and outer-approximation, were used in the aforementioned fundamental works,1,2 stochastic techniques have proven to be an appealing alternative to deterministic methods for HEN synthesis problems. These methods, commonly referred to as meta-heuristics, rely on randomized searches and heuristic rules for solutions guidance throughout the problem space. One of the key advances of meta-heuristics under the scope of HEN synthesis was the two-level HEN optimization established by Lewin and co-workers: These authors applied a scheme based on a Genetic Algorithm (GA) to models derived from superstructures with8 and without splits9. In these models, binary variables are optimized in the “upper-level” algorithm, which generates structures whose continuous variables are optimized in the “lower-level”. Ravagnani et al.10 employed a GA method along with Pinch Technology concepts. Luo et al.11 presented a hybrid GA which encompassed Simulated Annealing based strategies as well. This methodology was extended to the monogenetic algorithm of Fieg et al.12, which was able to handle larger problems by splitting them into smaller sub-problems. Silva et al.13 used Particle Swarm Optimization (PSO) for HEN synthesis with stream splits using the SWS. Pavão et al.14 presented a two-level GA/PSO method with parallel processing techniques, which was applied to the SWS with splits. Zhang et al.15 presented a “chessboard” HEN representation and used the Random Walking Algorithm of Compulsory Evolution (RWCE) meta-heuristic to solve the resulting model. Pavão et al.16 presented a two-level

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

SA/PSO method applied to a no-splits SWS, making use of parallel computing. Concatenating insights from GA/PSO14 and SA/PSO16 methods, Pavão et al.17 developed a two-level approach coupling SA and Rocket Fireworks Optimization (RFO). RFO was, in fact, a two-stage continuous optimization approach, which used a continuous variant of simulated annealing to find a promising solution prior to the application of PSO. The interested reader is also referred to the work of Toimil and Gómez,18 who presented a broad review on meta-heuristic methods applied to the HEN synthesis problem up to 2014, approaching particularities of each method such as superstructures used, optimization variables and solution schemes. As an essential component of industrial plants, HENs are subject to uncertain factors such as operating temperatures and flow-rates. A widely used approach to include such uncertainties in optimization models is that of the multiperiod HEN synthesis. In brief, such methodology entails designing a minimal-costs HEN able to operate feasibly (i.e., streams achieve their target temperatures) under a finite set of different operating conditions, each one characterizing a period with an estimated duration. Aaltola19 presented a simultaneous MINLP model derived from the SWS2 for synthesizing multiperiod HENs. In Aaltola’s model, the cost calculation associated to heat exchangers areas is based on an average of the areas required for all the conditions, instead of the more realistic assumption of employing the maximal area value. While this assumption avoided the disjunctions in the mathematical model arisen from the selection of the greater unit, it underestimates the real cost of the HEN since a heat exchanger able to operate under different conditions must always be sized according to the condition that demands the greater area. Verheyen and Zhang20 were able to debottleneck the simultaneous multiperiod HEN synthesis procedure, presenting a SWS-based model that could reach near-optimal solutions considering the greatest required areas. Pavão et al.21 applied a meta-heuristic approach to solve a multiperiod HEN model derived from the SWS, achieving promising results in terms of total annual costs. Jiang and Chang22 alleviated over-sized heat exchangers and unit idleness issues by presenting their so-called timesharing mechanism. Under that approach, a heat exchanger could perform different matches

ACS Paragon Plus Environment

Page 4 of 55

Page 5 of 55 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

according to the operating period. Miranda et al.23 presented an enhanced methodology based on that of Jiang and Chang22 for multiperiod HEN with timesharing mechanisms, achieving better solutions. Pavão et al.24 included a re-optimization meta-heuristic scheme after the heat exchangers rearrangement strategy of Jiang and Chang,22 leading to noteworthy multiperiod HEN results. In addition to operating conditions, economic parameters associated to energy sources, such as fuels and electricity costs, might also present variable behaviors throughout the plant lifetime. Such market volatilities might lead annual costs of a HEN to differ from those applicable by the time of the plant commissioning. Despite this, uncertainties in HEN economy-related factors have been approached in the literature only more recently. In the work of Mou and Qvale25, two scenarios were considered, with differences in utility costs and interest rate. A HEN was synthesized for the first scenario conditions, and further retrofitted to meet the conditions of the second one. The work of Nemet et al.26 addressed five scenarios for fuel costs, using pessimistic to optimistic predictions. Occurrence probabilities were assigned to each of the considered scenarios, which were examined simultaneously in a multiperiod optimization model. Later, additional economic risk assessment measures were included in the methodology.27 The methodology proposed by Pintarič and Kravanja28 considered several uncertain parameters. Heuristics were used to identify most critical scenarios, which guided their optimization procedure. Further, the obtained HEN configuration was tested regarding its flexibility via scenarios generation by Monte Carlo sampling (MCS) and by a flexibility index. Pavão et al.29 considered uncertain utility costs by modeling them as stochastic variables following a Geometric Brownian Motion (GBM) process30. The authors proposed a bicriteria optimization model considering the average TAC and the Downside Risk (DR) metric.31 Further, five bi-criteria optimization models were proposed by the authors,32 where the Expected Total Annual Costs (ETAC, average of the TAC of all scenarios) was always the first objective and the second objective was one of the five following financial risk metrics: Variance, Variability Index,33 Alternative Probabilistic Financial Risk (which was developed therein based on the

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

Probabilistic Financial Risk by Barbaro and Bagajewicz),34 Downside Risk31 and Worst Case optimization. Although Pavão et al.29 were the first to consider discretized stochastic variables into a large number of scenarios for addressing financial risk metrics in HEN synthesis, it is worth noting that such sort of approach was previously employed in several other fields in the Process Engineering literature. The use of such techniques for conducing robust process optimization was identified as an encouraging direction for future research by Sahinidis,35 who also presented a comprehensive state-of-the-art review up to the year of 2003. Investigations comprising uncertainties evaluation and stochastic modelling were conducted regarding case studies such as the natural gas commercialization in Asia,36 phosphoric acid production,37 absorption cooling systems,38 global supply chains planning,39 hydrogen supply chains,40 bioethanol and sugar infrastructures,41 chemical supply chains,42 environmental and economic performances of general industrial processes,43 biofuels supply chains,44 environmental efficiency for several electricity sources,45 and the reduction of emissions in the United States electricity production.46 To the best of the authors’ knowledge, no other work apart from those of Pavão et al.29,32 has attempted to handle uncertainties via stochastic discretization of uncertain variables and use of optimization models containing statistics-based risk metrics. This contribution aims to broaden these discussions by introducing a new optimization framework with four main novel features not yet employed in previous works. First, (i) different fuels — entailing independent cost fluctuations— can be used in multiple hot utility production units. Similarly, electricity cost fluctuations affecting the operating costs of auxiliary cooling units are also considered. Several future scenarios for these fluctuations are generated by simulating a GBM process30 via MCS. Secondly, (ii) we implement the aforementioned new features by adapting a recently developed HEN SWS-based mathematical model7 that incorporates new placement options for utilities. The model encompasses multiple utilities placement not only at stream ends, but also at intermediate stages in separate stream split branches. From that model, a bi-criteria formulation is

ACS Paragon Plus Environment

Page 6 of 55

Page 7 of 55 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

derived considering the expected total annual costs (i.e., the average of the annual costs of all scenarios) and the variability index (VI) risk metric. The Downside Risk (DR)31 and the Risk Area Ratio (RAR)36 metrics are considered as well. Since the implementation of the latter metric in the solution approach adopted here would be computationally burdening, (iii) we introduce a “pseudointersection” calculation scheme for providing a suitable approximation of that metric at reasonable computational time. Such new feature, which constitutes the third novelty of this contribution, led to the inclusion of new constraints and an additional variable to the problem. In order to evaluate different structural options, (iv) a flexibility study regarding utility sources is presented. In such experiment, in case that one utility might be produced in different sources, it is allowed that the cheaper is always used. It is thus possible to assess how such structural change would affect financial risk. The solution approach used is an enhanced version of SA-RFO,17 adapted for dealing with the bicriteria optimization of the new mathematical model, as well as with the aforementioned additional features of the pseudo-intersection RAR model. The bi-criteria optimization strategy was based on the popular ε-constraint method.47 It is worth noting that the use of a meta-heuristic method constitutes an interesting feature of this work, since previously reported uncertainty related investigations in HEN synthesis tackle their models with deterministic methods. From the technical point of view, some attractive features are present in meta-heuristics, such as being derivative-free methods, or allowing to deal with objective functions as “black-boxes”. These aspects yield some promising opportunities, such as using coding statements (if, else, for, while, etc.) within the objective functions. That facilitates, for example, the handling of constraints and disjunctions. Moreover, deterministic HEN optimization studies are, in their majority, implemented in commercial solving platforms, such as GAMS. On the other hand, this approach was coded in C++, a programming language that can be implemented in several free of charge integrated development interfaces. It is worth mentioning, though, that solutions’ optimality cannot be guaranteed. In fact,

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

to the best of the authors’ knowledge, no method, deterministic or stochastic, was able to reach global optimality for industrial-scale HEN synthesis. 2 2.1

Heat exchanger network synthesis under uncertainty Problem statement

A given set of |NH| hot and |NC| cold streams must be cooled/heated to target temperatures by exchanging heat among themselves and/or with one or more of the |NHU| types of hot utilities and/or |NCU| types of cold utilities. Heat capacities and temperatures for all streams and types of utilities are known. The cost function for heat exchangers, based on their areas, is available. Costs per energy unit are considered uncertain given the market fluctuations on the commodities used for the production of each utility. The energy costs are treated as stochastic parameters, which can be discretized with aid of sampling methods. Such tools can be used in forecasting future costs based on historical data. We aim to synthesize cost-efficient HEN’s considering also its associated financial risk evaluated via risk metrics. 2.2

Mathematical formulation

The mathematical model used here for HEN synthesis is an enhanced version of the SWS2 introduced by Pavão et al.7 As in the original SWS, this model uses the concept of stages to consider all the possibilities of heat exchange among process streams at each stage. To this end, a superstructure is built as follows. At each stage k, every stream is split into a set of branches, each containing a single unit (i.e., heat exchanger). In the original SWS2, heaters were always placed at cold streams ends, while coolers were enforced at the ends of hot streams. Conversely, in the model employed here,7 heaters/coolers, each for a different type of utility, are placed at every stream split branch at all the |Ns| superstructure stages. Figure 1 depicts the enhanced SWS in a hypothetical case with two hot and two cold streams (i.e., |NH| = |NC| = 2), illustrating, in addition, the multiple utility sources that are being considered. Specifically, there are two suppliers for hot utilities, each using a different fuel (i.e., |NHU| = 2), and two for cold utilities (|NCU| = 2). As a result, each of the

ACS Paragon Plus Environment

Page 8 of 55

Page 9 of 55 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

|NH| hot streams is split into |NC|+|NCU| branches at each stage k, while each of the |NC| cold streams is split into |NH|+|NHU| branches. This model is suitable to this work’s proposal since it comprises efficient options for the optimal allocation of multiple utilities.

Figure 1. Stage-wise superstructure representation with multiple utilities from different sources A fundamental concept for the development of the stochastic models is the use of scenarios. Commonly, in HEN synthesis models, each utility has an associated yearly cost per energy unit (e.g., $/kW·y). When those parameters are deemed as uncertain, they must be considered as stochastic variables. Formally, this means they can no longer be represented as scalars, but rather as distributions. These distributions propagate through all the calculations until they reach the objective function, which consequently, yields a distribution as well. A way of dealing with stochastic parameters is to discretize them into a set Nsc of scenarios s. That is, the stochastic parameter becomes a vector containing |Nsc| equally probable values, each resulting from a different realization of the associated uncertainty. Those values are obtained here via a Monte Carlo sampling methodology, which will be better elucidated during its application in the Case Study (Section 5). In the interim, the scenario concept and the fact that it is used to obtain a discretization of the uncertain parameters suffice for understanding the models. The total annual costs for each scenario (TACs) can thus be calculated as follows:

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 55

  , ℎ ,  +    , , , + 















  ℎ,  ℎ, ℎ +    ℎ, ℎ,, + 







   ,,  +  ∙ ,,  +    , ,  +  ∙ , ,  + 













   ℎ,,  +  ∙ ℎ,,  + 





(1)



  ,  +  ∙ ,  + 





  ℎ,  +  ∙ ℎ,  + 



 ∈ !" , # ∈ !"$ , % ∈ !& ,  ∈ !&$ , ' ∈ !( , ) ∈ !(* where zi,j,k is a binary variable indicating the existence of a match, Ai,j,k is the area of a unit, Qi,j,k is its heat load, and B, C and β are heat exchanger cost parameters. Terms with the cu and hu suffixes are associated respectively to cold and hot utilities and those with inter and outlet suffixes to intermediate and outlet stages of the SWS. For instance, Qhuinteri,n,k regards the heat load of heaters in intermediate stages. Qcuoutleti is the total heat that was not exchanged in intermediate stages and is still to be exhausted at the hot stream end via cold utilities. FhQcuoutleti,n is a fraction of such heat that will be exhausted by means of a given cold utility n. The definitions for Qhuoutletj and FcQhuoutletm,j are analogous, but for cold streams and hot utilities. Heat balances, area calculations, thermodynamic and feasibility constraints of the model can be found in Pavão et al.7 We aim to minimize the TAC average value among all scenarios which is referred here to as the expected TAC (ETAC), assuming equal chances of occurrence for all scenarios. The corresponding single-objective optimization (SOO) model for this is written as follows:

(,_.//) # ). .

1, =  23  4 

,6). (3) − (58) ?ã  >.A

(2)

where probs is the probability of occurrence of each scenario (considered equal for all scenarios) and ETAC stands for expected TAC. Note that all constraints are analogous to those in the work of

ACS Paragon Plus Environment

Page 11 of 55 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 model derived in Pavão et al.7 and can be found therein. Those are related to heat balances, temperatures and area calculations, as well as thermodynamic constraints for feasibility. 3

Risk metrics

The variance in TAC among all scenarios in ETAC-optimal solutions (i.e, those which yield the minimal average TAC of the scenarios) is typically high, which means that scenarios much more expensive as well as much cheaper than ETAC are possible. Therefore, ETAC-optimal solutions are commonly not those with lower financial risk. In this section, we present the concepts behind the three risk metrics used in the present work to evaluate the advantages attained by each HEN configuration under different risk-based points of view, and to aid in finding the most suitable solution for different investor profiles. In the methodology proposed here for analyzing solutions with different risks, firstly, a bi-criteria optimization of ETAC and a variance-related objective will be performed (see Section 4 for details on the solution approach and its adaptation for MOO). That provides a wide set of solutions each entailing a different risk. The Downside Risk and Risk Area Ratio metrics are then employed for evaluating the application of those metrics to the HEN problem. Specific comparisons among the metrics and designs attained will be presented as models are applied. For the sake of elucidation, the concepts behind the three metrics are illustrated in Figure 2, by means of cumulative distribution function (CDF) plots. These concepts are outlined throughout Sections 3.1-3.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

Figure 2. CDFs illustrating the fundamentals of the three metrics employed in this work 3.1

Variability index

While the variance is a useful metric for risks management, the variability index (VI), proposed by Ahmed and Sahinidis33 as an alternative to the formal variance calculation, will be used here instead for two reasons. First, because the use of the VI avoids the quadratic terms present in the variance calculation. Secondly, because in previous investigation,32 the solution approach used here (SARFO17) was in general more effective in reducing solutions’ variance by means of the VI than with the real variance itself. Considering a cumulative distribution function (CDF) as that in Figure 2a, the VI is defined as the area to the right of the average value (ETAC) and above the CDF curve. That is, the VI is given by

ACS Paragon Plus Environment

Page 12 of 55

Page 13 of 55 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 sum of all positive deviations between the average and the distribution curve (i.e., all the situations where the cost is higher than the average). Mathematically, the VI is defined as in Eqs. (3) and (4).

C =  − ,, <  ≥ , , ) ∈ !(* B  C = 0, ℎF)

(3)

GH =  23 C

(4)



where Ds is the aforementioned positive deviation from the average value. Note that Ds is null if such deviation is negative. In the HEN case, uncertainties arise from utility costs. Therefore, reducing utilities requirement in might be a cost-effective strategy. However, we can anticipate a conflictive behavior between variance and ETAC, since costs might increase impractically in solutions with excessive heat exchange among process streams due to the large units required to perform such tasks. Therefore, it is advantageous to minimize the variance without increasing the TAC average (ETAC), which can be achieved in a bi-criteria optimization model considering these two objectives (i.e., ETAC_VI_MOO), as described in Eq. (5).

(,_GH_I//) # ). .

J,, GHK ,6). (3) − (58) ?ã  >.A

(5)

In case that VI is the only objective being optimized, the optimization model is by analogy called VI_SOO.

3.2

Downside risk

The concept behind the downside risk (DR) is to some extent similar to that of the VI. In a CDF plot, the DR represents the area enclosed by the curve and a given target value (instead of the

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 14 of 55

average, as in VI). This concept is demonstrated in Figure 2b. Put in other words, the DR is the sum of positive deviations between the curve and the aforementioned target value, denoted by Ω. The DR value should be reasonably small, so as to reduce the sum of the costs of scenarios that are more expensive than Ω simultaneously with their expected risk of occurrence. For that reason, in a DRbased optimization model, Ω is generally fixed at a value contained in the reference (i.e., ETAC optimal) solution’s curve. This may provide a “candidate” solution where some scenarios are cheaper. A solution is referred here to as “candidate” when it is obtained with a financial risk management model aiming to reduce the risk in comparison to that with the optimal ETAC value (i.e., reference solution). The DR is calculated as follows:

B

L =  − Ω, <  ≥ Ω , ) ∈ !(* L = 0, ℎF)

(6)

CN =  23 L

(7)



where δS is the positive deviation between TAC of a scenario and Ω. The optimization problem considering the DR is finally formulated as follows:

(CN_.//(Ω)) # ). .

3.3

JCNK ,6). (3) − (58) ?ã  >.A

(8)

Risk Area Ratio and the pseudo-intersection models

Aseeri and Bagajewicz36 proposed the risk area ratio (RAR) as a form of evaluating solutions regarding opportunity loss and risk reduction trade-offs. The metric proposes to compare a candidate solution to a reference one (in this study, the ETAC-optimal) by means of the ratio of the “opportunity loss area” to the “risk reduction area”. In Figure 2c, the former corresponds to the area enclosed between the CDF curves at the left-hand side of their intersection (on the abscissa), while

ACS Paragon Plus Environment

Page 15 of 55 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 latter corresponds to the area enclosed at the right-hand side of it. As pointed out by the Aseeri and Bagajewicz,36 that ratio is always greater than one, which means that the risk reduction is never greater than the opportunity loss. However, RAR values that are closer to one entail a better opportunity loss versus risk reduction trade-off in comparison to the reference solution. The straightforward approach to calculate such areas from a given number of discretized scenarios is as follows. Firstly, for the risk reduction area, considering that TAC vectors (i.e., vectors containing TACs values for all scenarios |NSc|) for both the reference and candidate solutions are sorted in ascending order:

O

PQQ = (RST − (RSU , < (RST ≥ (RSU , ) ∈ !(* PQQ = 0, ℎF)

(9)

where Sol1 superscript indicates the reference solution (ETAC-optimal), while Sol2 refers to a candidate solution and ξsRR is an auxiliary variable that is non-zero only when, for a given scenario, the candidate solution has lower TAC than that of the reference one. In that situation, the value it receives is the difference between the TAC of a given scenario in Sol2 and Sol1. The risk reduction area is thus calculated with:

NN (V =  23 PQQ 

(10)

The SV superscript stands for “sorted vector”, which is a feature of this calculation method and is used to differentiate it from the method that will be presented later. For the opportunity loss area:

O

PWX = (RSU − (RST , < (RST < (RSU , ) ∈ !(* PWX = 0, ℎF)

(11)

where the role played by ξsOL is analogous to ξsRR. Here, though, such auxiliary variable (ξsOL) is non-zero for a given scenario in which the candidate solution is worse (in terms of TAC) than the

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 16 of 55

reference solution. Finally, opportunity loss area, and consequently the risk area ratio can be calculated with:

/Z(V =  23 PWX 

NN (V =

/Z(V NN (V

(12)

(13)

An alternative form of calculating these areas is to firstly allocate the coordinate of the intersection point in the abscissas and then proceed as follows:

[ QQ =  − Ψ, <  ≥ Ψ , ) ∈ !(* B  [QQ = 0, ℎF)

(14)

[WX = Ψ −  , <  < Ψ , ) ∈ !(* [WX = 0, ℎF)

(15)

O

where Ψ is abscissa of the pseudo-intersection of the candidate and the reference solutions, and γsRR and γsOL are the deviations of TAC from Ψ at each scenario. Note that Ψ is called here pseudointersection for two reasons. The first one is due to the fact that discrete scenarios, instead of continuous CDFs, are being dealt with, and it is unlikely that two TACs each from one of the two solutions are exactly equal, meeting the definition of intersection. The second issue is that, in some cases, given that scenarios are stochastically generated, curves (obtained by linking the discrete points with straight lines) might appear to intersect more than once. Then, the equation that yields the risk reduction (RR) is written as the approximated area enclosed between Ψ and the reference solution minus that enclosed between Ψ and the candidate one:

NN ]^ =  23 [QQ,(RST −  23 [QQ,(RSU 



The opportunity loss (OL) area is calculated analogously, as follows:

ACS Paragon Plus Environment

(16)

Page 17 of 55 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

/Z]^ =  23 [WX,(RST −  23 [WX,(RSU 



(17)

The PI superscript indicates the pseudo-intersection based calculations. Hence, RARPI can be calculated with:

NN ]^ =

/Z]^ NN ]^

(18)

These two methods for calculating the RAR have some drawbacks that hamper their implementation in an optimization algorithm. On the one hand, the calculation of RARSV (Eq. (13)), scenarios in both TAC vectors must be sorted in ascending order of TACs, so that a CDF is simulated and the nth value (regarding TAC) in the candidate solution is accordingly compared to its correspondent in the reference solution. That implies that throughout the execution of an optimization procedure, those vectors would have to be constantly ordered. In the case of a metaheuristic approach, a sorting routine would have to be executed at all iterations. Implementing such routine would be possible, but its recurrent execution would be a computationally onerous procedure. One the other hand, the calculation of RARPI (Eq. (18)) requires as well that TAC vectors are sorted so that an approximate intersection point can be selected visually or by interpolating points. That approach could work in a manual calculation yet it would not be suitable for tailoring an optimization model. Here we propose a systematic form of choosing a suitable Ψ based on finding a value for which the number of scenarios more expensive than it is equal for both the candidate and the reference solutions, while RARPI is minimal. The main advantage is that such approach does not require the TAC vector to be sorted. Figure 3a illustrates 1,000 unsorted values contained in the TAC vectors of two solutions to the case study of this work (see Section 5). It is worth noting that the depicted pseudo-intersection value found for minimal RARPI was found without the need of sorting the vectors with the approach that will be presented later. Both solutions

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

have 156 scenarios more expensive than Ψ. Figure 3b represents the same set of values after sorting, and Figure 3c depicts a zoom over the region of the pseudo-intersection.

Figure 3. For a reference and a candidate solution: (a) Unsorted TAC vectors; (b) Sorted TAC vectors; (c) Zoom over the pseudo-intersection region Note that there are cases to the right side of the intersection where TAC is lower for the reference solution, as well as cases to the left side that are cheaper in the candidate solution (see for instance subplot in Figure 3c). Such issues might lead to a marginal error when comparing values obtained via RARPI to those from RARSV, which will be quantitatively assessed in the case study (Section 5). With the afore-described concepts in mind, an optimization model for finding a solution with minimal RARPI is proposed here. The model treats Ψ as an additional continuous variable to be manipulated by the optimization approach simultaneously to those from the HEN model. In that manner, a solution with minimum RARPI as corresponding to the best Ψ value can be found in a one-step procedure. Feasibility constraints are also introduced as follows.

ACS Paragon Plus Environment

Page 18 of 55

Page 19 of 55 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 order to ensure that a pseudo-intersection always exists, both risk reduction and opportunity loss areas must be greater than zero, which always renders solutions with some advantageous scenarios. The feasibility constraints that follow are imposed to the model:

/Z]^ > 0

(19)

NN ]^ > 0

(20)

which ensure that the curves associated to Sol1 and Sol2 intersect and that Sol2’s one is steeper (i.e., has lower variance). Consequently, valid solutions are those which have both a risk reduction area and an opportunity loss area. Some additional constraints are implemented for aiding the method in finding a proper Ψ value so that a minimum RAR is achieved. These are as follows:

`QQ.(RST = 1, < (RST ≥ Ψ , ) ∈ !(* O `QQ,(RST = 0, ℎF)

(21)

`QQ.(RSU = 1, < (RSU ≥ Ψ , ) ∈ !(* O `QQ,(RSU = 0, ℎF)

(22)

 `QQ,(RST ≥  `QQ,(RSU 



(23)

where ηsRR is a scenario counting auxiliary variable, which is set to zero or one according to the conditions afore-described. It is assumed that the number of scenarios above (in the abscissas of Figure 4a and Figure 4b) the “real” intersection point (i.e., that which would be observed when the solution vectors are sorted) is equal for both candidate and reference solutions. The number of scenarios is evidently equal for those solutions below the “real” intersection as well. Implementing that as an equality constraint, however, would bottleneck the procedure since it would be difficult to find feasible values for Ψ.

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

The constraint concepts are illustrated in Figure 4, which aids in grasping them. In Figure 4a, we depict a situation where the number of scenarios more expensive than Ψ is greater in the reference solution. It is visually noticeable that in that situation the RRPI area is smaller or equal to the value that would be attained with the sorted vector method, leading to RAR ≥ 1. If the number of scenarios more expensive than Ψ is greater in the candidate solution, we can assume that the situation depicted in Figure 4b is occurring, and RR area value might be greater or lower than OL, and consequently RAR might be greater or lower than one. It is true that a RAR > 1 constraint could have been imposed. That was in fact implemented. However, the optimization scheme used here responded more efficiently when the afore-described method was activated, and was able to more easily find Ψ values close to the “real” intersection. For elucidative purposes, a numeric example is used. In Figure 4c the vector values depicted in Figure 3 are magnified to a narrower range. Figure 4d presents two parameters’ behaviors for evaluating constraints use: (i) the number of scenarios more expensive than Ψ in the reference solution minus that quantity in the candidate solution. Ideally, we expect that a solution is achieved where that number is zero; and (ii) the RAR value obtained for different values of Ψ (evaluated on TACs values). It is possible to observe that there are some rather small regions (in dark gray) that would satisfy the aforementioned scenarios number equality condition. One of them contains the pseudo-intersection value that is sought. Thus, an initial guess within that specific region would be advised. Otherwise, the optimization procedure might stagnate at non-optimal RAR values or not be able to even find a feasible solution. For those reasons the constraint was relaxed as follows. If the number of scenarios with a positive deviation from Ψ is greater in the candidate solution than in the reference solution, even though Ψ might not be in the real curves’ intersection, that is considered a valid solution (light gray areas in Figures 4c and Figure 4d). When that difference is negative, the RAR might be lower than one, which is an invalid situation. In that way, it is expected that the best RARPI value is going to be achieved for solutions where the aforementioned number of scenarios is equal for both curves at both sides of Ψ, and if the optimization method succeeds, when Ψ is at the

ACS Paragon Plus Environment

Page 20 of 55

Page 21 of 55 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

best pseudo-intersection value, i.e., the minimum RAR observed in Figure 4d. Note that, given the fact that data is stochastically generated, multiple intersections may be present (leading to the presence of the local minima observed in Figure 4d). Such local minima increase the complexity of the problem, which is already non-linear and non-convex. Hence, once a solution is obtained, it is important to assess its quality regarding the RAR value obtained via RARPI by manually computing RARSV. RARPI attempts to approximate RARSV, and thus some deviation is expected. However, a too large error may suggest local optimum stagnation. For that reason, an a posteriori screening of Ψ values is applicable to evaluate a solution’s quality. Moreover, with the relaxed constraint, an initial guess for Ψ is considerably easier to propose. For instance, there is a fair probability that a value around the midpoint between ETAC and the worst TAC in the reference solution is, at least, a feasible solution. Finally, the resulting optimization model is as follows:

(NN_.//) # ). .

JNN ]^ K ,6). (3) − (58) ?ã  >.A ,6). (19) − (23)

ACS Paragon Plus Environment

(24)

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

Figure 4. Possible situations with the RAR_SOO model: (a) with and (b) without imposing Eq. (23) constraint; (c) CDFs’ upper part showing the “pseudo-intersection”; (d) illustration of constraints’ effectiveness parameters In order to provide a better assessment on RAR, an additional RAR-based model is proposed. In it, one must define an acceptable RAR value a priori. The method then attempts to find solutions with

ACS Paragon Plus Environment

Page 22 of 55

Page 23 of 55 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

RAR lower or equal to that value with minimal variance. In other words, RARPI is transferred to a constraint, while VI is to be minimized. The constraints from the RAR_SOO model are maintained. In that form, VI (and consequently variance) is reduced while it is compulsory that the CDF associated to the achieved solution intersect that of the reference solution. The acceptable RAR model is described as follows:

(_NN(Φ)) # ). .

JGHK ,6). (3) − (58) ?ã  >.A ,6). (19) − (23) ]^ NN ≤ Φ

(25)

where Φ is the maximal value accepted for RARPI. 4

Solution methodology

The solution approach used is based on the SA-RFO strategy from Pavão et al.17 The methodology achieved promising solutions for large-scale HEN synthesis cases using Yee and Grossmann’s2 SWS without the isothermal assumption, as well as with the extended SWS model,7 which is used in this work. SA-RFO is a two-level optimization method. In the “upper” level, the SA strategy used has a simple move of adding a random heat exchanger at each iteration, creating a new slightly different topology. That structure is sent to RFO, which consists of a two-step continuous optimization method. The first step is a SA adaptation for continuous spaces (CSA), which, at each iteration, performs a small random move on a continuous variable (heat loads and stream split fractions). The best solution attained by CSA is included in a swarm of random solutions, which is then used by PSO to infer a promising search region. PSO then performs considerably better than with random solutions only, and it is very likely that the method will be able to find a better solution in the region of that provided by CSA. The ε-constraint based strategy described in Pavão et al.32 is also followed here for solving the ETAC_VI_MOO model (see Section 3.1). The ε-constraint method47 consists of firstly finding two anchor points (i.e., extreme solutions regarding the two objectives), then transferring one of the objectives (VI, in this study) to an auxiliary constraint. This is the so-called ε-constraint, where ε is an upper bound. Optimization runs are carried out for

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

several ε values, defined by equally dividing the space between the anchor points in the VI axis so that the desired number of Pareto solutions is obtained. In the methodology of Pavão et al.,32 when ε is changed, the topology of the best solution achieved for the previous ε is maintained as a starting point, which expedites the optimization procedure. For reaching solutions that are valid regarding the ε-constraint, the methodology proposes that the optimization is carried out first regarding the secondary objective. Once a feasible solution regarding the ε-constraint is found, then the optimization guidance is set to the primary objective. After reaching feasibility, if the meta-heuristic moves yield invalid solutions, these are not accepted by SA-RFO, which ensures that the εconstraint is not violated in final results. A more detailed description is contained in the work of Pavão et al.32 The method was implemented in C++ language in Microsoft Visual Studio 2015. A computer with an Intel® Core™ i5-4690 @ 3.50 GHz processor with 8.00 GB of RAM was used for carrying out the experiments. 5

Case study

In order to illustrate the concepts presented so far in this work, a large-scale HEN synthesis case is introduced. The process streams operating conditions were retrieved from the base case of Björk and Nordman’s48 investigation (temperatures, heat capacities and heat transfer coefficients can be found therein). This case study is a typical case from the HEN synthesis literature to which several HEN synthesis methods were applied.12,16,17 A new set of available utilities is suggested here according to the hypothetical situation that follows. Four heating utility sources are present in the industrial site where this plant is to be deployed. Two already aged boilers are available, the first being a coal-fired unit producing medium and low-pressure steam and the second burning pinewood for producing high and medium pressure steam. Two newer natural gas fired auxiliary heating sources are also present, the first one being a hot thermal oil supplier and the second one being a medium pressure steam boiler. Three types of cold utilities are present. Cooling water is available at 20 °C from cooling towers and cold water is available at 5 °C from a chiller. Air-cooling, the third cold utility, is available from air compression. The related data for these utilities can be found in

ACS Paragon Plus Environment

Page 24 of 55

Page 25 of 55 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 1. Deterministic costs and heat transfer coefficients are hypothetical values that were proposed inspired by values and proportions commonly found in the literature with some ad-hoc adjustments regarding, for instance, boiler efficiencies (e.g., older boilers’ operation is more expensive). Those monetary values reflect current operating costs per kW·y for each utility producing source. Note that, for concision, a labelling system is used for referring to each specific utility as follows (see the fourth column of Table 1). The first three characters refer to the utility producing unit (e.g., BO1 for boiler 1), the next two characters between parentheses represent the energy source (e.g., EL for electricity) and the three final characters that follow the hyphen refer to the utility produced (e.g., HPS for high pressure steam). Table 1. Available utilities set data Utility supplier

Energy

Utility

Label

source

Boiler 1

Boiler 2

Hot oil supplier Boiler 3 Air compressor Cooling tower Chiller

49

Coal

Pinewood

50

Operating

Heat

Temperature

coef.

transf.

Deterministic price ($/kWy)

2

(°C)

(kW/m K)

Medium pressure steam

BO1(CO)-MPS

200-200

1

80

Low-pressure steam

BO1(CO)-LPS

170-170

1

75

High pressure steam

BO2(PW)-HPS

255-255

1

60

Medium pressure steam

BO2(PW)-MPS

200-200

1

55

Natural Gas

51

Hot oil

HOS(NG)-OIL

330-250

0.5

65

Natural Gas

51

Medium pressure steam

BO3(NG)-MPS

200-200

1

55

52

Cooling air

COM(EL)-AIR

30-40

1

5

52

Cooling water 1

CTW(EL)-CW1

25-40

2

10

52

Cooling water 2

CHL(EL)-CW2

5-10

2

50

Electricity

Electricity

Electricity

The initial step in the methodology is to generate different forecasts for the utility prices, which is performed via MCS coupled with GBM. Utility prices from Table 1 are considered as current utility production costs (i.e., 2016, since those data were not available for 2017 by the time of the experiments). Historical data for fuels and electricity prices were obtained from different online sources free of charge. Those sources are the references that follow each of the energy sources in the second column of Table 1. It is supposed that the new plant is to be commissioned in the U.S. Hence, historical prices data (from 1997 up to 2016) for most energy sources were retrieved from

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

online databases of the United States government. However, American coal prices were not available for years prior to 2008 in such sources by the time of this investigation. Thus, Australian coal prices fluctuations, which are available for a wider time range, were used to simulate that commodity costs variation in the U.S. Although coming from different locations, prices for a single commodity were considered to typically vary in a similar manner. Since this is a conceptual case, this scenario can be considered satisfactory for research and illustration purposes. Prices for all utility sources were adjusted according to U.S. inflation. The methodology proceeds with the MCS, which is carried out for generating 1,000 costs variation scenarios for all available utilities throughout five years of plant operation. According to the GBM model, which is used with the MCS, the stochastic parameters follow a Normal distribution with higher uncertainty as we move away from the last known year. We then follow the approach by Pavao et al.32, where the multiperiod problem is transformed into a single period one by replacing the values for each year with the corresponding average along the time span. In other words, for each scenario, an average cost for that time span is calculated and stored in its respective utility index (m or n) and scenario position (s) in parameters Chum,s or Ccun,s. Figure 5 illustrates ten scenarios obtained via MCS for the costs of four of the utilities available in the case study. The figure also shows an additional axis at its right-hand side to illustrate the idea of using the average of the forecasted costs for the plant operation time. Therefore, the horizontal lines correspond to the average values of the predictions with matching the colors. These are the values that assigned to the corresponding s indexes of parameters Chum,s or Ccun,s.

ACS Paragon Plus Environment

Page 26 of 55

Page 27 of 55 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

Figure 5. Ten utility cost scenarios and their average costs obtained with Monte Carlo Simulation for four of the case study utilities Once the model parameters are known, the optimization approach is ready to be conducted. 5.1

Variance versus average TAC trade-offs

This stage of the methodology entails the evaluation of the trade-off between the expected TAC and variance. A bi-criteria optimization considering those two objectives is emulated by solving the ETAC_VI_MOO model (VI being used as a numerically simpler proxy for variance). Hence, an ETAC versus VI set of Pareto solutions is obtained via the ε-constraint based approach and variance can be manually calculated for those solutions a posteriori. Note that the term “Pareto solution” is used here for simplicity, since global optimality cannot be guaranteed. Firstly, a minimal ETAC solution must be found. The SA-RFO strategy is applied to the ETAC_SOO model, and the solution from Figure 6a is found after approximately three hours of processing time. It is worth noting that such solution uses two utilities (HPS and MPS) produced in Boiler 2, which burns pinewood for producing steam. It can be noted in Figure 5 that pinewood prices present an historical descending

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

trend. Hence, it can be expected that low ETAC solutions might apply utilities produced by using such commodity. The ETAC-optimal solution topology has great potential for heat exchange and is maintained to be used as initial topology (i.e., starting point for binary variables) for finding a solution with minimal variability index (i.e., solving VI_SOO model). The optimal configuration for VI is depicted in Figure 6b and was achieved after two hours of processing time. Henceforth, the optimal solutions found for ETAC_SOO and VI_SOO models will be referred to as ETAC-optimal and VI-optimal.

Figure 6. (a) ETAC-optimal and (b) VI-optimal designs After those two solutions are found, we consider them as anchor points in the ε-constraint approach. Eight intermediate solutions to the ETAC_VI_MOO model trading off ETAC and VI are found with the ε-constraint-based method described in Section 4, which leads to a ten-point Pareto front. This set of Pareto solutions is illustrated in Figure 7. Since VI was used as proxy for variance, the formal variance for all solutions was calculated a posteriori, so that the Pareto front representation

ACS Paragon Plus Environment

Page 28 of 55

Page 29 of 55 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

could be finally rendered in terms of ETAC and Variance (Figure 7a). Note that the distances between Pareto solutions in the abscissa are fairly uniform, and that the curve is monotonically decreasing. That indicates, as also seen in Pavão et al.,32 that the VI-based model can reasonably emulate an ETAC versus variance bi-criteria model. Design features of each solution, depicted in a double-ordinate (area and utility heat loads) chart, are presented as well (Figure 7b). Note that both graphs’ abscissas are the same and refer to the solutions’ variance. The conflictive behavior of variance and ETAC can be noticed as variance is reduced, with the increase in ETAC. That growth becomes impractical at a point. Note that the variance-optimal solution has variance 46.5% lower than the reference solution (ETAC-optimal) at expense of a 74.9% ETAC increase. From Figure 7b, it is possible to notice that in general variance is lower on solutions with less utility requirements. That is an expected behavior, since the only parameters deemed as uncertain are utility costs. However, there is an interesting feature in the solution with optimal variance. Its total utility requirements are greater than those from the next solution in the Pareto front (9265 kW/y versus 9194 kW/y for total hot utilities and 6890 kW/y versus 6819 kW/y for total cold utilities). It should be noted that since reducing variance was the only objective sought when obtaining that extreme solution (VI_SOO model), more utilities which yield lower variance values were employed, even though their total heat load and expected costs were higher. Regarding utility types, solutions with lower ETAC rely mostly on utilities from BO2, which burns pinewood, as expected given that fuel’s cost decreasing trend. Concerning cold utilities, it is worth noting air-cooling was used in the lower ETAC range, while water from cooling tower was present in solutions with lower variances. Figure 7c presents the TAC CDF curves corresponding to each ETAC_VI_MOO Pareto solution. That figure can provide a complementary analysis of the afore-presented solutions. It can be observed from the CDFs that, despite the substantial reduction in variance, the solution with optimal variance has no real advantages regarding risks. Not only its average TAC is prohibitively higher (which could already be observed in the Pareto front), but in all possible scenarios TACs is higher than in the reference CDF (ETAC-optimal). From the Pareto front, ε(1) solution can draw some

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

attention with a considerable variance reduction of 41.9% and an ETAC increase of only 13.4%, which is a better trade-off than that achieved in the Variance-optimal solution and might suggest that it constitutes an interesting choice. Extending the analysis to the CDFs, it is possible to notice in Figure 7c, and in details in Figure 7d, that the ε(1) solution CDF does not intersect the reference CDF at any point. That means that, as in the Variance-optimal solution, the solution for the smallest ε brings no risk reduction at any of the scenarios evaluated. Other solutions’ curves indeed intersect the reference curve, which is advantageous. However, it should be noted also that below the intersections, scenarios’ TAC increase. For a more complete discussion, the investigation is extended with aid of the Downside risk and the Risk area ratio metrics.

Figure 7. (a) ETAC_VI_MOO Pareto solutions projected in ETAC versus Variance plot; (b) design aspects of each solution in the Pareto front; (c) TAC CDFs associated to those solutions; (d) zoom over the upper area of the TAC CDFs

ACS Paragon Plus Environment

Page 30 of 55

Page 31 of 55 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

5.2

Risk management with Downside Risk and Risk Area Ratio

In order to better evaluate the financial risk of the solutions obtained with the ETAC_VI_MOO model and to put the potential of the models from Sections 3.2 and 3.3 into perspective, we firstly assess the Pareto solutions obtained in Section 5.1 regarding their Downside Risk (Ω = 1.25M$/y, 1.30M$/y, 1.35M$/y and 1.4M$/y) and Risk Area Ratio. Those results are presented in Table 2. Table 2. DR and RAR for the optimal solutions obtaining with the ETAC_VI_MOO model (best results in bold) Var.-opt.

ε(1)

ε(2)

ε(3)

ε(4)

ε(5)

ε(6)

ε(7)

ε(8)

ETAC-opt.

Variance (109$2/y2)

6.34

6.90

7.48

8.03

8.62

9.24

9.92

10.49

11.17

11.87

ETAC ($/y)

2,100,534

1,361,319

1,323,956

1,298,581

1,264,876

1,233,047

1,221,940

1,211,221

1,204,738

1,200,959

DR(1.25M$/y)

850.53

113.16

80.90

62.98

43.37

29.87

27.25

25.09

24.30

24.35

800.53

70.32

45.95

33.99

22.25

15.12

13.99

12.99

12.76

13.01

750.53

37.76

23.01

16.61

10.67

7.26

6.87

6.47

6.50

6.76

700.53

17.74

10.57

7.61

4.97

3.45

3.40

3.38

3.47

3.68

RARSV

-

-

18137.0

1581.0

433.1

89.0

64.0

34.2

14.9

-

RRSV ($/y)

0.00

0.00

6.78

61.79

147.91

364.43

333.00

309.45

272.19

-

OLSV ($/y)

899,574

160,359

123,003

97,684

64,064

32,452

21,313

10,571

4,051

-

(k$/y) DR(1.30M$/y) (k$/y) DR(1.35M$/y) (k$/y) DR(1.40M$/y) (k$/y)

Under the DR metric, it can be noted that ε(8) solution performed better for the two lower targets (1.25 and 1.30 M$/y), while the best results for the two more conservative targets (1.35 and 1.4 M$/y) were achieved with the ε(7) solution. We next aim to solve the ETAC versus DR model for different Ω targets, where we will use the information obtained from Table 2 in order to facilitate the optimization procedure. Specifically, the ε(8) solution is used as initial solution for problems DR_SOO(1.25 M$/y) and DR_SOO(1.30 M$/y), while ε(7) solution is employed on DR_SOO(1.35 M$/y) and DR_SOO(1.40 M$/y). TAC

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

CDF curves associated to each of those solutions are presented in Figure 8, where the larger plots regard the upper part of the curves, which is the region where the DR_SOO is driven by the reduction in scenarios’ costs, and the smaller plots present an overview of the entire CDFs. The influence of the target (Ω) choice can be clearly inferred from the curves presented, where it can be noted that the higher Ω, the lower the variance. Solutions for higher targets had their costs reduced mostly for the most expensive scenarios (i.e., the higher parts of the curve), which means that choosing higher targets might lead to solutions appealing for risk-averse investors. By comparing the solutions achieved with the DR_SOO model to those from Table 2, it could be noted that DR improved in all cases, which is not surprising giving that it became the optimized objective in DR_SOO. Regarding solutions’ design features, it is possible to note a behavior similar to that found in the variance optimization. As DR targets increase, DR and variance become lower, reducing as well utilities use and increasing in turn the heat exchanger areas. It could also be observed that two types of hot utilities (BO2(PW)-MPS and BO3(NG)-HPS) are used in all solutions. However, the employment of BO2(PW)-MPS in the DR(1.4 M$/y) optimal solution contributed to lowering the downside risk.

ACS Paragon Plus Environment

Page 32 of 55

Page 33 of 55 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

Figure 8. Solutions obtained with DR_SOO model for different targets

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

Figure 9. DR_SOO(1.3 M$/y) solution Although DR is an interesting metric for reducing risks of expensive scenarios, a drawback is that the opportunity loss (i.e., scenarios that become more expensive in comparison to the reference solution) is not taken into account while assessing a candidate solution. When observing the lower part of the CDFs in Figure 8, it becomes evident that more scenarios become more expensive as more conservative targets are chosen. Figure 9 presents the solution to DR_SOO(1.3 M$/y). Despite employing hot utilities from two boilers, this configuration uses two less units than the ETACoptimal solution (four for hot plus two for cold vs five plus three in ETAC-optimal, see Figure 6a), constituting an attractive option. Moreover, among solutions obtained to all Ω targets using DR_SOO, the solution to the DR_SOO(1.3 M$/y) yields the greatest absolute downside risk reduction (9.47×102 $/y) in comparison with the ETAC-optimal solution evaluated with the same target. As an interesting alternative, the RAR metric can be used for tackling the issue that a solution may become too expensive in several scenarios in conservative solutions, providing a better risk evaluation. The method directly compares candidate solutions to a reference one (in this study, ETAC-optimal) regarding the trade-off between risks reduction and opportunities loss (see Section 3.3). At a first attempt of applying the metric, the RAR_SOO model was employed, being ε(8) the initial solution used. The solution with the best RARPI that could be achieved with the RAR_SOO model had an objective function value of 1.08. When compared to the value obtained with Eqs. (9)-

ACS Paragon Plus Environment

Page 34 of 55

Page 35 of 55 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

(13), which is the calculation form that does not require intersection values, the difference between them is in the order of 10-7. Such small deviation suggests that using the pseudo-intersection concept in the optimization models proposed here yields a reasonable approximation for RARSV. The aforementioned solution, though, was practically the same as the ETAC-optimal: ETAC was 1.5 $/y higher, with nearly negligible variations in the design. That configuration was thus discarded, and the ACC_RAR model, which also uses the pseudo-intersection concept, was used for achieving better possible solutions at pre-determined RARPI upper boundary values. The values used were 2.0, 5.0 and 10.0. Figure 10 depicts the TAC CDF curves of each solution, which analogous format with those in Figure 8. The value found for the pseudo-intersection (ψ) is also represented in each graph.

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

Figure 10. Solutions obtained with ACC_RAR model for different targets The solutions achieved for ACC_RAR(2.0) and ACC_RAR(5.0) are structurally very similar to the ETAC-optimal ones. For comparison, see Figure 11, which depicts the ACC_RAR(5.0) solution, and Figure 6a. The topologies are equivalent, and marginal changes in heat loads yield slightly lower utility requirements, reducing also the variance. Moreover, regarding risks, ACC_RAR(5.0) solution can be considered a middle-ground solution among those obtained with ACC_RAR.

ACS Paragon Plus Environment

Page 36 of 55

Page 37 of 55 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

ACC_RAR(2.0) configuration is the closest to the risk-takers’ option (ETAC-optimal), while ACC_RAR(10.0) solution presents the higher risk reduction, but also the greatest opportunity loss, constituting a conservative alternative.

Figure 11. ACC_RAR(5.0) solution For ACC_RAR(10.0), BO3(NG)-MPS was used, leading to lower variance but higher ETAC. The pseudo-intersection based calculation performed well as a RARSV approximation. The greatest deviation (-0.6) occurred with the target of RARPI = 10.0. We discarded the hypothesis of local minimum stagnation by manually screening Ψ values, similarly to Figure 4d, and observing that the value found for Ψ led to the lowest RARPI. Thus, inaccuracy is due to the existence of multiple intersections. In order to evaluate RR versus OL trade-offs obtained using the DR metric model, solutions obtained with that approach were evaluated under the RAR concept as well. Table 3 presents those results. Table 3. Comparison of solutions to RAR and DR models regarding risk reduction, opportunity loss and RAR (best values in bold) RAR best

RAR (2.0)

RAR (5.0)

RAR (10.0)

DR

(1.25

DR

(1.3

DR(1.35

DR

M$/y)

M$/y)

M$/y)

M$/y)

Variance (109$2/y2)

11.85

11.67

11.43

10.24

10.42

9.82

9.71

9.56

ETAC ($/y)

1,200,961

1,201,211

1,202,312

1,209,546

1,208,063

1,212,817

1,213,803

1,224,193

RARSV

1.1

2.0

5.0

9.4

7.8

12.2

13.2

50.1

RRSV ($/y)

20

252

338

1,028

1,043

1,061

1,054

474

ACS Paragon Plus Environment

(1.4

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

OLSV ($/y)

21

503

1,691

9,615

8,147

Page 38 of 55

12,918

13,897

23,707

It can be observed that, with both metrics, as risks are reduced (RR), opportunity loss (OL) increases at a greater proportion, which justifies the application of a metric such as RAR for assessing those trade-offs. However, it can be noted in Table 2 and in the last three columns of Table 3 that, after a point, as variance decreases, RR also decreases. In nearly optimal variance values, there is no advantage at all (as observed for ε(1) and Variance-optimal solutions in Table 2). It is notable that RAR can be used as an efficient form of evaluating solutions that present some advantage regarding risks reduction, or in simpler words, solutions whose TAC CDF intersect a reference solution’s one. Including the RAR into an optimization model can efficiently lead to more promising solutions among those advantageous solutions, and the pseudo-intersection RAR model is an efficient way of doing so. If such model is not at hand, a target-based metric such as DR can be applied with different targets and evaluated regarding RAR a posteriori as well, providing interesting results. 5.3

Utility source flexibility study

The case presented so far considers fixed sources for utilities in a given solution. For instance, if BO1(CO)-MPS is selected for use in a heater, that unit uses MPS only from BO1 during the plant operation. However, there are three different sources for medium pressure steam (MPS): boilers 1, 2 and three, each burning a different fuel. Therefore, assuming continuous operation of the three boilers and in the absence of capacity or structure limitations in the plant, additional piping and valves could be implemented so as to ensure that the cheaper source is always used. Since piping costs are neglected in this work as we focus on the HEN, the structural costs related to the flexible configuration are also left out of the scope. Those costs must, evidently, be evaluated at the plant design stage. In order to explore such flexible operation scheme, we compare the costs of MPS production among the three sources for the 1,000 scenarios along the five years of plant lifetime,

ACS Paragon Plus Environment

Page 39 of 55 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

selecting and storing always the lowest MPS cost for each year. That is, the three MPS sources are lumped into a single MPS source, thus simulating the situation of using always the cheapest source. This affects the m dimension of the Chum,s parameter, which now presents only a general MPS utility, originated always at the cheapest source, instead of three MPS sources as in the original parameter. The probabilities of sources yielding the cheapest MPS costs for a given year, according to a Monte Carlo simulation, are 7.1% for BO1(CO), 44.4% for BO2(PW) and 48.5% for BO3(NG). With the yearly values for the lumped MPS source at hand, the mean value for the five years is then computed, in the same spirit as before (see Figure 5, for instance). With the new MPS utility source matrix, the ETAC_SOO, DR_SOO, RAR_SOO and ACC_RAR models are solved again. Results are presented in Table 4. Table 4. Comparison of solutions to ETAC_SOO, DR_SOO, RAR_SOO and ACC_RAR models for the source flexibility study regarding RR, OL and RAR MP

ETAC best

RAR best

RAR (1.5)

RAR (2.0)

RAR (5.0)

RAR (10.0)

DR (1.25 DR M$/y) M$/y)

Variance (109$2/y2)

10.29

10.27

7.10

6.74

6.23

5.33

6.21

5.91

5.58

5.66

ETAC ($/y)

1,193,915

1,193,920

1,196,720

1,199,453

1,209,895

1,224,068

1,199,992

1,204,127

1,210,435

1,210,790

-

1.1

1.5

2.0

5.0

9.9

1.9

2.7

4.0

4.2

RAR RR

SV

SV

(1.3 DR(1.35 M$/y)

DR M$/y)

($/y)

-

46

5,610

5,541

4,001

3,370

6,541

6,110

5,421

5,215

OLSV ($/y)

-

51

8,415

11,078

19,980

33,523

12,618

16,322

21,940

22,090

As expected, ETAC, as well as variance values, are slightly reduced on average in comparison to those obtained with fixed utility sources. The RAR values achieved with the DR model with different Ω values are considerably lower with the flexible condition. As in the base case, the solution to RAR_SOO is of little use, since it is again only marginally different from the ETAC_SOO solution. The DR_SOO(1.25M$/y) model yields a RAR of 1.9, which is the best RAR value for this case found without using a RAR-based model. For that reason, the ACC_RAR model was also solved with maximal value of 1.5, demonstrating that it can outperform the DR_SOO(1.25M$/y) model results regarding the risk reduction versus opportunity loss trade-off

ACS Paragon Plus Environment

(1.4

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

(see second column in Table 4). Another interesting aspect of the flexible case is that, in absolute values, risk reduction is considerably greater than in the base case. This fact can be better observed in Figure 12, which presents the best two solutions regarding RAR (DR_SOO(1.25M$/y) and ACC_RAR(1.5)). In such illustration, values where the areas enclosed by the curves are more noticeable even by visual inspection. The figure also presents the optimal design for the ACC_RAR(1.5) model. Note that there are sequential heaters placed at cold stream 7, a configuration that can only be found thanks to the employed enhanced stage-wise superstructure.7 H5 uses MPS, which operates at 200 ºC. This imposes a thermodynamic constraint to that unit, which is not able to heat cold stream 7 up to its target temperature. A heater using HPS is used in subsequently. If the MPS heater were removed and the HPS heater were used alone in the same position performing the whole heating service, RAR would be of 11.0, ETAC would increase by 1,258 $/y and variance by 35%. Note that the remaining required heat is received further in the cold stream via heat exchange with hot stream 2 via HE2. On the other hand, if the MPS heater were again removed, but the HPS heater were placed at the end of the cold stream (instead of in the same position), RAR would be 136 (nearly no risk reduction), ETAC would raise by 7,639 $/yr and variance by 35%. Placing the MPS/HPS heaters sequence after the heat exchanger yields an infeasible solution. Maintaining the MPS heater at its position and placing the HPS one at the end of the stream yields the same variance, but RAR of 3.9 and ETAC 7,525 $/yr higher.

ACS Paragon Plus Environment

Page 40 of 55

Page 41 of 55 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

Figure 12. (a) Flexible design achieved with ACC_RAR(1.5) model and (b) TAC CDFs for DR(Ω=1.25M$/y) and ACC_RAR(1.5) 5.4

Benchmarking

In order to put the advances presented in this work into perspective, the case study from Pavão et al.32 is here revisited so that the new methodology can be applied and the results compared to those presented therein. The case is now analyzed with the use of the enhanced SWS,7 which enables new structural options and thus it is expected that better solutions could be found. Moreover, the RAR metric is applied, making it possible to compare the trade-offs between opportunity loss and risk reduction among the solutions reported in the aforementioned work and those identified here. Firstly, the best solution found here for ETAC_SOO has lower ETAC (2,923,070 versus 2,930,842 $/yr). An analysis is then performed using the RAR perspective to evaluate the solutions reported in Pavão et al.32 Taking as reference solution the one with the best ETAC there reported, RAR was computed to all other solutions identified in that work. The one with the best RAR (2.43) was the solution to the TETAC_VAR model (presented therein). That model consists of fixing an ETAC

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

value (2.95M$/y in that case) and finding the minimal variance solution. As previously stated, the ETAC-optimal solution found in this work is better than that from the referred work. Thus, it is reasonable to use it as reference solution for RAR. In that case, the solution to the TETAC_VAR model presented by Pavão et al.32 yields RAR = 2.08. We thus applied the ACC_RAR(1.5) model in order to find a solution with lower RAR value. Figure 13 (a) presents the optimal design. The CDF for ETAC-optimal solutions from the present work and from the literature are compared in Figure 13 (b). Figure 13(c) presents the CDF for the ETAC-optimal solution from this work (used as reference), for the solution here obtained with the ACC_RAR(1.5) model and for the TETAC_VAR model obtained in the literature32. It is worth noting that a heater is placed at a stream split branch in the presented design, which means that the new design possibilities led to a broader HEN evaluation regarding their risks. Moreover, a better risk reduction versus opportunity loss trade-off evaluation could be achieved with the proposed ACC_RAR model.

ACS Paragon Plus Environment

Page 42 of 55

Page 43 of 55 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

Figure 13. (a) ACC_RAR(1.5) design; (b) CDFs for ETAC-optimal solution from this work and from Ref. 32; (c) CDFs for ETAC-optimal solution from this work, for the solution with best RAR from Ref. 32 and for the ACC_RAR(1.5) solution

6

Conclusions

A methodology for analyzing and reducing financial risk in heat exchanger networks synthesis comprising multiple utilities from different sources was presented. Optimization models based on three risk metrics (VI, DR and RAR) were developed. A novel “pseudo-intersection” calculation

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

method for RAR enabled the RAR metric to be efficiently applied in an optimization model at low computational burden. A large-scale case study comprising several utilities and energy sources was investigated. The VI bi-criteria model provided an interesting evaluation on variance versus average TAC. For instance, an appealing solution entailing a variance reduction of 41.9% and an average cost increase of only 13.4% (in relation to the ETAC-optimal solution) was identified with ε(1) (see Section 5.1). However, it could be observed in TAC CDFs of low-variance configurations, that reducing such metric (which is performed by means of the VI in this work) presents advantages only up to a certain point. Low-variance solutions, such as the aforementioned ε(1) one, presented few or no advantageous scenarios at all. The downside risk based model demonstrated to be able to provide more advantageous solutions regarding general risks reduction in relation to the ETAC-optimal solution. However, such metric does not assess how solutions that had their risks reduced in some scenarios may become more expensive in others. Conversely, the RAR metric is able to capture such trade-off, and therefore we advocate to use such metric for a more complete financial risk evaluation. When calculating RAR for solutions obtained with other metrics, it could be observed that, although some might yield high risk reduction, the opportunity loss was also large (for instance, DR_SOO(1.3 M$/y) solution has the higher RR, but OL is 12.2 times greater). The acceptable RAR (ACC_RAR) model was able to efficiently lead to solutions with pre-determined RARPI values. With it, an investor can directly set how conservative a solution can be regarding the risk reduction versus opportunity loss trade-offs and achieve a solution with the desired profile. The pseudointersection RAR (RARPI) calculation demonstrated to be a reliable approximation to the RAR calculated with the sorted vectors method (RARSV). That metric can simultaneously evaluate risks reduction and opportunity loss. From the design point of view, it was observed that low-variance solutions had lower utility requirements, as all uncertain parameters are utility-related. It is interesting to note that in lower ETAC solutions, most of the hot utilities come from Boiler 2, which burns pinewood. A historical

ACS Paragon Plus Environment

Page 44 of 55

Page 45 of 55 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

decreasing trend can be observed for that commodity’s price in Figure 5. In solutions with reduced variance, though, other utilities are used. For cold utilities, most advantageous solutions used aircooling. Water from cooling towers were present only in low-variance solutions. In order to consider that different sources might be used for producing the same utility during the plant operation and that the cheapest one would always be selected, a flexibility study was also performed. Specifically, we consider that the less costly MPS was always used. As expected, ETAC of solutions to all studied metrics was slightly lower than the ETAC found to these same metrics in the fixed utility sources case. Moreover, it was observed that, for similar values of RAR, the risk reduction values were higher in the flexibility study than in the fixed-sources case. For instance, for RAR = 5, RR = 4,001 $/y for the flexible case and RR = 338 $/y for the fixed-sources case. The techniques presented in this work were applied to the case study evaluated in a previous investigation. The present framework achieved a better solution for risk-takers (lower ETC value in the ETAC-optimal solution) as well as a better-balanced solution regarding the trade-off between risk reduction and opportunity loss. Finally, it can be concluded that the inclusion of risk metrics to an enhanced HEN SWS-based superstructure model solved with an efficient meta-heuristic optimization approach can yield more robust HEN designs regarding the financial risk to be taken. Future works might extend this investigations by evaluating, for instance, the retrofit of HEN in plants whose utilities had already became too expensive. Uncertain utility costs may be considered based on different utility costs forecasts when performing such retrofitting task. 7

Acknowledgements

L.V. Pavão, C.B.B. Costa and M.A.S.S. Ravagnani gratefully acknowledge the financial support from the Coordination for the Improvement of Higher Education Personnel (CAPES, Brazil) and the National Council for Scientific and Technological Development (CNPq, Brazil). L. Jiménez acknowledges the financial support from Spanish Ministry of Education and Competitiveness (CTQ2016-77968, MINECO/FEDER).

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

8

Nomenclature

Models Single-objective optimization for variability index at an acceptable risk area

ACC_RAR

ratio DR_SOO

Single objective optimization model for downside risks

ETAC_SOO

Single objective optimization model for expected total annual costs

ETAC_VI_MOO

Bi-criteria optimization model for expected total annual costs and variability index

RAR_SOO

Single objective optimization model for risk area ratio

VI_SOO

Single objective optimization model for variability index

Parameters B

[$/y]

Equipment fixed costs

C

[$/(m2βy)]

Annualized capital cost factor

Ccu

[$/kWy]

Cold utility costs

Chu

[$/kWy]

Hot utility costs

prob

[-]

Probability of occurrence of a given scenario

β

[-]

Capital cost exponent



[$/y]

Target costs

Variables A

[m2]

Heat exchanger area

Acuinter

[m2]

Intermediate stage heater area

Acuoutlet

[m2]

Outlet stage cooler area

Ahuinter

[m2]

Intermediate stage cooler area

Ahuoutlet

[m2]

Intermediate stage cooler area

D

[$/y]

Positive deviation from ETAC in variability index

ACS Paragon Plus Environment

Page 46 of 55

Page 47 of 55 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

ETAC

[$/y]

FcQhuoutlet [-]

Expected total annual costs Fraction of the required heat in outlet stage of a cold stream provided with a given hot utility

FhQcuoutlet [-]

Fraction of the available heat in outlet stage of a hot stream exhausted with a given cold utility

OLPI

[$/y]

Opportunity loss area calculated with the pseudo-intersection method

OLSV

[$/y]

Opportunity loss calculated with the “sorted vector” method

Qcuinter

kW

Heat load of an intermediate stage cooler

Qcuoutlet

kW

Available heat in outlet stage of a hot stream

Qhuinter

kW

Heat load of an intermediate stage heater

Qhuoutlet

kW

Required heat in outlet stage of a cold stream

RARPI

[-]

Risk area ratio calculated with the “pseudo-intersection” method

RARSV

[-]

Risk area ratio calculated with the “sorted vector” method

RRPI

[$/y]

Risk reduction area calculated with the “pseudo-intersection” method

RRSV

[$/y]

Risk reduction calculated with the “sorted vector” method

TAC

[$/y]

Total annual costs for a scenario

VI

[$/y]

Variability index

z

[-]

Heat exchanger existence binary variable

zcuinter

[-]

Intermediate stage cooler existence binary variable

zcuoutlet

[-]

Outlet stage cooler existence binary variable

zhuinter

[-]

Intermediate stage heater existence binary variable

zhuoutlet

[-]

Outlet stage heater existence binary variable

γOL

[$/y]

Auxiliary variable for calculating OLPI

γRR

[$/y]

Auxiliary variable for calculating RRPI

δ

[$/y]

Positive deviation from Ω in downside risk metric

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

ηRR

[-]

Scenarios counting auxiliary variable

ξOL

[$/y]

Auxiliary variable for calculating OL

ξRR

[$/y]

Auxiliary variable for calculating RR

Φ

[-]

Maximal value accepted for RARPI

Ψ

[$/y]

Pseudo-intersection

Indexes i

Hot stream

j

Cold stream

k

Stage

m

Hot utility

n

Cold utility

s

Scenario

Sets NC

Cold streams

NCU

Cold utilities

NH

Hot streams

NHU

Hot utilities

NS

Stages

NSc

Scenarios

9 (1)

References Floudas, C. A.; Ciric, A. R. Strategies for Overcoming Uncertainties in Heat Exchanger Network Synthesis. Comput. Chem. Eng. 1989, 13 (10), 1133–1152.

(2)

Yee, T. F.; Grossmann, I. E. Simultaneous Optimization Models for Heat integration—II. Heat Exchanger Network Synthesis. Comput. Chem. Eng. 1990, 14 (10), 1165–1184.

(3)

Ponce-Ortega, J. M.; Serna-González, M.; Jiménez-Gutiérrez, A. Synthesis of Heat

ACS Paragon Plus Environment

Page 48 of 55

Page 49 of 55 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

Exchanger Networks with Optimal Placement of Multiple Utilities. Ind. Eng. Chem. Res. 2010, 49 (6), 2849–2856. (4)

Huo, Z.; Zhao, L.; Yin, H.; Ye, J. Simultaneous Synthesis of Structural-Constrained Heat Exchanger Networks with and without Stream Splits. Can. J. Chem. Eng. 2013, 91 (5), 830– 842.

(5)

Huang, K. F.; Karimi, I. A. Simultaneous Synthesis Approaches for Cost-Effective Heat Exchanger Networks. Chem. Eng. Sci. 2013, 98, 231–245.

(6)

Kim, S. Y.; Jongsuwat, P.; Suriyapraphadilok, U.; Bagajewicz, M. Global Optimization of Heat Exchanger Networks. Part 1: Stages/Substages Superstructure. Ind. Eng. Chem. Res. 2017, 56 (20), 5944–5957.

(7)

Pavão, L. V.; Costa, C. B. B.; Ravagnani, M. A. S. S. An Enhanced Stage-Wise Superstructure for Heat Exchanger Networks Synthesis with New Options for Heaters and Coolers Placement. Ind. Eng. Chem. Res. 2018, 57 (7), 2560–2573.

(8)

Lewin, D. R. A Generalized Method for HEN Synthesis Using Stochastic Optimization - II. The Synthesis of Cost-Optimal Networks. Comput. Chem. Eng. 1998, 22 (10), 1387–1405.

(9)

Lewin, D. R.; Wang, H.; Shalev, O. A Generalized Method for HEN Synthesis Using Stochastic Optimization – I. General Framework and MER Optimal Synthesis. Comput. Chem. Eng. 1998, 22 (10), 1503–1513.

(10)

Ravagnani, M. A. S. S.; Silva, A. P.; Arroyo, P. A.; Constantino, A. A. Heat Exchanger Network Synthesis and Optimisation Using Genetic Algorithm. Appl. Therm. Eng. 2005, 25 (7), 1003–1017.

(11)

Luo, X.; Wen, Q.-Y.; Fieg, G. A Hybrid Genetic Algorithm for Synthesis of Heat Exchanger Networks. Comput. Chem. Eng. 2009, 33 (6), 1169–1181.

(12)

Fieg, G.; Luo, X.; Jeżowski, J. A Monogenetic Algorithm for Optimal Design of LargeScale Heat Exchanger Networks. Chem. Eng. Process. Process Intensif. 2009, 48 (11–12), 1506–1516.

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

(13)

Silva, A. P.; Ravagnani, M. A. S. S.; Biscaia, E. C.; Caballero, J. A. Optimal Heat Exchanger Network Synthesis Using Particle Swarm Optimization. Optim. Eng. 2010, 11 (3), 459–470.

(14)

Pavão, L. V.; Costa, C. B. B.; Ravagnani, M. A. S. S. Automated Heat Exchanger Network Synthesis by Using Hybrid Natural Algorithms and Parallel Processing. Comput. Chem. Eng. 2016, 94, 370–386.

(15)

Zhang, H.; Cui, G.; Xiao, Y.; Chen, J. A Novel Simultaneous Optimization Model with Efficient Stream Arrangement for Heat Exchanger Network Synthesis. Appl. Therm. Eng. 2017, 110, 1659–1673.

(16)

Pavão, L. V.; Costa, C. B. B.; Ravagnani, M. A. S. S. Heat Exchanger Network Synthesis without Stream Splits Using Parallelized and Simplified Simulated Annealing and Particle Swarm Optimization. Chem. Eng. Sci. 2017, 158, 96–107.

(17)

Pavão, L. V.; Costa, C. B. B.; Ravagnani, M. A. S. S.; Jiménez, L. Large-Scale Heat Exchanger Networks Synthesis Using Simulated Annealing and the Novel Rocket Fireworks Optimization. AIChE J. 2017, 63 (5), 1582–1601.

(18)

Toimil, D.; Gómez, A. Review of Metaheuristics Applied to Heat Exchanger Network Design. Int. Trans. Oper. Res. 2017, 24 (1–2), 7–26.

(19)

Aaltola, J. Simultaneous Synthesis of Flexible Heat Exchanger Network. Appl. Therm. Eng. 2002, 22 (8), 907–918.

(20)

Verheyen, W.; Zhang, N. Design of Flexible Heat Exchanger Network for Multi-Period Operation. Chem. Eng. Sci. 2006, 61 (23), 7730–7753.

(21)

Pavão, L. V.; Miranda, C. B.; Costa, C. B. B.; Ravagnani, M. A. S. S. Efficient Multiperiod Heat Exchanger Network Synthesis Using a Meta-Heuristic Approach. Energy 2018, 142, 356–372.

(22)

Jiang, D.; Chang, C. T. A New Approach to Generate Flexible Multiperiod Heat Exchanger Network Designs with Timesharing Mechanisms. Ind. Eng. Chem. Res. 2013, 52 (10), 3794–

ACS Paragon Plus Environment

Page 50 of 55

Page 51 of 55 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

3804. (23)

Miranda, C. B.; Costa, C. B. B.; Caballero, J. A.; Ravagnani, M. A. S. S. Heat Exchanger Network Optimization for Multiple Period Operations. Ind. Eng. Chem. Res. 2016, 55 (39), 10301–10315.

(24)

Pavão, L. V.; Miranda, C. B.; Costa, C. B. B.; A.S.S. Ravagnani, M. Synthesis of Multiperiod Heat Exchanger Networks with Timesharing Mechanisms Using MetaHeuristics. Appl. Therm. Eng. 2018, 128, 637–652.

(25)

Mou, C.; Qvale, B. The Consequences of Unpredictable Development of Economic Conditions on Heat Exchanger Network Configurations and Economic Results. Energy Convers. Manag. 2002, 43 (9), 1549–1562.

(26)

Nemet, A.; Klemeš, J. J.; Kravanja, Z. Minimisation of a Heat Exchanger Networks’ Cost over Its Lifetime. Energy 2012, 45 (1), 264–276.

(27)

Nemet, A.; Klemeš, J. J.; Kravanja, Z. Optimising Entire Lifetime Economy of Heat Exchanger Networks. Energy 2013, 57, 222–235.

(28)

Novak Pintarič, Z.; Kravanja, Z. A Methodology for the Synthesis of Heat Exchanger Networks Having Large Numbers of Uncertain Parameters. Energy 2015, 92, 373–382.

(29)

Pavão, L. V.; Pozo, C.; Costa, C. B. B.; Ravagnani, M. A. S. S.; Jiménez, L. A MetaHeuristic Approach for Financial Risks Management in Heat Exchanger Networks. Comput. Aided Chem. Eng. 2017, 40, 955–960.

(30)

Hull, J. C. Options, Futures and Other Derivatives, 8th ed.; Prentice Hall, 2012.

(31)

Eppen, G. D.; Martin, R. K.; Schrage, L. OR Practice—A Scenario Approach to Capacity Planning. Oper. Res. 1989, 37 (4), 517–527.

(32)

Pavão, L. V.; Pozo, C.; Costa, C. B. B.; Ravagnani, M. A. S. S.; Jiménez, L. Financial Risks Management of Heat Exchanger Networks under Uncertain Utility Costs via MultiObjective Optimization. Energy 2017, 139, 98–117.

(33)

Ahmed, S.; Sahinidis, N. V. Robust Process Planning under Uncertainty. Ind. Eng. Chem.

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

Res. 1998, 37 (5), 1883–1892. (34)

Barbaro, A.; Bagajewicz, M. J. Managing Financial Risk in Planning under Uncertainty. AIChE J. 2004, 50 (5), 963–989.

(35)

Sahinidis, N. V. Optimization under Uncertainty: State-of-the-Art and Opportunities. Comput. Chem. Eng. 2004, 28 (6–7), 971–983.

(36)

Aseeri, A.; Bagajewicz, M. J. New Measures and Procedures to Manage Financial Risk with Applications to the Planning of Gas Commercialization in Asia. Comput. Chem. Eng. 2004, 28 (12), 2791–2821.

(37)

Bojarski, A. D.; Guillén-Gosálbez, G.; Jiménez, L.; Espuña, A.; Puigjaner, L. Life Cycle Assessment Coupled with Process Simulation under Uncertainty for Reduced Environmental Impact: Application to Phosphoric Acid Production. Ind. Eng. Chem. Res. 2008, 47 (21), 8286–8300.

(38)

Gebreslassie, B. H.; Guillén-Gosálbez, G.; Jiménez, L.; Boer, D. Economic Performance Optimization of an Absorption Cooling System under Uncertainty. Appl. Therm. Eng. 2009, 29 (17), 3491–3500.

(39)

You, F.; Wassick, J. M.; Grossmann, I. E. Risk Management for a Global Supply Chain Planning under Uncertainty: Models and Algorithms. AIChE J. 2009, 55 (4), 931–946.

(40)

Sabio, N.; Gadalla, M.; Guillén-Gosálbez, G.; Jiménez, L. Strategic Planning with Risk Control of Hydrogen Supply Chains for Vehicle Use under Uncertainty in Operating Costs: A Case Study of Spain. Int. J. Hydrogen Energy 2010, 35 (13), 6836–6852.

(41)

Kostin, A. M.; Guillén-Gosálbez, G.; Mele, F. D.; Bagajewicz, M. J.; Jiménez, L. Design and Planning of Infrastructures for Bioethanol and Sugar Production under Demand Uncertainty. Chem. Eng. Res. Des. 2012, 90 (3), 359–376.

(42)

Ruiz-Femenia, R.; Guillén-Gosálbez, G.; Jiménez, L.; Caballero, J. A. Multi-Objective Optimization of Environmentally Conscious Chemical Supply Chains under Demand Uncertainty. Chem. Eng. Sci. 2013, 95, 1–11.

ACS Paragon Plus Environment

Page 52 of 55

Page 53 of 55 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

(43)

Sabio, N.; Pozo, C.; Guillén-Gosálbez, G.; Jiménez, L.; Karuppiah, R.; Vasudevan, V.; Sawaya, N.; Farrell, J. T. Multiobjective Optimization under Uncertainty of the Economic and Life-Cycle Environmental Performance of Industrial Processes. AIChE J. 2014, 60 (6), 2098–2121.

(44)

Santibañez-Aguilar, J. E.; Guillen-Gosálbez, G.; Morales-Rodriguez, R.; Jiménez-Esteller, L.; Castro-Montoya, A. J.; Ponce-Ortega, J. M. Financial Risk Assessment and Optimal Planning of Biofuels Supply Chains under Uncertainty. BioEnergy Res. 2016, 1–17.

(45)

Ewertowska, A.; Pozo, C.; Gavaldà, J.; Jiménez, L.; Guillén-Gosálbez, G. Combined Use of Life Cycle Assessment, Data Envelopment Analysis and Monte Carlo Simulation for Quantifying Environmental Efficiencies under Uncertainty. J. Clean. Prod. 2017, 166, 771– 783.

(46)

Galán-Martín, A.; Pozo, C.; Azapagic, A.; Grossmann, I. E.; Mac Dowell, N.; GuillénGosálbez, G. Time for Global Action: An Optimised Cooperative Approach towards Effective Climate Change Mitigation. Energy Environ. Sci. 2018, 11, 572–581.

(47)

Ehrgott, M. Multicriteria Optimization, 2nd ed.; Springer Berlin Heidelberg, 2005.

(48)

Björk, K.-M.; Nordman, R. Solving Large-Scale Retrofit Heat Exchanger Network Synthesis Problems with Mathematical Optimization Methods. Chem. Eng. Process. Process Intensif. 2005, 44 (8), 869–876.

(49)

Index

Mundi

-

Coal,

Australian

thermal

coal

Monthly

Price

http://www.indexmundi.com/commodities/?commodity=coal-australian&months=360 (accessed Nov 1, 2017). (50)

Texas

A&M

Forest

Service

-

Landowner

Assistance:

Timber

Price

Trends

http://tfsweb.tamu.edu/timberpricetrends/ (accessed Nov 1, 2017). (51)

U.S. Energy Information Administration - Henry Hub Natural Gas Spot Price https://www.eia.gov/dnav/ng/hist/rngwhhdd.htm (accessed Nov 1, 2017).

(52)

U.S.

Energy

Information

Administration

ACS Paragon Plus Environment

-

Electricity

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

https://www.eia.gov/electricity/data.php#sales (accessed Nov 1, 2017).

ACS Paragon Plus Environment

Page 54 of 55

Page 55 of 55 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

For Table of Contents Only

ACS Paragon Plus Environment