A Mixed-Integer Linear Programming Formulation for Optimizing Multi-Scale Material and Energy Integration

1. Introduction

Increasing demand for energy and raw materials necessitate commensurate improvements in material and energy efficiency in major process industries to slow the depletion of natural resources. Process integration techniques, often proposed out of necessity from industries to respond to scarcity or market conditions, have been developed to analyse process efficiency and propose solutions. Such methods are becoming increasingly important and complex due to additional objectives such as emission reductions and job creation/retention as well as additional economic, political, or environmental constraints.

Mixed-integer linear programming (MILP) is often used for system analysis and optimization as it presents a flexible and powerful method for solving large, complex problems such as the case with industrial symbiosis and process integration. Several limitations have also lead researchers toward bi-level optimization (Aviso et al., 2010), non-linear deterministic optimization (Kantor et al., 2015), or stochastic methods (Gonela et al., 2015) to address certain aspects of processing. Such extensions often require other simplifications and the development of heuristic methods for limiting the solution space and determining solution globality.

Process integration methods can be extended from the unit process or plant level to include exchange possibilities between plants to discover new opportunities for process efficiency improvements which address modern requirements for meeting environmental and social constraints while maintaining economic competitiveness. Such opportunities have traditionally been considered on an individual basis without decision-support from optimization tools or advanced mathematical programming methods. The opportunity to optimize potentially interconnected, complex processing systems by applying modern optimization methods is a burgeoning research area, though a holistic framework to identify optimal industrial symbiosis solutions with a mixed-integer programming approach is still lacking. Identifying solutions for a given situation is often completed manually by professionals in the area of industrial symbiosis but solutions are likely to be sub-optimal, not taking advantage of many opportunities which may exist for integrating material and energy streams between plants. The opportunities for symbiosis solutions have been identified as promising by several authors (Jacobsen, 2006; Hashimoto et al., 2010; Kantor et al., 2015), though deep analysis of the promising solutions is required to ensure that site specificities are considered and that detailed design aspects such as materials and safety are also considered.

Connections between plants, typical symbiosis solutions, can also be augmented by a superstructure of new units which provide services within or between industries, notable examples being polymerization processes to provide low-temperature heat, organic Rankine cycles to convert waste heat into electricity, or separation processes to split mixed streams into useful fractions for one or more industrial applications. Using a comprehensive MILP optimization framework identifies the best connections between plants, utilities and modifications to achieve maximum resource efficiency, minimum environmental impact, minimum total cost, or other objectives.

This work addresses the need for comprehensive optimization frameworks. Section 2 presents a review of the current state and previous work to provide comprehensive treatment of large-scale integration problems and provides motivation for this work. Section 3 presents the MILP formulation of the problem to provide such a framework and section 4 presents an example case study and its solution. The future improvements of the formulation are briefly discussed in section 5 and conclusions of this work are drawn in section 6.

2. Background

Process Integration (PI) is a field of research which progresses beyond heat integration to identify processes and utilities which can be optimally integrated to achieve higher efficiencies at lower cost than conventional processes by also including material and energy flows. In classical heat integration, following the concepts of Pinch Analysis (PA), the system is divided into two subsystems, namely heat source and heat sink (Klemeš and Kravanja, 2013). The reduction in heating utility requirement results in an equal cooling utility requirement reduction. Subsequently, CO₂ emissions are also reduced, which is a primary target of modern Heat Integration (HI) and PI (Dhole and Linnhoff, 1993).

The earliest work on HI is based on PA. Edward Hohmann (1971) introduced the concept of recovery pinch, which is accepted as the first systematic approach to obtain heating and cooling targets by using the feasibility table. Linnhoff (1972) and Linnhoff and Hindmarsh (1983) proposed PA using the problem table method, which has been widely used by other researchers (Maréchal and Kalitventzeff, 1996; Klemeš and Kravanja, 2013). In PA, results are typically represented using temperature—enthalpy profiles called Composite Curve (CC) and Grand Composite Curve (GCC). An example of such curves is given in Figure 1. Constructing the CCs involves combining the heat loads of hot streams (those that need cooling) by temperature interval to form the hot composite curve and the same for cold streams (those that need heating) to obtain the cold composite curve. When plotted together, the area between the two CCs represents the potential for heat recovery, where the heating requirement of the cold streams can be satisfied by the hot streams of the system. The distance between the composites at the top of the curves gives the minimum heating requirement and the one at the bottom gives the minimum cooling requirement. This information can also be combined to form the GCC, in which all streams are combined within temperature intervals to give a profile of process thermal energy requirements. From the GCC, one can read the temperatures at which heating and cooling are required and proper utilities can be proposed to satisfy those requirements.

FIGURE 1

Figure 1. Example of a composite curve (Left) and a grand composite curve (Right), using data from Papoulias and Grossmann (1983b).

PA can also be adapted to other flows, provided that a stream is defined with quality (temperature, concentration, material properties, etc.) and quantity (heat load, flowrate etc.) and thus has applications beyond targeting of heat recovery. The connection between mass and heat transfer led to the extension of the method to mass pinch analysis. Mass pinch analysis was first applied for reuse of solvents by El-Halwagi and Manousiouthakis (1989) and later for other material flows such as cooling water (Kim and Smith, 2001), industrial water (Wang and Smith, 1994; Olesen and Polley, 1997), and hydrogen (Towler et al., 1996; Alves and Towler, 2002; Government of Canada, 2002). Related work in the domain by Kermani et al. and Wallerand et al. also address specific problems such as optimal integration of organic Rankine cycles (Kermani et al., 2018) and heat pumps (Wallerand et al., 2018a) with industrial processes, simultaneous optimization of water and energy (Kermani et al., 2019a), and holistic industrial system design (Kermani et al., 2019b). Industrial retrofit and planning strategies using similar methods have been proposed by Bütün et al. (2018, 2019), and approaches to heat exchanger network design with utility systems (Ciric and Floudas, 1991) including new solution strategies (Mian et al., 2016) are also of interest within this domain. Concepts of solar utilities and various storage types, together with utilities like heat pumps, complicate the targeting formulation as described by Maréchal and Kalitventzeff (2003), but yield potential benefits in both the industrial and urban settings, as investigated by several researchers (Becker, 2012; Becker and Maréchal, 2012; Suciu et al., 2018; Wallerand et al., 2018b).

A similar basis of PA and mathematical programming was suggested by Kastner et al. (2015) in a review of quantitative tools for exchanges between actors in industrial parks. They provided a review of tools and existing industrial park applications of exchanges, classified by type. Methods for optimizing the exchanges were also reviewed by the authors, but without recommending specific methods or frameworks for completing such studies. The authors of the review suggest that methods need to be more comprehensive and to include aspects which have been developed by different authors and research group in the literature; thereby, many aspects can be treated by a multi-layer and integrated approach, such as that proposed in this work. Various aspects must still be added, but initial steps toward integrating disparate aspects are taken herein.

While providing information on the minimum heating and cooling requirements of a given system, PA also gives indications on the utilities which should be integrated. However, for complex systems in which multiple utilities can satisfy the demand, determining the optimal configuration of utilities is difficult based solely on classical PA. Mathematical Programming (MP) techniques have been adapted in the field to fill this gap (Papoulias and Grossmann, 1983a,b). MP transforms the problem into a purely mathematical form to obtain the optimal solution by closing the energy balance, considering both thermal and electrical energy. Since energy flows are directly linked to mass flows entering and leaving the system, mass balance constraints are also added to the formulation.

Maréchal and Kalitventzeff proposed a methodology taking advantage of both PA (targeting) and MP (synthesis) (Maréchal and Boursier, 1989). The method is extended by calculating the optimal utility configuration for a given system in three steps, named AGE (Analyse-Generate-Evaluate) (Maréchal and Kalitventzeff, 1996). In the Analyse step, PA was used to determine the maximum heat recovery and the minimum heating and cooling requirements. The Generate step used Mixed Integer Linear Programming (MILP) optimization to select and size utilities from a superstructure of potential technologies. The Evaluate step used a new composite curve definition, called Integrated Composite Curve (ICC), to determine the integration of a sub-system with the rest of the system. This work was further improved by including a steam network in the mathematical formulation to optimize co-generation of heat and electricity (Maréchal and Kalitventzeff, 1997, 1999). In Maréchal and Kalitventzeff (2003), the MILP is extended to a multi-time formulation to account for different operating modes at different time steps. Kermani et al. (2019b) proposed a holistic design approach based on these concepts, also including many additional aspects in a holistic design approach but without reiterating the basis problem formulation.

When the thermodynamic models or compositions of material streams are considered, the PI problem becomes non-linear. Both deterministic and stochastic optimization methods have been used in the literature to solve the resulting Mixed Integer Non-Linear Programming (MINLP) problem. Kantor et al. (2015) formulated the problem with a single objective, combining environmental, and economic objectives by using weight factors and simplifying several operations. Li et al. (2018) proposed a block-superstructure framework using MINLP to solve a system of balanced blocks for mass, energy, and other properties. Li et al. proposed that such an approach could be used for a variety of applications but also noted that large problems become intractable without simplifying assumptions. Gassner and Maréchal (2009) proposed a solution strategy including two steps: master optimization and slave optimization. In the slave optimization, the MILP problem was solved using the formulation of Maréchal and Kalitventzeff (2003) for each iteration of the master optimization. The master optimization used the Genetic Algorithm (GA) developed by Leyland (2002), with multiple objectives: Synthetic Natural Gas (SNG) production (economic) and energy efficiency (environmental).

