**1. Introduction**

More now than ever, industrial processes are integrated to increase efficiency [1]. Process integration (PI) dates back to 1970 when PI was the response to the oil crisis, and, since then, this field has been a hot topic as it can be utilised to minimise energy and water usage [2], waste as well as emissions [3]. Since then, tools developed for Heat and Energy Integration have become essential elements of Process Design [4], including Energy Storage Systems [5]. Meanwhile, PI increases the efficiency of the technologies, and the increased complexity also complicates operations. Heat Integration also often results in more complex but less operable technologies [6]. Thus, a methodology to support the integrated system design must have a good qualitative and quantitative model that highlights the characteristics of the processes [7].

Although many works deal with the controllability [8] and observability [9] of complex systems, the connection between the complexity of the system and the difficulty of operations has not been examined in details. However, the analysis of the structural motifs of the building elements of complex systems and their models can highlight potentially useful information that can support the design and operation. The structure-relevant analysis of Heat Exchanger Networks (HENs) has already been proven to be beneficial, e.g., the resilience index (RI) quantifies the ability of a HEN to deal with disturbances [10], the controllability index measures the controllability of the HENs [11], and structural analysis can be applied to determine the locations of additional sensors that can reveal otherwise indistinguishable faults [12].

The connection between the structural properties of HENs designed by different methodologies and the operability of the system has ye<sup>t</sup> to be examined in details. The optimisation of HENs can be formalised in several ways. Algorithms can focus on the minimum number of matches [13], the Maximum Energy Recovery (MER) [14], the Minimum Energy-Capital cost [15], the minimum Total Annualised Cost (TAC) [16], the minimum number of exchangers [17], or, in the case of retrofit design, the minimum number of additional exchangers and the additional area of the exchangers or piping costs [18]. The goals listed above can be achieved by different algorithms, such as the Pinch methodology [19], dual-temperature approach method [20], pseudo-pinch [21], Supertargeting [14], State-Space approach [16], branch- and bound-based algorithms [22], or the application of Genetic Algorithm and Simulated Annealing (GA and SA, respectively) [23].

Although there are no systematic studies related to how optimisation algorithms affect the structural properties of HENs with regard to operability, there are some well-known relations. The pinch methodology determines the minimum temperature difference, Δ *Tmin*, which divides the HEN into two sub-networks: above the Pinch and below the pinch [19]. The method enables an MER design to be created, i.e., it minimises the energy required. The price of the MER design is the increased number/area and installation costs of heat exchangers. In contrast, optimisation based on the minimum number of matches decreases the number of heat exchangers in addition to the installation cost [24], and makes the operation easier [17]. As this approach has a negative impact on TAC and utility costs, other approaches aim to minimise them [25]. The tradeoff between the targets is continuously changing from one study to the next, and it is unequivocal that the tradeoff has a significant effect on the structure of the designed HEN, along with the controllability and observability of the system.

The dynamical characteristics of the HENs have increasingly become the focus of attention, such as controllability [26,27], observability [28], flexibility [29] or operability [30,31]. The trends mentioned above highlight that it is more critical to automatically qualify the operability-related dynamical properties of HENs based on structural information. For this purpose, HENs can be transformed into the networks of state variables [32], streams and matches [33], or networks of state-space representations [34].

Over the last five years, the system analysis based on the network has spread quickly, as the number and location of necessary actuators [8] and sensors [9] can be determined efficiently by utilising the structure-based maximum matching algorithm. The maximum matching algorithm is a hot topic in other fields such as pattern recognition [35] or machine vision [36] and method to determine maximum matching in large graphs also proposed [37]. The maximum matching algorithm is widely applied also in the fields of pattern recognition [35], machine vision [36], and it can scale up to handle large graphs [37]. With this methodology, the systems are analysed in terms of how does the correlation degree affect the controllability [38], robustness [39] and energy demand [40]. How structural motifs in the network influence the controllability of a transcription network [41], a human protein-protein interaction network [42] or cancer metabolic networks [43] has also been studied.

These promising applications sugges<sup>t</sup> that it could be beneficial to utilise this methodology to reveal the effect of the complexity of HENs in terms of their operability and to address the question concerning what kind of measures of network science can be applied to compare the complexity of HENs obtained by different design methodologies.

