Optimal Design of Heat-Integrated Water Allocation Networks

: Industrial operations consume energy and water in large quantities without accounting for potential economic and environmental burdens on future generations. Consumption of energy (mainly in the form of high pressure steam) and water (in the form of process water and cooling water) are essential to all processes and are strongly correlated, which requires development of systematic methodologies to address their interconnectivity. To this end, the subject of heat-integrated water allocation network design has received considerable attention within the research community in recent decades while further growth is expected due to imposed national and global regulations within the context of sustainable development. The overall mathematical model of these networks has a mixed-integer nonlinear programming formulation. As discussed in this work, proposed models in the literature have two main difﬁculties dealing with heat–water speciﬁcities, which result in complex formulations. These difﬁculties are addressed in this work through proposition of a novel nonlinear hyperstructure and a sequential solution strategy. The solution strategy is to solve three sub-problems sequentially and iteratively generate a set of potential solutions through the implementation of integer cut constraints. The novel mathematical approach also lends itself to an additional innovation for proposing multiple solutions balancing various performance indicators. This is exempliﬁed with both a literature test case and an industrial-scale problem. The proposed solutions address a variety of performance indicators which guides decision-makers toward selecting the most appropriate conﬁguration(s) among a large number of potential possibilities. Results exhibit that despite having a sequential solution strategy, better performance can be reached compared to previous approaches with the additional beneﬁt of providing many potential solutions for further consideration by decision-makers to select the best case-speciﬁc solution.