The fundamental problem formulation in this work, which evolved from the formulation of Maréchal (Maréchal and Boursier, 1989; Maréchal and Kalitventzeff, 1996, 1997, 1999, 2003) uses PI in the form of an MILP. The framework has been used to analyse and optimize energy systems at large scale (Moret et al., 2014; Girones et al., 2015) and smaller scale (Girardin, 2012; Henchoz, 2016). Moreover, it has proven to be an effective tool in urban system design (Rager, 2015), biomass conversion systems (Peduzzi, 2015), and several industrial sectors such as chemicals (Pouransari, 2015; Bungener, 2016), metals (Hubert et al., 2009), food (Muller, 2007; Becker et al., 2012), and pulp and paper (Kermani et al., 2019a,b). (Liu et al., 2011) also provided insight into a generic problem formulation, using the term polygeneration to denote the multi-product nature of such arrangements. They propose a generic mathematical formulation to cover various systems such as commercial buildings, biomass and biological systems, and infrastructure planning. Liu et al. (2011) also mention future perspectives for the domain from a high level and some representative results; however, a detailed mathematical formulation or practical framework is not provided for engineers to follow such methods in practice. The current work complements the work of Liu et al. (2011) by also providing such a practical formulation, allowing engineers and practitioners to complete studies for real cases. Ng et al. (2015) Proposed optimization applications for industrial symbiosis with biorefinery/bioenergy applications, allocating various flows from decomposition of biomass feedstock to units yielding the best final arrangement for maximizing value from the biomass. The optimization approach differed from other studies in the domain by using a two-sided fuzzy optimization approach for the specific application explored in the work. Despite several reviews, such as Boix et al. (2015) noting the necessity of multiple objectives, Ng et al. (2015) omitted environmental and social factors in the work but noted the capacity of the model to treat uncertainty in feedstock quality and delivery scheduling constraints.

Recent scientific contributions in literature have also cautioned practitioners and researchers from relying on a single, global, optimum as there are inherent differences between solutions, which are not (or cannot) be captured by the mathematical formulation. Recent developments in this area extend to the use of parallel coordinate plots to display multiple solutions with their key settings and performance indicators. Effective use of such tools in the domain of process integration are notably demonstrated by Kermani et al. (2019b), who generated multiple industrial system designs using a holistic and systematic approach. This elucidates the need for advanced graphing and/or decision-making support tools to also be linked to computational frameworks, adding another element required to attain practicable methods for process integration solutions.

There have been several attempts to describe a comprehensive process integration framework. Bolliger (2010) described an integrated toolbox and explained the methodology behind it. Moret et al. (2016) summarized an MILP formulation in a compact form, covering the main sets, parameters, and variables. Computational burdens for large problems have led to increasingly more efficient implementations of the framework in various languages, such as Lua (Miller et al., 2009). Yoo et al. (2015) introduced such a Lua-based concept and gave a global overview; however, the available literature has not fully represented the problem and optimization structure in a comprehensive way. The previous efforts in the field, mentioned here, are often lacking a full problem formulation to focus on the novel aspect of the specific work; therefore, omitting the basis problem formulation. Kuznetsova et al. (2016) proposed a framework for application in eco-industrial park design and optimization, including many commonalities with the proposed work; however, no detailed mathematical framework is actually presented in the work, therefore inhibiting others from adopting or applying the approach. The authors proposed a source-sink mathematical model of interconnected black-box units and specific applications, but such a framework naturally omits many exchanges and improvements to be made within plants. Tools for optimizing process design for cleaner production were reviewed by Fan et al. (2020), delving into the basis methods such as PA and P-graph approaches in addition to novel techniques such as artificial intelligence and modern first-principles modeling tools. They introduced perspectives related to modern problems and noted findings of various active researchers in the field of process optimization, specifically with respect to the environmental impact of various scenarios and technologies. They proposed a variety of applications and tools under exploration by various authors as having large potentials in the field, including advances modeling and optimization approaches such as that proposed in the current work. An optimization-based framework for industrial symbiosis was developed by Valenzuela-Venegas et al. (2020) with application to well-known, existing, eco-industrial networks. They showed that, depending on the objective function, the implemented exchanges in existing arrangements are not necessarily optimal, highlighting the importance of methodological developments to find and implement the best solutions and not to approach the problem ad hoc. The mathematical formulation of the problem is similar to that proposed by Kantor et al. (2015), using mixed-integer non-linear programming to describe complex relationships where linear functions often introduce error. Additional effort was made to correctly model the flow of exchanged materials including detailed flow and piping calculations, though necessitating the use of non-linear constraints for correctly modeling the detailed correlations.

Development in optimization of industrial exchange networks, termed eco-industrial parks, eco-industrial networks, or simply eco-parks is also a burgeoning research field but has not been succinctly and sufficientl y addressed by recent work. Boix et al. (2015) reviewed recent work in the field, but used excessively specific search terms to identify the recent research in the field and therefore excluded many relevant publications. Regardless, this review drew attention to a variety of relevant topics in the domain and classified them by their main focus, being water, material (more generally), and energy. The authors concluded that multi-objective optimization methods are clearly required to address the complexity and potential benefits of performing such studies on potential exchange networks. The need for treating uncertainty, flexibility, operability and retrofit of energy systems was also discussed by Andiappan (2017), who reviewed the state of the art in optimization of energy systems, also relating the domain of process systems engineering to the study of industrial parks, and classifying studies based on the approach taken. Andiappan (2017) also noted various approaches for treating multi-objective optimization problems, uncertainties and design/operation decoupling, but also did not elaborate a mathematical formulation to structure and solve the model. Boix et al. (2015) also pointed to the interconnected nature of many aspects in such industrial networks and that properly formulating the mathematical problem to simultaneously optimize water, material, and energy is a key for future development. Additional gaps such as uncertainty and data collection are the subject of other publications in the field, such as Moret et al. (2017).

Such gaps as methods for simultaneous optimization of resources and energy were identified in the published literature, exhibiting a lack of thorough optimization frameworks and their formulations to consider complex industrial symbiosis problems and find novel solutions for process integration within and between industrial sites. This work, therefore, combines and builds on previous work in the field to present a comprehensive and holistic process integration framework which forms a powerful basis for solving system integration and industrial symbiosis problems in the MILP domain. This work encompasses the method and optimization framework, independent of language choice and implementation; therefore, presenting an approach and not software, as such. The extensible framework proposed herein can be facilitated by implementation in any programming language but requires coupling with an optimization tool for solving the optimization problem.

3. Methodology

The framework developed for this approach uses a mixed-integer linear programming formulation, balancing imprecision from linearization with benefits of global optima and rapid resolution time when compared to mixed-integer non-linear programming formulations. MILP formulations have been suggested and effectively applied by other researchers in the domain, as the problem can be effectively formulated using such an approach. Other problems in the domain of process integration, such as heat exchanger network design, determination of the minimum approach temperature for heat exchange, necessitate non-linear approaches to treat their specificities. Inclusion of complex units/operations are treated by connecting with external flowsheeting software or by the use of surrogate models, which is discussed briefly in section 3.6. Rapid resolution time is important when considering multiple solutions, as elucidated in the previous section to be increasingly relevant in this field, as long solution times inhibit this systematic generation of solutions and the associated feedback from decision-makers. Li et al. reported solution times between 1 and 24 h for several academic case studies (Li et al., 2018) using an MINLP solver, whereas larger problem sizes can yield solution times from 7 to 10 days (Kantor et al., 2015). Generating multiple solutions for large problems can therefore become computationally prohibitive using an MINLP framework, highlighting the importance of rapid solution generation using MILP. The problem formulation in this chapter is principally inspired by Maréchal and Kalitventzeff (2003), with learning from other research in the domain, summarized in section 2, including the extension of the targeting approach to include material flows, which are especially relevant for industrial symbiosis applications. Kermani et al. (2019b) proposed a holistic design approach based on these principles, exploring additional options for utility integration and a solution strategy with a genetic algorithm guiding the MILP formulation to explore a large and complex solution space with many objectives. Such a MILP sub-problem is explicitly formulated here to provide the basis for such additional work.

3.1. Definition of Sets

The MILP problem presented in this article is essentially a simultaneous sizing and scheduling optimization problem, where the objective is to minimize costs and/or environmental impact, and where the system to be optimized is represented by a number of units belonging to the set U. Units can be seen as nodes that supply, demand or convert material and energy streams (∈ S) associated with different layers (∈ L). In each layer, a constraint is imposed to close the mass or energy balance.

The set of units is further divided into process units (subset PU ⊂ U) and utility units (subset UU ⊂ U). A process unit has fixed size and operation, which means that it cannot be optimized (i.e., the associated decision variables are fixed). They typically represent a process demand (e.g., chemical plant) or service (e.g., space heating in a building) that has to be fulfilled. In contrast, utility units are not fixed and can therefore be optimized in terms of existence, size, and operation.

The different sets are represented in Figure 2 using an illustrative example of an energy system. In this example, five layers, corresponding to different resources and forms of energy, are defined (represented by different colors on the figure): Heat, Electricity, Gas, Wood, and SNG¹. Each unit (represented by boxes) has at least one associated stream (represented by arrows) and each must belong to a layer (e.g., the stream EG_elec is associated to ElecGrid and is of layer Electricity).

FIGURE 2

Figure 2. Illustrative example showing the main sets of the MILP problem.

Two process units, ProcessA and ProcessB, impose a demand of heat and electricity to be fulfilled by a set of utility units representing energy conversion technologies: Gasif , CHP², SCHP³, HP⁴. The latter require inputs of resources that are provided by additional utility units (WoodRes, GasGrid, ElecGrid, and WasteHeat). Finally, excess heat that has not been used can be disposed of in the AeroCond utility unit to close the system energy balance.

As mentioned, each stream s ∈ S is associated to a unit u ∈ U and layer l ∈ L, such that a subset SoU_u ⊂ S containing all streams associated to unit u and a subset SoL_l ⊂ S containing all streams associated to layer l can be defined. From these subsets, it is possible to deduce the subsets UoS_s ⊂ U (Equation 1) and UoL_l ⊂ U (Equation 2) representing the unit to which a stream s is associated, and the units that have streams of layer l associated to them, respectively.

$\begin{array}{ll}{\text{UoS}}_s=\left\{u\in \text{U}:s\in {\text{SoU}}_u\right\}\forall s\in \text{S}& (1)\end{array}$

$\begin{array}{ll}{\text{UoL}}_l={\displaystyle \bigcup _{s\in {\text{SoL}}_l}}{\text{UoS}}_s\forall l\in \text{L}& (2)\end{array}$