Section 2 presents three different network-based representations of HENs to explore motifs that have a significant effect on structural controllability and observability. To evaluate the structural controllability and observability of HENs, the minimum sets of driver and sensor nodes were generated. Additionally, an extended and more exhaustive version of the maximum matching-based approach was used [8] since it can underestimate the number of controllers required [44]. The generated set of driver nodes Was evaluated structurally through the number and location of driver nodes to determine control-relevant installation costs. The relative degree of the resulting dynamical system was also analysed to assess the sluggishness [45] and difficulty of the operation. A new methodology to decrease the relative degree is introduced in Section 2. Section 3 presents how the proposed approach can be used in the systematic analysis of almost 50 HEN design problems and how the developed measures can be used to compare different design methodologies.

#### **2. Network-Based Evaluation of the Complexity, Controllability and Observability of HENs**

The workflow of the proposed method is depicted in Figure 1. The method handles several goal-oriented representations of HENs to analyse their specific properties. As this work focuses on graph-based approaches, approaches based on temperature-enthalpy, heat-content or temperature-interval diagrams [46] are not discussed. As our goal was to study dynamic structural characteristics, this section presents the three most relevant models that can be used for such a purpose. After the introduction of the network-based representations, the systematic analysis of the extracted network is presented (see Figure 1). For this purpose, a MATLAB toolbox was developed, which is available at www.abonyilab.com.

**Figure 1.** Workflow of the proposed methodology.

#### *2.1. Network Representations of HENs*

Besides the widely applied Process Flow Diagrams (PFDs), HENs of hot and cold streams are represented in three different ways, as illustrated in Figure 2.

The first classical representation is the state-space (SS) approach (Figure 2a), which consists of the Distribution Network (DN) and the superstructure operator [34]. The DN determines how units are located on the streams, while the superstructure operator shows the interactions between the streams. The SS-based network representation (Figure 2b) can be defined as *GSS*(*VSS*, *ESS*) based on the *VSS* set of vertices that contain the inlet and outlet of the streams and the *ESS* edges that present the connections between the heat exchangers.

The second classical representation shows how streams are matched according to the units, namely heat exchangers, utility heaters and utility coolers (Figure 2c) [33]. The second SM-based network representation is almost the same as the classical one (Figure 2d), and it can be defined as *GSM*(*VSM*, *ESM*), where the nodes ( *VSM* = {*Vhs*, *Vcs*, *Vhp*, *Vcw*}) denote the sets of hot streams ( *Vhs*)d an cold streams ( *Vcs*), the high pressure steam ( *Vhp*) and the cold water ( *Vcw*). The edges (*ESM* = {*Ehe*, *Euh*, *Euc*} ⊆ *VSM* × *VSM*) represent the heat exchangers (*Ehe* ⊆ *Vhs* × *Vcs*), the utility heaters (*Euh* ⊆ *Vhp* × *Vcs*) and the utility coolers (*Euc* ⊆ *Vhs* × *Vcw*). The two partitions of the nodes are the hot streams (Process and utility, {*Vhs* ∪ *Vhp*}) and the cold streams (Process and utility, {*Vcs* ∪ *Vcw*}). The weights of the nodes represent the heat capacity of the streams, while the weights of the edges represent the heat load on the units.

**Figure 2.** Classical representations of HENs and their graph-based models. The representations are as follows: (**a**) state-space approach (SS); (**b**) network of state-space approach; (**c**) streams and matches (SM); (**d**) bipartite network representation of streams and matches; (**e**) grid diagram; and (**f**) network representation of state variables.

The third and most common classical representation of a HEN is the Grid diagram (GD) [6,47], where the streams are presented as they are connected through the heat exchangers and influenced by utility units (Figure 2e). The dynamics of the units, i.e., the heat exchangers, utility coolers and utility heaters can be described by a differential-algebraic system of equations (DAE) as [48]:

$$\frac{dT\_{ho}}{dt} = \frac{\upsilon\_{h}}{V\_{h}}(T\_{hi} - T\_{ho}) + \frac{lIA}{c\_{ph}\rho\_{h}V\_{h}}(T\_{co} - T\_{ho})\,\,\,\tag{1}$$

$$\frac{dT\_{co}}{dt} = \frac{v\_c}{V\_c}(T\_{ci} - T\_{co}) + \frac{lIA}{c\_{pc}\rho\_c V\_c}(T\_{ho} - T\_{co})\,\,\,\tag{2}$$

where *Thi*, *Tci*, *Tho* and *Tco* denote the temperatures of the hot input, cold input, hot output and cold output streams, respectively; *vh* and *vc* represent the flow rates of the cold and hot streams; *Vh* and *Vc* stand for the volumes of the hot and cold side tanks of the heat exchanger; *U* is the heat transfer coefficient; *A* denotes the heat transfer area of the heat exchanger; and *cp* and *ρ* are the specific heat and density of the streams.