Introduction
Heat-Integrated Water Allocation Network (HIWAN) design has been extensively studied in the literature in recent decades [1][2][3][4][5]. Multiple water-and energy-related specifications of such systems have been addressed by the literature [5], including multi-contaminant problems [6] and interplant operations [7][8][9][10]. Several researchers have also proposed methodologies to incorporate other aspects of an industrial process in the pursuit of exploiting potential synergies between them. Ahmetović et al. [4] and Kermani et al. [5] have carried out systematic review of the available methodologies and have identified several gaps that have received no/little attention within the research community. These include: 1. providing systematic holistic methodologies addressing non-water processes [10][11][12][13], thermal utilities [10,14,15], and other energy conversion and heat recovery technologies [15], 2. considering multi-period operations and uncertainty analysis of HIWANs addressing variations of operating conditions, 3. proposing Multi-Criteria Decision Making (MCDM) tools addressing conflicting set of objectives and providing a set of optimal solutions to guide decision-makers in selecting the most appropriate configuration(s) [15], and 4. proposing efficient solution strategies to address the increasing difficulty of solving large-scale cases [16,17].
Kermani et al. [5] further provided a comprehensive analysis of the proposed HIWAN methodologies with particular focus on mathematical approaches. Mathematical methodologies are capable of addressing a variety of water specifications by constructing thorough superstructures that encompass all possible interconnections within the system. The overall problem comprises binary variables and nonlinear (including bilinear) constraints and hence is a non-convex Mixed-Integer Nonlinear Programming (MINLP) problem. To address the aforementioned gaps in a systematic manner using mathematical programming techniques, one of the main difficulties to consider is the interaction between the water network and the Heat Exchanger Network (HEN), that has to be well-understood and investigated. As discussed by Kermani et al. [5], this difficulty can be classified into two distinct problems of: • selecting water streams that participate in heat exchange, and • selecting the thermal state of these water thermal streams (i.e., whether they are hot or cold).
Depending on the origin or the destination of a water flow, a water stream can be categorized as freshwater, wastewater, inlet (to a process), outlet (of a process), recycling (among processes), or reuse (within a process). Any of these flows can be subject to thermal duties; however, considering all possibilities within the mathematical formulation introduces computational complexities. Therefore, a subset of these streams can be integrated with HEN. Kermani et al. [5] has shown that all proposed methodologies consider freshwater and wastewater interactions with HEN design while other types may or may not be addressed. Having set the water thermal streams, the second issue is to choose whether the stream is hot or cold as the HEN formulation is generally formulated by knowing these sets in advance. Within the literature, freshwater and wastewater streams are typically modeled as a succession of heat exchangers (modeled as cold streams) and splitters, and heat exchangers (modeled as hot streams) and mixers, respectively. For other types of streams, this procedure is carried out by defining two thermal streams and assigning binary variables and imposing constraints to ensure the activation of only one state. These interactions increase the complexity of the mathematical formulations which necessitates development of robust and efficient solution strategies. Kermani et al. [5] provided a detailed classification of the proposed solution strategies. The main classes can be categorized as follows: • Simultaneous: These approaches consider all variables and constraints in a single formulation, which are mainly formulated as MINLP models. Ahmetović and Kravanja [18] proposed such an approach by placing the potential heat exchangers on freshwater, wastewater, inlet (after mixer), and outlet (before splitter) streams. They solved the water network model to provide a good initialization point for the subsequent MINLP model. Ahmetović and Kravanja [19] further extended this approach by adding potential heat exchangers on recycling and reuse water streams and essentially extending the heat exchange possibilities to all water streams. Nonetheless, they have only modelled freshwater thermal streams as cold thermal streams and wastewater as hot thermal streams. • Decomposition: One technique to reduce the complexity of solving large MINLP models is via decomposing the model into two steps. Papalexandri and Pistikopoulos [20] and Leewongtanawit and Kim [21] proposed solution strategies by decomposing the overall MINLP model into a Mixed-Integer Linear Programming (MILP) and a Nonlinear Programming (NLP) model and solved the problem iteratively until no improvement was reached. All binary decision variables, i.e., the network configuration, were allocated to the MILP model. Despite the iterative nature of the solution strategies, only one solution was presented as the optimum one. • Sequential: Unlike decomposition approaches, sequential approaches lack iterations and hence are considered as uni-directional solution strategies. The problem can be formulated in two or more steps, where each step potentially solves one aspect of the problem. Liao et al. [22] proposed a two-step sequential solution strategy, by optimizing the potential thermal matches and the freshwater consumption in the first step and solving the MINLP model of HEN design in the second step.
Kermani et al. [14] used a more decomposed approach by proposing a three-step solution strategy, i.e., by minimizing the utility consumption in the first step and minimizing the number of thermal matches in the second step. The third step of their solution strategy is to apply the Pinch Design Method for HEN synthesis. Hong et al. [23] also proposed a three-step solution strategy by formulating an NLP model as the first step and then solving two subsequent MINLP models. The first step minimizes the freshwater consumption, however, lacks consideration of heat integration. Overall, none of the proposed sequential solution strategies addressed the possibility of having multiple solutions to the problem. • Transformation: Another technique to reduce the complexity of solving an MINLP model is via transforming and reformulating the problem. Yan et al. [24] proposed such an approach by essentially relaxing the integrality of binary variables and hence providing an NLP formulation. Despite this transformation, the proposed superstructure still possesses the same limitation of other ones, namely, the choice of water thermal streams and their state.
Furthermore, Ibrić et al. [16] and Wang et al. [17] have proposed several heuristic to simplify the overall superstructure via eliminating impractical and/or infeasible interconnections and have shown that these techniques can lessen the burdens of solving large problems. Jeżowski [1] have also highlighted the importance of sequential, decomposition, and combination of meta-heuristic techniques to facilitate finding a global or "good" local optimum.
In the design of HIWANs, there exist several sets of conflicting objective functions (e.g., operating cost vs total investment cost). Due to presence of these types of trade-offs, no single solution can be regarded as the best alternative with respect to all objective values. This necessitates developing multi-criteria optimization methodologies to provide a set of "good" solutions to the problem. Chen and Wang [2] and Kermani et al. [5] have highlighted the importance and benefits of multi-objective optimization techniques in HIWAN synthesis problems. Very few research works have considered these techniques in their methodologies [15,25,26]. Boix et al. [26] selected four objective functions defined as minimizing the freshwater consumption, energy consumption, number of heat exchangers, and number of water interconnections. In the first stage, they solved the problem of water allocation and thermal streams selection by combining the first two objectives using -constraint while fixing the value of the latter two objectives. Although they proposed a Pareto frontier of a set of solutions, they applied the classical TOPSIS method [27,28] to choose the best alternative among all in order to proceed to the second stage. In the second stage, having fixed the water allocation network and extracting the thermal streams, they formulated an MINLP model for minimizing the total HEN cost. Their approach, however, neglects the possibility of mixing and splitting among various water streams.
The innovative aspects included in this work can therefore be summarized as the novel decomposition and sequential solution strategy introduced, a systematic generation and presentation of multiple solutions for MCDM, and the application of these aspects on an industrial case study in addition to validation using a literature test case. The proposed method thus addresses the aforementioned difficulties noted as points 1, 3, and 4, above. The overall HIWAN problem is decomposed into three sub-models. The first two models were introduced in previous work, namely: targeting (problem P1) [14] and Heat Load Distribution (HLD) (problem P2) models [29,30]. Problem P1 is formulated as MILP and its solution provides minimum targets in thermal utilities and freshwater consumption together with potential water thermal streams and their states. Adding integer variables allows all water flows to participate in heat exchange duties. Having solved problem P1, problem P2 provides a feasible set of heat exchange matches among these thermal streams. This work adapts the NLP model of HEN design [31] by including a water network model to solve the total HIWAN synthesis problem as the third model (problem P3). This model is formulated using NLP with the objective of minimizing the annualized capital investment subject to utility targets (results of problem P1, set as upper bound) and heat exchange matches (results of problem P2). To address the trade-offs between resource consumption and investment cost and to further consider non-objective performance indicators, a novel solution strategy is proposed by which the three sub-problems are solved sequentially and iteratively to generate a set of potential solutions through the implementation of integer cut constraints. Furthermore, due to the presence of many other non-objective performance indicators in the system and due to their intertwined interactions, it would be ill-advised to present only a single solution while ignoring other "promising" ones. To this end, MCDM tools which allow decision-makers and practitioners to explore the set of solutions and select the best with respect to their own criteria are promising alternatives to the classical MCDM techniques. They further allow practitioners to understand why they are selecting a specific solution and how this solution's performance is compared to the rest. Among these tools, parallel coordinates visualization tool [32] is a powerful and practical way for evaluation and decision-making in high-dimensional spaces with applications in urban planning [33,34] and process integration [15]. This tool has been selected to visualize and further select promising alternatives similar to the approach presented by Kermani et al. [15]. The governing equations of problem P3 and the proposed solution strategy are described in Sections 3 and 4, respectively. The methodology is validated by well-established literature case studies in Section 5.