The system is divided into sub-systems called clusters (set C) which imply the closure of certain types of balances. Hence, each unit is located inside a cluster, and the subset UoC_c ⊂ U contains all units located inside cluster c ∈ C. This concept is illustrated graphically in Figure 3.

FIGURE 3

Figure 3. Example of cluster set.

Different layer types, on which different types of constraints may apply, are defined in the set LT. The different layer types defined here are mass balance (mb), resource balance (rb), and heat cascade (hc). Each layer has an associated layer type, and the subset LoT_lt ⊂ L contains all layers of type lt ∈ LT. The subsets MBU_{l, c} ⊂ U and RBU_l ⊂ U, containing all units with at least one stream of type mb and rb associated to them, respectively, are obtained using (3) and (5).

$\begin{array}{ll}{\text{MBU}}_{l,c}={\text{UoL}}_l\cap {\text{UoC}}_c\forall l\in {\text{LoT}}_{mb},\forall c\in \text{C}& (3)\end{array}$

$\begin{array}{ll}{\text{RBU}}_l={\text{UoL}}_l\forall l\in {\text{LoT}}_{rb}& (4)\end{array}$

3.2. Objective

Considering multiple objectives such as economics, environmental impact and social obligations, the process integration MILP problem could have many objectives, two of which are often cost and environmental impact. The goals can be achieved by minimizing one of the objectives shown here as Equation (5) for operating cost, Equation (6) for investment cost, or Equation (7) for environmental impact. The objectives mentioned here are typical selections for many problems in process integration but other objectives could be considered with an appropriate mathematical formulation. Objective selection must be completed as required for the specific user and application.

$\begin{array}{l}{\text{C}}^{op}={\displaystyle \sum _{u\in \text{U}}}{\displaystyle \sum _{p\in \text{P}}}{\displaystyle \sum _{t\in \text{T}}}\left({\text{c}}_{u,p,t}^{op,fix}·y_{u,p,t}^{\prime }+{\text{c}}_{u,p,t}^{op,var}·f_{u,p,t}^{\prime }\right)·\Delta {\text{t}}_t^{op}\\ \forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (5)\end{array}$

where ${\mathrm{\text{c}}}_{u,p,t}^{op,fix}$ and ${\mathrm{\text{c}}}_{u,p,t}^{op,var}$ are the fixed and variable operating costs of unit u ∈ U in period p ∈ P and time t ∈ ToP_p, respectively, and $\Delta {\mathrm{\text{t}}}_t^{op}$ is the operating time of t.

$\begin{array}{l}{\text{C}}^{inv}={\displaystyle \sum _{u\in \text{U}}}{\displaystyle \sum _{p\in \text{P}}}\left({\text{c}}_u^{inv,fix}·y_{u,p}+{\text{c}}_u^{inv,var}·f_{u,p}\right)·\Delta {\text{t}}_t^{op}\\ \forall u\in \text{U},\forall p\in \text{P}& (6)\end{array}$

where ${\mathrm{\text{c}}}_u^{inv,fix}$ and ${\mathrm{\text{c}}}_u^{inv,var}$ are the fixed and variable investment cost of unit u ∈ U in period p ∈ P, respectively.

$\begin{array}{l}\text{W}={\displaystyle \sum _{u\in \text{U}}}{\displaystyle \sum _{p\in \text{P}}}{\displaystyle \sum _{t\in \text{T}}}\left({\text{e}}_u^{fix}·y_{u,p}+{\text{e}}_u^{var}·f_{u,p}\right)·\Delta {\text{t}}_t^{op}\\ \forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (7)\end{array}$

where ${\mathrm{\text{e}}}_u^{fix}$ and ${\mathrm{\text{e}}}_u^{var}$ are the fixed and variable impact of unit u ∈ U in period p ∈ P and time t ∈ ToP_p, respectively. The impact can be defined as any quantifiable environmental impact, such as CO₂ emissions.

If environmental costs can be monetized, Equation (7) can be included in the operating and investment cost functions by adding the cost of the emission, c^e, as shown in Equations (8) and (9) for operating cost and investment cost with emissions, respectively.

$\begin{array}{l}{{\text{C}}^{op}}^{+e}={\displaystyle \sum _{u\in \text{U}}}{\displaystyle \sum _{p\in \text{P}}}{\displaystyle \sum _{t\in \text{T}}}(\left({\text{c}}_{u,p,t}^{op,fix}+{\text{c}}^e·{\text{e}}_{u,p,t}^{fix}\right)·y_{u,p,t}^{\prime }\\ +\left({\text{c}}_{u,p,t}^{op,var}+{\text{c}}^e·{\text{e}}_{u,p,t}^{var}\right)·f_{u,p,t}^{\prime })·\Delta {\text{t}}_t^{op}\\ \forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p,\forall e\in \text{E}& (8)\end{array}$

$\begin{array}{l}{{\text{C}}^{inv}}^{+e}={\displaystyle \sum _{u\in \text{U}}}{\displaystyle \sum _{p\in \text{P}}}{\displaystyle \sum _{t\in \text{T}}}(\left({\text{c}}_u^{inv,fix}+{\text{c}}^e·{\text{e}}_u^{fix}\right)·y_{u,p,t}^{\prime }\\ +\left({\text{c}}_u^{inv,var}+{\text{c}}^e·{\text{e}}_u^{var}\right)·f_{u,p,t}^{\prime })·\Delta {\text{t}}_t^{op}\\ \forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p,\forall e\in \text{E}& (9)\end{array}$

The total cost can also be calculated which finds the minimum cost solution by combining the benefits from reducing the operating cost and associated capital expenditure. This total cost is calculated as shown by Equation (10) or, including the monetized environmental impact, Equation (11).

$\begin{array}{ll}{C}^{total}={\text{C}}^{op}+\text{a}·{\text{C}}^{inv}& (10)\end{array}$

$\begin{array}{ll}{C}^{total+e}=C^{op+e}+\text{a}·C^{inv+e}& (11)\end{array}$

Where a is the annualization factor calculated by Equation (12) which is based on the interest rate d imposed by the industry or project investor and the project lifetime z, shown in Equation (12).

$\begin{array}{ll}\text{a}=\left(\frac{\text{d}{(1+\text{d})}^{\text{z}}}{{(1+\text{d})}^{\text{z}}-1}\right)& (12)\end{array}$

3.3. Sizing and Scheduling

The sizing of each unit is governed by Equation (13), where f_{u, p} is a continuous variable representing the capacity factor, y_{u, p} is a binary variable representing the existence of the unit, and ${\mathrm{\text{f}}}_u^{\mathrm{\text{min}}}$ and ${\mathrm{\text{f}}}_u^{\mathrm{\text{max}}}$ are parameters representing the minimum and maximum capacity of unit u. Sizing is defined for each period p, representing high level time scales (e.g., years) and used, for instance, to account for investment stages. Equation (13) ensures that, if a unit exists in period p, its capacity is within given boundaries which are set for each unit.

$\begin{array}{ll}{\text{f}}_u^{\text{min}}·y_{u,p}\le f_{u,p}\le {\text{f}}_u^{\text{max}}·y_{u,p}\forall u\in \text{U},\forall p\in \text{P}& (13)\end{array}$

The scheduling of each unit is governed by Equations (14–16), where $f_{u,p,t}^{\prime }$ is a continuous variable representing the usage factor, $y_{u,p,t}^{\prime }$ is a binary variable representing the activation, and ${\mathrm{\text{l}}}_{u,p}^{\mathrm{\text{min}}}$ and ${\mathrm{\text{l}}}_{u,p}^{\mathrm{\text{max}}}$ are parameters representing the minimum and maximum load factor of unit u. Scheduling is defined for each time t, representing lower-level time scales when compared to periods (e.g., hours) and used for operational scheduling⁵. Equation (14) ensures that, if a unit is activated in time t, its usage factor is within given boundaries defined as a fraction of the capacity f_{u, p} of the unit. Moreover, if a unit is not activated in a given time t, its usage factor must also equal 0. This is imposed by the combination of: (a) multiplying the lower bound by the binary variable $y_{u,p,t}^{\prime }$ and, (b) imposing the constraint in Equation (15). Finally, Equation (16) ensures that a unit can only be activated if it exists in period p.

$\begin{array}{ll}{\text{l}}_{u,p}^{\text{min}}·y_{u,p,t}^{\prime }·f_{u,p}\le f_{u,p,t}^{\prime }\le {\text{l}}_{u,p}^{\text{max}}·f_{u,p}\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (14)\end{array}$

$\begin{array}{ll}{f}_{u,p,t}^{\prime }\le {\text{l}}_{u,p}^{\text{max}}·{\text{f}}_u^{\text{max}}·y_{u,p,t}^{\prime }\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (15)\end{array}$

$\begin{array}{ll}{y}_{u,p,t}^{\prime }\le y_{u,p}\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (16)\end{array}$

The multiplication of $y_{u,p,t}^{\prime }$ and f_{u, p} in Equation (14) introduces a non-linearity which must be linearized to maintain the formulation within the MILP domain; therefore, the lower bound in this equation is linearized by a standard technique which requires creating a bounded auxiliary variable. For this purpose, a new decision variable $v_{u,p,t}^1$ is created and Equation (14) is replaced by Equations (17–20), which are an equivalent set of linear constraints.

$\begin{array}{ll}{v}_{u,p,t}^1\le f_{u,p,t}^{\prime }\le {\text{l}}_{u,p}^{\text{max}}·f_{u,p}\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (17)\end{array}$

$\begin{array}{ll}{v}_{u,p,t}^1\ge {\text{l}}_{u,p}^{\text{min}}·f_{u,p}-(1-y_{u,p,t}^{\prime })·{\text{l}}_{u,p}^{\text{min}}·{\text{f}}_u^{\text{max}}\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (18)\end{array}$

$\begin{array}{ll}{v}_{u,p,t}^1\ge 0\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (19)\end{array}$

$\begin{array}{ll}{v}_{u,p,t}^1\le {\text{l}}_{u,p}^{\text{min}}·f_{u,p}\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (20)\end{array}$