In the above representation, the state variables are the temperatures of the outlet streams of the heat exchanger *<sup>x</sup>*(*t*)=[*Tho*, *Tco*]*<sup>T</sup>*, the temperature of the inlet streams are regarded as disturbances *d*(*t*)=[*Thi*, *Tci*]*<sup>T</sup>*, and *vc*(*t*) and *vh*(*t*) are time-varying parameters [48]. The dynamics of the utility units are similar. Utility coolers can be described by Equation (1) by considering the temperature of the cold water stream *Tco* as a controlled input variable. Analogously, utility heaters are modeled by Equation (2) where *Tho* denotes the controlled inputs. HENs based on these building blocks can be represented by linear state-space models in the general form of Equations (3) and (4).

$$
\dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}u + \Gamma d\tag{3}
$$

$$y = \mathbb{C}x + \mathbb{D}u\tag{4}$$

Recently, Liu et al. introduced a network science-based representation of dynamical systems and a methodology to determine the minimal input configuration that ensures the controllability of the system [8]. The third, state-space model-based network representation (Figure 2f) easily lends itself to the application of this methodology. The resultant network can be described as a graph *GDAE* = (*VDAE*, *EDAE*), where vertices represent the state variables *x*, while the edges are derived from the structure matrix of the state-transition matrix **A** as if **<sup>A</sup>***ij* = 0, then an edge from *xj* to *xi* exists. As a result of Equations (1) and (2), the simplest building blocks of HENs can be defined as the heat exchanger cells, the utility coolers and the utility heaters. Their grid diagram and state-space-based network representations can be seen in Figure 3.

**Figure 3.** The grid diagram and state-space network representations of: a heat exchanger cell (**<sup>a</sup>**,**d**); a utility cooler (**b**,**<sup>e</sup>**); and a utility heater (**<sup>c</sup>**,**f**). Red edges denote loops in the network representation, which create a strongly connected component in each elementary building block of HENs. In the case of the heat exchanger cell, the hot and cold sides of the exchanger belong to the same component even though they also have their own intrinsic dynamics, as in the case of utility units.

The SS-based network representation can be easily applied to HENs even when the streams are mixed and split. In the other two network representations, the handling of mixers and bypasses is less straightforward.

The number of vertices and edges of the SM-based network representation directly represents the streams and units, respectively. Loops and paths can be identified more easily in this network representation than in any other representations. This property is essential when the goal is not the MER design but to minimise the number of units that can be readily determined by the equation *Uun* = *Ns* + *L* − *S*, where *Uun* denotes the number of units, *Ns* the number of streams (|*VGD*|), *L* the number of independent loops and *S* the number of separate components in the network [33]. Loops and paths should be excluded from the designs, where the target is to minimise the number of units. The presence of splitters and mixers does not influence the property *Uun* = *Ns* + *L* − *S* when Pinch Partitioning is excluded [49]. Furthermore, with this representation the degree of freedom of the heat loads can also be determined as *H* − *Ns* + *S*, where *H* denotes the number of heat exchangers (|*Ehe*|) in the network.

The DAE-based network representation is suitable to analyse the operability of HENs. A system is controllable if it can be derived from an initial state to any desired final state over a finite period of time, while it is observable if any of the internal states can be reproduced by the knowledge of the initial state in addition to all the inputs and outputs [50]. The system is structurally controllable (observable) if the controllability (or observability) matrix, C = [**<sup>B</sup>**, **AB**, ... , **<sup>A</sup>***<sup>n</sup>*−1**B**] (or O = [**C***T*,(**CA**)*<sup>T</sup>*, ... ,(**CA***<sup>n</sup>*−<sup>1</sup>)*<sup>T</sup>*]*<sup>T</sup>*) is of full rank, *rank*(C) = *n* (or *rank*(O) = *n*) [51].

Liu et al. utilised the maximum matching algorithm on the DAE-based network representation to determine the minimum number of actuators (or sensors) and create a controllable (or observable) input (or output) configuration [8]. Nevertheless, how the network is represented is critical. If the adjacency matrix of the network is identical to the structure matrix of the state-transition matrix of the linear dynamical equation seen in Equation (3), then observability can be analysed, while in the reversed direction the controllability of the system can be investigated.