Problem Statement
Two sets of water unit operations (P WAN in , P WAN out ) are given. Each unit j in P WAN in requires water at temperature T j , with a maximum allowed inlet contamination level of c k,max j . Each unit i in P WAN out supplies water at temperature T i , with maximum allowed outlet contamination level of c k,max i . In addition, for mass-transfer water unit operations, a mass load of L u must be removed from unit u. Furthermore, sets of hot and cold non-water process streams (P H ,P C , respectively) are also considered. Thermal hot and cold utilities (U H ,U C , respectively) are also available within the system to close the energy balance. In addition, multiple freshwater sources and wastewater sinks (U WAN out ,U WAN in ) are given. The objective is to design a HIWAN that exhibits minimum Total Annualized Cost (TAC).

Mathematical Formulation
As introduced in Section 1 the overall MINLP formulation of HIWAN is decomposed in three sub-problems denoted P1, P2 and P3. In problem P1, a discretization technique is applied in which any water stream can be heated or cooled from its initial temperature to reach any other temperature in the water network [14]. By assigning binary variables to the existence of each thermal stream, the solution to problem P1 yields the minimum operating cost of the network by selecting potential thermal streams within the water network. This allows cold and hot streams to be selected for both source and sink units. Problem P2 minimizes the number of thermal matches among this set of streams. Appendix A and Appendix B provide the governing equations of problems P1 and P2, respectively. The HIWAN hyperstructure proposed for the third step (Figure 1), and the focus of this paper, is based on the NLP hyperstructure of Floudas and Ciric [31] (Figure 1a) and has been adapted to the HIWAN problem by applying the following modifications:

•
Final mixers are removed from superstructures of streams associated with water source units. Hence, the modified superstructure allows for stream splitting from the outlet of any water source unit heat exchanger. (Figure 1b) • Initial splitters are removed from superstructures of streams associated with water sink units. Hence, the modified superstructure allows for stream mixing at the heat exchanger inlet of sink water units. (Figure 1c) • Bypass streams are added between different superstructures associated to one source or one sink (shown as dashed arrows in Figure 1b,c). This mimics f B,k l,l in the HEN hyperstructure (Figure 1a), which characterizes the flow from the splitter following the match between k and l to the mixer preceding the match between k and l within the superstructure of stream k.

•
As opposed to the HEN hyperstructure for which heat loads of each match (Q k l ) are fixed at their optimal values from the HLD model, the proposed HIWAN hyperstructure relaxes this constraint; hence,Q k l can become zero indicating that the match is no longer necessary. For non-water thermal streams (cold or hot), the sum over hot or cold streams, respectively, is fixed to the stream heat load. In essence,Q k l for non-water thermal streams are allowed to deviate from their optimal values determined by the HLD model. These modifications produce the following novelties: • Freshwater and wastewater streams can have hot or cold thermal duties. This addition treats situations in which freshwater cooling or wastewater heating is required.

•
In most superstructures proposed in the literature, only one thermal stream is assigned to the outlet of a source water unit operation. As a result, no water stream can be split at the outlet of a heat exchanger for mixing with other streams. This issue is handled in the proposed HIWAN hyperstructure by removing the final mixers of source units. This issue is also treated for sink units by removing their initial splitters.   In order to present the constraints in a coherent way, similar set definitions of Floudas and Ciric [31] have been adapted to represent stream superstructures and matches: (1) Set MA is the set of thermal matches (k, l): MA = (k, l) | k ∈ HS and l ∈ CS and y hld k,l = 1 (1) (2) Set HCS is the set of all thermal streams: where HS and CS are the sets of hot and cold streams, respectively. (3) Set SM k is the set of streams, k , matching with stream k ∈ HCS [31], i.e., set of all streams, k , so that either of the two matches (k, k ) or (k, k ) is in the set of thermal matches, MA: (4) Set SM i,k is the set of streams matching with stream k ∈ S i of water unit i ∈ WAN, where, S i is the set of all water thermal streams related to water unit i: For mass flow variables, superscripts I, O, and B denote heat exchanger inlet, outlet, and bypass streams, respectively.

Objective Function
The objective function is the total network investment cost, C inv , which is the sum of investment costs over all heat exchangers: where, Q i,j is the heat load of match (i, j) which, as stated before, is not fixed at its optimal value from problem P2 but rather defined as a variable. This makes the objective function non-convex. Parameters α 1 , α 2 , and β are fixed, area, and exponential cost parameters for heat exchangers. The LMTD i,j is the logarithmic mean temperature difference of match (i, j) defined as: where, ∆T1 i,j and ∆T2 i,j are the temperature differences at the two ends of heat exchanger of match (i, j). Equation (6) is simplified using Chen's approximation [35]: Figure 1 illustrates the basics of the proposed hyperstructure for a water source, a water sink, and a non-water thermal stream which are governed by the equations presented in this sub-section.