Some circumstances require the formulation to be flexible in allowing forced activation of specific units. Parameter ${\mathrm{\text{y}}}_{u,p,t}^{force}$ and Equation (21) have been added to force the activation of a unit in a given period and time.

$\begin{array}{ll}{y}_{u,p,t}^{\prime }\ge {\text{y}}_{u,p,t}^{force}\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (21)\end{array}$

Given that they have a fixed size and operation, process units (u ∈ PU), have some parameters which are fixed. Hence, the following sizing parameters used in Equations (13–21) must all be set to 1: ${\mathrm{\text{f}}}_u^{\mathrm{\text{min}}}$, ${\mathrm{\text{f}}}_u^{\mathrm{\text{max}}}$, ${\mathrm{\text{l}}}_{u,p}^{\mathrm{\text{min}}}$, ${\mathrm{\text{l}}}_{u,p}^{\mathrm{\text{max}}}$, ${\mathrm{\text{y}}}_{u,p,t}^{force}$.

3.4. Mass and Energy Balance

For each unit u, supply ${\dot{M}}_{l,u,p,t}^{+}$ and demand ${\dot{M}}_{l,u,p,t}^{-}$ of a specific layer l ∈ L are calculated with Equations (22) and (23), respectively, where ${\stackrel{\cdot }{\mathrm{\text{m}}}}_{l,u,p,t}^{out}$/ ${\stackrel{\cdot }{\mathrm{\text{m}}}}_{l,u,p,t}^{in}$ represent the reference output/input flow⁶. Note that these are only defined for layers of type mb and rb. It should also be mentioned that the reference flows are originally properties of the streams. Indeed, each mb and rb stream is defined with a layer, a direction (in/out) and a positive reference flow. For these layer types, there can only be one stream of a given layer defined per unit and the reference flow value of the unit for that layer corresponds to the value of the associated stream (and based on its direction).

$\begin{array}{l}{\dot{M}}_{l,u,p,t}^{+}={\stackrel{\cdot }{\text{m}}}_{l,u,p,t}^{out}·f_{u,p,t}^{\prime }\\ \forall l\in {\text{LoT}}_{lt}:lt\in \left\{mb,rb\right\},\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (22)\end{array}$

$\begin{array}{l}{\dot{M}}_{l,u,p,t}^{-}={\stackrel{\cdot }{\text{m}}}_{l,u,p,t}^{in}·f_{u,p,t}^{\prime }\\ \forall l\in {\text{LoT}}_{lt}:lt\in \left\{mb,rb\right\},\forall u\in \text{U},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (23)\end{array}$

For each layer of type mb, the balance is closed inside each cluster c, meaning that streams of this layer type cannot exchange across clusters. This is enforced by Equation (24), which ensures the equalization of supply and demand for layer l, and thus the supply of all units equals the demand of all units located inside the cluster. Conversely, for layers of type rb, the balance is closed for the whole system (Equation 25), in which case the streams can exchange across clusters. The distinction between whether a particular flow is defined in a mb or rb layer is not strictly defined and depends on the specific situation of whether a flow could be shared between clusters or not. The best usage of rb layers is for resources which can be easily transported or defined as a network such as electricity, natural gas, steam, industrial gases such as hydrogen or chloride, among others. For considering industrial symbiosis, the most flexibility is provided by defining all streams which could be exchanged as belonging to rb layers which allows the greatest number of possible connections between disparate industries, sectors or sites. Streams which are difficult or impossible to exchange should be defined in mb layers to prevent the exchange between clusters.

$\begin{array}{l}{\displaystyle \sum _{u\in {\text{MBU}}_{l,c}}}{\dot{M}}_{l,u,p,t}^{+}={\displaystyle \sum _{u\in {\text{MBU}}_{l,c}}}{\dot{M}}_{l,u,p,t}^{-}\\ \forall l\in {\text{LoT}}_{mb},\forall c\in \text{C},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (24)\end{array}$

$\begin{array}{l}{\displaystyle \sum _{u\in {\text{RBU}}_l}}{\dot{M}}_{l,u,p,t}^{+}={\displaystyle \sum _{u\in {\text{RBU}}_l}}{\dot{M}}_{l,u,p,t}^{-}\\ \forall l\in {\text{LoT}}_{rb},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (25)\end{array}$

Two additional variables are defined to represent the quantity of a given layer transferred from one unit to another. ${x^{mb}}_{l,c,i,j,p,t}$ is the quantity of layer l ∈ LoT_mb that is transferred from unit i ∈ MBU_{l, c} to unit j ∈ MBU_{l, c} (both inside cluster c ∈ C), in period p and time t. x^rb follows the same principles but applies to rb layers and their corresponding units and therefore considers transfer between units belonging to different clusters. x^mb and x^rb are determined by Equations (26) and (27), respectively. To ensure that units only transfer what they supply (i.e., to avoid transfer nodes), Equations (28) and (29) are added.

$\begin{array}{l}{\dot{M}}_{l,u,p,t}^{+}+{\displaystyle \sum _{\begin{array}{c}i\in {\text{MBU}}_{l,c}\\ i\ne u\end{array}}}{x}_{l,c,i,u,p,t}^{mb}={\dot{M}}_{l,u,p,t}^{-}+{\displaystyle \sum _{\begin{array}{c}j\in {\text{MBU}}_{l,c}\\ j\ne u\end{array}}}{x}_{l,c,u,j,p,t}^{mb}\\ \forall l\in {\text{LoT}}_{mb},\forall c\in \text{C},\forall u\in {\text{MBU}}_{l,c},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (26)\end{array}$

$\begin{array}{l}{\dot{M}}_{l,u,p,t}^{+}+{\displaystyle \sum _{\begin{array}{c}i\in {\text{RBU}}_l\\ i\ne u\end{array}}}{x}_{l,i,u,p,t}^{rb}={\dot{M}}_{l,u,p,t}^{-}+{\displaystyle \sum _{\begin{array}{c}j\in {\text{RBU}}_l\\ j\ne u\end{array}}}{x}_{l,u,j,p,t}^{rb}\\ \forall l\in {\text{LoT}}_{rb},\forall u\in {\text{RBU}}_l,\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (27)\end{array}$

$\begin{array}{l}{\displaystyle \sum _{\begin{array}{c}i\in {\text{MBU}}_{l,c}\\ i\ne u\end{array}}}{x}_{l,c,u,i,p,t}^{mb}\le {\dot{M}}_{l,u,p,t}^{+}\\ \forall l\in {\text{LoT}}_{mb},\forall c\in \text{C},\forall u\in {\text{MBU}}_{l,c},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (28)\end{array}$

$\begin{array}{l}{\displaystyle \sum _{\begin{array}{c}i\in {\text{RBU}}_l\\ i\ne u\end{array}}}{x}_{l,u,i,p,t}^{rb}\le {\dot{M}}_{l,u,p,t}^{+}\\ \forall l\in {\text{LoT}}_{rb},\forall u\in {\text{RBU}}_l,\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (29)\end{array}$

An illustrative example of a layer balance is given in Figure 4. It shows the mass balance for the Gas layer of type rb, in a given period and time. In this example, the unit GasGrid supplies 21 kg/s of Gas to the CHPa and CHPb units that have respective demands of 5 and 16 kg/s, hence closing the balance. ${\stackrel{\cdot }{\mathrm{\text{m}}}}_{Gas,GasGrid}^{out}$, ${\stackrel{\cdot }{\mathrm{\text{m}}}}_{Gas,CHPa}^{out}$ and ${\stackrel{\cdot }{\mathrm{\text{m}}}}_{Gas,CHPb}^{out}$ are fixed parameters of the MILP problem, while the variables are decisions of the optimization. It should be noted that as this is an rb layer; therefore, the balance is closed over all clusters and not within clusters.

FIGURE 4

Figure 4. Illustrative example of a rb layer balance.

3.5. Heat Cascade

The heat cascade constraints are applied to layers of type hc. They enforce the first and the second law of thermodynamics, ensuring that the energy balance in each cluster is closed (first law) and heat is only transferred from higher temperature intervals to lower ones (second law). Implicitly, this formulation includes concepts of heat integration by allowing the cascade of heat from higher temperatures to lower temperatures within a cluster.

The streams of a heat cascade layer are divided into two subsets, HS_{l, c} and CS_{l, c}, for hot and cold streams, respectively and are defined according to Equations (31) and (32). Each of these streams is characterized by an inlet (${\mathrm{\text{T}}}_s^{in}$) and outlet temperature⁷ (${\mathrm{\text{T}}}_s^{out}$), and an inlet (${\stackrel{\cdot }{\mathrm{\text{h}}}}_s^{in}$) and outlet enthalpy (${\stackrel{\cdot }{\mathrm{\text{h}}}}_s^{out}$). A hot stream is cooled down and therefore delivers energy to the system, while a cold stream is heated up and takes energy from the system.

$\begin{array}{ll}{\text{HS}}_{l,c}=\left\{s\in {\text{SoL}}_l\cap {\text{SoC}}_c:{\text{T}}_s^{out}<{\text{T}}_s^{in}\right\}\forall l\in {\text{LoT}}_{hc},\forall c\in \text{C}& (30)\end{array}$

$\begin{array}{ll}{\text{CS}}_{l,c}=\left\{s\in {\text{SoL}}_l\cap {\text{SoC}}_c:{\text{T}}_s^{out}>{\text{T}}_s^{in}\right\}\forall l\in {\text{LoT}}_{hc},\forall c\in \text{C}& (31)\end{array}$

where

$\begin{array}{l}{\text{SoC}}_c={\displaystyle \bigcup _{u\in {\text{UoC}}_c}}{\text{SoU}}_u\forall c\in \text{C}\end{array}$

For each heat cascade layer and each cluster, the list of all stream inlet and outlet temperatures is extracted and ordered to generate the set of temperature intervals⁸ K. Hence, T_k represents the lower temperature of interval k ∈ K. Equation (32) ensures that the energy balance is satisfied in each temperature interval⁹. This is achieved by residual heat Ṙ_{l, c, k, p, t} which transfers excess heat from higher temperature intervals (k) to lower temperature intervals (k − 1). As the summation is applied to the streams, a usage factor for the stream is used, which is equal to that of the unit to which the stream belongs.