As each motif in the DAE-based network representation possesses a self-loop (coloured in red in Figure 3), i.e., the diagonal elements are non-zero values in the state-transition matrix, maximum matching selects these edges and leaves no unmatched node that can be a driver (or sensor) node. Since each HEN is constructed from these motifs, the approach is unusable in the case of HENs. Nevertheless, this method can be extended to be suitable for dynamical systems that exhibit this kind of behaviour. This method is referred to as the path-finding method [44] as it is also based on maximum matching, but following the maximum matching it searches for circles which can be cut off into Hamiltonian paths.

Since only the DAE network representation is detailed enough to determine the location of actuators and sensors, the relative degree can only be interpreted by this representation. To evaluate the "physical closeness" or "direct effect" of the control configuration of a nonlinear system with state equation *x*˙ = *f*(*x*) + ∑*m j*=1 *gj*(*x*)*uj* + ∑*<sup>p</sup> <sup>κ</sup>*=1 *<sup>w</sup>κ*(*x*)*dκ* and output equation *yi* = *hi*(*x*), relative order is introduced as a structural measure of the initial sluggishness of the response [45]. The relative degree can be determined by the standard Lie derivative. Therefore, for scalar field *hi*(*x*) and vector field *f*(*x*) it is defined as L*f hi*(*x*) = ∑*n <sup>l</sup>*=<sup>1</sup>(*∂h*(*x*)/*∂xl*)*fl*(*x*)), where *fl*(*x*) denotes the row element *l* of *f*(*x*). Higher order Lie derivatives are defined as L*<sup>k</sup> f hi*(*x*) = <sup>L</sup>*f*L*k*−<sup>1</sup> *f hi*(*x*), while mixed Lie derivatives <sup>L</sup>*gj* L*k*−<sup>1</sup> *f hi*(*x*) in an obvious way [45]. Then, relative degree *ri* is defined for output *yi* as the smallest integer for which [L*g*<sup>1</sup>L*ri*−<sup>1</sup> *f hi*(*x*)···L*gm* L*ri*−<sup>1</sup> *f hi*(*x*)] ≡ [0 ··· 0] with respect to input vector **u** if exists, otherwise *ri* = ∞. The relative order *rij* of output *yi* with respect to input *uj* is the smallest integer for which <sup>L</sup>*gj* L*rij*−<sup>1</sup> *f hi*(*x*) ≡ 0 if exists, otherwise *rij* = ∞. The relationship between *ri* and *rij* can be defined as *ri* = *min*(*ri*1,*ri*2, ... ,*rim*). Analogously to *rij*, the relative order *ρiκ* of output *yi* with respect to disturbance *dκ* can be defined as the smallest integer <sup>L</sup>*<sup>w</sup>κ*L*ρi<sup>κ</sup> f hi*(*x*) ≡ 0 if exists, otherwise *ρiκ* = ∞. For linear systems that are defined above, *rij* can be defined as the smallest integer for which *ci***A***rij*−<sup>1</sup>*bj* = 0, where *ci* is row *i* of **C** and *bj* is column *j* of **B**. The relative order *ρiκ* with respect to the disturbance *d* is the smallest integer for which *ci***A***ρiκ*<sup>−</sup>1*γκ* = 0, where *γκ* denotes column *κ* of Γ. By utilizing the DAE network representation, the relative degree *rij* can be defined as the shortest path *ij* from input *uj* to output *yi* minus one, *rij* = *ij* − 1. Analogously, the relative degree *ρiκ* can be defined as the geodesic path *iκ* from disturbance *dκ* to output *yi* minus one (*ρiκ* = *iκ* − 1). The visual representation of the relative degree can be seen in Figure 4.

As can be seen in Figure 4, the shaded area represents clusters of state variables that can be achieved by a relative degree that is smaller than *ri*. The operability of HENs can be improved by minimising the highest relative degree with the addition of actuators or sensors [52] as will be presented in the following subsection.

The analysis of the previously introduced network representations can highlight the structural properties of HENs. The developed network-based measures are summarised in Tables 1 and 2.


**Table 1.** Building elements of the studied network-based representations of HENs.


*Energies* **2019**, *12*,

**Figure 4.** Visual representation of relative degree in the DAE-based network representation. (**a**) The relative degree is equal to one if three actuators were placed in the network, while (**b**) it is three if only two actuators were assigned to the system.

#### *2.2. Operability-Focused Sensor and Controller Placement in HENs*

Two methods are used to extend the minimal input and output configurations [58]. The first approach utilises the set-covering method. Firstly, the allowed maximal relative degree *rmax* is defined. Secondly, the set of nodes *Wi* reachable from node *i* over a maximum of *rmax* steps was determined. In the formulation of the algorithm, U denotes the set of all the state variables, *C* the set of the actuators necessary to ensure structural controllability, and *O* the set of sensors necessary for structural observability. In the case of input (or output) configuration, *J* represents the set of necessary driver (or sensor) nodes, such that *P* is the set of state variables that is covered by *J*, namely *P* = <sup>∪</sup>*j*∈*<sup>J</sup>Wj*. The goal is to minimise *J* such that *P* = U and *C* ⊂ *J* (or *O* ⊂ *J*); moreover, *ru* ≤ *rmax*, and ∀*u* ∈ U. When initial input or output configurations are not given, then the method yields a global optimum for the problem. A greedy algorithm was applied to solve the set-covering problem [59].

The second approach uses two network-specific measures, namely the closeness and the node betweenness centrality measures. The closeness centrality of node *i* can be calculated using Equation (5).

$$\mathbb{C}\mathcal{L}(i) = \frac{N-1}{\sum\_{\vec{j}\neq\vec{i}\ i\neq\vec{j}}^{\prime}\ell\_{i,\vec{j}}} \tag{5}$$

Betweenness centrality calculates the number of geodesic paths that intercept node *i* (*<sup>σ</sup>st*(*i*)) and divides this by the number of all the geodesic paths (*<sup>σ</sup>st*) for each start *s* and target *t* node, such that *s* = *i* = *t*. The betweenness centrality of node *i* is shown in Equation (6).

$$B\boldsymbol{\varepsilon}(\boldsymbol{i}) = \sum\_{s:\boldsymbol{i}\neq\boldsymbol{i}\neq\boldsymbol{t}} \frac{\sigma\_{\rm st}(\boldsymbol{i})}{\sigma\_{\rm st}} \tag{6}$$

For each component that exceeds *rmax*, a node with the highest centrality measure is selected as an additional actuator (or sensor), i.e., *C* = *C* ∪ {*i* : max (*Cc*(*i*)*Bc*(*i*)), *i* ∈/ *C*} in the case of the input configuration and *O* = *O* ∪ {*i* : max (*Cc*(*i*)*Bc*(*i*)), *i* ∈/ *O*} in the case of the output configuration. The steps are repeated iteratively until all the components exhibit an order not in excessive of *rmax*.

For all HENs presented, the number of additional actuators and sensors required to manage the system with *rmax* was calculated by the set-covering method (global optimum), the retrofit set-covering method (using existing configurations) and the retrofit centrality-based method (using current configurations). In this analysis, the maximal order was determined as *rmax* = *dmax*/4, where *dmax* denotes the diameter of the network.

#### **3. Systematic Analysis of HENs**

#### *3.1. The Studied Benchmark Problems*

The well-known benchmark sets of Furman and Sahinidis (2004) [60], Chen et al. (2015) [61,62] and Grossmann (2017) [63] were used to study how the proposed measures can be used to evaluate HENs and compare different design methodologies of HENs. The benchmark problems and the applied methods are summarised in Table 3, while the notation of the utilised methods and their objective functions are shown in Table 4. These problem sets contain 48 problems and 23 methods, as well as 639 different HENs, of which 539 are unique. Fifty-three different measures, as summarised in Table 5, were generated for all the HENs.

**Table 3.** Benchmark HENs and their optimisation methods that are used during the analysis. As the heuristic methods, CP, CRR, CSH, FLPR, GP, GSH, LFM, LHM, LHM-LP, LRR, SS, WFG and WFM were applied to each problem. This set of methods is denoted as HEU referring to heuristic methods. The methods are introduced in Table 4, where abbreviations are also presented.






#### *3.2. Analysis of the 9sp-al1 Problem*

The 9sp-al1 benchmark problem was studied in detail to demonstrate the applicability of the proposed methodology. As can be seen in Figure 5, the solutions to the problem generated by the E, MER, EC and ECR methods significantly differed [14]. The first two solutions were based on the pinch point analysis to ensure Minimal Energy Requirement design, while the second two solutions were obtained by minimising energy and capital costs in the case of a new (EC) and retrofit (ECR) design.

The networks extracted from the HEN 9sp-al1-E solution (Figure 5a) can be seen in Figure 6. In DAE-based network representation (Figure 6b), the driver and sensor nodes are denoted by greenand orange-coloured symbols, respectively, while in SM-based network representation, a critical path is denoted by the colour green (Figure 6c).

**Figure 5.** Four designed HEN for the 9sp-al1 problem [14]: (**a**) existing network of an aromatic complexes in Europe; (**b**) maximum Energy Recovery design of the aromatic complexes; (**c**) grassroo<sup>t</sup> (new) design with optimisation for the Energy-Capital tradeoff; and (**d**) retrofit design with optimisation for the Energy-Capital tradeoff.

**Figure 6.** Networks extracted from the HEN 9sp-al1-E: (**a**) grid diagram of the HEN seen in Figure 5a; (**b**) DAE-based network representation, in which the green symbols represent the nodes where the input signal should be shared, and the orange symbols stand for the outputs to gran<sup>t</sup> structural controllability and observability; and (**c**) SM-based network representation, in which green edges show paths between high-pressure steam and cold water, i.e., between utility units which should be broken in order to reach the minimal number of units. In this HEN, no loop appears.

The clusters (communities) in the DAE-based network representation were also analysed by community detection algorithms [75]. The clusters of the HEN 9sp-al1-MER can be seen in Figure 7. Such clustering of the network can highlight useful information, e.g., it can show how the HEN is integrated and how the pinch isolates the clusters.

**Figure 7.** Communities of the DAE-based network representation of HEN 9sp-al1-MER were also identified. The colours of the nodes and their labels denote the community IDs. The highlighted edges are intercept Pinch.

All measures introduced in Table 5 were calculated for the presented HENs, and the results of the analysis are presented in Table 6.

As our analysis also aimed to identify structural parameters that have an impact on the dynamical properties of the HEN, firstly, the correlation between the proposed measures in these four examples was studied (see Figure 8). The results are discussed in the following section together with the analysis of all 639 HENs.

**Figure 8.** Correlations between the proposed network measures according to four HENs of the 9sp-al1 problem. The rows and columns were ordered based on similarities calculated based on the Spearman's rank correlation. The colours of the map represent the rank values.


**Table 6.** Calculated measures and their short descriptions.

#### *3.3. Results and Discussion of Systematic Correlation Analysis*

Figure 9 shows the similarity-based ordering of the 539 unique networks (on the rows) and the measures (on the columns). The colours of the map represent rankings of the networks based on the given measures. As can be seen, the map is clustered and on the right-hand side measures that are negatively correlated to the measures on the left-hand side are presented. This phenomenon originates from the calculation of the measures (some of these are calculated mostly as the reciprocal of the measures on the left-hand side). The similarities between the most important measures are visualised by a dendrogram in Figure 10.

Based on the analysis of the results, the following conclusions can be made:

Firstly, it is visible that the studied HENs were clustered and the members of such clusters exhibited similar dynamical properties. To confirm this, the Spearman distance for each network based on the measures was generated and visualised on a heat map, as can be seen in Figure 11. On the heat map, three groups can be easily distinguished, and these groups consisted mostly of problems similar in size.

Secondly, the number of unmatched nodes that were generated by maximum matching was equal to zero for all HENs. This was caused by the SCCs and the intrinsic dynamics introduced in Figure 3. Thus, additional driver and sensor nodes should be determined by methods such as path-finding [44].

**Figure 9.** Similarity-based ordering of the 539 HENs (rows) and developed measures (columns). The networks are ranked according to each measure, the similarities were calculated by Spearman's rank correlation and the colours of the map represent the rank values.

**Figure 10.** Correlations between the measures. Dendrogram presenting the distances between the measures for all unique HENs.

**Figure 11.** Spearman's rank correlation-based clustering of 539 unique HENs. The clusters of HENs can be easily detected. The method ordered the networks according to their similarities which were mainly determined by the complexity and the applied design methodology.

Thirdly, the number of the additional driver and sensor nodes correlated with the number of components in DAE-based representation, i.e., the number of sub-networks in the HEN. Thus, more interconnected streams demanded fewer driver and sensor nodes.

Fourthly, a smaller number of driver and sensor nodes resulted in an increased relative order, thus more complex operability.

Fifthly, the number of independent loops in the SM-based network representation correlated with the number of units in the HEN, and the average degree in SS- and SM-based network representations. The number of loops also correlated negatively with connectivity indexes. As the number of independent loops resulted in enhanced heat recovery, it attracted an increase in the number of heat exchangers too. Thus, a larger average degree meant that energy integration was greater, and more heat exchangers were utilised.

Finally, the comparison of the proposed methods showed that the results significantly varied from problem to problem. Only one correlation could be determined: the CRR, LHM and LRR algorithms solved the problem using more exchangers than other methods.