Water Network Mass and Energy Balances
(1) Mass balance at the initial splitter of source unit i ∈ WAN out : where,ṁ i is the flow rate of source unit i,ṁ k,I i,l is the flow rate from source unit i towards the inlet of heat match (k, l) in the superstructure of stream k,ṁ i,j is the flow rate between source unit i and sink unit j, andṁ k,I i,j,l is the flow rate from source unit i towards the inlet of heat match (k, l) in the superstructure of stream k of sink unit j.
(2) Mass balance at the splitter following a match of source unit i ∈ WAN out (k ∈ S i , l ∈ SM i,k ): where,ṁ ex,k i,l is the flow rate through heat exchanger of match (k, l) in the superstructure of stream k of source unit i,ṁ B,k ,k i,l ,l is the flow rate from the outlet of heat exchanger of match (k, l) in the superstructure of stream k towards the heat exchanger inlet of match (k , l ) in the superstructure of stream k where both heat exchangers belong to source unit i,ṁ k,O i,l,j is the flow rate towards sink unit j, andṁ k,O,k ,I i,l,j,l is the flow rate towards the inlet of heat exchanger of match (k , l ) in the superstructure of sink unit j.
(3) Mass balance at the mixer preceding a match of source unit i ∈ WAN out (k ∈ S i , l ∈ SM i,k ): where,ṁ ex,k i,l is the flow rate through heat exchanger of match (k, l) in the superstructure of stream k of source unit i,ṁ B,k,k i,l,l is the flow rate from the outlet of heat exchanger of match (k , l ) in the superstructure of stream k towards the heat exchanger inlet of match (k, l) in the superstructure of stream k where both heat exchangers belong to source unit i, andṁ k,I i,l is the flow rate from source unit i towards the inlet of heat match (k, l) in the superstructure of stream k. (4) Energy balance at the mixer preceding a match of source unit i ∈ WAN out (k ∈ S i , l ∈ SM i,k ): where t ex,k,I is the temperature at the inlet (outlet) of heat exchanger of match (k, l) (k , l ) in the superstructure of stream k(k ), and T i is the temperature of source unit i.
where,ṁ j is the flow rate of sink unit j,ṁ k,O j,l is the flow rate from the outlet of heat match (k, l) in the superstructure of stream k towards sink unit j,ṁ i,j is the flow rate between source unit i and sink unit j, andṁ k,O i,l,j is the flow rate from the outlet of heat match (k, l) in the superstructure of stream k of source unit i towards sink unit j. (6) Mass balance at the splitter following a match of sink unit j ∈ WAN out (k ∈ S j , l ∈ SM j,k ): where,ṁ ex,k j,l is the flow rate through heat exchanger of match (k, l) in the superstructure of stream k of sink unit j andṁ B,k ,k j,l ,l is the flow rate from the outlet of heat exchanger of match (k, l) in the superstructure of stream k towards the heat exchanger inlet of match (k , l ) in the superstructure of stream k where both heat exchangers belong to sink unit j. (7) Mass balance at the mixer preceding a match of sink unit j ∈ WAN out (k ∈ S j , l ∈ SM j,k ): where,ṁ B,k,k j,l,l is the flow rate from the outlet of heat exchanger of match (k , l ) in the superstructure of stream k towards the heat exchanger inlet of match (k, l) in the superstructure of stream k where both heat exchangers belong to sink unit j,ṁ k,I i,j,l is the flow rate from source unit i towards the inlet of heat match (k, l) in the superstructure of stream k of sink unit j, andṁ k ,O,k,I i,l ,j,l is the flow rate from the outlet of heat match (k , l ) in the superstructure of stream k of source unit i towards the inlet of heat exchanger of match (k, l) in the superstructure of stream k of sink unit j. (8) Energy balance at the mixer preceding a match of sink unit j ∈ WAN out (k ∈ S j , l ∈ SM j,k ): is the temperature at the inlet (outlet) of heat exchanger of match (k, l) ((k , l )) in the superstructure of stream k (k ) and T j is the temperature of sink unit j. (9) Energy balance in matches of stream k ∈ S u of water unit u ∈ WAN:

Contamination Constraints
(10) Contamination equality constraint at the final mixer of sink unit j ∈ WAN in for contaminant k ∈ C: where c k,in j is the contamination level of contaminant k at the inlet of sink unit j, and c k,out i is the contamination level of contaminant k at the outlet of source unit i. (11) Contamination at the inlet of unit j ∈ WAN in should be less than a maximum allowed value: whereṁ j =ṁ i and WAN u,in (WAN u,out ) is the set of sinks (sources) associated to mass transfer unit u. (13) Outlet contamination of unit i ∈ WAN out should be less than a maximum allowed value:

Mass and Energy Balances in Non-Water Streams Superstructure
The mass and energy balances of non-water thermal streams follow the same constraints of the HEN synthesis problem of Floudas and Ciric [31] (Figure 1a introduces the full set of related variables and parameters): (14) Mass balance at the initial splitter of stream k ∈ HCS: where, F k is the heat capacity flow rate of stream k and f k,I l is the heat capacity flow rate of stream k towards the inlet of heat exchanger of match (k, l)/(l, k) in the superstructure of stream k. (15) Mass balance at the mixer preceding a heat exchanger of stream k ∈ HCS: (16) Mass balance at the splitter following a heat exchange match of stream k ∈ HCS: (17) Energy balance at the mixer preceding a heat exchange match of stream k ∈ HCS: (18) Energy balance of heat exchange match (k, l) ∈ MA (as mentioned previously, heat load is defined as a variable in the HIWAN model): Energy balance of hot and cold non-water thermal streams: This constraint holds for all thermal streams k ∈ HCS that do not belong to any water unit operation. Q k is the heat load of non-water thermal stream k.