$\begin{array}{l}{\displaystyle \sum _{\begin{array}{c}s\in {\text{HS}}_{l,c}\\ {\text{T}}_s^{out}\ge {\text{T}}_{l,c,k}+ϵ\end{array}}}{\stackrel{\cdot }{\text{M}}{C}_P}_s·{\text{f}}_{s,p,t}^{\prime \prime }·({\text{T}}_s^{in}-{\text{T}}_s^{out})\\ +{\displaystyle \sum _{\begin{array}{c}s\in {\text{HS}}_{l,c}\\ {\text{T}}_s^{out}\le {\text{T}}_{l,c,k}\le {\text{T}}_s^{in}\end{array}}}{\stackrel{\cdot }{\text{M}}{C}_P}_s·{\text{f}}_{s,p,t}^{\prime \prime }·({\text{T}}_s^{in}-{\text{T}}_{l,c,k})\\ -{\displaystyle \sum _{\begin{array}{c}s\in {\text{CS}}_{l,c}\\ {\text{T}}_s^{in}\ge {\text{T}}_{l,c,k}+ϵ\end{array}}}{\stackrel{\cdot }{\text{M}}{C}_P}_s·{\text{f}}_{s,p,t}^{\prime \prime }·({\text{T}}_s^{out}-{\text{T}}_s^{in})\\ -{\displaystyle \sum _{\begin{array}{c}s\in {\text{CS}}_{l,c}\\ {\text{T}}_s^{in}\le {\text{T}}_{l,c,k}\le {\text{T}}_s^{out}\end{array}}}{\stackrel{\cdot }{\text{M}}{C}_P}_s·{\text{f}}_{s,p,t}^{\prime \prime }·({\text{T}}_s^{out}-{\text{T}}_{l,c,k})\\ -{\dot{R}}_{l,c,k,p,t}=0\\ \forall l\in {\text{LoT}}_{hc},\forall c\in \text{C},\forall k\in \text{K},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (32)\end{array}$

where

$\begin{array}{l}{\stackrel{\cdot }{\text{M}}{C}_P}_s=\frac{{\stackrel{\cdot }{\text{h}}}_s^{in}-{\stackrel{\cdot }{\text{h}}}_s^{out}}{{\text{T}}_s^{in}-{\text{T}}_s^{out}}\forall s\in {\text{HS}}_{l,c}\cup {\text{CS}}_{l,c}\end{array}$

and

$\begin{array}{l}{\text{f}}_{s,p,t}^{\prime \prime }=f_{u,p,t}^{\prime }\forall s\in \text{S},\forall u\in {\text{UoS}}_s,\forall p\in \text{P},\forall t\in {\text{ToP}}_p\end{array}$

The residual heat can be a minimum of zero when no heat is transferred from the corresponding temperature interval to the lower ones. This corresponds to the pinch temperature as defined for the eponymous method. Equation (33) sets the lower bound of Ṙ_{l, c, k, p, t}.

$\begin{array}{ll}{\dot{R}}_{l,c,k,p,t}⩾0\forall l\in {\text{LoT}}_{hc},\forall c\in \text{C},\forall k\in \text{K},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (33)\end{array}$

Residual heat in the first interval Ṙ_{l, c, 1, p, t} is also zero since lower temperature intervals do not exist to accept the residual heat. Similarly, since heat cannot be delivered to the n_k^th interval from a higher temperature interval, Ṙ_{l, c,nk + 1, p, t} is zero. Equation (34) imposes these constraints.

$\begin{array}{l}{\dot{R}}_{l,c,1,p,t}=0,{\dot{R}}_{l,c,n_k+1,p,t}=0\\ \forall l\in {\text{LoT}}_{hc},\forall c\in \text{C},\forall p\in \text{P},\forall t\in {\text{ToP}}_p& (34)\end{array}$

To illustrate this point, Figure 5 represents an example of a heat cascade inside a cluster PlantA containing a process unit and four utility units, for a given period and time. In this example, each unit has at least one hot or cold stream¹⁰ associated with it, all of which are characterized by an inlet and outlet temperature and enthalpy parameter. Thus, it is possible to draw the hot (red) and cold (blue) composite curves that show the enthalpy-temperature characteristics of the aggregated hot and cold streams, respectively, as shown in Figure 5.

FIGURE 5

Figure 5. Illustrative example of a heat cascade “balance”.

Moreover, the temperature intervals, which include all temperatures of all streams, are represented on this figure. Focusing on the upper portion of the curve shows the heat cascade for a specific interval k. It can be observed that for this interval, the constraint from Equation (32) imposes that the sum of the enthalpy difference of all hot streams situated above temperature T_k is equal to the sum of the enthalpy difference of all cold streams situated above temperature T_k plus the heat Ṙ_k that is cascaded to the next interval k − 1. One can therefore understand that the pinch temperature, also shown on the graph, is where the lower bound of the constraint shown in Equation (33) is reached, and where no heat is cascaded to a lower interval.

3.6. Integration With External Software

Implementation of the formulation presented in this paper could be completed in a variety of ways. Thus, the purpose of this effort is not to provide software, but rather a method which can be implemented flexibly in any programming language. For optimization, either an MILP solver must be constructed for this purpose or an input file prepared and executed within an existing optimization platform, which is the recommended approach for researchers or practitioners who are not experienced in developing optimization algorithms. A minimalist scripting language is recommended for structuring the problem and interacting between different components such as databases, flowsheeting software, plotting, thermodynamics, reporting, and the mentioned optimization languages. The optimization problem can then be written using any language suited to the purpose and linked to the appropriate solvers. The authors have used Lua (Miller et al., 2009) as a scripting language, with the optimization problem solved using AMPL (Fourer et al., 1990) or GLPK (Free Software Foundation Inc., 2019) with existing MILP solvers provided therein.

The MILP domain allows for a certain degree of modeling flexibility but inherently is limited by the requirement of linearity. Non-linear and highly complex processing systems are often simulated using flowsheeting software or purpose-built industry-specific tools. When non-linearities are integral within a process and functions cannot be linearized by well-defined techniques, these processes can be included by connecting flowsheeting or non-linear equation solvers as black-box calculation tools. Input data are prepared and output data are received from the external software without interacting directly with the reaction scheme, systems of differential equations or other non-linearities which traditionally have required simplification in MILP models or construction of surrogate models. Most software packages allow interaction from a terminal execution command, database trigger or other method which can then be linked with the optimization. Flowsheeting software which have been tested for such communications are Belsim Vali (Belsim, 2018), Aspen Plus (Al-Malah, 2016), and gPROMS [Process Systems Enterprise, PSE]; thus, any flowsheet model constructed in these software can be solved externally to provide the reference case for the MILP. This strategy obviates the need for surrogate model construction and allows a large degree of flexibility, though poses convergence risks for complex systems. The choice of whether to employ the external simulation or use of surrogate models should be evaluated for each case depending on the problem. Thermodynamic property databases such as Refprop (NIST, 2013) or Coolprop (Bell et al., 2014) can also be connected in this way to calculate and return the appropriate properties of pure components or mixtures for a range of conditions. In each example mentioned here, external software outputs must simply be connected to appropriate variables defined in the optimization framework for seamless integration.

4. Case Study

To illustrate how the model can be applied to industrial symbiosis studies, it has been tested on a fictional case study represented in Figure 6. As shown on this figure, the case study is composed of five process units, each representing an industrial process from a different sector. The processes are simple black-box models with basic inputs of raw materials and energy, and outputs of products, co-products and waste. It should be Highlighted that the models used for the case study are purely fictional, and included for the sole purpose of illustrating the potential of the model. Despite the simplistic models used for the illustrative case, the approach remains the same as for real applications with a different level of detail/complexity.

FIGURE 6

Figure 6. Superstructure representation used for the case study.

Each process unit is in a cluster of its own, which also contains a list of utility units that are common to all industrial processes. Many are relevant for locally producing or transforming heat (heat cascade layer), but also electricity production, steam production, and waste water treatment units are also included. Cluster A also contains an additional utility unit that purifies a waste stream into a usable material stream. The raw materials used as inputs by the process units are each supplied by a utility unit, all of which are regrouped in the Resources cluster. Electricity, gas, and water are provided by utility network units included in the Networks cluster. Products of the different processes are sold to utility units in the Market cluster. Finally, unsaleable co-products and waste are sent to utility units in the Waste cluster. Except the heat cascade layer to which the hot and cold streams belong and the Biomass layer which is a mass balance layer, the layers considered in this example are resource balance layers, meaning that streams of this type can exchange between clusters.

For illustrating discovery/optimization of industrial symbiosis, some products, co-products or waste from process units can be used by others, either directly (e.g., Process B outputs Material 7, which can be used by Process C) or indirectly with an intermediate utility unit (e.g., Waste 1 produced by Process C can be upgraded to Material 3 using the Purification unit, which can then be used by Process A). Moreover, excess heat (i.e., hot streams) which is not directly valorized in the process itself can be sent to another cluster via the Steam supply unit which converts a cold stream into a resource balance layer stream (Steam), which can be sent to another cluster where it can be converted into a hot stream (with Steam intake unit) to be used by the process unit. The steam network is considered to be a closed loop which implies that all condensate is recovered and recycled.

The list of utility units and associated parameter values are given in the Appendix (Table 4). For all utility units, the minimum load factor is set to zero (${\mathrm{\text{l}}}_{u,p}^{\mathrm{\text{min}}}=0$) and maximum load factor to one (${\mathrm{\text{l}}}_{u,p}^{\mathrm{\text{max}}}=1$); therefore, they are not included in the table. For process units, all cost and emissions factors are set to zero in this example and their other sizing parameters are set to unity (see section 3). The list of streams and their characteristics are also given in the Appendix (Tables 5, 6 for hc and rs/ms streams, respectively) and are ordered by units to which they belong.

