Deterministic global optimization of steam cycles using the IAPWS-IF97 model

The IAPWS-IF97 (Wagner et al. (2000) J Eng Gas Turbines Power 122:150) is the state-of-the-art model for the thermodynamic properties of water and steam for industrial applications and is routinely used for simulations of steam power cycles and utility systems. Its use in optimization-based design, however, has been limited because of its complexity. In particular, deterministic global optimization of problems with the IAPWS-IF97 is challenging because general-purpose methods lead to rather weak convex and concave relaxations, thus resulting in slow convergence. Furthermore, the original domains of many functions from the IAPWS-IF97 are nonconvex, while common global solvers construct relaxations over rectangular domains. Outside the original domains, however, many of the functions take very large values that lead to even weaker relaxations. Therefore, we develop tighter relaxations of relevant functions from the IAPWS-IF97 on the basis of an analysis of their monotonicity and convexity properties. We modify the functions outside their original domains to enable tighter relaxations, while we keep them unchanged on their original domains where they have physical meaning. We discuss the benefit of the relaxations for three case studies on the design of bottoming cycles of combined cycle power plants using our open-source deterministic global solver MAiNGO. The derived relaxations result in drastic reductions in computational time compared with McCormick relaxations and can make design problems tractable for global optimization.


Introduction
Water is one of the most important substances in energy conversion systems. Today, around three quarters of the global power generation relies on steam power cycles (IEA 2019), either in pure steam power plants or as part of combined-cycle plants. Furthermore, water is the most common heat transfer fluid for applications ranging from domestic heating to combined heat and power supply in utility systems for industrial production sites (Podolski et al. 2008). Given the widespread use of water in energy conversion systems, even small improvements in the design and operation of such systems can have a significant impact.
Despite the long industrial experience in the design of steam-based energy conversion systems, the optimal design of such systems (e.g., with respect to energy efficiency or cost) for given boundary conditions is still an active field of research. Steam cycle design is tackled with a variety of model-based approaches (Wang et al. 2019) ranging from simulation-based evaluation of designs manually derived from engineering experience through advanced exergy-based analyses for identifying promising points for improvement to optimization-based methods.
The state-of-the-art model for the required thermodynamic properties of water and steam is the IAPWS-IF97 (Wagner et al. 2000) in its revised version (IAPWS 2007a). It is recommended by the International Association for the Properties of Water and Steam (IAPWS) "for industrial use (primarily the steam power industry) for the calculation of thermodynamic properties of ordinary water in its fluid phases" (IAPWS 2007a). It was developed on the basis of the IAPWS-95 model (Wagner and Pruss 2002) with the goal of providing explicit expressions for all common calculations without iterative solution and with low computational effort.
While the IAPWS-IF97 is in fact the standard model for the simulation of steam cycles used in many simulators, its use in optimization-based approaches has been more limited since it is often considered too complex (Wang et al. 2019). Many existing optimization studies on steam power cycle or utility system design have either directly opted for a simpler model (Bruno et al. 1998;Bongartz and Mitsos 2017) or replaced the IAPWS-IF97 with a simplified surrogate model, typically polynomials of lower degree (Ahadi-Oskui et al. 2010). The latter approaches have sometimes been combined with smoothing techniques to remedy the nondifferentiabilities that occur at phase boundaries (Tică et al. 2012;Åberg et al. 2017). Those studies using the IAPWS-IF97 itself have done so employing different methods: stochastic optimizers (Nadir and Ghenaiet 2015), often coupled with a simulation software implementing the IAPWS-IF97 (Koch et al. 2007;Luo et al. 2011;Wang et al. 2014); local NLP solvers (Zebian et al. 2012); convex MINLP solvers (Savola et al. 2007;Manassaldi et al. 2011Manassaldi et al. , 2016; or a combination of local and stochastic optimizers (Wang et al. 2015(Wang et al. , 2016. The existing optimization approaches using the IAPWS-IF97 share the limitation that the local or stochastic optimizers employed cannot guarantee to find a global optimum since the resulting optimization problems are nonconvex and can have multiple local solutions even for the simplest cycles (and even with much simpler models) (Bongartz and Mitsos 2017). Therefore, deterministic global optimization is desirable for solving steam cycle design problems with the IAPWS-IF97. The most rigorous existing approach that we are aware of is the work of Ahadi-Oskui et al. (2006), who use a branch-and-cut approach for optimizing the design of combined cycle power plants modeled with the IAPWS-IF97. However, their approach is also partially heuristic since it relies on sampling of black box model functions (Nowak and Vigerske 2008).
In deterministic global optimization, a key challenge is the construction of tight convex and concave relaxations and range bounds of the functions used in the model (Locatelli and Schoen 2013). Although the form of the IAPWS-IF97 functions is such that they could in principle be addressed with general purpose methods for deriving relaxations such as BB (Androulakis et al. 1995), the auxiliary variable method (AVM) (Smith and Pantelides 1997;Tawarmalani and Sahinidis 2002), or the McCormick technique (McCormick 1976), this is challenging for the following reasons: First, the resulting relaxations are expected to be rather weak given the complexity of the functions. In global optimization, weak relaxations result in slow convergence that makes larger problems intractable. Second, when using the AVM, the subproblems for computing lower bounds will get very large because many auxiliary variables will be added, again because of the complexity of the functions. To obtain tighter relaxations and reduce computational time in global optimization, Schweidtmann et al. (2019) recently proposed replacing complex thermodynamic models with artificial neural networks as surrogate models in the design of organic Rankine cycles. However, given the very high accuracy of the IAPWS-IF97 as well as its role as an established and trusted model implemented in most simulators, we would like to avoid replacing it with a surrogate model.
In this work, we retain the original model functions from the IAPWS-IF97 in those regions where they have physical meaning and construct tighter convex and concave relaxations and range bounds to speed up the global optimization of steam cycle design problems. Our approach is based on the observation that functions in thermodynamic models often exhibit useful monotonicity and convexity properties that can be exploited to derive tighter relaxations than those obtained by applying general purpose methods to their complex functional forms (Najman et al. 2019a, b).
In the following, we first provide some background on the methods employed for constructing relaxations. Next, we analyze the relevant functions from the IAPWS-IF97 model and describe the derived relaxations. Finally, we discuss the application of the relaxations to three case studies of combined cycle power plants to demonstrate the computational speedup compared with a plain application of the McCormick technique to the IAPWS-IF97. The results demonstrate that deterministic global optimization of problems with the IAPWS-IF97 can get intractable for larger problems when using McCormick relaxations as a general purpose technique, whereas larger problems can be solved with the proposed tailored relaxations. The proposed relaxations are implemented in the MC++ library (Chachuat et al. 2015) used by our open-source global optimizer MAiNGO (Bongartz et al. 2018). 1

Preliminaries and methods used
In the present work, we aim at enabling the use of IAPWS-IF97 functions as intrinsic functions in factorable programming techniques, in particular the multivariate McCormick technique. In the following, we thus first summarize the multivariate McCormick technique before discussing two methods that are used in Sect. 3.4 to derive relaxations of the IAPWS-IF97 functions, namely methods for relaxations of componentwise convex or concave functions and variants of the BB method.
Throughout this article, we represent scalars as lowercase or uppercase letters (e.g., p or T), vectors in boldface (e.g., ), and sets in calligraphic typeface (e.g., X ). ℝ n denotes the set of all nonempty compact interval subsets of ℝ n . The superscripts L and U denote the left and right end points of an interval, respectively (e.g.,

Multivariate McCormick relaxations
Established methods for deterministic global optimization are based on spatial branch-and-bound (B&B) (Falk and Soland 1969) and rely on the availability of valid convex and concave relaxations (Locatelli and Schoen 2013).

Definition 1 (Convex and concave relaxation)
In Sects. 2.3 and 3.4, we make use of non-standard under-and overestimators to derive tighter relaxations.

Definition 2 (Concave underestimator and convex overestimator)
The original McCormick technique (McCormick 1976) provides rules for computing relaxations of so-called factorable functions consisting of finite compositions of binary sums, binary products, and a library of univariate functions. The latter are called intrinsic functions, and for these functions both convex and concave relaxations and range bounds need to be known.
Tsoukalas and Mitsos (2014) extended the McCormick technique to allow for multivariate intrinsic functions. Note that such intrinsic functions are also used in the context of the AVM, such that the relaxations and range bounds developed herein could also be applied in that context. Unlike the AVM (Smith and Pantelides 1997;Tawarmalani and Sahinidis 2002), the McCormick technique provides convex and concave relaxations in the original variable space. A challenge in the application of the multivariate McCormick method is that the evaluation of the relaxation requires the solution of a convex but possibly nonlinear and nonsmooth problem (Theorem 2, Tsoukalas and Mitsos (2014)). To this end, it is highly desirable for the developed relaxations to have specific monotonicity properties that allow us to either determine closed-form solutions of the problem described by Theorem 2 of Tsoukalas and Mitsos (2014) analytically, or compute it numerically with relatively little effort, e.g., by solving a one-dimensional nonlinear equation, which is solved via Newton's method.

Componentwise convex functions
Componentwise convexity of multivariate functions enables the construction of tight relaxations.

Definition 4 (Componentwise convexity)
For twice continuously differentiable functions, a convenient way of confirming componentwise convexity or concavity with respect to a variable consists in examining the second derivative of the univariate function f in (1) (and hence the second partial derivative with respect to the variable of interest) and proving that it is non-negative or non-positive over the considered set (cf. Rockafellar 1970), either analytically or through global maximization or minimization. However, since some of the functions considered herein are only piecewise continuously differentiable, we use a specific procedure that enables an analysis for functions where a partial derivative need not exist at all points. The functions we consider herein are of the form where X 1 , X 2 ∈ ℝ , f 1 , f 2 are both smooth on X 1 × X 2 , f is continuous on X 1 × X 2 , and x 2 (x 1 ) is a strictly monotonic function. To determine whether functions of form (2) are componentwise convex, we proceed as described in the following. First, we solve the optimization problems to global optimality. Note that these problems are much smaller and less complex than the design problems within which the functions are intended to be used and can be solved with general purpose methods. Furthermore, since we only consider certain known functions, the problems need to be solved only once in order to construct tighter relaxations for later use in design problems. If the solution values of both problems (3) and (4) are non-negative, we know that f 1 is componentwise convex with respect to x 2 on S 1 ∶= {(x 1 , x 2 ) ∈ X 1 × X 2 |x 2 ≥x 2 (x 1 )} and f 2 is componentwise convex with respect to x 2 on S 2 ∶= {(x 1 , x 2 ) ∈ X 1 × X 2 |x 2 ≤x 2 (x 1 )} . To show that f is componentwise convex with respect to x 2 on X 1 × X 2 , it remains to show that the difference between the right and left derivatives of f with respect to x 2 at the boundary between the subdomains S 1 and S 2 is non-negative. To this end, we solve a third optimization problem given as and examine the sign of the optimal objective value. When applying the above procedure for testing componentwise convexity with respect to x 1 , the order of the derivatives of f 1 and f 2 in (5) needs to be exchanged in case x 2 (x 1 ) is increasing. An analogous procedure can be used for testing for componentwise concavity.
Once it has been established that a function is componentwise convex (concave), its concave (convex) envelope is known to be vertex polyhedral (Tardella 2004). Since the multivariate functions considered herein are only bivariate, this envelope can be computed efficiently using the method of Meyer and Floudas (2005). A convex (concave) relaxation, on the other hand, can be computed using the method of Najman et al. (2019a) if the derivatives of the function fulfill min certain additional requirements. However, since some of the functions herein are nonsmooth, we need to adapt the latter method. First, we define partial subderivatives in analogy to the definition of subgradients for convex functions (cf. Rockafellar 1970).
Definition 5 (Subgradient and subdifferential) For a convex and concave function f cv , f cc ∶ X → ℝ , where X ∈ ℝ n , we call cv , cc ∈ ℝ n a convex and a concave subgradient of f cv , f cc at ̂ , respectively, iff The convex and concave subdifferential of f at ̂ denoted by cv f (̂ ), cc f (̂ ) , respectively, are the sets of all convex and concave subgradients of f at ̂ , respectively.
Definition 6 (Partial subderivative and subdifferential) Let f ∶ X → ℝ , where X ∈ ℝ n , be componentwise convex (concave) with respect to some x i , i ∈ {1, … , n} . We call s ∈ ℝ a partial subderivative of f with respect to x i at ̂ ∈ X iff it is a subgradient of the univariate function f as defined in (1) at x i . The set of all partial subderivatives of f with respect to x i at ̂ ∈ X is called the partial subdifferential of f with respect to x i at ̂ ∈ X and denoted by i f (̂ ).
Using this definition of the partial subdifferential, we introduce the following adapted version of Theorem 1 from Najman et al. (2019a): If it holds that for any fixed x 2 ∈ X 2 , then f cv,u sum is a convex relaxation of f.

Proof
The proof is analogous to that of Theorem 1 of Najman et al. (2019a) when replacing the partial derivatives therein with the suitable minimum and maximum elements of the partial subdifferentials as in (7), and replacing the derivatives in Lemma 1 of Najman et al. (2019a) with the subgradients for convex (rather than differentiable) functions. ◻ In particular, Theorem 1 can be applied to componentwise convex functions of form (2). Similar to Najman et al. (2019a), up to 2 valid relaxations can be obtained from Theorem 1 by reordering the variables and changing the sign of the coordinates. Furthermore, an analogous way for constructing concave relaxations of componentwise concave functions is given by considering −f instead.
In Corollary 1 of Najman et al. (2019a), a convenient procedure based on mixed second-order partial derivatives is provided for confirming that the conditions of Theorem 1 therein are satisfied. This procedure is not directly applicable for confirming the validity of (7) because functions of form (2) are possibly not twice differentiable, despite the fact that f 1 and f 2 are both twice differentiable. Therefore, we describe an alternative way of confirming the validity of (7) for functions considered of form (2), which is analogous to the test for componentwise convexity described above: First, we solve the two optimization problems to global optimality. If the solution values of both problems (8) and (9) are non-negative, Corollary 1 of Najman et al. (2019a) suggests that the corner c = (x L 1 , x L 2 ) can be used in Theorem 1 in order to construct the convex relaxation. Assume w.l.o.g. that f equals f 2 at c . In order to guarantee that the relaxation resulting at c is also valid for f, a sufficient condition is that is non-negative for all x 1 ∈ X 1 if x 2 (x 1 ) is decreasing and non-positive if x 2 (x 1 ) is increasing. Since this condition is equivalent to the one for confirming componentwise convexity via problem (5) or concavity via the equivalent maximization problem, respectively, the condition is automatically satisfied for componentwise concave functions in case x 2 (x 1 ) is increasing and componentwise convex functions in case x 2 (x 1 ) is decreasing. An analogous procedure can be used for testing whether the corner point c = (x L 1 , x U 2 ) can be used in Theorem 1.

Variants of ˛BB
The BB method was introduced as a general-purpose method for constructing relaxations of twice continuously differentiable functions (Androulakis et al. 1995).
The basic idea is to add a quadratic function with parameters to construct convex underestimators, where the parameters are chosen such that the quadratic function offsets the smallest eigenvalue of the Hessian of the function over the considered domain. These parameters can either be precomputed and thus be independent of (8) min min the interval, or they can be computed using interval methods when constructing the relaxation for a specific interval (Adjiman et al. 1998).
Hasan (2018) recently introduced a variant of BB that does not use the parameters to construct convex underestimators directly, but rather to construct a componentwise concave underestimator, and then uses the vertex polyhedral convex envelope of the latter (cf. Sect. 2.2) as convex underestimator. Herein, we use a similar approach in the sense that we only aim at achieving componentwise convexity or concavity and not convexity or concavity of the function directly.

Proof
The result for f cv,u, BB follows when applying the original BB method (Androulakis et al. 1995) to the univariate function f in (1). For f cv,o, BB , it follows from the results of Hasan (2018) applied to −f . ◻

Proof
The proof is analogous to that of Lemma 1. ◻ Since both the original BB method and the variant of Hasan (2018) require the function to be twice continuously differentiable, we consider the following procedure for two dimensional functions of form (2): If the solution of problem (5) is nonnegative, we use the minimum of the solution values of problems (3) and (4) to construct under-and overestimators as in Lemma 1 that are componentwise convex with respect to x 1 . Instead, if the solution value of the maximization problem analogous to (5) is non-positive, we solve the maximization problems analogous to problems (3) and (4) and use the maximum of their solution values to construct under-and overestimators as in Lemma 2 that are componentwise concave with respect to x 1 . The same procedure is used to construct under-and overestimators that are componentwise convex or concave with respect to x 2 . In the present work, we only use BB to compute relaxations of intrinsic functions within larger factorable functions, which are then relaxed using the multivariate McCormick theorem (cf. Sect. 2.1). Therefore, we need to be able to compute minima of convex relaxations and maxima of concave relaxations, respectively, over boxes, either analytically or with little effort. However, the BB-type relaxations (10)-(13) may not be monotonic even when the original functions are, thus complicating the application of the multivariate McCormick theorem. We therefore add a linear term to functions described in (10)-(13) to make the final BB relaxation monotonic. While this makes the final relaxations slightly weaker, it greatly simplifies the application of the multivariate McCormick theorem.
In order to apply Theorem 1, it may additionally be necessary to alter the mixed second-order partial derivatives of a function f, which we achieve with mixed BB terms. For example, to obtain an underestimator that has a non-negative mixed second-order partial derivative with respect to x i x j in analogy to Lemma 1. This procedure can be used analogously to achieve a non-positive mixed second-order partial derivative and for overestimators.

Construction of relaxations of functions from the IAPWS-IF97
The IAPWS-IF97 is divided into five regions (cf. Fig. 1). Region 1 represents subcooled liquid up to a temperature of T max 1 = 623.15 K , Region 2 superheated vapor, and Region 4 the two-phase region. Of the latter, we only consider the part up to T min 3 = T max 1 (and, accordingly, p min 3 ≈ 16.5292 MPa ) in which the saturated liquid and vapor states lie in Regions 1 and 2, respectively. This part will be denoted as Region 4-1/2 in the following. Additionally, we restrict Region 2 to pressures p ≥ p min 2 ∶= p min 1 = 611.2127 Pa instead of the original definition p ≥ 0 (IAPWS 2007a) to avoid the difficulty that entropy goes to infinity as pressure goes to zero. Since such low pressures are not typically encountered in steam cycles, this restriction of Region 2 is not a practical limitation for the intended application. We do not consider Region 5, since it corresponds to much higher temperatures than those encountered in steam cycles. Region 3, which represents the region around and above the critical point, can be relevant for steam cycles and could be included as future work.
Within the considered regions, we are interested in those functions that relate temperature T, pressure p, specific enthalpy h, specific entropy s, and vapor fraction y. Specifically, we consider the following functions taken directly from the IAPWS-IF97: In Region 2, we do not consider the functions T 2 (p, h) and T 2 (p, s) from the IAPWS-IF97 although they would be useful to eliminate optimization variables in reduced-space optimization formulations (Bongartz and Mitsos 2017). The reason is that these functions have piecewise definitions over subdomains in the p-h or p-s space. Since the functions were independently fitted to the original data for each of these subdomains during the development of the IAPWS-IF97, the resulting functions T 2 (p, h) and T 2 (p, s) are not continuous over their entire domains. Although the differences between the one-sided limits at the interface between different subdomains is within the tolerated error for model development (Wagner et al. 2000), this discontinuity does complicate the analysis of the functions for deriving relaxations. Similar complications arise in Region 3, for which useful supplementary functions with piecewise definitions have been introduced (IAPWS 2007b).
In addition to the functions taken directly from the IAPWS-IF97, we define the following auxiliary functions for calculating specific enthalpy and entropy of saturated vapor and liquid states: Using (14)- (17), we define the following functions for computing specific enthalpy and entropy as well as vapor fraction in the two-phase region: For the bivariate functions, we define two types of domains that are useful for the analysis of the functions.
Definition 7 (Physical domain) Given a function f (x 1 , x 2 ) , we denote by P f (x 1 ,x 2 ) ⊂ ℝ 2 the domain specified in the IAPWS-IF97 (IAPWS 2007a) and we call In the following, the lower and upper variable bounds of the box domains are denoted with the superscripts min and max, respectively, i.e., . The physical domains, which are often nonconvex (cf. Fig. 1), can be represented in terms of the box domains and additional inequalities. This representation is useful because in global optimization the functions typically need to be evaluated on rectangular subsets of the box domains (called nodes of the B&B tree) while the inequalities are enforced as constraints. To distinguish the bounds of nodes from those of the box domains, we denote the former by the superscripts L and U, e.g., .
All functions of the IAPWS-IF97 and by extension also those in (14)-(21) are given as explicit expressions, in most cases compositions of linear and signomial functions, and can be written as factorable functions. Therefore, general purpose methods for constructing relaxations of factorable functions such as the (multivariate) McCormick technique (McCormick 1976;Tsoukalas and Mitsos 2014) or the auxiliary variable method (Smith and Pantelides 1997;Tawarmalani and Sahinidis 2002) could be applied to the original, factorable representation of the functions. However, this is problematic for two reasons. First, the functional forms of the expressions are rather long and complex. Therefore, general purpose methods often result in rather weak relaxations. Second, during optimization, boxes need to be considered that contain regions of the box domain that lie outside the physical domain. Even though such regions are made infeasible using a suitable constraint, the functions may still need to be evaluated at points in such regions and relaxations need to be constructed over boxes containing these regions as well. However, the functions were never intended to be used outside their physical domain, and many functions take values of extremely large magnitude when evaluated in such regions. This can lead to very weak relaxations over the physical domain when considering boxes containing these regions, as well as numerical problems when solving the linear and nonlinear subproblems.
To address these challenges, we consider the IAPWS-IF97 functions directly as intrinsic functions instead of using their factorable representations. To derive the required information to treat them as intrinsic functions, we conduct the following steps that are explained in detail in the following subsections: 1. We determine the physical domains and box domains for bivariate functions. 2. We modify the bivariate functions outside their physical domains to improve the mathematical properties of the functions in these regions. 3. We analyze the functions for monotonicity properties to derive range bounds. 4. We analyze the functions for convexity properties to derive convex and concave relaxations.

Determination of physical domains and box domains
For the bivariate functions, we first determine the physical domains and box domains according to Definitions 7 and 8. For the functions h 1 (p, T) , s 1 (p, T) , h 2 (p, T) , and s 2 (p, T) , which are obtained as derivatives of the basic equations of the IAPWS-IF97, the domains can be taken directly from the model definition (IAPWS 2007a). For most other functions, the physical domains are not explicitly given in the IAPWS-IF97. In this case, the physical domains need to be determined through an analysis of the monotonicity properties of related functions from the IAPWS-IF97 and the box domains need to be determined via globally maximizing and minimizing each variable over the physical domain. Examples for this procedure both for a straightforward case and a more involved case can be found in "Appendix 1" along with a summary of all box and physical domains in Tables 5 and 6.

Model modification outside the physical domain
We analyze the model functions and their derivatives to determine whether the functions exhibit excessive peaks or other undesired behavior in those regions of their box domains that are outside the physical domain. In case of undesired behavior, we replace the functions with a suitable extrapolation outside (a relaxation of) the physical domain. The resulting piecewise defined functions are of form (2)  Beyond these intermediate modifications, we restrict the range of each function to that achieved over its physical domain. This helps to avoid domain violations when considering compositions of functions from the IAPWS-IF97. The resulting functions will be called final modified functions in the following and denoted with the superscript mod.
The final modified functions are the ones that will be used in the optimization problems. The relaxations, however, will be constructed based on the intermediate modified functions that have more useful convexity properties (cf. Sect. 3.4). Based on these relaxations of the intermediate modified functions, relaxations for the final modified functions are obtained using the rules for composition with the max and min functions (Tsoukalas and Mitsos 2014). Figure 2 shows the application of the above procedure to the function h 2 (p, T) . A more detailed description for this function can be found in "Appendix 2". Similar modifications are conducted for the functions h 1 (p, T) , s 1 (p, T) , T 1 (p, h) , and s 2 (p, T) . The intermediate modifications are summarized in Table 7 in "Appendix 2". The remaining bivariate functions are only cut at the minimum and maximum values occurring on their physical domains but not extended otherwise since their properties outside their physical domains are already satisfactory.
Note that the modifications of the model functions themselves are only conducted in parts of the box domain where the model has no physical meaning and that are excluded via suitable constraints when using the functions in an optimization problem. They thus do not alter the solutions of meaningful process models using these functions. Note also that the modifications do make the functions nonsmooth outside of the physical domains (both the replacement of the function with an extrapolation for some of the intermediate modified function and the restriction of the range of the functions to that over their physical domains when constructing the final modified functions). However, we did not find this to be an issue in practice so far (cf. discussion in Sect. 4).

Range bounds
To derive range bounds, we analyze the intermediate modified functions for monotonicity properties by globally maximizing and minimizing their first partial derivatives over the corresponding box domains using our open-source global solver MAiNGO (Bongartz et al. 2018). For functions that are replaced by an extrapolation outside their physical domains and that are hence of form (2), the maximization and minimization is done separately within the relaxed physical domain and the domain of the extrapolation. In both subdomains, the functions are differentiable by construction. Even though the functions may be nonsmooth at the boundary between the subdomains, monotonicity results within the subdomains can still be translated Since this boundary is nonsmooth, we introduce a relaxed physical domain delimited by the smooth function p lim 2 (T) instead. b Original function h 2 (p, T) with undesired peaks outside P h 2 (p,T) that go beyond 4.5 × 10 9 kJ/kg and −6.1 × 10 55 kJ/kg and were cut off to improve readability. c Intermediate modified function h int 2 (p, T) that was modified (green) for p > p lim 2 (T) to have favorable properties for deriving relaxations. d Final modified function h mod 2 (p, T) = max(h int 2 (p, T), h min 2 ) , where h min 2 is the minimum of h 2 (p, T) over P h 2 (p,T) . The region where this lower bound is active is shown in orange. Note the different scale of the z-axes compared with b and c into global monotonicity properties since the functions are continuous. These properties also hold for the final modified functions since taking the maximum or minimum with a constant does not change the monotonicity properties.
For functions with monotonicity properties over the entire (box) domain, we immediately obtain exact range bounds. For other functions, exact range bounds can only be obtained in certain cases. In case a function is only monotonic with respect to one variable, an exact upper or lower bound can sometimes still be obtained in case of componentwise convexity or concavity with respect to the other variable (cf. Sect. 3.4). As a last resort, we obtain natural interval extensions from the FILIB++ library (Lerch et al. 2011) by evaluating the factorable representation of the function using the interval datatypes from FILIB++ either with respect to one or both variables over a suitable part of its domain. Examples for both a straightforward case and a more involved case can be found in "Appendix 3", and all exploited monotonicity properties are summarized in Table 8.

Convex and concave relaxations
Similar to the monotonicity analysis described above, we analyze the intermediate modified functions for (componentwise) convexity by globally maximizing and minimizing their second partial derivatives over their box domains using MAiNGO. For bivariate functions with piecewise definition of form (2), we proceed as described in Sect. 2.2. The identified properties are summarized in Table 9 in "Appendix 4". These properties are used to derive relaxations as described in the following.

Univariate functions
For univariate functions that are convex or concave over their entire domain (cf .  Table 9), we trivially obtain the convex and concave envelopes over any interval as the function itself and the secant between the end points of the interval. The functions h liq 4-1/2 (p) , s liq 4-1/2 (p) and s vap 4-1/2 (p) , on the other hand, are not convex or concave on their entire domains. However, they are either mostly convex or mostly concave, while the remaining parts are essentially linear. We can therefore obtain rather tight relaxations via the BB variants (10)-(13) with fixed values for the parameters. These values are obtained as described in Lemmata 1 and 2 from the maximum or minimum values (whichever is smaller in magnitude) of the second derivatives over the entire domain that were precomputed through global optimization in MAiNGO. When constructing relaxations on a given interval during the B&B procedure, the BB terms in (10)-(13) are only used if the node is not fully in a region where the function is convex or concave (cf. Table 9). Since the functions are univariate, convex and concave envelopes could instead also be obtained using the technique described by McCormick (1976, Section 4). However, in general this method requires iterative solution of a nonlinear equation using, e.g., Newton's method, for each evaluation of the relaxations, and we would like to avoid this computational effort.
As an example, Fig. 3c shows the proposed relaxations for h liq 4-1/2 (p) , which are discussed in more detail in "Appendix 4". The relaxations are orders of magnitude tighter than the McCormick relaxations (cf. Fig. 3a). Quantitatively, the maximum difference between the value of h liq 4-1/2 (p) and that of the convex relaxation is approximately 1.3 × 10 7 for the McCormick relaxation whereas it is less than 750 for the proposed relaxation. The proposed convex relaxation, which is based on the secant of the concave underestimator (12) and hence the maximum of the second derivative of h liq 4-1/2 (p) , is also much tighter than the regular BB relaxation (10) based on the minimum of the second derivative of h liq 4-1/2 (p) . For this regular BB relaxation, the maximum difference between the values of the function and that of the relaxation is approximately 1.5 × 10 4 (cf. Fig. 3b).
Beyond the tightness of the relaxations, another factor that can impact the performance of a global solver is the computational cost for computing the relaxations. To quantify this computational cost, we evaluated both the McCormick relaxations and the proposed relaxations (using the implementation described in Sect. 4) for all considered univariate functions on 100 evenly spaced points on each of 1000 randomly generated intervals within their domains and measured the CPU time.

Bivariate functions
The considered bivariate intermediate modified functions are neither convex nor concave on most parts of their domains. However, they are often componentwise convex or concave with respect to one or both variables on large parts of their domain (cf . Table 9). Additionally, they are often almost componentwise convex or concave over the entire domain in the sense that either the maximum of the second partial derivative is much larger in magnitude than the minimum or vice versa. This property is exploited to construct relaxations as described in the following, where we consider the function h mod 2 (p, T) as an example. First, in addition to the identification of regions where the intermediate modified function (e.g., h int 2 (p, T) ) is componentwise convex or concave with respect to some variable, we identify the minimum and maximum values of its second partial derivatives via global optimization using MAiNGO (separately over the subdomains in For those regions where the function is not componentwise convex with respect to a variable, we apply componentwise BB with respect to that variable using either (10) and (11) or (12) and (13), depending on whether the function is almost componentwise convex or concave with respect to that variable in the above sense (cf. Fig. 4b). If necessary, we add a linear term to (10)-(13) to ensure that the relaxation is monotonic, and we include a mixed BB term to ensure a constant sign of the mixed second-order partial derivatives (cf. Sect. 2.3).
In most cases, the resulting under-and overestimators are already either componentwise convex or concave (with respect to both variables). If, however, an underestimator is componentwise convex with respect to one variable and componentwise concave with respect to the other, we construct a componentwise convex underestimator by taking the secant with respect to the concave variable (Najman et al. 2019a). For an overestimator, we instead construct a componentwise concave overestimator by taking the secant with respect to the convex variable. For componentwise convex underestimators and componentwise concave overestimators, we use Theorem 1 to obtain convex underestimators and concave overestimators, respectively. For componentwise concave underestimators and componentwise convex overestimators, we use the method of Meyer and Floudas (2005) instead (cf. Fig. 4c).
Up to this point, the relaxations were constructed based on the intermediate modified functions (e.g., h int 2 (p, T) ) before cutting at the minimum and maximum values over the physical domain. We therefore apply the rules for composition with max and min functions (Tsoukalas and Mitsos 2014) to derive valid relaxations for the final modified functions. For h mod 2 (p, T) , the relaxations are shown in Fig. 4d. They are significantly tighter than the McCormick relaxations of h 2 (p, T) (cf. Fig. 4a).
In some cases, however, the resulting relaxations are still somewhat weak when considering large boxes, in particular when requiring large values for achieving componentwise convexity (cf. Fig. 5a). In these cases, we manually construct  Fig. 4 are somewhat weak over larger boxes. b We therefore add ad hoc relaxations that do not converge but are valid over the entire domain and provide tighter relaxations on large boxes additional (linear or nonlinear) ad hoc relaxations. These ad-hoc relaxations are valid over the entire box domain and are independent of the subset of the box domain over which they are evaluated in a B&B procedure. Therefore, they do not converge to the function when considering boxes of decreasing size, but they help tighten the relaxations for large boxes (cf. Fig. 5b). To construct the ad-hoc relaxations, we visually inspect the graph of the respective functions by plotting them in MATLAB. We then simultaneously plot linear or convex nonlinear functions, the potential ad-hoc relaxations, with parameters that we adjust through trial and error until the functions visually appear to be valid relaxations (e.g., in case of a convex relaxation, they appear to be below the original function on the entire box domain). We then confirm the validity of these potential ad-hoc relaxations by globally minimizing the difference between them and the original function in MAiNGO and examining the sign of the optimal objective value.
Finally, to evaluate the proposed relaxations at a given point, we select the maximum among the described convex relaxations and the minimum among the concave relaxations, including the multiple facets that are obtained from the application of both Theorem 1 and the method of Meyer and Floudas (2005). Additionally, we cut the relaxations off at the lower and upper range bounds.
The full procedure described above is used for the functions h 1 (p, T) , s 1 (p, T) , T 1 (p, h) , h 2 (p, T) , s 2 (p, T) , and y 4−1∕2 (p, h) . For the functions T 1 (p, s) , h 4−1∕2 (p, y) , s 4−1∕2 (p, y) , and y 4−1∕2 (p, s) , the McCormick relaxations of the factorable representations or the compositions using the relaxations of the univariate functions discussed in Sect. 3.4.1 are relatively good already so that these are used instead. For T 1 (p, s) and y 4−1∕2 (p, s) , we merely add ad-hoc relaxations for large boxes as described above.
While the proposed relaxations of the considered bivariate functions are at least as tight as and typically significantly tighter than the McCormick relaxations, unlike in the univariate case, they are not always cheaper to evaluate. When evaluating them on 100 evenly spaced points on each of 1000 randomly generated boxes within their box domains, only the proposed relaxations of h 2 (p, T) , s 2 (p, T) , h 4−1∕2 (p, y) , s 4−1∕2 (p, y) , and y 4−1∕2 (p, s) are faster to evaluate (80-99%) than the McCormick relaxations. The time for evaluating the proposed relaxation for T 1 (p, s) is virtually the same as that for evaluating the McCormick relaxation, since we merely add an ad-hoc relaxation, which is very cheap to compute. For the functions h 1 (p, T) , s 1 (p, T) , T 1 (p, h) , and y 4−1∕2 (p, h) , however, the proposed relaxations take 160-6500% longer to evaluate than the McCormick relaxations. The reason for this higher computational cost is that for these functions, iterative solution of onedimensional nonlinear equations via Newton's method is required to determine the correct point to be used in the multivariate McCormick composition theorem (cf. Sect. 2.1). Nevertheless, computational experiments confirmed that this additional effort for computing the relaxations typically pays off in global optimization thanks to the much better tightness than the McCormick relaxations.

Case studies
To demonstrate the benefit of the derived relaxations, we use the design problems from our previous work (Bongartz and Mitsos 2017) that consider bottoming cycles for combined power plants in three levels of complexity: • Case Study 1 Basic Rankine cycle with fixed temperature of the hot gas turbine exhaust at the outlet of the heat recovery steam generator (HRSG). • Case Study 2 Basic Rankine cycle but with variable gas outlet temperature and with a turbine bleed to an open feedwater heater that also serves as deaerator. • Case Study 3 Dual-pressure cycle where the outlet of the high-pressure (HP) turbine is mixed with the outlet of the low-pressure (LP) superheater before entering the LP turbine. The latter also has a bleed stream to the deaerator.
Two objective functions are considered for each flowsheet, namely maximization of the net power output ( Ẇ net ) and minimization of the levelized cost of electricity (LCOE). The latter objective is more complex because it includes investment cost and hence, the models contain sizing and cost correlations for the process units.
In our previous work (Bongartz and Mitsos 2017), we demonstrated the benefit of modeling the design problems in a reduced-space optimization formulation where only the design variables and very few model variables remain in the optimization problem and most other model variables (and equations) are collapsed in a sequential evaluation of the model going through the cycle. In this previous work, we used very simple models for the thermodynamic properties of water, mostly ideal gas and liquid with constant heat capacities and constant heat of vaporization at a reference temperature. In the present work, we consider the same design problems and replace the simple thermodynamic models with the IAPWS-IF97 and solve the resulting problems using MAiNGO.

Modeling and implementation
Compared to the original reduced-space formulation (Bongartz and Mitsos 2017), we need to make some slight changes to the model. These changes are required because the functions T 2 (p, h) and T 2 (p, s) for computing temperatures for given pressure and enthalpy or entropy in the gas phase are not available (cf. discussion in Sect. 3). Specifically, we need to use steam temperatures at the superheater outlets as optimization variables, either instead of the corresponding enthalpies as degrees of freedom (for Case Studies 2 and 3) or as an additional variable with a corresponding equality constraint (for Case Study 1). For Case Study 3, we also need to add the steam temperature at the inlet of the LP turbine, i.e., after mixing the HP turbine outlet with the LP superheater outlet, as well as the temperature of the hypothetical isentropic state at the outlet of the HP turbine. Additionally, we use mass flow rates of all branches of the flowsheet as design variables instead of the overall mass flow rate and split fractions, which we found to improve the relaxations.
We implement the models in C++ using two different ways of handling the IAPWS-IF97 functions: using their factorable representation, and considering them as intrinsic functions. In the former case, the factorable representation of each IAPWS-IF97 function is incorporated into the directed acyclic graph (DAG) built in MAiNGO by the MC++ library (Chachuat et al. 2015). Range bounds are then obtained via natural interval extensions by FILIB++ (Lerch et al. 2011), relaxations and their subgradients via the multivariate McCormick technique by MC++ , and gradients via automatic differentiation by FADBAD++ (Bendtsen and Stauning 2012). In the latter case, each IAPWS-IF97 function occurs as a single node in the DAG. In this case, the range bounds and relaxations derived in the previous sections are used instead. For gradients, we use automatic differentiation via FAD-BAD++ for each of the piecewise defined regions of the functions and arbitrarily use one of the one-sided limits at the nonsmooth kinks (the same is done for the max and min functions). While we are aware that this could potentially cause difficulties for the local subsolvers that rely on gradients, we did not experience such difficulties in the case studies (cf. Sect. 4.2). This is likely due to the fact that the performance of the local subsolvers for upper bounding is often not as crucial for the performance of global solvers. Nevertheless, in the future, it would be desirable to use recently published methods for generalized derivatives (Khan and Barton 2015) instead. Finally, we use custom relaxations for the equipment cost correlations and pressure factors for heat exchangers and deaerator vessels (Najman et al. 2019b), and for the logarithmic mean temperature difference in heat exchangers (Mistry and Misener 2016;Najman and Mitsos 2016). The process models are available via our website, 2 while the functions from the IAPWS-IF97 and the proposed relaxations are available in the MC++ library used by MAiNGO.
All problems are solved with MAiNGO v0.2.0 (Bongartz et al. 2018) using CLP 1.17.0 (Forrest et al. 2019) for the linear lower bounding problems constructed on the basis of the subgradients of the relaxations, Ipopt 3.12.12 (Wächter and Biegler 2006) as local NLP solver during pre-processing, and no local solver (pure function evaluation of the lower-bounding solution point) during the B&B. For range reduction, we use optimization-based bound tightening with trivial filtering (Gleixner et al. 2017), duality-based bound-tightening (Ryoo and Sahinidis 1995), and basic constraint propagation. The feasibility-tolerances are set to 10 −6 and the relative optimality tolerance to 10 −2 . All calculations are conducted on an Intel ® Core ™ i5-3570 CPU with 3.4 GHz running Fedora Linux 30.

Numerical results
The optimal solution points and objective values of the three case studies are summarized in Tables 1, 2 and 3. To avoid confusion with the symbols introduced in the previous sections, the stream indices are prefixed by S in the tables. For comparison, the results with the very simple ideal thermodynamic models of Bongartz and Mitsos (2017) are also shown. Note that these deviate slightly from the original values (Bongartz and Mitsos 2017) when minimizing LCOE because the cooling water temperature in the condenser needed to be modified for the chosen condenser pressure to remain feasible when using the IAPWS-IF97. While the optimal objective values differ by less than 6% between the results with the IAPWS-IF97 model and the results with the ideal model, the optimal solution points and hence the optimal cycle designs differ significantly. For example, for Case Study 1, at the solution for maximum power output the cycle pressure p S2 is at its upper bound (note that the bound p S2 ≤ 10 MPa was chosen for consistency with the original problem of Bongartz and Mitsos (2017); therein, it was chosen because of the simplistic model that was expected to perform best at moderate pressures) when using the IAPWS-IF97 model, whereas with the simple model, it is approximately 50% lower (cf. Table 1). This highlights the importance of using detailed thermodynamic models for the design of steam power cycles. When using McCormick relaxations for the IAPWS-IF97 functions, only Case Study 1 can be solved (rather quickly) within the given time limit of 24 h (see Table 4). For Case Studies 2 and 3, the global solution and thus the correct final upper bound (UBD; all problems are cast as minimization problems when implementing them for MAiNGO) is found during pre-processing as well, but the lower bound (LBD) barely improves from the values for the root node, resulting in large or even undefined (for Ẇ net in Case Study 3; the lower bound remains at minus infinity, i.e., no valid overall lower bound was identified) relative optimality gaps, defined as (UBD − LBD)∕UBD . Additionally, numerous numerical difficulties are encountered by CLP when solving the linear lower bounding problems (e.g., false infeasibility claims or solution points that are not actually feasible). These are likely due to the extremely large function values encountered outside the physical domains as well as the very weak relaxations (cf. Sect. 3). When using the proposed relaxations, the solution of Case Study 1 takes almost 97% less B&B iterations and 86% (for Ẇ net ) to 95% (for LCOE) less CPU time than the solution with McCormick relaxations. Unlike with McCormick relaxations, Case Study 2 can be solved quickly as well with either objective. For Case Study 3, only the problem of maximizing Ẇ net can be solved to the desired relative optimality tolerance of 1% within a few minutes, while the minimization of LCOE terminates at the CPU time limit of 24 h with a remaining relative optimality gap of 3.7%. Nevertheless, the optimality gap is closed much faster with the proposed relaxations than with McCormick relaxations (cf. Fig. 6). Furthermore, this problem can not be solved to the desired accuracy with the current version of MAiNGO when using the very simple ideal model of Bongartz and Mitsos (2017) either (cf. Table 4 and Fig. 6). This indicates that the difficulty with this problem is not purely due to the complexity or relaxations of the IAPWS-IF97.

Conclusion
We have derived relaxations for relevant functions from the IAPWS-IF97 that are orders of magnitude tighter than those obtained from general purpose methods like the McCormick technique. To derive the relaxations, the functions were modified  (Bongartz and Mitsos 2017) outside their original domains in regions where the model has no physical meaning but where evaluation is required during global optimization. The functions were then analyzed for monotonicity properties to construct tight range bounds and (componentwise) convexity properties to construct tight convex and concave relaxations using variants of the BB method as well as methods for relaxation of componentwise convex or concave functions.
The relaxations were tested on three bottoming cycles for combined cycle power plants of increasing complexity. For all but the simplest example, global optimization of the cycle design for either power output or levelized cost of electricity was not possible within reasonable computational time with McCormick relaxations but only with the relaxations developed herein. For the largest cycle, the minimization of the levelized cost of electricity could not be solved to the desired accuracy with the proposed relaxations either, although the optimality gap was closed much faster than with McCormick relaxations.
Future work could aim at improving the relaxations even further, which is in principle possible because the proposed relaxations are no envelopes, except for some of the univariate functions. Beyond better relaxations of the functions considered herein, compositions could be considered (e.g., h 1 (p, s) ∶= h 1 (p, T 1 (p, s)) ), given the fact that good relaxations for intrinsic functions do not always lead to good relaxations of the composite function . Hence, tighter relaxations could be achieved by considering these composite functions as intrinsic functions themselves.
In addition to having tight relaxations for the thermodynamic models, computational advantages can also result from suitable modeling of the process flowsheets, especially in the context of reduced-space optimization formulations as considered herein (Bongartz and Mitsos 2017) that aim at a sequential evaluation of large parts of the process model. To this end, it would be beneficial to enable the use of the backward equations T 2 (p, h) and T 2 (p, s) from the IAPWS-IF97. These would allow for more freedom in the way the flowsheet is modeled, and for example allow to eliminate more optimization variables and equality constraints from the optimization problem and move them into the flowsheet evaluation, thus leading to a smaller problem and potentially further reduced runtime. When developing relaxations for T 2 (p, h) and T 2 (p, s) , care needs to be taken to handle the discontinuities induced by the piecewise definition of the functions. Similar difficulties arise when extending the present approach to Region 3 of the IAPWS-IF97 (Wagner et al. 2000;IAPWS 2007b), which would enable optimization of transcritical and supercritical cycles.
Finally, from a modeling perspective, it would be desirable to avoid having to specify which point in the flowsheet lies in which region of the IAPWS-IF97. This could be achieved either by using integer variables to let the optimizer choose between subregions, or by considering functions with piecewise definition over the regions, e.g., a function T(p, h) that consists of the respective functions in the different regions. Similar approaches have already been used for local dynamic optimization (Tică et al. 2012;Åberg et al. 2017), but in conjunction with simplifications (and smoothing) of the IAPWS-IF97. In analogy to the present approach, the functions could also be kept unchanged where they have physical meaning and instead be analyzed to derive tighter relaxations.
The parametric optimization problems in (22) and (23) h 1 (p, T max 1 ), otherwise.  (T) (cf. Fig. 2b). These peaks destroy the monotonicity properties that h 2 (p, T) exhibits on its physical domain  with parameters k 5 and k 6 . The extension for p > p lim 2 (T) is chosen such that h int 2 (p, T) is continuous, it is increasing with respect to p and decreasing with respect to T, it is componentwise concave with respect to p, componentwise concave with respect to T < 0 ). Since � PT ⊂ P h 1 (p,T) , we have h mod 1 (p,T) = h 1 (p,T) for every (p,T) ∈ � PT and thus h mod 1 (p,T) has a factorable representation over P T and we can use natural interval extensions from FILIB++ to obtain an underestimation of the minimum function value. Another, potentially tighter, lower bound can be obtained by exploiting the monotonicity with respect to T that implies that a lower bound over P × {T} is a valid lower bound over P × {T L } for every T ≤ T L . In particular, this holds for T = 510 K , for which we know that the lower bound is attained at p L .

Appendix 4: Convexity and relaxations
As an example for a univariate function, we consider the function h liq 4-1/2 (p) , which is defined on P h liq 4-1/2 (p) = [611.2127 × 10 −6 , 16.5292] MPa but is convex only on [611.2127 × 10 −6 , 14.48] MPa (cf. Table 9). By globally maximizing the second derivative of h liq 4-1/2 (p) , we obtain ∶= 0.5 × 1.0592301 kJ∕(kgMPa 2 ) ≥ 0.5 × max p∈P h liq 4-1/2 (p) d 2 h liq 4-1/2 dp 2 . Given a nondegenerate interval [p L , p U ] and range bounds [h L , h U ] , we construct a convex relaxation as the secant of a concave underestimator based on (12) as and a concave relaxation based on (13) as The max and min functions in (28) and (29) potentially tighten the relaxation in case we do not have an envelope anyway. The resulting relaxations are orders of magnitude tighter than those obtained by applying standard McCormick relaxations to the factorable representation of h liq 4-1/2 (p) (cf. Fig. 3a, c). Furthermore, the convex relaxation (28) obtained from the BB variant (12) by Hasan (2018) is significantly tighter than that obtained from the regular BB version (10) (cf. Fig. 3b, c). This is due to the fact that the function itself is almost concave in the sense that the maximum of the second derivative is much smaller in magnitude than the minimum.