Temperature Difference Constraints
(20) Temperature difference at the two ends of heat exchanger (k, l) ∈ MA: (21) Minimum approach temperature difference:

Solution Strategy
As discussed in Section 3, the HIWAN synthesis problem is decomposed into three sub-problems: Step 1-Targeting Model (problem P1): At this step, a value of HRAT is assumed for water thermal streams. For non-water thermal streams, this value depends on the nature of the stream and problem requirements or can be assigned based on general assumptions. A solution of this model provides the utility targets of the HIWAN synthesis problem and a list of water thermal streams with their corresponding heat loads. Moreover, a water allocation network will be generated satisfying the temperature and contamination constraints (Appendix A) [14].
Step 2-Heat Load Distribution Model (problem P2): The input of this model is the list of thermal streams with their corresponding heat loads from Step 1. A solution of this model provides a set of potential thermal matches with their corresponding heat loads which exhibits the minimum number of heat exchangers. Several solutions may exist with the minimum number of heat exchangers. This is addressed in the proposed solution strategy by applying Integer-Cut Constraints (ICCs). The water allocation network is fixed at its solution from Step 1 (Equation (A2)) for this stage.
Step 3-HIWAN Hyperstructure Model (problem P3-NLP): The model takes the water data and utility targets from Step 1, while the thermal matches and their heat loads are provided from Step 2. Since the water allocation network superstructure is embedded in this model, the water allocation network from Step 1 can be changed but provides an initial, feasible solution. Additionally, heat loads of thermal matches are relaxed at this step; hence, heat exchanger flows can reach zero, indicating that the optimized water network flows and temperatures do not require the existence of this exchanger. The variables of this problem, i.e., flows, temperatures, and heat loads of all heat exchangers, are initialized by solving the HEN hyperstructure of Floudas and Ciric [31] (model P3 init -NLP). It should be highlighted that a solution to P3 init -NLP is a feasible solution to the HIWAN synthesis problem.
Solving the above-mentioned steps in sequence provides, at most, two solutions to the HIWAN synthesis problem: the solutions of P3 init -NLP and P3-NLP. Compared to simultaneous approaches, sequential approaches are easier to solve due to reduced problem size at each step. Nevertheless, similar drawbacks observed for sequential solution strategies that are applied in HEN synthesis problems [31] are also valid in this case: 1. The trade-off between operating cost and investment cost (i.e., number of matches and area) cannot be fully captured in sequential solution strategies. 2. The nature of the HIWAN synthesis problem is non-convex and hence global optimality cannot be guaranteed. 3. Problem P2-MILP provides only one feasible solution providing the minimum number of heat exchangers. Therefore, sequentially minimizing the number of heat exchangers followed by HEN synthesis does not guarantee a HIWAN exhibiting the lowest total cost.
To address these drawbacks, an iterative sequential solution strategy is proposed by applying an integer cut constraint on problems P1 and P2. In addition, the problem is solved for different values of HRAT to address the pinch point decomposition of temperature intervals and to get different utility targets with the goal of achieving lower investment costs. Applying integer cut constraints at Steps 1 and 2 provides a double benefit: By assuming the overall synthesis problem encompasses all water thermal streams and states (hot and cold), problems P1 and P2 limit the search space to a specific region, i.e., a specific set of potential thermal streams with their corresponding states and a specific set of thermal matches. This allows for a smaller problem size at Step 3 and hence provides guidance towards a good solution. Applying an integer cut constraint generates several "reduced-size" problems and increases the opportunities for identifying good solutions. By applying, for example, N P1 icc and N P2 icc integer cuts at Steps 1 and 2, respectively, (N P1 icc + 1)(N P2 icc + 1) solutions to the HIWAN synthesis problem can be generated. These solutions can be ranked based on different key performance indicators and thus guide the decision-making process toward a "good" solution.
This approach can be regarded as iteratively optimizing the overall MINLP model within different regions of the global search space (as imposed by the formulations of problems P1 and P2 at each iteration). Figure 2 illustrates the proposed solution strategy. It should be noted that, similar to other methodologies, the choice of solver and its options affects the path toward a solution.