Two single objective optimizations were completed, using the case study presented herein. The objective function for the first case was to minimize the total annualized cost (including operating, maintenance and investment costs). For the second, the objective function was to minimize the total annual direct CO₂ emissions. In both cases, the time horizon considered was 1 year ($\Delta {\mathrm{\text{t}}}_t^{op}$ = 8,760 h), and a constant operating condition was assumed over the full period. Additionally, a total cost optimization was carried out on each of the five process unit clusters independently, yielding a reference without industrial symbiosis as a baseline for comparison.

Tables 1–3 give the resulting utility unit sizes in each cluster for each optimization. The utility units can be divided into four categories; the first category includes the utilities that are required to meet the process heat requirements (boilers and gas engine) or to evacuate excess heat (cooling tower). The second category includes the units with streams that are integrated into the heat cascade of the process. These recover heat from the process for either heat pumping or electricity generation (Rankine cycle, ORC). The third category relates to transfer of heat between clusters (steam supply and intake). The fourth category deals with material streams, which includes treatment, upgrading and disposal of waste streams and co-products.

TABLE 1

Table 1. Utility unit sizing in each cluster (reference optimization).

TABLE 2

Table 2. Utility unit sizing in each cluster (total cost optimization).

TABLE 3

Table 3. Utility unit sizing in each cluster (emissions optimization).

The first category of utilities are always present to provide the heat required for process operation. For economic objectives in the optimization, the gas engine is typically selected as the priority provider of heat. Excepting the maximum capacity parameter, the sizing is limited by the maximum electricity or heating demand (whichever is reached first). The reference case, where each cluster is optimized independently and heat cannot be transferred from one cluster to another, is typical of conventional plant operation. In that case, extra heat is provided by the gas boiler. Conversely, when all clusters are optimized simultaneously, heat can be exchanged via the steam network, which results in higher utilization of the gas engine power and less of the gas boiler, since heating and electricity supply/demand are combined from all clusters. Optimization on emissions reduces the gas consumption to zero, and all gas boilers and gas engines are replaced with biomass boilers at their maximum size to maximize electricity production. Furthermore, biomass boilers, which had a prohibitive investment cost in the total cost optimization, become the most attractive option as costs are no longer considered in the objective function.

The second category of utilities (heat pumps, Rankine cycles, and ORCs) may be selected if integration with process streams is possible, and if there is an advantage to do so for the given objective function. Such opportunities can be visualized on the composite curves, which are represented for each cluster of the combined optimization in Figure 7. On Figure 7C, one can see that a heat pump is integrated with process C, while on Figure 7D, an ORC is integrated with process E. Conventional Rankine cycles are never present in cost-optimal solutions because the process hot stream temperatures are too low compared to the evaporation temperature of the Rankine cycle. Hence, it would require heat from a hot utility (first category), and additional investment for this extra capacity cannot be justified. Conversely, the environmental impact optimization maximizes the use of Rankine cycles to produce as much local electricity as possible, which results in reducing the CO₂ emissions from imported electricity. The heat source for this case is biomass-derived and is thus neutral with respect to CO₂ emissions. As such, the load factors of biomass boilers are maximized, but are insufficient provide the entire electricity demand.

FIGURE 7

Figure 7. Composite curves of the case study processes. (A) Process A, (B) Process B, (C) Process C and (D) Process E.

The integration options (second category of utilities) are in competition with the third category: the units that allow transfer of heat from one cluster to another (steam network in this case). This leads to different integration options in the independent cluster optimization and the combined cluster optimization. In the latter, greater cost savings are achieved for the overall system if excess heat from some clusters is sent to other clusters via the steam network, rather than valorized locally. The composite curves show how the steam network (supply and intake units) is integrated with the process, and how it competes with the other options. On Figures 7B,D, it can be observed that the steam network recovers heat from the process, and is distributed to process A (see Figure 7A), reducing the required size of the gas engine. It should be noted that the utilities in the second and third categories do not have an associated investment cost and that the results could be very different had they been included.

The fourth and last category of utilities deal with the treatment of material streams. Only two units of this type are considered: waste water treatment plant (WWTP) and purification. The WWTP converts a wastewater stream into a (clean) water stream using heat and electricity while also producing some biowaste. The purification unit upgrades a waste stream (Waste1) to a input material stream (Material3) and also consumes heat end electricity. Optimizing either cost or emissions, the WWTP is never selected from the superstructure. With the parameters imposed for the WWTP, it is always more beneficial to export wastewater and have it treated externally due to capital cost of the treatment unit¹¹. The purification unit is selected in cost optimization, but not for emissions optimization. In the reference case, the sizing of purification is limited by the quantity of Waste1 produced by process A, and therefore some import of Material3 is still required. When clusters are optimized together, a symbiosis opportunity appears, where Waste1 produced by process A is sent to the purification unit in cluster A. In that case, import of Material3 is no longer required. This can be observed in the Sankey diagram of Figure 8.

FIGURE 8

Figure 8. Sankey diagram showing the resulting flows for total cost optimization.

The flow of material streams and industrial symbiosis opportunities are represented using here using Sankey diagrams to show the flow of resources between clusters¹², and are shown in Figures 8, 9. These show, for instance, that Material5 from process E and Material7 from process B are sent to process C. It is also observed that biomass boilers in the emissions optimization are primarily supplied by biowaste resulting from process A. In reality, all three sources of biomass are equivalent in the emissions optimization, as none have associated CO₂ emissions (directly or indirectly).

FIGURE 9

Figure 9. Sankey diagram showing the resulting flows for environmental impact optimization.

The optimization was carried out using an open source MILP solver named GLPSOL, present within the GLPK software¹³. Writing the files for the solver and the optimization itself used less than one second of computational time on a laptop computer running Windows 7 with a dual-core Intel core i7 2.80 GHz processor with 12 GB of RAM. The default solver options (e.g., integrality tolerance, optimality gap, etc.) of GLPK/GLPSOL were used, which are described in the solver documentation (Free Software Foundation Inc., 2019).

While the case study presented in this section is fictitious, the approach and framework described herein is broadly applicable in realistic scenarios to find good solutions for process system integration at different scales and settings. The approach is applicable at multiple levels of detail and temporal resolution while considering numerous energy and material streams. This provides a flexible framework to solve real system integration problems in many domains using an optimization approach to provide valuable solutions for complex problems.

5. Recent Developments and Future Improvements

A comprehensive and succinct baseline modeling structure for site optimization and industrial symbiosis has been presented; however, additional aspects should be added to provide more flexibility in system analysis. Further improvements and additions to the methodology have already been completed or should be foreseen as additions to the basis formulation:

• Mass and thermal storage: A new unit type which will allow energy or material storage between times or periods will be formulated to allow flexibility over short or long time horizons and in batch processes. The storage unit will act as a sink with a variable but bounded capacity for energy or material from a specific layer, and act as a source in a future time or period. This requires sequentiality of times, as well as a cyclic constraint to ensure that the storage level returns to its initial state in the last time/period.

• Heat exchanger network design: The current model solves the heat cascade, thus enabling the identification of the best possible integration target, but without specifying the heat exchangers required to achieve it. Several methods have been proposed to design optimal heat exchanger networks, such as that proposed by Ciric and Floudas (1991), and one of the most recent should be added which will accurately account for the capital costs and configuration of the network.

• Steam networks: A steam network model will be added which allows optimization of heat loads, pressure levels and compression/expansion stages considering production of electricity through turbine expanders.

• Material pinch: Pinch analysis has been successfully applied to different material streams (Government of Canada, 2002) (hydrogen, water, compressed air, etc.), where the heat load is replaced by flowrate, and temperature by quality (such as impurities hydrogen, pressure for compressed air). Such a method could also be applied in the formulation presented in this paper for common resource streams.

• Transportation options and costs: Including transportation media and their associated cost would provide more accurate solutions for real symbiosis solutions considering plants which are not geographical co-located. Various transportation media would be considered (road, rail, ship, pipeline, conveyor belt), each of which would be associated to different costs and emissions as in Kantor et al. (2015). The cost and emissions would also be included in the calculation of the objective function. In the formulation presented here, x^rb and x^mb would be indexed over the set of transportation options with additional constraints considering the type of materials and its transportability by the different media.

• Coupling with GIS: Geographical coordinates should be introduced to assign real locations with units. The location of new units could also be optimized, as this would have an impact on the transport costs/piping networks involved. Coupling the optimization structure to a GIS database is closely linked to the consideration of transportation options.

• Investment planning: Modifications on the framework presented here will be proposed to allow the optimal selection of not only what technologies to invest in from the superstructure, but also which period to make the investment and how to cascade investments to prevent short-sighted planning or to account for availability of capital investment funds.

• Additional KPIs: Adding Key Performance Indicators (KPIs) that can either be calculated or added to the objective function. This can include indicators such as human health impact or job creation. Theoretically, any impact associated with technology selection from the superstructure could be included, provided that it can be linearized for application in the MILP formulation.

• Optimization under uncertainty: Though often proposed to be fixed values, parameters used in the optimization are subject to uncertainty. Taking this uncertainty into account in the optimization model would not only allow the identification of the optimal solution, but the most probable optimal solution given a variety of uncertain events. For industrial symbiosis, consideration of random unit failures can also assist in understanding the risk of inter-plant dependencies.

• Parametric optimization: This is a well-known method to carry out multi-objective optimization for a limited number of objectives which cannot be characterized by the same physical units, normalized or otherwise combined into a single objective. While one of the objectives is optimized, a constraint with varying bounds is applied on each of the other objectives. This technique permits the identification of many solutions (Pareto frontier) which are beneficial to some degree in unrelated areas and provides additional options for decision-makers to select from a set of optimal solutions.

• Flexibility optimization: In order to take into account the variable electricity availability of renewable sources like sun or wind, a methodology to assess the robustness of a system to real time changes and variations should be developed.

6. Conclusion

In this paper, a methodology implemented using MILP was presented, which optimizes process integration opportunities on different scales, providing solutions for integration within a plant and potential industrial symbiosis options. The presented model builds from many previous contributions in process integration and industrial process optimization. It deals with both material and energy streams of different types and incorporates pinch analysis fundamentals to obey the second law of thermodynamics. Given fixed requirements by processes, the optimization solves the sizing and scheduling problem simultaneously, resulting in the optimal set of technologies selected from a superstructure to fulfil the requirements (heating, waste treatment, etc.), their size and operation.