Illustrative Example
The four-water-unit-operation, single-contaminant case study of Savulescu and Smith [36] was selected to demonstrate the generation of the proposed HIWAN hyperstructure. First, problem P1 is solved providing minimum freshwater, hot, and cold utility consumption of 90 kg/s, 3780 kW, and 0 kW, respectively. Problem P1 also provides a feasible water allocation scheme. Given the list of thermal streams, problem P2 is solved which results in seven matches (Figure 3a). There are several solutions to problem P2 with an objective value of seven. For the sake of illustration, the fourth solution (by applying three integer cuts) is shown in Figure 3a and will be discussed hereafter. The CPLEX solver [37] was used for both MILP problems (Table 1). Figure 3b presents the HIWAN hyperstructure constructed based on the results provided in Figure 3a. The HIWAN hyperstructure generation follows the same approach as the HEN hyperstructure [31], as shown in Figure 3a. The inlet of unit u 2 has one cold thermal stream which has one match with hot utility u H , one match with wastewater unit u ww and one match with its outlet; hence, one superstructure with three matches is defined for this stream. The Wastewater unit u ww has two cold streams with one and three matches, respectively. This results in a hyperstructure consisting of two stream superstructures as illustrated in Figure 3b. As presented in Figure 1 for the proposed hyperstructure, mass flows can be split from source unit heat exchanger outlets and mixed at sink unit heat exchanger inlets. This has been illustrated in Figure 3b with thicker, colored arrows.
HIWAN synthesis is solved for problem P3 as well as for the HEN model of Floudas and Ciric [31]. SNOPT solver [38] is used to solve both models (Table 1). Its algorithm is based on sequential quadratic programming (SQP) and its solutions are locally optimal. It should be noted that problem P3 is non-convex and hence the global optimality of the solutions cannot be guaranteed. The solution of the NLP hyperstructure of Floudas and Ciric [31] is shown in Figure 4a. Since this solution is based on the results of HLD, all seven heat exchangers are active. The area of heat exchangers 3 (55.4 m 2 ) and 7 (196.4 m 2 ) are relatively small compared to the others ( Table 2). This approach, however, does not account for the possibility of stream mixing and splitting. The HIWAN hyperstructure addresses these possibilities by embedding the water allocation network within the HEN hyperstructure. Figure 4b presents the results of the proposed hyperstructure. By considering additional possibilities for stream mixing and splitting, hence increased possibility of non-isothermal mixing, three out of the seven targeted heat exchangers are determined not to be required. The flows in bold illustrate the new possibilities embedded in the proposed HIWAN hyperstructure. Consequently the HEN cost is reduced by 22.47 % (290,799 compared to 375,089 kUSD/yr). Table 2 provides the temperatures, flows and heat loads associated to heat exchangers in both cases.

Validation and Discussion
Two test cases from the literature were selected to validate the mathematical formulation and the proposed iterative sequential solution strategy. The cost data are presented in Table 3. The maximum number of integer cuts was limited to 50, though one sub-problem (P1 or P2) may become infeasible before reaching its respective limit, at which point the algorithm terminates the corresponding integer cut loop. In combination with the application of integer cut constraints, small minimum sizes were introduced for all water units in problem P1 (0.05 kg/s). Similarly, for problem P2, the minimum allowed heat exchange between hot stream i and cold stream j in pinch interval l is set at a small positive number (Q min i,j,l = 0.1 kW in Equation (A6)). These two additions prevent generation of multiple identical solutions. Furthermore, to investigate the trade-off between investment cost and operating cost, the HRAT in problem P1 was varied between 1 • C and 10 • C with step of 1 • C (N HRAT = 10). These values are passed to problem P3 to be used as the value of ∆T min . Test case I is a threshold problem and hence the value of HRAT does not affect the utility consumption. This, however, becomes a decisive factor for test case II where both hot and cold utilities are required. In this section, only solutions of the HIWAN hyperstructure are presented. A parallel coordinate visualization tool is used to illustrate all solutions and selected key performance indicators. It should be noted that test case I was intentionally designed by the scientific community to validate their mathematical formulations and solution strategies. Due to the complexity of solving non-convex MINLP problems that arise in HIWAN synthesis, and solver tendency toward local optimality in such problems, it can be argued that these test cases are not easily solved without carefully manipulating solvers and their intrinsic options to reach a single,´´good" solution. Thus, the main goal and contribution of the proposed solution strategy in this work is to solve the HIWAN problem by first reducing the search space, then generating many solutions to the problem. To this end, many solutions can be generated for problem P3 which can be seen as approaching an optimal point from different starting points (solutions of problems P1 and P2).

Test Case I-Single-Contaminant Problem
This test case in its optimal configuration [14,40,41] exhibits 90 kg/s of freshwater consumption and 3780 kW of hot utility consumption (minimum targets, 2397.0 kUSD/yr in operating cost). Table A1 provides the assumptions and solver options for this test case. Out of 26,010 potential solutions, 6462 solutions exist to problem P2-this arises from the fact that problem P1 converged for all integer cuts, while problem P2 was found to be infeasible for low HRAT values after relatively few integer cuts. 3611 solutions exist to problem P3. They exhibit a wide spectrum of solutions with HEN cost (the objective function in problem P3) ranging from 255.1 kUSD/yr (the minimum value within the literature [5]) to 660.9 kUSD/yr. Among these, 3506 solutions exhibit HEN cost in the range reported in the literature [5] (the maximum observed cost is 413.0 kUSD/yr). To visualize this set of solutions, several filters were imposed to further reduce the number of solutions: • C HEN ≤ 300 kUSD/yr; • N m s ≤ 15, maximum value reported in the literature for this test case [5]; • N mixer ≤ 13, maximum value reported in the literature for this test case [5]; • N N I M mixer ≤ 10, maximum value reported in the literature for this test case in [5]; Figure 5 presents 2247 solutions which remain after filtering (the color spectrum is based on HEN cost). Selected results from the literature [5] with HEN cost lower than 300 kUSD/yr are also plotted in this figure in black bold. The number of heat exchangers (N HE ) varies between two and five, where two is the lowest number of heat exchangers reported for this test case by the proposed solution strategy. Figure 5a illustrates an inverse correlation between the number of heat exchangers and their total area. The total area (A total HEN ) varies between 3350 m 2 and 7090.7 m 2 . The total area of heat exchangers reported in the literature varies between 3111 m 2 (Wan Alwi et al. [42]) and 6300 m 2 . However, it should be noted that the HEN design reported by Wan Alwi et al. [42] is infeasible due to a violation of ∆T min .