Several objective functions were proposed for use within this model, namely: operating costs, investment costs and environmental impacts of different kinds (e.g., CO₂ emissions) but other objectives could also be formulated. MILP formulations typically optimize a single objective and thus modification of the objective functions (e.g., by monetizing environmental impacts) also provides an option for assessing differing objectives within this framework. Additional methods, such as parametric optimization, are proposed in future work to consider multiple objectives without aggregation.

To illustrate the methodology, a fictional case study consisting of five industrial processes was constructed with potential for industrial symbiosis through the exchange of energy and materials. The methodology identifies the technologies necessary to fulfil the requirements of the processes in an optimal way, with respect to different objectives. It shows which utilities should supply heat to the different processes (boilers, CHPs), which utilities can be integrated with the heating profiles of the processes (heat pumps, Rankine cycles), which exchanges are best between processes and the associated technologies required to transfer and/or treat the recovered waste or co-products.

Numerous future additions to the model have been identified to improve the methodology, to better represent the systems and handle a wide variety of cases in industrial symbiosis studies. These include the addition of new types of technologies (e.g., storage, detailed steam network, heat exchanger network), the inclusion of new methodologies (e.g. material pinch), increasing temporal (investment planning) and spatial (coupling with GIS) granularity of the model and extending the possibilities of the model with multi-objective optimization (using parametric optimization) and optimization under uncertainty.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Author Contributions

IK and J-LR contributed equally to the conceptualizing, writing, and editing of the article. HB also contributed to allaspects of the article. FM contributed to the conceptualization and discussion.

Conflict 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.

Acknowledgments

The authors gratefully acknowledge the financial support from the EPOS and SCCER-EIP projects. The EPOS project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 679386. This research project is part of the Swiss Competence Center for Energy Research SCCER EIP of the Swiss Innovation Agency Innosuisse. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The authors also gratefully acknowledge the expert advice of Anna Sophia Wallerand and Maziar Kermani who provided critical review and contributions to improve the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenrg.2020.00049/full#supplementary-material

Abbreviations

LHV, lower heating value; HHV, Higher Heating Value; CEPCI, Chemical Engineering Plant Cost Index; HI, Heat Integration; PI, Process Integration; PA, Pinch Analysis; MILP, Mixed Integer Linear Programming; MINLP, Mixed Integer Non-Linear Programming; CC, Composite Curve; GCC, Grand Composite Curve; ICC, Integrated Composite Curve; CCC, Carnot Composite Curve; MP, Mathematical Programming; GA, Genetic Algorithm; SNG, Synthetic Natural Gas.

Footnotes

1. ^SNG: Synthetic natural gas.

2. ^CHP: Combined heat and power (natural gas).

3. ^SCHP: Combined heat and power (synthetic natural gas).

4. ^HP: Heat pump.

5. ^set ToP_p.

6. ^${\stackrel{\cdot }{\mathrm{\text{m}}}}_{l,u,p,t}^{out}\ge 0$ and ${\stackrel{\cdot }{\mathrm{\text{m}}}}_{l,u,p,t}^{in}\ge 0$.

7. ^Corrected temperatures used for pinch analysis (Linnhoff and Hindmarsh, 1983).

8. ^K is a set of ordered integers {1,2,3,…} (but not the temperature values themselves).

9. ^Equation (32) is expressed in such a way that the heat balance is closed above the lower temperature of each interval. The second law is respected since the balance is applied to all intervals starting from the top.

10. ^Streams are assumed to belong to the same hc layer.

11. ^With the emissions optimization, the comparison is unfair since no CO₂ emissions are associated with the WasteWaterEx unit.

12. ^They are, in fact, representations of the x^rb and x^mb variables.

13. ^https://www.gnu.org/software/glpk/

References

Al-Malah, K. I. (2016). “Introducing aspen plus,” in Aspen PlusⓇ (Hoboken, NJ: John Wiley & Sons, Inc), 1–47.

Google Scholar

Alves, J. J., and Towler, G. P. (2002). Analysis of refinery hydrogen distribution systems. Indus. Eng. Chem. Res. 41, 5759–5769. doi: 10.1021/ie010558v

CrossRef Full Text | Google Scholar

Andiappan, V. (2017). State-of-the-art review of mathematical optimisation approaches for synthesis of energy systems. Process. Integr. Optim. Sustain. 1, 165–188. doi: 10.1007/s41660-017-0013-2

CrossRef Full Text | Google Scholar

Aviso, K. B., Tan, R. R., Culaba, A. B., and Cruz, J. B. (2010). Bi-level fuzzy optimization approach for water exchange in eco-industrial parks. Process Safe. Environ. Prot. 88, 31–40. doi: 10.1016/j.psep.2009.11.003

CrossRef Full Text | Google Scholar

Becker, H., and Maréchal, F. (2012). Targeting industrial heat pump integration in multi-period problems. Comput. Aid. Chem. Eng. 31, 415–419. doi: 10.1016/B978-0-444-59507-2.50075-5

CrossRef Full Text | Google Scholar

Becker, H., Vuillermoz, A., and Maréchal, F. (2012). Heat pump integration in a cheese factory. Appl. Therm. Eng. 43, 118–127. doi: 10.1016/j.applthermaleng.2011.11.050

CrossRef Full Text | Google Scholar

Becker, H. C. (2012). Methodology and thermo-Economic optimization for integration of industrial heat pumps (Ph.D. thesis). EPFL, Lausanne, Switzerland.

Google Scholar

Bell, I. H., Wronski, J., Quoilin, S., and Lemort, V. (2014). Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library coolProp. Ind. Eng. Chem. Res. 53, 2498–2508. doi: 10.1021/ie4033999

PubMed Abstract | CrossRef Full Text | Google Scholar

Belsim (2018). Software – Belsim Vali v4.

Google Scholar

Boix, M., Montastruc, L., Azzaro-Pantel, C., and Domenech, S. (2015). Optimization methods applied to the design of eco-industrial parks: a literature review. J. Clean. Prod. 87, 303–317. doi: 10.1016/j.jclepro.2014.09.032

CrossRef Full Text | Google Scholar