Test Case II-Simplified Industrial Case Study
The simplified industrial case study introduced by Kermani et al. [14] is revisited here. Similar to previous case studies, problems P1 and P2 are solved for 50 integer cuts.  (Table A2); Overall, 2053 solutions satisfy the above criteria which are visualized in Figure 6a (the color spectrum is based on HEN cost). The number of heat exchangers is limited to being 8 or 9, while the area varies between 2100 and 2750 m 2 . It can be observed from this figure that, for a given number of mixers and mass streams, the HEN cost varies over the entire range which enunciates that minimizing HEN cost and providing a single solution would not satisfy other important criteria. The lowest total cost reached by the proposed solution strategy is 960.6 kUSD/yr which is less than the lowest reported in the literature (971.4 kUSD/yr [13]). It should also be highlighted that the solution proposed by Ibrić et al. [13] considered the optimization of tank temperatures, while these temperatures are fixed for the analysis presented herein. Moreover, all proposed solutions in this work exhibit 11% lower freshwater consumption than the solutions reported in the literature. Figure 6b illustrates all of the solutions (240) proposed by the sequential solution strategy which exhibit lower total cost than 971.4 kUSD/yr. Figure 6b shows that for a given number of heat exchangers (8 for these solutions), increasing the number of non-isothermal mixing points lowers the total cost. Figure 7 presents one of the optimal HIWAN designs. Table 4 also provides a comparison between problems P2 and P3 in terms of thermal matches and their heat loads for the selected solutions. As described in Section 5, the HLD model and consequently HEN hyperstructure model do not consider stream mixing. The proposed HIWAN model takes these possibilities into account while optimizing the temperatures and flows of heat exchangers. This increases the possibility of replacing a heat exchanger by non-isothermal mixing which led to the elimination of the two smallest heat exchangers as highlighted in Table 4. In both cases, the total heat load of each non-water thermal stream remained constant, while the flow and temperature of water thermal streams adapted to address the shift in the heat load of non-water thermal streams.

Conclusions
Two main complexities have been observed within the HIWAN methodologies in the literature. One difficulty is selecting the state of water thermal streams in the HEN design, as this choice affects the overall mathematical formulation in MINLP models. This, in turn, introduces the second difficulty of selecting thermal matches with undefined states of water thermal streams; thus, resulting in a complex MINLP formulation. These two complexities have been addressed by solving problems P1 (targeting model) and P2 (HLD) in sequence prior to the final HIWAN design which identify the states and potential thermal matches, respectively.
A novel HIWAN hyperstructure with NLP formulation has been proposed in this work for the final design of HIWANs which exhibits higher degrees of non-isothermal mixing and lower investment costs. In addition, an iterative sequential solution strategy is proposed in this work by solving the sub-problems with integer cut constraints applied in both problems P1 and P2. As discussed, this approach can be viewed as a technique of reducing the search space at each iteration to a set of potential water thermal streams and a set of potential matches followed by solving the HIWAN synthesis problem for the reduced space.
Applying the proposed methodology on two test cases from the literature illustrated that it is not only able to reach the minimum total system cost (the only objective function used in the literature), but also generates a set of alternative "good" solutions exhibiting various results with respect to other key performance indicators. The first test case, a well-known problem in the literature, yielded 375 solutions with the same (minimum) value of the single objective used. This highlights the plethora of solutions which could be selected by decision-makers to achieve cost targets but with distinguishing features in non-objective or non-quantifiable goals. The second test case yielded 240 alternatives with lower cost than the minimum reported in the literature with an additional benefit of reducing the freshwater consumption by 11%. Similar to the first test case, this reinforces the approach by proposing many "good" solutions which can then be evaluated by engineers or decision-makers based on their other merits. Complexities related to application of the proposed solution strategy and associated MCDM tool on real industrial cases, as the one in a kraft pulp process [15], is limited to characterizing the industrial plant, its specifications, and collecting relevant data.
The proposed approach is intended for the grassroot design of HIWANs and cannot be directly applied to retrofit problems; however, as shown for problem P3, the operating conditions of a water allocation network (i.e., temperatures, flow rates, and heat exchange areas) can be optimized for a given set of thermal matches. Thus, for retrofit cases in which additional heat exchangers are not envisaged, the proposed superstructure can be applied at the current level. For investigating the addition of new heat exchangers, problem P2 can be reformulated by adding new constraints. Multi-period operation of HIWANs (by considering daily/seasonal variations in process operating conditions, temperatures of freshwater, etc.) and uncertainty analysis of HIWANs (to find the most resilient networks given system uncertainties, including costs and operating conditions) should be addressed in future work.
The choice of HEN cost as the main objective in the literature cases should not be the sole measure of optimality, as the basic functions for heat exchanger cost reflect average, observed correlations for overall networks and hence should not be used as the decisive factor in proposing a final design. The proposed approach in this work overcomes this limitation by generating and analyzing a set of solutions and applying multi-criteria optimization to propose competing solutions considering multiple decision criteria. Acknowledgments: This research project is financially supported by the Swiss Innovation Agency, Innosuisse, and is part of the Swiss Competence Center for Energy Research SCCER EIP and has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 818011.