Bolliger, R. (2010). Méthodologie de la synthése des systémes énergétiques industriels, n 4867 (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Google Scholar

Bungener, L. G. S. (2016). Energy efficiency and integration in the refining and petrochemical industries, n 7038 (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Google Scholar

Bütün, H., Kantor, I., and Maréchal, F. (2019). An optimisation approach for long-term industrial investment planning. Energies 12:4076. doi: 10.3390/en12214076

CrossRef Full Text | Google Scholar

Bütün, H., Kantor, I., Mian, A., and Maréchal, F. (2018). “A heat load distribution method for retrofitting heat exchanger networks,” in Computer Aided Chemical Engineering, volume 43 of 28 European Symposium on Computer Aided Process Engineering, eds A. Friedl, J. J. Klemeš, S. Radl, P. S. Varbanov, and T. Wallek (Graz: Elsevier), 1395–1400.

Google Scholar

Ciric, A. R., and Floudas, C. A. (1991). Heat exchanger network synthesis without decomposition. Comput. Chem. Eng. 15, 385–396. doi: 10.1016/0098-1354(91)87017-4

CrossRef Full Text | Google Scholar

Dhole, V. R., and Linnhoff, B. (1993). Total site targets for fuel, co-generation, emissions, and cooling. Comput. Chem. Eng. 17, S101–S109.

Google Scholar

Edward Hohmann (1971). Optimum networks for heat exchange (Ph.D. thesis). University of Southern California, Los Angeles, CA, United States.

Google Scholar

El-Halwagi, M. M., and Manousiouthakis, V. (1989). Synthesis of mass exchange networks. AIChE J. 35, 1233–1244.

Google Scholar

Fan, Y. V., Chin, H. H., Klemeš, J. J., Varbanov, P. S., and Liu, X. (2020). Optimisation and process design tools for cleaner production. J. Clean. Prod. 247:119181. doi: 10.1016/j.jclepro.2019.119181

CrossRef Full Text | Google Scholar

Fourer, R., Gay, D. M., and Kernighan, B. W. (1990). A modeling language for mathematical programming. Manage. Sci. 36, 519–554. doi: 10.1287/mnsc.36.5.519

CrossRef Full Text | Google Scholar

Free Software Foundation Inc. (2019). GLPK - GNU Project - Free Software Foundation (FSF) (Boston, MA).

Google Scholar

Gassner, M., and Maréchal, F. (2009). Methodology for the optimal thermo-economic, multi-objective design of thermochemical fuel production from biomass. Comput. Chem. Eng. 33, 769–781. doi: 10.1016/j.compchemeng.2008.09.017

CrossRef Full Text | Google Scholar

Girardin, L. (2012). A gis-based methodology for the evaluation of integrated energy systems in urban area, n 5287 (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Google Scholar

Girones, V. C., Moret, S., Maréchal, F., and Favrat, D. (2015). Strategic energy planning for large-scale energy systems: a modelling framework to aid decision-making. Energy 90, 173–186. doi: 10.1016/j.energy.2015.06.008

CrossRef Full Text | Google Scholar

Gonela, V., Zhang, J., and Osmani, A. (2015). Stochastic optimization of sustainable industrial symbiosis based hybrid generation bioethanol supply chains. Comput. Indus. Eng. 87, 40–65. doi: 10.1016/j.cie.2015.04.025

CrossRef Full Text | Google Scholar

Government of Canada (2002). Pinch Analysis - For the Efficient Use of Energy, Water and Hydrogen. Technical Report M154-60/1-2012E, Government of Canada.

Google Scholar

Hashimoto, S., Fujita, T., Geng, Y., and Nagasawa, E. (2010). Realizing CO2 emission reduction through industrial symbiosis: a cement production case study for Kawasaki. Resour. Conserv. Recycl. 54, 704–710. doi: 10.1016/j.resconrec.2009.11.013

CrossRef Full Text | Google Scholar

Henchoz, S. (2016). Potential of refrigerant based district heating and cooling networks, n 6935 (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Google Scholar

Hubert, T., Maréchal, F., Viart, M., Griffay, G., and Jumel, S. (2009). “Process integration analysis to save energy in large integrated iron and steel plants,” in Energy Efficiency and CO2 Reduction in the Steel Industry – Event Presentation (Budapest).

Google Scholar

Jacobsen, N. (2006). Industrial symbiosis in Kalundborg, Denmark: a quantitative assessment of economic and environmental aspects. J. Indust. Ecol. 10, 239–255. doi: 10.1162/108819806775545411

CrossRef Full Text | Google Scholar

Kantor, I., Betancourt, A., Elkamel, A., Fowler, M., and Almansoori, A. (2015). Generalized mixed-integer nonlinear programming modeling of eco-industrial networks to reduce cost and emissions. J. Clean. Prod. 99, 160–176. doi: 10.1016/j.jclepro.2015.03.017

CrossRef Full Text | Google Scholar

Kastner, C. A., Lau, R., and Kraft, M. (2015). Quantitative tools for cultivating symbiosis in industrial parks; a literature review. Appl. Energy 155, 599–612. doi: 10.1016/j.apenergy.2015.05.037

CrossRef Full Text | Google Scholar

Kermani, M., Kantor, I. D., and Maréchal, F. (2019a). Optimal design of heat-integrated water allocation networks. Energies 12:2174. doi: 10.3390/en12112174

CrossRef Full Text | Google Scholar

Kermani, M., Kantor, I. D., Wallerand, A. S., Granacher, J., Ensinas, A. V., and Maréchal, F. (2019b). A holistic methodology for optimizing industrial resource efficiency. Energies 12:1315. doi: 10.3390/en12071315

CrossRef Full Text | Google Scholar

Kermani, M., Wallerand, A. S., Kantor, I. D., and Maréchal, F. (2018). Generic superstructure synthesis of organic Rankine cycles for waste heat recovery in industrial processes. Appl. Energy 212, 1203–1225. doi: 10.1016/j.apenergy.2017.12.094

CrossRef Full Text | Google Scholar

Kim, J.-K., and Smith, R. (2001). Cooling water system design. Chem. Eng. Sci. 56, 3641–3658. doi: 10.1016/S0009-2509(01)00091-4

CrossRef Full Text | Google Scholar

Klemeš, J. J., and Kravanja, Z. (2013). Forty years of heat integration: pinch analysis (PA) and mathematical programming (MP). Curr. Opin. Chem. Eng. 2, 461–474. doi: 10.1016/j.coche.2013.10.003

CrossRef Full Text | Google Scholar

Kuznetsova, E., Zio, E., and Farel, R. (2016). A methodological framework for Eco-Industrial Park design and optimization. J. Clean. Prod. 126, 308–324. doi: 10.1016/j.jclepro.2016.03.025

CrossRef Full Text | Google Scholar

Leyland, G. B. (2002). Multi-objective optimisation applied to industrial energy problems, n 2572 (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Google Scholar

Li, J., Demirel, S. E., and Hasan, M. M. F. (2018). Process integration using block superstructure. Ind. Eng. Chem. Res. 57, 4377–4398. doi: 10.1021/acs.iecr.7b05180

CrossRef Full Text | Google Scholar

Linnhoff, B. (1972). Thermodynamic analysis of the cement burning process (Ph.D. thesis). Dissertation for the Doctoral Degree. ETH Zurich (in German), Zürich, Switzerland.

Google Scholar

Linnhoff, B., and Hindmarsh, E. (1983). The pinch design method for heat exchanger networks. Chem. Eng. Sci. 38, 745–763.

Google Scholar

Liu, P., Georgiadis, M. C., and Pistikopoulos, E. N. (2011). Advances in energy systems engineering. Ind. Eng. Chem. Res. 50, 4915–4926. doi: 10.1021/ie101383h

CrossRef Full Text | Google Scholar

Maréchal, F., and Boursier, I. (1989). Synep1:a methodology for energy integration and optimal heat exchanger network synthesis. Comput. Chem. Eng. 13, 603–610.

Google Scholar

Maréchal, F., and Kalitventzeff, B. (1996). Targeting the minimum cost of energy requirements: a new graphical technique for evaluating the integration of utility systems. Comput. Chem. Eng. 20, S225–S230.

Google Scholar

Maréchal, F., and Kalitventzeff, B. (1997). Identification of the optimal pressure levels in steam networks using integrated combined heat and power method. Comput. Chem. Eng. 52, 2977–2989.

Google Scholar

Maréchal, F., and Kalitventzeff, B. (1999). Targeting the optimal integration of steam networks: Mathematical tools and methodology. Comput. Chem. Eng. 23, 133–136.

Google Scholar

Maréchal, F., and Kalitventzeff, B. (2003). Targeting the integration of multi-period utility systems for site scale process integration. Appl. Therm. Eng. 23, 1763–1784. doi: 10.1016/S1359-4311(03)00142-X

CrossRef Full Text | Google Scholar

Mian, A., Martelli, E., and Maréchal, F. (2016). Framework for the multiperiod sequential synthesis of heat exchanger networks with selection, design, and scheduling of multiple utilities. Ind. Eng. Chem. Res. 55, 168–186. doi: 10.1021/acs.iecr.5b02104

CrossRef Full Text | Google Scholar

Miller, F. P., Vandome, A. F., and McBrewster, J. (2009). Lua (Programming Language). Orlando, FL: Alpha Press

Google Scholar

Moret, S., Codina Gironès, V., Bierlaire, M., and Maréchal, F. (2017). Characterization of input uncertainties in strategic energy planning models. Appl. Energy 202(Suppl. C), 597–617. doi: 10.1016/j.apenergy.2017.05.106

CrossRef Full Text | Google Scholar

Moret, S., Codina Gironès, V., Maréchal, F., and Favrat, D. (2014). “Swiss-energyscope. ch: a platform to widely spread energy literacy and aid decision-making,” in Pres 2014, 17th Conference on Process Integration, Modelling and Optimisation for Energy Saving and Pollution Reduction, Pts 1-3, Vol.39 (Prague: Aidic Servizi Srl), 877–882.

Google Scholar

Moret, S., Peduzzi, E., Gerber, L., and Maréchal, F. (2016). Integration of deep geothermal energy and woody biomass conversion pathways in urban systems. Energy Conver. Manage. 129, 305–318. doi: 10.1016/j.enconman.2016.09.079

CrossRef Full Text | Google Scholar

Muller, D. (2007). Web-based tools for energy management in large companies applied to food industry, n 3785 (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Google Scholar

Ng, R. T. L., Ng, D. K. S., and Tan, R. R. (2015). Optimal planning, design and synthesis of symbiotic bioenergy parks. J. Clean. Prod. 87, 291–302. doi: 10.1016/j.jclepro.2014.10.045

CrossRef Full Text | Google Scholar

NIST (2013). REFPROP. database/software. Gaithersburg, MD.

Google Scholar

Olesen, S. G., and Polley, G. T. (1997). A simple methodology for the design of water networks handling single contaminants. Chem. Eng. Res. Design 75, 420–426.

Google Scholar

Papoulias, S. A., and Grossmann, I. E. (1983a). A structural optimization approach in process synthesis-I: utility systems. Comput. Chem. Eng. 7, 695–706.

Google Scholar

Papoulias, S. A., and Grossmann, I. E. (1983b). A structural optimization approach in process synthesis-II: heat recovery networks. Comput. Chem. Eng. 7, 707–721.

Google Scholar

Peduzzi, E. (2015). Biomass to liquids: thermo-economic analysis and multi-objective optimisation, n 6529 (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Google Scholar

Pouransari, N. (2015). Towards practical solutions for energy efficiency of large-scale industrial sites, n 6390 (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Google Scholar

Process Systems Enterprise (PSE) (1997). gPROMS. London: Software – Process Systems Enterprise gProms.

Google Scholar

Rager, J. M. F. (2015). Urban energy system design from the heat perspective using mathematical programming including thermal storage, n 6731 (Ph.D. thesis). École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Google Scholar

Suciu, R., Girardin, L., and Maréchal, F. (2018). Energy integration of CO2 networks and power to gas for emerging energy autonomous cities in Europe. Energy 157, 830–842. doi: 10.1016/j.energy.2018.05.083

CrossRef Full Text | Google Scholar

Towler, G. P., Mann, R., Serriere, A. J., and Gabaude, C. M. (1996). Refinery hydrogen management: cost analysis of chemically-integrated facilities. Indus. Eng. Chem. Res. 35, 2378–2388.

Google Scholar

Valenzuela-Venegas, G., Vera-Hofmann, G., and Dìaz-Alvarado, F. A. (2020). Design of sustainable and resilient eco-industrial parks: planning the flows integration network through multi-objective optimization. J. Clean. Prod. 243:118610. doi: 10.1016/j.jclepro.2019.118610

CrossRef Full Text | Google Scholar

Wallerand, A. S., Kermani, M., Kantor, I., and Maréchal, F. (2018a). Optimal heat pump integration in industrial processes. Appl. Energy 219, 68–92. doi: 10.1016/j.apenergy.2018.02.114

CrossRef Full Text | Google Scholar

Wallerand, A. S., Kermani, M., Voillat, R., Kantor, I., and Maréchal, F. (2018b). Optimal design of solar-assisted industrial processes considering heat pumping: case study of a dairy. Renew. Energy 128, 565–585. doi: 10.1016/j.renene.2017.07.027

CrossRef Full Text | Google Scholar

Wang, Y. P., and Smith, R. (1994). Wastewater minimisation. Chem. Eng. Sci. 49, 981–1006.

Google Scholar

Yoo, M.-J., Lessard, L., Kermani, M., and Maréchal, F. (2015). “OsmoseLua-an integrated approach to energy systems integration with LCIA and GIS,” in 12th International Symposium on Process Systems Engineering and 25th European Symposium on Computer Aided Process Engineering, Vol. 37 (Copenhagen: Elsevier), 587–592.

Google Scholar

Nomenclature