Conflicts of Interest:
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Appendix A. Problem P1-Targeting
Problem P1 is formulated as an MILP model following the methodology proposed by [14]. The objective function is defined as minimizing the total annualized cost, C tac , of the system and is given by whereṁ f w andṁ ww are the flow rates of freshwater source f w and wastewater sink ww, respectively, C f w and C ww are the cost of freshwater f w and wastewater ww (for treatment), respectively,ḟ H u , q H u , and C H u (ḟ C u , q C u , and C C u ) are the flow rate, unit heat load and cost of hot (cold) utility u, respectively. The investment cost, C inv is annualized given the interest rate, irr and operating lifetime, n year , of the system, I WAN f ,i,T and I WAN p,i,T are the linearized coefficients of the penalizing factor. The optimization is subject to mass, energy, and contamination constraints. For details of the calculations of the penalizing cost and the constraints refer to Kermani et al. [14], Maréchal and Kalitventzeff [43], Kermani et al. [44]. The results provide water and energy targets as well as a feasible water network with a set of water exchanges and heat loads of all thermal streams. Having solved problem P1, the second problem can be formulated to find the promising thermal connections.

Appendix B. Problem P2-Heat Load Distribution
The set of all hot (HS hld ) and cold (CS hld ) thermal streams in the network can be listed at this stage knowing the solution of problem P1. In addition, several pinch points will be identified that divide the temperature range into several subnetworks (SN). At this stage, the HLD model will be applied to each subnetwork to find the potential thermal matches. HLD is formulated as an MILP model, based on minimizing the total number of heat exchange matches among all hot and cold streams [29,30]. Similar notations introduced by Papoulias and Grossmann [29] are used here for familiarity purposes.
It is important to note that the result of HLD is one feasible solution among many possible solutions that exhibit the minimum number of matches. Moreover, to keep the targeted thermal utilities intact, no thermal exchange is allowed crossing the pinch point, i.e., matches are limited between thermal streams within each subnetwork. For any subnetwork l ∈ SN, one can readily define the set of hot and cold thermal streams within each subnetwork (HS hld l ,CS hld l , respectively). HS hld l,r is defined as the set of hot streams within subnetwork l that are present in temperature interval r and above, while CS hld l,r is defined as the set of cold streams within subnetwork l that are present in temperature interval r. Q i,j,r ≥ 0 is defined as the heat exchanged between hot stream i and cold stream j in interval r ∈ SN l . Upper and lower bounds can be defined for each match (Q min i,j,l ,Q max i,j,l ) to better incorporate practical infeasibilities. To avoid violating thermal utility targets, one should follow the same constraints as in problem P1. A heat cascade formulation should be completed for each hot stream (R i,r ≥ 0). By defining binary variables y i,j,l ∈ {0, 1} to denote the existence of a match between thermal streams within each subnetwork l, the HLD problem will minimize the number of heat exchangers. Ranking factors (p i,j,l ) can be defined to favor or avoid certain matches:

Appendix C. Assumptions and Solver Options for Test Cases
Both test cases are run on a Windows machine with Intel(R) Core(TM) i7-4600U CPU 2.10CHz 2.7 GHz and 8.00 GB RAM. The models are written in AMPL programming language. Table A1 provides the assumptions and solver options for test cases. Relative tolerance for optimizing integer variables: stop if abs((best bound) -(best integer)) < mipgap * (1 + abs(best bound)). integrality 10 −10 10 −4 Amount by which an integer variable can differ from the nearest integer and still be considered feasible. P2 (CPLEX) mipgap 10 −7 10 −5 See above integrality 10 −10 10 −8 See above P3 (initialization) (SNOPT) Problem P3 is solved by minimizing the sum of flows between heat exchangers meminc 1 1 Increment to minimum memory allocation iterations 10 6 10 6 Minor iteration limit feas_tol 10 −6 10 −4 Minor feasibility tolerance for all variables and linear constraints P3 and P3 HIWAN hyperstructure (SNOPT) meminc = 1 iterations = 10 6 feas_tol = 10 −4 (see above) HRAT = ∆T min for each iteration. Enthalpy values of thermal stream are rounded for HEN design to 10 −6 . Optimality gap of heat balance constraints (HIWAN) at mixing points is set at 0.01%. For problem P1, sum of mass exchanges is added to the objective function (Equation (A1)) to reduce the number of mass exchanges while minimizing the total cost.