# Flexible networks to promote the development of industrial symbioses: A new optimization procedure

Marianne Boix, Stéphane Negny, Ludovic Montastruc, Florent Mousqué

► To cite this version:

Marianne Boix, Stéphane Negny, Ludovic Montastruc, Florent Mousqué. Flexible networks to promote the development of industrial symbioses: A new optimization procedure. Computers & Chemical Engineering, 2023, 169, pp.108082. 10.1016/j.compchemeng.2022.108082. hal-03997208

HAL Id: hal-03997208

https://hal.science/hal-03997208v1

Submitted on 6 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

Flexible networks to promote the development of industrial

2

symbioses: a new optimization procedure

3

Marianne Boix*, Stéphane Négny, Ludovic Montastruc, Florent Mousqué

4

Laboratoire de Génie Chimique, UMR 5503 CNRS, Toulouse INP, UPS, 4 Allée Emile Monso, 31432 Toulouse

5

Cedex 4, France

6

Corresponding Author : marianne.boix@ensiacet.fr

7

8

Abstract

The development of eco-industrial parks is now proven a local efficient solution to reduce environmental impacts of industrial areas as well as increasing economic benefits. Furthermore, they are a good way to produce goods while reducing environmental footprints by exchanging energies, water and raw materials instead of consuming too much resources; that way it improves circularity of resources with no disruption of production. However, in an uncertain context caused by different factors such as environmental problems, sanitary crises or geopolitical conflicts, the design of an interconnected network must take into account possible deviations of nominal conditions. In this study, an integrated procedure is developed so that an optimal flexible network is designed with a given flexibility threshold. The flexibility is defined as the minimum of the maximal deviations for the nominal inlet parameters that can be handle by the network. Impact of such a flexibility is visible on the total investment cost of the network. Indeed, when adding some flexibility to the network, the investment cost is increased. The integrated procedure developed is able to find a global optimum in a very large computational time (up to 500 hours); in this article, two other procedures are then tested over a large-scale case study (15 companies) to validate the approach. One of them consists in oversizing the network so that it can handle large deviations but the increase of the cost is important. The other procedure consists in a decomposition of the integrated approach in two steps: the first one aims at defining the structure of the nominal solution (minimum of the cost) and the second one evaluates several deviation scenarios to take into account disruption of inlet parameters. The results show that the two-steps procedure obtains results close to the integrated procedure with reasonable computational time (up to 40 minutes). This study helps to provide solutions to design optimal flexible networks, also adapted to solve large-scale eco-industrial parks; it can also be used as an optimization-decision aid tool for the territorial development of industrial ecology.
31

# 1. Introduction

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

Industrial symbioses are now well-known to be organizations able to change our linear production systems into circular ones. Eco-industrial parks (EIP), which are a typical application of industrial ecology, propose to gather different companies in a same geographic area so that they can exchange flows of water, energies, raw material or services, equipment and infrastructures (Hartley et al., 2020). The global warming and more generally the current environmental context or other punctual crises (such as COVID pandemic) can lead to several uncertain phenomena that can compromise availability of resources, disturbance in supply chains or sudden increase/decrease of production for example (Alkahtani et al., 2021). In this context, when several companies are interconnected, the question of adaptability of these networks in a changing and uncertain environment becomes relevant (Xu et al., 2021). The design of eco-industrial parks, also called industrial symbioses in the literature, needs to use optimization methods in order to achieve environmental and economic benefits for the different participants at a meso-scale (Kristensen and Moosgard, 2020). However, complex systems such as inter-industry networks are uncertain, so they often face unplanned disruption. Therefore, the design of flexible networks can be crucial to ensure the operational efficiency of an exchange network deviating from these nominal operating conditions (Grossmann and Morari, 1983). Hence, the development of procedures using mathematical optimization to design flexible networks needs to be developed according to the review of Andiappan (2017).

49

50

51

52

53

54

55

56

57

58

The recent literature has handled this field and the emergence of discussions about definition of concepts arise. The resilience is commonly defined as the ability to recover from and/or adjust to a disturbance or change (Meerow and Newell, 2015). The flexibility is more an indicator of resilience such as adaptability, diversity or efficiency. Flexibility commonly refers to the ability of a given system to cope with a fluctuating context, either actively or passively (that is, with or without modification of the system). In the case of an industrial system, it is the capacity of a system to adapt to a set of uncertain parameters (Hoch and Eliceche, 1996). The measure of the flexibility of a system is therefore the measure of the maximum value that its input parameters can reach while ensuring that its output variables remain within their operating range (Figure 1), delimited by the dimensioning of the system.

<figure><img src="image_2.png" /></figure>

59

60

<caption>Figure 1. Illustration of the flexibility</caption>
The range in which the input variables can vary is then called the feasibility domain of the system. The exact determination of this domain is not a trivial problem, because of the large number of constraints, parameters, variables, and interrelationships within the system.

To go further, operational efficiency, also known as operability, corresponds to the ability of a system to cope with disturbance scenarios while meeting the performance criteria defined during the design phase (Sharifzadeh, 2013). Operability therefore calls for different scientific concepts: flexibility, controllability, reliability and security of industrial systems (Grossmann and Morari, 1983). However, the aspects of controllability, reliability and safety of the systems will not be detailed in this study, as they are generally the subject of a detailed design on a smaller scale. In an eco-industrial park, the notion of operational efficiency will be considered as the ability of a system to be flexible in a real operational environment while optimizing the various design criteria adopted. The development of EIP's also suffers from the fear of companies to be interconnected each other. This interdependency between companies is more important when the flow crossing the interconnection is high. In order to evaluate the interdependence of industries throughout the EIP, the network interconnections are studied according to the number of interconnections but also the distribution of flows exchanged by the different interconnections (Mousqué et al., 2020).

By considering these various definitions, the aim of this study consists in integrating a flexibility indicator at the design stage of industrial symbioses. The necessity of developing optimization approaches to design optimal resources network in eco-industrial parks has been widely emphasized in previous studies. However, the optimal networks obtained are often not flexible (Montastruc et al., 2013); this study proposes, after a literature review, the development of an optimization procedure allowing the design of flexible resources exchanges network for the development of eco-industrial parks. In our previous work (Mousqué et al., 2020) an optimization model has been developed to handle the coupling between the utility system with the hybrid power system to offer different possibilities of energy production in an eco-industrial park. In that previous study, a focus has been made on an interdependence indicator. This tool helps to take multicriteria decision making when a design must be chosen. In this study, another main bottleneck is addressed: the flexibility, that remains a barrier when an eco-industrial park is studied. Flexibility is typically an optimization application challenge and concentrates a lot of difficulties in terms of modelling, formulation and numerical resolution. The development of an optimization procedure has been tested and validated on a large-scale case study consisting in an eco-industrial park with 15 companies. The aim is to develop a generic approach able to evaluate the cost of the flexibility to serve at the decision stage.

93
## 2. Bibliographic Analysis

The literature search on the terms: "design", "eco-industrial parks" and "flexibility" led to only a small number of publications that actually address the subject of flexible eco-industrial park design.

First, the article of Aguilar et al (2008) that deals with the availability and reliability of industrial utility systems through the design of flexible networks. The method used to make the utility system flexible consists in optimizing the system throughout several operating scenario. The authors developed an integrated approach where sizing and operational variables are optimized simultaneously. The optimal network is then determined among the different possibilities obtained in an integrated step through a MILP (Mixed Integer Linear Programming) approach.

Next, Valenzuela-Venegas et al. (2018) introduce in a first paper an indicator to assess the resilience of EIPs. This indicator is based on two main characteristics of EIP: the number of connections between participants, and the ability of each flow to change in magnitude when a participant suddenly stops sharing with other industries. In a second paper, Valenzuela et al. (2020) propose a multi-objective optimization model for the design of resilient industrial networks using this indicator. The optimization is then performed on the minimization of the economic criterion, the environmental impact and the maximization of the flexibility criteria developed in the previous article.

Among the techniques that exist to make an energy network flexible in an active way, we can mention the Demand-Side Management techniques, which consist in piloting a part of the demand to face production shortages or demand peaks. However, these techniques are not yet widely used, as they are only applicable to a limited number of demands (storage, controllable processes and non-contracts).

Regarding the methods commonly used to design a flexible network, a classically employed approach is to conceive a network capable of meeting a scenario where all the demands are simultaneously at their maximum value (Bandyopadhyay 2011). This approach usually leads to a sub-optimal and over-designed solution, because most of the time, the input parameters vary independently of each other, and are sometimes antagonist. The scenario with all the disturbances at their maximum value is therefore not likely to occur.

Another method is to identify the limiting facility. This facility is then considered as the blocking element that reduces the overall flexibility. To find this limiting element, the flexibility is evaluated at different thresholds and the facilities working at higher conditions than their feasible value are observed. To increase flexibility, the system is then re-designed in an iterative process (Schneider,
127

1997). Given that this approach is iterative, based on “trial and error”, and may involve an important number of iterations, we can expect the optimality of the final solutions to be far from the global optimum allowed by an integrated method.

128

129

130

Another previous study conducted by Montastruc et al., (2013) proposed a design method to manage the flexibility of industrial water networks. The results indicated that optimal design methods without considering flexibility did not lead to flexible solutions. The flexibility studied for the water system was the maximum deviation of the pollutant flow from the different processes. The solution to increase this flexibility was to increase the number of interconnections. This was done at the expense of the overall cost, the number of interconnections is then increased in an iterative process, until a certain profitability threshold is reached for the ecopark.

131

132

133

134

135

136

137

The first publication presenting an indicator of flexibility was written by Swaney and Grossmann (1985). This indicator, named Swaney and Grossman Flexibility Index (FSG), evaluates the smallest maximum deviation that the fluctuating parameters can bear. For this purpose, all possible combinations of parameter deviations are evaluated one by one. For each combination, the maximum positive and negative deviation of the input parameters that ensures the feasibility of the network is determined. Among all these values, the smallest recorded is thus equal to the value of the FSG indicator. Since, for each combination, the maximum positive and negative deviation is studied, the number of deviations to be studied is equal to $2^n$, with $n$ the number of parameters considered.

138

139

140

141

142

143

144

145

146

Shortly after, Saboo et al (1985) proposed an indicator called the resilience index (RI). They define this indicator as a measure of the capacity of a system to adapt to an unexpected event or change. Although called resilience, this indicator evaluates this capacity of adaptation nor by modifying the properties of the designed system but by looking for the smallest maximum deviation of the input parameters subject to variations. This indicator mathematically evaluates the same criterion as the one introduced previously (FSG), the difference residing in the method used to vary the different input parameters. Indeed, the RI considers the deviation of only one parameter at a time, which results in a total of $2^n$ deviations to be measured. Compared to FSG, the computational efforts of RI are therefore significantly reduced, especially for a system with a large number of input parameters.

147

148

149

150

151

152

153

154

155

The Dynamic Flexibility (DF) indicator (Pistikopoulos et al., 1995) was developed to evaluate the flexibility of dynamic systems. In fact, the previous indicators were dedicated to the study of steady state systems. They do not take into account the transient functioning between these states. The DF indicator is identical in principle to the FSG, its particularity is to integrate the time variable to its input parameters, adapted to dynamic models. Since the time variable is continuous, it can be

156

157

158

159
160

considered the number of decision variables is infinite. To solve this issue, another approach
consisting in reducing the dimension of the problem has been developed.

161

162

The previous indicators do not represent reality because they consider that each failure has the same
probability of occurrence. Therefore, the stochastic flexibility (SF) indicator of Pistikopoulos and
Mazzuchi (1990) has been developed and it allows to associate to each deviation a probability. While
the previous indicators were dedicated to deterministic models, the SF applies to stochastic models.

163

164

165

166

A probability distribution is then associated to each uncertain input parameter. With this indicator, a weighted estimate of the flexibility is then provided, the deviations having a small chance of appearing will then have only little influence in the measurement of the indicator. The main drawback of this indicator compared to the others is the need to know the system in depth and to be able to evaluate the probability of occurrence of each deviation.

167

168

169

170

171

More recently, the flexibility of an industrial system has been studied by Di Pretoro et al. (2019) who
proposed an extensive comparison of the different indexes applied to the design of distillation
columns. The authors pointed out the main advantages and drawbacks of the different indexes.

172

173

174

Regarding resilience, the term has been more largely employed for larger scale systems such as
supply chains for instance. Chrisandina et al. (2022), in a review article, proposed to define the
resilience according to several definitions of several authors. For them, resilience is not an index but
"a set of capabilities a supplychain needs to demonstrate during a disruptive event". It consists then
to prepare and anticipate disruptions, absorb the impact of disruption, adjust operations and recover
operations.

175

176

177

178

179

180

In the industrial ecology field, the levels of details are needed at the process scale although the
adaptability to fluctuations is modelled and evaluated at a larger scale. Regarding all these definitions
and points of view, in this study, we propose to talk about flexibility for the eco-industrial park
network. Flexibility of an EIP is defined as the ability of a design to handle inlet parameters
fluctuations without changing its structure. That is to say, the connections already defined can
overcome an over-production or other events without changing the connection itself or any other
component of the network. A classical approach to take into account flexibility as it has been defined
is to design a larger network what leads to a great increase of the overall cost. Here it is proposed to
develop an approach able to estimate the cost of the design for a flexible network.

181

182

183

184

185

186

187

188

189
190

## 3. Scientific challenges and problem statement

191

The scientific challenge addressed in this study is to develop approaches to design a flexible network that remains optimal with respect to the design criteria in a changing operational context. The multiobjective optimization procedure should therefore minimize the interdependence between industries while optimizing the other criteria targeted by the EIP actors such as the cost for example.

192

193

194

195

However, the mathematical program can be difficult to solve using a classical multiobjective method, especially when the EIP considered is large and complex. Indeed, during the optimization stage, despite a large computation time, the size of the possibility tree exceeded often the available computational memory. Moreover, the binary variables used to define the existence of interconnections between companies have been identified as super variables that can prevent to find an optimal solution.

196

197

198

199

200

201

To achieve this, the first objective of this study is to develop a procedure for the design of an optimal and flexible eco-industrial network. The second goal is to propose a repository based on the different procedures imagined, so that the flexibility of the different solutions obtained will then be analysed and compared.

202

203

204

205

The objective of the developed procedures is to design a network able to supply the industries despite a deviation of their demand, while remaining optimal with regard to the various optimization criteria. This flexibility study therefore focuses on the maximum deviation in input parameters that a designed industrial network is able to withstand. The reader can note that in this study, the expression "maximum deviation" refers to a parameter, whereas the notion of "flexibility" refers more broadly to the whole system. The optimization problem is stated as follows:

206

207

208

209

210

211

Given a set of input parameters and industrial demands to satisfy, the aim is to find the design variables of the eco-industrial network constituted by production, conversion and transport facilities of the resource and energy flows, by minimizing the overall cost and maximizing the flexibility of the network over different time periods. The input parameters studied are the variations in resource and energy demands of the different companies in the EIP. Several procedures are developed in order to provide a benchmark, to compare the approaches and the results obtained with the best procedure are detailed. The mathematical model including energy, water and mass balances is based on the superstructure described in figure 2.

212

213

214

215

216

217

218

219
<figure><img src="image_3.png" /><figcaption>Figure 2. Superstructure of the network</figcaption></figure>

This research work aims at developing a generic procedure, i.e. a procedure that is widely applicable to inter-industrial physical flow networks (resources, energy, water, waste, etc.). The case study used is presented in section 4.3, and the different results are shown and the procedures are discussed in section5. The discussion proposes an analysis of the procedures developed and of the results obtained.

## 4. Optimization procedures

### 4.1. Mathematical formulation

The mathematical model for the design of the exchanges in the eco-industrial park is formulated as follows:

#### a. Cost objective function

A Net Present Cost is calculated to evaluate the overall cost of the network; it is constituted by actualized investment and operational costs:

$$ \min \text{ NPC} = (Cost^{Raw} + Cost^{Boilers} + Cost^{Turbines} + Cost^{Pipes} + Cost^{Elec\_sources}) \quad (1) $$

Each component of the Net Present Cost is detailed below:

The cost of raw materials includes fuel, electricity and water:

$$ Cost^{Raw} = (Cost^{Fuel} + Cost^{Elec} + Cost^{Water}) \times DFA \quad (2) $$
The detailed costs composing the cost of raw materials are detailed below:

$$240 \qquad Cost^{Fuel} = \sum_{b,f,t,c} F_{b,f,t,c}^{purchase} \times Pri_{f}^{fuel} \qquad (3)$$

Regarding the cost of consumed electricity, it is possible to purchase or to sell electricity:

$$242 \qquad Cost^{Elec} = (\sum_{t,c} E_{t,c}^{purchase} \times PriElec^{purchase} - \sum_{t} E_{t}^{prodSell} \times PriElec^{sales}) \qquad (4)$$

Water is purchased for the needs of the steam network and some feedings can be done during the project to compensate water leaks or losses:

$$245 \qquad Cost^{Water} = (\sum_{t,c} W_{t,c}^{purchase} + \sum_{t,c} W_{t,c}^{feeding}) \times Pri^{water} \qquad (5)$$

Boilers and turbines represent important parts of the total cost with both the investment needed (i.e. CAPEX) and the operational cost (i.e. OPEX).

$$249 \qquad Cost^{Boilers} = (\sum_{b,c} y_{b,c}^{sel} \times CAPEX_{b} \times DFA) + (\sum_{b,t,c} y_{b,t,c}^{prod} \times OPEX_{b} \times DFA) \qquad (6)$$

$$250 \qquad Cost^{Turbines} = \sum_{tu,c} y_{tu,c}^{sel} \times (CAPEX_{tu} \times DFI_{tu} + OPEX_{tu} \times DFA) \qquad (7)$$

The discount factor for allocation of annual operational expenditure and resources over the life of the project (DFA) is applied to the OPEX (Eq. 8), with n the project duration.

$$253 \qquad DFA = \frac{(1+i)^{n-1}}{i(1+i)^n} \qquad (8)$$

The diameter of the pipe is determined by the highest flowrate associated to the connection over the different periods.

$$255 \qquad Cost^{Pipes} = \sum_{h,c,c'} y_{h,c,c'}^{pipe} \times CAPEX_{pipe}^{fix} + \sum_{h,c,c'} Interc_{h,c,c'}^{Diameter} \times CAPEX_{pipe}^{var} \times DFA + \sum_{h,c,c'} Interc_{h,c,c'}^{Diameter} \times OPEX_{pipe} \times DFA \qquad (9)$$

The cost of electricity production is composed of investment and operational costs depending on the nominal installed capacity $P_{sr}^{rated}$ (9).

$$257 \qquad Cost^{Elec\_sources} = \sum_{sr} (P_{sr}^{rated} \times DFA \times (CAPEX_{sr} + OPEX_{sr})) \qquad (10)$$

b. Mathematical formulation

For the sake of clarity, the whole mathematical model is proposed in Appendix 2. The constraints are based on mass and energy balances for each element of the superstructure (Figure 2). Physical and thermodynamic constraints are imposed to the different equipment such as turbines, boilers or steam headers. Figure 3 proposes to sum up this formulation. while the corresponding equations can be found in detail in our previous article The aim of this paper is to present a procedure to solve this MILP (Mixed Integer Linear Programming) model by taking into account the overall flexibility of the network at the optimal solution. During the development of the procedure, the different constraints
cited in Figure 3 will be changed and completed to evaluate the interdependency between companies and flexibility.

<figure><img src="image_4.png" /><figcaption>Figure 3: Formulation of the optimization problem as an MILP</figcaption></figure>

## 4.2. Evaluating the flexibility over time periods

Among the indicators described in section 2, the resilience index (RI) has been chosen for taking into account the flexibility of the solutions obtained. Indeed, this indicator does not require prior knowledge of the system behaviour (unlike the SF indicator), so in a search for genericity it is applicable to a larger number of cases.

The resilience index RI is defined as the smallest maximum deviation (DevMax) of the input parameters i that a system is able to withstand without becoming infeasible, the deviations of the parameters are studied independently of each other. Its mathematical expression is presented in equation (11):

$$RI = \min_{i}(DevMax_{i}) \qquad (11)$$

To obtain the maximum deviation of one parameter, the network sizing variables as well as all input parameters are fixed except the one under study. In the case of a positive deviation, the value of the studied parameter is then maximized by varying the operational variables (generated flows, converted flows, flowrates in the networks). In the case of a negative deviation, the operation is
290 identical by minimizing the value of the parameter. DevMax is calculated by taking the ratio 291 between the admissible deviation and the value of the parameter at nominal conditions as described 292 in equation (12).

$$293 \quad DevMax = \frac{\text{Allowable deviation}}{\text{Nominal Value}} \qquad (12)$$

294 Evaluating RI and DevMax relative to the nominal values, allows each parameter and more broadly 295 each company included in the EIP to be considered with the same relative importance.

296 Compare to the FSG, the use of RI is also better in terms of computational times and more suitable 297 for the procedure developed. For example, for n parameters studied, there are 2n possible 298 deviations with the RI (one positive and one negative deviation, with respect to the nominal, for each 299 parameter n). Whereas with the FSG indicator, 2n deviations need to be considered, because this 300 indicator takes into account all the possible combinations of variation of the parameters. For a large 301 number of parameters studied, such as the numerous demands of the large network of an EIP, the RI 302 index is therefore the most suitable. Finally, since the FSG indicator studies all the combinations, it 303 provides a more reductive and conservative value of flexibility than the RI indicator.

304 The nominal demand is discretized over 4 periods representing the 4 seasons, so that each scenario 305 (nominal and disturbance) includes 4 periods. In addition, the study parameter considered is the 306 steam demands of one company including the demand for very high pressure (VHP), high pressure 307 (HP), medium pressure (MP) and low pressure (LP) steam of this company. The parameter evaluated 308 corresponds to the smallest maximum deviation among the demands of a company. The demands i 309 of the company considered are left variable and the smallest variation of the demands is maximized, 310 as described in equation (3). In fact, in this evaluation step, the demands of the other companies and 311 the facilities are fixed.

$$312 \quad \text{Flexibility}_{\text{Company } c} = \min_{i}[(\text{DevMax}_{i})] \qquad (13)$$

313 The flows circulating in the installations are variable over the time periods. The sizing of the system 314 consists in designing the facilities that will meet the demand of all the time periods. At each time 315 step, the planning of the system corresponds to the determination of the values of the operational 316 variables. In this model and for a given system, its sizing variable, Vdim, (not time-dependent) is 317 equal to the largest value of the operational variable, Vope(t), as presented in equation (4).

$$318 \quad Vdim \geq Vope(t) \qquad (14)$$

319 This relationship is only true if Vdim is minimized, achieved by minimizing the infrastructure cost 320 related to the sizing variable.
321

322

## 4.3. Case study

323

The case study chosen has been already studied in previous studies (Kim et al., 2010; Mousqué et al., 2020); it consists in an industrial park including 15 companies, the Yeosu park located in South Korea.

325

This large-scale park is composed of various companies of different processes operating over four seasons. Demands of processes represent four processes per company (60 demands) over four time periods (240 demands), they are presented in Appendix 1. The demands to satisfy are various and can be either electricity, low, medium, high, very high-pressure levels of steam and water. Steam is produced by boilers fed by natural gas in the different companies and electricity can be produced by renewable sources (wind turbines, solar PV panels) or can be taken to the grid. This case has been chosen because it constitutes a large-scale case study with a good representation of eco-industrial parks. Indeed, the corresponding mathematical model includes more than 7000 constraints, 7200 variables and among them, more than 1300 are binary variables. In terms of modelling and resolution, the large size of this case study and especially, the large number of binary variables make the resolution and the flexibility study particularly difficult to handle.

326

327

328

329

330

331

332

333

334

335

336

337

## 4.4. An integrated approach

338

For the sake of developing a systemic approach, the procedure called "Integrated procedure" is developed in a first intention (Figure 4). It consists in several steps able to find a flexible network: the first step finds the network with the cost minimum. This "initial network" constitutes a basis and the optimal network that takes into account flexibility is found starting from this point. In the integrated procedure, the structure of the network at the cost minimum is set and constraints based on demand variations are added to the design stage in order to find the optimal network. The final optimal solution must handle nominal demands and also all the deviation scenarios tested.

339

340

341

342

343

344

<figure><img src="image_5.png" /><figcaption>Figure 4. Integrated procedure</figcaption></figure>

345

346
In this procedure, the principle is to design a network to meet the nominal demand as well as all expected deviations for each demand at the first design stage. Indeed, an integrated approach can provide significant benefits depending on the optimization criteria considered, as the approach simultaneously considers a larger number of operating parameters as well as the relationships between these operating states (Andiappan, 2017). The network is optimized according to the design objectives. In addition, to achieve a flexible network, the facilities will have to be able to meet the nominal demand as well as the scenario demand. A scenario is constructed for each parameter (i.e. for each demand). For each scenario, one parameter is at its flexibility threshold value and the other parameters are at their nominal values. Thus, scenario after scenario and therefore parameter after parameter, the solution becomes flexible. This flexibility is obtained on the whole network, because once optimized it must be able to operate for all the scenarios with the same infrastructure. Since this approach integrates the nominal and fluctuating operating parameters, it is expected that the computation time will be important and even tremendous. The goal is to obtain overall better solutions because the model is optimized with the fluctuations that have already been taken into account.

The calculations over the case study presented before were performed on an Intel® Core™ i5-7440HQ processor running at 2.80 GHz. The optimization model is solved within ILOG environment with the help of CPLEX solver. The initial solution over four time periods is obtained in two minutes. The solution of the Integrated procedure, on the other hand, is performed over 64 time steps, i.e., over the scenario at nominal conditions and 15 disruption scenarios, i.e., 16 scenarios of 4 time steps. It requires 483 hours to obtain a solution able to overcome 10% of overall flexibility and for other flexibility thresholds (20 or 30%), no solution is found by the procedure even when the number maximal of iteration is enlarged.

### 4.5. Proposal of new procedures

Some declinations of the integrated procedure have been developed and tested in order to either decrease computational times or help the convergence of the optimization model including many binary variables. Two other procedures are studied: "Two-step procedure" and "Peak-demand procedure", the characteristics of these procedures are summarized in Figure 5. During this study, different procedures have been developed and studied. The procedures differ by the ability of the optimal solution to handle different operating conditions so that several scenario are tested according to the *Rl* presented in the previous subsection.
Peak-demand procedure

<figure><img src="image_6.png" /></figure>

a.

378

Two-steps procedure

<figure><img src="image_7.png" /></figure>

b.

379

<caption>Figure 5. Formulation of the procedures developed</caption>

380

(a. Peak-demand procedure and b. Two-steps procedure).

381

In order to avoid high computational times observed with the integrated procedure, the "Peak-demand procedure" has been also tested (Figure 5a). This procedure has the advantage of being an integrated method where all the fluctuations are taken into account at the same design stage. This approach consists in sizing the system to meet the nominal demand as well as the peak demand (i.e. maximum demand) that each process can achieve independently of the others. Although this is a common method of sizing networks, significant oversizing can be expected compared to previous methods, as the scenario where all demands are at their maximum level at the same time is a scenario that should not occur.

382

383

384

385

386

387

388

389

A third procedure, called "Two-steps procedure" has been developed in order to overcome the different drawbacks of the previous ones. According to this approach, during a first step, the network is designed according to the nominal demand. In the second step, this network is re-designed (retrofit) so that the nominal demands and all the demands of the deviation scenarios are integrated in this optimal design step. The principle is to start with an initial network, although not flexible, and

390

391

392

393

394
395

to re-design it in one-step by expanding its infrastructure and network. This situation can be
compared to the case of a retrofit design of an existing ecopark. It is expected that the computation
time should be reduced compared to the integrated procedure since in the re-design step a set of
variables is already established in the first step and becomes constraints in the second step. This
method is thus comparable to a compromise between the two previous methods.

396

397

398

399

400

401

4.6. Numerical resolution

402

With the same computational characteristics previously enounced, while for the integrated
procedure, it requires 483 hours, the integrated retrofit solutions are obtained in a range of 9 to 45
minutes. Indeed, the multiplication of the number of scenarios exponentially increases the number
of time-related variables. This time includes the 2 minutes to obtain the nominal solution over 4
periods plus the integrated retrofit time over 64 periods. The solution for the peak procedure is
obtained within 8 and 11 minutes over 8 time periods: the 4 nominal periods and the 4 periods of the
scenario with the peak demand.

403

404

405

406

407

408

409

410

5. Results and analyses

411

5.1. Analysis of the nominal solution

412

The nominal solution is the network configuration representing the minimum of the cost without considering any flexibility indicator. Figure 6 shows the flexibility of each company included in the park over the four periods. Let us recall that the flexibility announced represents the smallest maximal deviation of studied parameters. Companies 2, 9, 12 and 13 are the only companies able to be flexible to all the four periods considered. They are flexible over the different time periods even if the flexibility for the first period is very limited (less than 10%), whereas is very huge for the third period of company 12 (more than 300%). The other companies have a flexibility equal to zero at least at one period so their overall flexibility is null when the indicator RI is considered because the overall flexibility is the minimum value over the four time periods.

413

414

415

416

417

418

419

420
<figure><img src="image_8.png" /><figcaption>Figure 6. Flexibility of industries in the EIP over time periods</figcaption></figure>

On the first time period companies are globally less flexible. Indeed, in the model the four periods represent the four seasons of a year and the first one simulates demands in winter. In this period, the vapour demand is very high so the network is designed to provide resources according to the more constrained period. On the contrary, in the third period, more flexibility is observed because in the summer the demand in energy is lower so the network becomes more flexible. The network can handle and assume a large deviation of the demand in this period. However, some companies like 6 and 15 remain not flexible at all; it means that the designed network only handle the initial demand and cannot support any change in flows exchanged. When the cost is minimized, the equipment are sized to satisfy the constraints imposed by demands, not more. The capacity production is then found exactly to be able to answer the nominal demand for the most restricting companies. In companies where a little flexibility can be obtained, it is attributed to steam turbines: when a demand fluctuation occurs, it is possible to stop the turbines and the save steam directly supplies processes. The steam excess is then equal to the vapour lost with condensation process in steam turbines. Most particularly, the flexibility is driven by energy demands; for instance, when a company requires a great amount of very high pressure vapour, the system is constrained and the flexibility decreases. It is also very dependent of the equipment set (which is also a variable), if the equipment is able to provide larger amount of steam, the related pipes will be flexible but if the equipment already produces its maximal value, the flexibility will be very low. More generally, in a network, flexibility can be anticipated at the design stage by adding different equipment such as storage or controlled elements that can be stopped on demand.
445

## 5.2. Cost of flexibility

446

447

448

449

Following the different procedures, the flexibility of the solutions found have been analysed. Figure 7 shows the cost differences obtained by the different procedures developed. These cost differences are relative to the nominal solution and represent the additional total cost (CAPEX and OPEX) needed to reach the required flexibility.

<figure><img src="image_9.png" /></figure>

450

451

<caption>Figure 7. Additional overall cost in function of the flexibility required with the different procedures</caption>

452

453

454

455

456

457

As it was expected, adding flexibility to networks induces an increase of the investment cost for the different companies included in the EIP. The flexibility reported is the minimum change accepted on the overall network considering each company over each time-period. Without any change, the nominal solution is not flexible at all and does not adapt to any fluctuation of inlet parameters. For obtaining a flexibility minimum of 10%, the three procedures converge to three different optimal solutions:

458

459

460

461

462

463

464

465

- **_The peak-demand procedure_**: the solution obtained consists in a network oversized with more interconnections, higher boilers production and larger turbines compared to nominal solution. Even if it can accept 10% of changes in its inlet parameters, the cost has an increase of 20 million €, representing an increase of 4.7% of the nominal CAPEX. Then, when the minimum flexibility increases of 10%, the cost doubles. This observation remains case study-dependant and here, it can be explained by the data on costs equipment's. Indeed, when more flexibility is required, another equipment (turbine or boiler) is implanted, so that it doubles the additional cost.
- The two-steps procedure: the optimal solution found is close to the global optimal solution found with the integrated procedure. In terms of cost, this solution is interesting, as it represents an additional cost of 10 million € accounting for an increase of 2.2% of the nominal CAPEX.

- Integrated procedure: The integrated procedure offers the advantage to find a global optimum to the optimization problem. The main drawback of this methodology is the computation time and the ability of the algorithm to converge to the global optimum. For a flexibility of 10%, an optimal point has been found in 483 hours of computation time with a result very close to the two-steps procedure developed in this study (Table 1).

<table><caption>Table 1. Computational times for obtaining optimal solutions</caption><thead><tr><th>Flexibility</th><th>Procedure</th><th>CPU time</th></tr></thead><tbody><tr><td>0%</td><td>Nominal solution</td><td>2 min</td></tr><tr><td rowspan="3">10%</td><td>Peak-demand</td><td>9 min</td></tr><tr><td>Two-steps</td><td>9 min</td></tr><tr><td>Integrated</td><td>483 hr</td></tr><tr><td rowspan="3">20%</td><td>Peak-demand</td><td>11 min</td></tr><tr><td>Two-steps</td><td>45 min</td></tr><tr><td>Integrated</td><td>-</td></tr><tr><td rowspan="3">30%</td><td>Peak-demand</td><td>8 min</td></tr><tr><td>Two-steps</td><td>45 min</td></tr><tr><td>Integrated</td><td>-</td></tr></tbody></table>

Table 1 shows the different CPU times for the solutions obtained in order to compare the procedures developed. Several trends are important to notice: for a limited flexibility (10%), the constraints over binary variables are not too important and the two-procedures developed, (peak-demand and two-steps) are very close in terms of computational times. However, when flexibility imposed increases, the time for the two-step procedure also increases mainly due to important constraints on binary variables (existence of connections and of equipment's). Regarding the integrated procedure, the time is very long but a solution has finally been found. For other flexibility thresholds (20 and 30%) the solver does not reach to find a global optimum.

### 5.3. Interconnections of the companies in the EIP

Another important view to analyse solutions is to understand the link between interconnections in the Eco-industrial Park and its flexibility (Figure 8).
<figure><img src="image_10.png" /><figcaption>Figure 8. Cost related to additional connections for each optimal solutions</figcaption></figure>

The nominal solution involves a network with 15 interconnections between the different companies of the park. The peak-demand procedure finds the most expensive solution but the additional charge is not attributed to connections and is mainly due to the capacity production of turbines and boilers.

For the two-step procedure, interconnections are added to provide flexibility over time periods and are generated a cost higher than in the peak-demand solution. Indeed, by increasing the possibilities to exchange flows between companies, it is easier to adapt to fluctuations of demands.

<figure><img src="image_11.png" /><figcaption>Figure 9: Network designed with 20% of flexibility</figcaption></figure>
Another observation to notice is related to the interdependency of companies. hen many interconnections are involved between companies, it can be more difficult to handle the leave of one company and it increases the dependency between the participants of the park (Mousqué et al., 2020). The final solution with 20% of flexibility for the network studied is represented on Figure 9. If we consider other solutions obtained by the other procedures the network obtained are different even if the criteria are close. The changes can be inside a particular company or related to a connection between two companies, it is variable.

## 5.4. Discussion

Among the procedures developed in this study (Table 2), although the integrated procedure should be the most suitable to find global optimal flexible solutions, the CPU time and the capacity to converge towards a solution constitute serious bottlenecks. The solutions obtained by the peak-demand procedure is interesting in terms of CPU time but it allows to design larger capacity network, able to adapt to fluctuations of demands with a raise up of the overall cost of the network. Then, we suggest that it can be used in a first step to approximate some cost calculations. A good compromise remains the two-steps procedure developed in this study because it allows obtaining a flexible network with a cost near from the solution with the integrated procedure.

Although longer than the peak-demand procedure, the CPU time for the two-step procedure remains largely reasonable, and the most important, the algorithm converges towards a solution in all the cases tested.

<table><caption>Table 2. Summary of the main advantages/drawbacks of the procedures developed</caption><thead><tr><th>Procedure</th><th>CPU Time</th><th>Capacity of convergence</th><th>Flexibility of the designed network</th><th>Cost of the designed network</th><th>Interdependency of companies</th></tr></thead><tbody><tr><td>Classical monobjective optimization</td><td>+++</td><td>++</td><td>--</td><td>++</td><td>+</td></tr><tr><td>Peak demand procedure</td><td>++</td><td>++</td><td>-</td><td>-</td><td>++</td></tr><tr><td>Two-step procedure</td><td>+</td><td>++</td><td>++</td><td>+</td><td>-</td></tr><tr><td>Integrated procedure</td><td>-</td><td>--</td><td>++</td><td>+</td><td>-</td></tr></tbody></table>
517
518

Even if the different procedures have several advantages and drawbacks, the two-step procedure remains a good compromise to obtain economic and flexible networks in reasonable CPU time.

519
520
521
522
523
524
525
526

Another interesting approach could be to use it in order to have recommendation of which company should be included in a park without taking an important risk. At a tactical scale, when companies are already set up in an industrial zone, this approach could also serve to determine the limits of the EIP to be constructed so that it will be able to overcome the fluctuations in a changing world. Referring to the case study of this paper, some companies such as n°6 or n°15 can threaten the whole flexibility of the EIP if they become participants of the EIP. To avoid this kind of situation, a particular attention must be paid to these companies so that they can hold a storage equipment or gas turbines to overcome fluctuations.

527
528

Finally, this kind of numerical optimization-based approach can also give some positive arguments to promote the development of industrial symbioses at a territorial scale.

530

## 6. Conclusions

531
532
533
534
535
536

Eco-industrial parks are now well-known to be a solution to transform our classical linear industrial zones into circular areas where production level and environmental preservation can go towards the same direction. However, their development in local territories remains limited mainly due to the existence of bottlenecks such as interdependencies or flexibility of the companies included in the network. In this study, a generic methodology has been developed for the optimal design of flexible networks in eco-industrial parks.

537
538
539
540
541
542
543
544
545
546
547
548

To evaluate the flexibility of the energies network designed, the Resilience Indicator RI (RI) has been chosen. By evaluating independently the deviations of the parameters, the procedure can be applied to complex and large-scale problems with a linear formulation for the optimization model. In the approach, both the inlet parameters and demands can vary and the chosen indicator doesn't need to anticipate any occurrence probability of the fluctuations to come. The RI represents the maximal deviation from the nominal conditions that a network can handle without changing its configuration, whereas the other companies carry on to run at their nominal conditions. The developed procedures are based on a flexibility threshold fixed in an optimization-based approach, which can be fixed either by experts or by decision-makers. In this study, it is observed that flexibility is expensive with the different procedures: if the flexibility threshold is high, the network will be adaptable to fluctuations but its overall cost will be higher than the nominal network around 5% (more than the nominal CAPEX) for a flexibility of 10%.
On a numerical point of view, a two-step procedure is presented as a good solution to handle complex and large-scale problems by facilitating the convergence to a feasible solution. The decomposition of the method in two steps allows looking for the nominal network and starting from its structure, flexibility is evaluated over different thresholds for different scenarios. In this procedure, the network is designed so that production must satisfy the demand of the nominal situation as well as the situation of several scenarios of deviation.

Future works should analyse the link between dependencies of the companies in the EIP and the flexibility of the overall network. Another important focus for future research could be the importance of the level of decision (tactical, strategic) on the design of flexible networks. Indeed, even if the procedure developed in this study can bring solutions at a tactical scale, it could be interesting to go further, at the operational scale, to understand the balance of the solutions proposed between the different scales.

## References

Aguilar, O., Kim, J.K., Perry, S., Smith, R., 2008. Availability and reliability considerations in the design and optimisation of flexible utility systems. Chem. Eng. Sci. 63, 3569–3584. https://doi.org/10.1016/j.ces.2008.04.010

Alkahtani, M., Omair, M., Khalid, Q.S., Hussain, G., Ahmad, I., Pruncu, C., 2021. A COVID-19 Supply Chain Management Strategy Based on Variable Production under Uncertain Environment Conditions. Internat. J. Environ. Res. Publ. Health, 18 (4) https://doi.org/10.3390/ijerph18041662

Andiappan, V., 2017. State-Of-The-Art Review of Mathematical Optimisation Approaches for Synthesis of Energy Systems. Process Integr. Optim. Sustain. 1, 165–188. https://doi.org/10.1007/s41660-017-0013-2

Bandyopadhyay, S., 2011. Design and optimization of isolated energy systems through pinch analysis, in: Asia-Pacific Journal of Chemical Engineering. pp. 518–526. https://doi.org/10.1002/apj.551

Di Pretoro, A., Montastruc, L., Manenti, F., Joulia, X., 2019. Flexibility analysis of a distillation column: Indexes comparison and economic assessment. Comput. Chem. Eng. 124, 93–108. https://doi.org/10.1016/j.compchemeng.2019.02.004

Grossmann, I.E., Morari, M., 1983. Operability, Resiliency, and Flexibility: process design objectives for a changing world. https://doi.org/10.1184/R1/6467234.V1

Hartley, K., van Santen, R., Kirchherr, J., 2020. Policies for transitioning towards a circular economy: Expectations from the European Union (EU). Resour. Conserv. Recycl., 155.

Hoch, P.M., Eliceche, A.M., 1996. Flexibility analysis leads to a sizing strategy in distillation columns. Comput. Chem. Eng. https://doi.org/10.1016/0098-1354(96)00034-8
583 Kim, S.H., Yoon, S.G., Chae, S.H., Park, S., 2010. Economic and environmental optimization of a multi-site utility network for an industrial complex. J. Environ. Manage. 91, 690-705.
584 https://doi.org/10.1016/j.jenvman.2009.09.033
585

586 Kristensen, H.S., Mosgaard, M.A., 2020. A review of micro level indicators for a circular economy - moving away from the three dimensions of sustainability? J. Clean. Prod., 243.
587

588 Meerow, S., Newell, J.P., 2015. Resilience and Complexity: A Bibliometric Review and Prospects for Industrial Ecology. J. Ind. Ecol. 19, 236-251. https://doi.org/10.1111/jiec.12252
589

590 Montastruc, L., Boix, M., Pibouleau, L., Azzaro-Pantel, C., Domenech, S., 2013. On the flexibility of an eco-industrial park (EIP) for managing industrial water. J. Clean. Prod. 43, 1-11.
591 https://doi.org/10.1016/j.jclepro.2012.12.039
592

593 Mousqué, F., Boix, M., Montastruc, L., Domenech, S., Negny, S., 2020. Optimal design of eco-industrial parks with coupled energy networks addressing complexity bottleneck through an interdependance analysis. Comp. Chem Eng, 138, 106859.
594
595

596 Pistikopoulos, E.N., Mazzuchi, T.A., 1990. A novel flexibility analysis approach for processes with stochastic parameters. Comp. Chem. Eng. https://doi.org/10.1016/0098-1354(90)87055-T
597

598 Pistikopoulos, E.N., 1995. Uncertainty in process design and operations. Comp. Chem. Eng., 19, 553-563.
599

600 Saboo, A.K., Morari, M., Woodcock, D.C., 1985. Design of resilient processing plants-VIII. A resilience index for heat exchanger networks. Chem. Eng. Sci. https://doi.org/10.1016/0009-2509(85)80097-X
601
602

603 Schneider, D.F, 1997. Debottlenecking Options and Optimization. Report, Stratus Engineering.

604 Sharifzadeh, M., 2013. Integration of process design and control: A review. Chem. Eng. Res. Des. 9, 2515-2549. https://doi.org/10.1016/j.cherd.2013.05.007
605

606 Swaney, R.E., Grossmann, I.E., 1985. An index for operational flexibility in chemical process design. Part I: Formulation and theory. AIChE J. 31, 621-630. https://doi.org/10.1002/aic.690310412 1
607

608 Valenzuelas-Venegas, G., Henriquez-Henriquez, F., Boix, M., Montastruc, L., Arenas-Araya, F., Miranda-Perez, J., Diaz-Alvarado, F.A., 2018. A resilience indicator for Eco-industrial Parks. J. Clean. Prod. 174, 807-820. https://doi.org/10.1016/j.jclepro.2017.11.025
609
610

611 Valenzuela-Venegas, G., Vera-Hofmann, G., 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. https://doi.org/10.1016/j.jclepro.2019.118610
612
613

614 Xu, X.D., Sun, W.Q., Abeysekera, M., Qadrdan, M., 2021. Quantifying the flexibility from industrial steam systems for supporting the power grid. IEEE Transac. Pow. Systems, 36 (1), pp. 313-322.
615 https://doi.org/10.1109/TPWRS.2020.3007720
616

617
618
619 Appendix 1 – Mathematical formulation of the optimization model

620 Nomenclature

621 Sets

622 b: boiler

623 f: fuel

624 tu: turbine

625 so: photovoltaic (PV) solar panel

626 wi: wind turbine

627 sr: electrical sources (so; wi ∈ sr)

628 sr<sub>AC</sub>: Alternative Current (AC) electrical power sources (wi ∈ sr<sub>AC</sub>)

629 sr<sub>DC</sub>: Direct Current (DC) electrical power sources (so ∈ sr<sub>DC</sub>)

630 C: company (c ∈ C, c' ∈ C, c ≠ c')

631 H: steam header (h ∈ H, h' ∈ H, h < h')

632 T: time period (t ∈ T, tm ∈ T, ∀ tm ≤ t)

633

634 Parameters

635 M: large positive value

636 Nbtime: number of step time

637

638 i: discount rate

639 n: project lifespan (years)

640 SL<sup>h</sup>: percentage of steam losses by trap at steam header

641 WL<sup>wu</sup>: percentage of water losses at water unit

642

643

644 Cal<sub>f</sub>: fuel calorific power (kWh/t)

645 Eff<sub>b,f</sub>: boiler efficiency for each fuel

646 Eff<sub>tu</sub>: turbine efficiency

647 ConP: condensate parameter for turbine

648 h<sub>h</sub>: steam enthalpy depending on pressure level (kWh/t)

649 h<sub>w</sub>: water enthalpy before boiler heating (kWh/t)

650

651 D<sup>s</sup><sub>h,t,c</sub>: steam demand (ton/h)

652 D<sup>e</sup><sub>tc</sub>: electricity demand (kWh)

653

654 ECons<sup>fix</sup><sub>b</sub>: fixed electricity consumption for boiler (kWh)

655 ECons<sup>var</sup><sub>b</sub>: variable electricity consumption for boiler (kWh)

656

657 P<sup>Max</sup><sub>b,h</sub>: maximum steam production at boiler b (ton/h)

658 P<sup>min</sup><sub>b,h</sub>: minimum steam production (ton/h)

659 S<sup>Max</sup><sub>h'tu</sub>: maximum steam in turbine tu (ton/h)

660 S<sup>min</sup><sub>h'tu</sub>: maximum steam in turbine tu (ton/h)

661

662 LoadFactor<sub>so,t</sub>: load factor for the solar PV production

663 LoadFactor<sub>wit</sub>: load factor for the wind turbines production

664

665 v<sub>wit</sub>: wind speed in the time interval t (m/s)

666 v<sup>rated</sup><sub>wi</sub>: rated wind speed (m/s)

667

668 DFA: Discount factor for allocation of annual operational expenditure and resources over the life of the project.

670

671 Pri<sup>fuel</sup><sub>f</sub>: fuel purchase price (USD/ton)

672 Pri<sup>water</sup>: water purchase price (USD/ton)

673 PriElec<sup>purchase</sup>: electricity purchase price (USD/kWh)

674 PriElec<sup>sales</sup>: electricity sales price (USD/kWh)
675

676 F<sub>b,f,t,c</sub><sup>stock_safety</sup>: fuel safety stock (ton)

677

678

679 CAPEX<sub>b</sub>: capital expenditure for boiler (USD)

680 OPEX<sub>b</sub>: operational expenditure for boiler (USD per year)

681 DFI<sub>b</sub>: discount factor for investment and replacement cost of boiler

682

683 CAPEX<sub>tu</sub>: capital expenditure for boiler (USD)

684 OPEX<sub>tu</sub>: operational expenditure for turbine (USD per year)

685 DFI<sub>tu</sub>: discount factor for investment and replacement cost for turbine

686

687 CAPEX<sub>pipe</sub><sup>fix</sup>: fixed capital expenditure for pipeline (USD)

688 CAPEX<sub>pipe</sub><sup>var</sup>: variable capital expenditure for pipeline (USD/ton/h)

689 OPEX<sub>pipe</sub>: operational expenditure for pipeline (USD/ton/h))

690 DFI<sub>pipe</sub>: discount factor for investment and replacement cost of boiler

691

692 minInterc<sub>Threshold</sub><sup>Diameter</sup>: minimum diameter of interconnection pipes (ton/h)

693 MaxInterc<sub>Threshold</sub><sup>Diameter</sup>: Maximum diameter of interconnection pipes (ton/h)

694

695 CAPEX<sub>sr</sub>: capital expenditure for electric power sources (USD/kW)

696 OPEX<sub>sr</sub>: operational expenditure for electric power sources (USD/kW per year)

697 DFI<sub>sr</sub>: discount factor for investment and replacement cost of electric power sources

698 Effinv: inverter efficiency for DC sources

699

700 MaxEprodTur: maximum production for turbines (kWh)

701 MaxInstall<sub>sr</sub>: maximum power installed for electrical source (kW)

702

703 Variables

704 NPV: Net Present Value of the overall project

705 Cost<sup>Raw</sup>: raw material cost for the project lifespan (USD)

706 Cost<sup>Boilers</sup>: boilers cost for the project lifespan (USD)

707 Cost<sup>Pipes</sup>: pipes cost for the project lifespan (USD)

708 Cost<sup>Tubines</sup>: turbines cost for the project lifespan (USD)

709 Cost<sup>Elec.sources</sup>: electrical sources cost for the project lifespan (USD)

710 Cost<sup>Elec_cons</sup>: electricity cost for the project lifespan (USD)

711 Cost<sup>Fuel</sup>: fuel cost for the project lifespan (USD)

712 Cost<sup>Water</sup>: water cost for the project lifespan (USD)

713

714 Nb<sub>interc</sub>: number of interconnection in the utility system

715

716 T<sub>h,t,c,c'</sub><sup>s</sup>: steam transfer from company c to company c' (ton/h)

717 T<sub>t,c,c</sub><sup>w</sup>: water transfer from company c to company c' (ton/h)

718

719 y<sub>b,c</sub><sup>sel</sup>: binary variable to select boiler technologies

720 y<sub>b,t,c</sub><sup>prod</sup>: binary variable to regulate boiler steam production

721 y<sub>tu,c</sub><sup>sel</sup>: binary variable to select turbine technologies

722 y<sub>h,c,c'</sub><sup>pipe</sup>: binary variable to select interconnection pipes from company c to company c'

723

724 P<sub>sr</sub><sup>rated</sup>: rated installed power of electrical power source (kW) (so; wi ∈ sr)

725

726 Interc<sub>h,c,c'</sub><sup>Diameter</sup>: diameter of an interconnection pipe (ton/h)

727

728 F<sub>b,f,t,c</sub>: fuel consumed by boiler (ton)

729
730 $F_{b,f,t,c}^{\text{purchase}}$: fuel purchase by boiler (ton)

731 $F_{b,f,t,c}^{\text{stock}}$: fuel stock by boiler (ton)

732

733 $S_{b,h,t,c}^{\text{prod}}$: steam production from boiler to steam header (ton/h)

734 $S_{h',h,t,c}$: steam flow from header h' to steam header h (ton/h)

735 $S_{h,t,u,t,c}^{\text{In}}$: steam flow from steam header to turbine (ton/h)

736 $S_{h,t,u,t,c}^{\text{Out}}$: steam flow from turbine to steam header (ton/h)

737 $S_{h,t,c}^{\text{trap}}$: steam losses by condensate trap at steam header (ton/h)

738 $S_{h,t,c}^{\text{vent}}$: steam vented at steam header (ton/h)

739

740 $W_{h,t,c}$: water flow from process use to water unit (ton/h)

741 $W_{b,t,c}$: water flow from water unit to boiler (ton/h)

742 $W_{tu,t,c}^{\text{Out}}$: water flow from turbine to water unit (ton/h)

743 $W_{t,c}^{\text{stock}}$: quantity of water in the utility network during time interval t (ton)

744 $W_{t,c}^{\text{purchase}}$: water purchase for water unit feeding (ton)

745 $W_{t,c}^{\text{feeding}}$: water purchase for water unit feeding (ton)

746 $W_{t,c}^{\text{losses}}$: Water losses in water feed boiler (ton)

747 $W_{t,c}^{\text{dump}}$: water dumped from the utility network (ton)

748

749

750 $E_t^{\text{prodSources}}$: produced electric power in the HPS (kWh)

751 $E_{sr,t}^{\text{prod}}$: produced electric power depending on the electrical source (kWh) (so; $w_i \in sr$)

752 $E_t^{\text{prodCons}}$: produced electric power consumed (kWh)

753 $E_t^{\text{prodSell}}$: electric power sold per step time (USD)

754 $E_{t,c}^{\text{purchase}}$: electric power purchased per step time by a company (USD)

755 $E_{b,t,c}^{\text{boiler}}$: consumed electric power by boiler (kWh)

756 $E_{tu,t,c}^{\text{turbine}}$: electric power produced by turbine (kWh)

757 $E_{sr,t}^{\text{sources}}$: electric power produced by sources in the HPS (kWh)

758

### 759 Mass and energy balances for equipment, based on the superstructure of Figure 2:

760 **Steam headers:** A steam header can take steam produced by the boilers, released by turbines or by the pressure release valve. Then it allows feeding other turbines and the demands of the processes within or between companies.

$$ \sum_b S_{b,h,t,c}^{\text{prod}} + \sum_{h'} S_{h',h,t,c} + \sum_{tu} S_{h,t,u,t,c}^{\text{Out}} + \sum_{c'} T_{h,t,c',c}^{s} = D_{h,t,c}^{s} + \sum_{h'} S_{h,h',t,c} + S_{h,t,c}^{\text{trap}} + S_{h,t,c}^{\text{vent}} + \sum_{tu} S_{h,t,u,t,c}^{\text{In}} + \sum_{c'} T_{h,t,c,c'}^{s} \quad (1) $$

$$ (\forall c' \neq c, h < h') $$

763

$$ \text{With, } S_{h,t,c}^{\text{trap}} = SL^h \times \left( \sum_b S_{b,h,t,c}^{\text{prod}} + \sum_{h'} S_{h',h,t,c} + \sum_{tu} S_{h,t,u,t,c}^{\text{Out}} + \sum_{c'} T_{h,t,c',c}^{s} \right) \quad (2) $$

### 764 Water unit:

765 The water flowing into the water unit can come from other companies, turbines, processes, and from water feeding to compensate water losses.
$$ \Sigma_{c'} T_{t,c',c}^{w} + \Sigma_{tu} W_{tu,t,c}^{\text{Out}} + \Sigma_{h} W_{h,t,c} + W_{t,c}^{\text{feeding}} = \Sigma_{b} W_{b,t,c} + W_{t,c}^{\text{losses}} + \Sigma_{c} T_{t,c,c'}^{w} \quad (3) $$

767 Water losses are a part of the water flowing into the water unit (16).

$$ W_{t,c}^{\text{losses}} = WL^{wu} \times (\Sigma_{tu} W_{tu,t,c}^{\text{Out}} + \Sigma_{h} W_{h,t,c}) \quad (4) $$

768

$$ W_{t,c}^{\text{feeding}} = W_{t,c}^{\text{losses}} + \sum_{h} (S_{h,t,c}^{\text{trap}} + S_{h,t,c}^{\text{vent}}) \quad (5) $$

769 The quantity of water in the system $W_{t,c}^{\text{stock}}$ corresponds to the quantity consumed by the boilers during a time step t (18).

770

$$ W_{t,c}^{\text{stock}} = \sum_{b} W_{b,t,c} \quad (6) $$

771 **Multi-period constraints:**

772 Water stock:

$$ W_{t,c}^{\text{stock}} = \sum_{t=0}^{\text{tm}} (W_{t,c}^{\text{purchase}} - W_{t,c}^{\text{dump}}) \quad (\forall \text{ tm} \le t) \quad (7) $$

773

774 Fuel stock:

775

$$ F_{b,f,t,c}^{\text{stock}} = F_{b,f,t-1,c}^{\text{stock}} + F_{b,f,t,c}^{\text{purchase}} - F_{b,f,t,c} \quad (\forall t>0) \quad (8) $$

776

$$ F_{b,f,t,c}^{\text{stock}} \ge F_{b,f,t,c}^{\text{stock\_safety}} \quad (9) $$

777

$$ F_{b,f,t=0,c}^{\text{stock}} = F_{b,f,t=N\text{btime},c}^{\text{stock}} \quad (10) $$

778 Turbines:

779 Inlet steam is constrained by a minimum and a maximum value, if it exists:

$$ S_{h',tu}^{\min} \times y_{tu,c}^{\text{sel}} \le S_{h',tu,t,c}^{\text{In}} \le S_{h',tu}^{\text{Max}} \times y_{tu,c}^{\text{sel}} \quad (11) $$

780

781 An upper limit threshold to the amount of steam leaving the

$$ S_{h',tu}^{\min} \times y_{tu,c}^{\text{sel}} \le S_{h',tu,t,c}^{\text{Out}} \le S_{h',tu}^{\text{Max}} \times y_{tu,c}^{\text{sel}} \quad (12) $$

782

783 In the turbines, the steam entering through higher pressure inlet h' is distributed at the outlets

784 between the different lower pressure levels h and the condensed water (25).
$$S_{h',tu,t,c}^{\text{In}} = \sum_{h} S_{h,tu,t,c}^{\text{Out}} + W_{tu,t,c}^{\text{Out}} \quad (h < h') \qquad (13)$$

785

786 The electricity produced depends on the efficiency of the turbines but also on the pressure of the
787 steam at the inlet and its pressure at the outlets (26) (Agha et al., 2010).

$$E_{tu,t,c}^{\text{turbine}} = \text{Eff}_{tu} \times \left( S_{h',tu,t,c}^{\text{In}} \times h_{h'} - \sum_{h} S_{h,tu,t,c}^{\text{Out}} \times h_{h} - W_{tu,t,c}^{\text{Out}} \times h_{w} \right) \quad (h < h') \qquad (14)$$

788 A proportion corresponding to ConP of the steam entering the turbine is condensed (27).
$$W_{tu,t,c}^{\text{Out}} = \text{ConP} \times S_{h',tu,t,c}^{\text{In}} \qquad (15)$$

789 **Boilers:**

790 The steam production of the boilers is between its minimum power production $P_{b,h}^{\min}$ and its
791 maximum capacity $P_{b,h}^{\max}$ (28).

$$P_{b,h}^{\min} \times y_{b,t,c}^{\text{prod}} \leq S_{b,h,t,c}^{\text{prod}} \leq P_{b,h}^{\max} \times y_{b,t,c}^{\text{prod}} \qquad (16)$$

792

793 Steam production depends on the enthalpy difference between the water at the inlet and the steam
794 at the outlet, on the quantity of fuel consumed, on the calorific value of the fuel, and also on the
795 efficiency of the boiler.

$$S_{b,h,t,c}^{\text{prod}} = \frac{\sum_f \text{Cal}_f \times F_{b,f,t,c} \times \text{Eff}_{b,f}}{(h_h - h_w)} \qquad (17)$$

796 The amount of water entering the boilers corresponds to the amount of steam produced (30).

$$S_{b,h,t,c}^{\text{prod}} = W_{b,t,c} \qquad (18)$$

797 The electricity consumed by the boilers is described by equation (31).

$$E_{b,t,c}^{\text{boiler}} = E\text{Cons}_{b}^{\text{fix}} \times y_{b,c}^{\text{sel}} + E\text{Cons}_{b}^{\text{var}} \times S_{b,h,t,c}^{\text{prod}} \qquad (19)$$

798 If a boiler is started at least once then it is selected, that means purchased (32).

$$y_{b,c}^{\text{sel}} \geq y_{b,t,c}^{\text{prod}} \qquad (20)$$

799

800 **Hybrid Power System:**

801 The electricity produced by the energy sources in the HPS is equal to the electricity produced by the
802 turbines, by the AC (i.e. wind turbines) and DC (i.e. solar PV panels) sources multiplied by the
803 efficiency of the converter:

$$E_{t}^{\text{prodSources}} = \sum_{sr \in sr_{AC}} E_{sr,t}^{\text{prod}} + \sum_{tu,c} E_{tu,t,c}^{\text{turbine}} + \sum_{sr \in sr_{DC}} E_{sr,t}^{\text{prod}} \times \text{Eff}_{inv} \qquad (21)$$

804

$$E_{t}^{\text{prodSources}} = E_{t}^{\text{prodCons}} + E_{t}^{\text{prodSell}} \qquad (22)$$
805

$$ \sum_{c} D_{t,c}^{e} = E_{t}^{prodCons} + E_{t}^{purchase} \qquad (23) $$

806

The production of solar PV panels and wind turbines depends on the installed capacity and the load factor during the time period (36).

807

$$ E_{so,t}^{prod} = P_{so}^{rated} \times LoadFactor_{so,t} \qquad (24) $$

808

$$ E_{wi,t}^{prod} = P_{wi}^{rated} \times LoadFactor_{wi,t} \qquad (25) $$

809

The maximum power that can be installed for an electrical source is limited; this constraint depends on the resource available on the site (38).

810

$$ P_{sr}^{rated} \leq MaxInstall_{sr} \qquad (26) $$

811

The electric power production of the turbines is also limited depending on the possibility of installation on site (39).

812

$$ E_{tu,n,t,c}^{turbine} \leq MaxEproTur \qquad (27) $$

813

**Interconnection pipes:**

814

An exchange of steam between industries is only possible if the interconnection exists ($y_{h,c,c'}^{pipe} = 1$)

815

(40). M is a big positive value used to satisfy the existence constraint of a pipe.

$$ y_{h,c,c'}^{pipe}/M \leq T_{h,t,c,c'}^{s} \leq y_{h,c,c'}^{pipe} \times M \quad (\forall c' \neq c) \qquad (28) $$

816

817

If an interconnection exists, its diameter must be larger than the threshold, that represents the minimum technically feasible section of pipes (41).

818

$$ \text{Interc}_{h,c,c'}^{\text{Diameter}} \geq y_{h,c,c'}^{\text{pipe}} \times \text{minInterc}_{\text{Threshold}}^{\text{Diameter}} \quad (\forall c' \neq c) \qquad (29) $$

819

820
<table><caption>822 Steam (ton/hr) and electricity (kWh) demands for the 15 companies of Yeousou industrial park</caption><thead><tr><th>Period</th><th>Company</th><th>VHP</th><th>HP</th><th>MP</th><th>LP</th><th>Elec</th><th>Period</th><th>Company</th><th>VHP</th><th>HP</th><th>MP</th><th>LP</th><th>Elec</th></tr></thead><tbody><tr><td rowspan="15">t1</td><td>1</td><td>0</td><td>0</td><td>32</td><td>0</td><td>45</td><td rowspan="15">t3</td><td>1</td><td>0</td><td>0</td><td>23</td><td>0</td><td>45</td></tr><tr><td>2</td><td>0</td><td>0</td><td>84</td><td>127</td><td>91</td><td>2</td><td>0</td><td>0</td><td>18</td><td>36</td><td>73</td></tr><tr><td>3</td><td>0</td><td>0</td><td>54</td><td>45</td><td>73</td><td>3</td><td>0</td><td>0</td><td>50</td><td>41</td><td>73</td></tr><tr><td>4</td><td>0</td><td>0</td><td>13</td><td>18</td><td>45</td><td>4</td><td>0</td><td>0</td><td>9</td><td>16</td><td>45</td></tr><tr><td>5</td><td>0</td><td>0</td><td>0</td><td>5</td><td>18</td><td>5</td><td>0</td><td>0</td><td>0</td><td>5</td><td>18</td></tr><tr><td>6</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>6</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>7</td><td>0</td><td>0</td><td>0</td><td>23</td><td>45</td><td>7</td><td>0</td><td>0</td><td>0</td><td>18</td><td>45</td></tr><tr><td>8</td><td>0</td><td>59</td><td>42</td><td>6</td><td>91</td><td>8</td><td>0</td><td>24</td><td>21</td><td>0</td><td>91</td></tr><tr><td>9</td><td>0</td><td>0</td><td>153</td><td>87</td><td>63</td><td>9</td><td>0</td><td>0</td><td>61</td><td>0</td><td>63</td></tr><tr><td>10</td><td>408</td><td>408</td><td>0</td><td>0</td><td>109</td><td>10</td><td>393</td><td>145</td><td>0</td><td>0</td><td>109</td></tr><tr><td>11</td><td>408</td><td>494</td><td>156</td><td>105</td><td>109</td><td>11</td><td>377</td><td>313</td><td>156</td><td>105</td><td>109</td></tr><tr><td>12</td><td>0</td><td>0</td><td>28</td><td>166</td><td>45</td><td>12</td><td>0</td><td>0</td><td>9</td><td>50</td><td>45</td></tr><tr><td>13</td><td>0</td><td>105</td><td>31</td><td>13</td><td>63</td><td>13</td><td>0</td><td>100</td><td>23</td><td>13</td><td>63</td></tr><tr><td>14</td><td>553</td><td>203</td><td>191</td><td>104</td><td>181</td><td>14</td><td>532</td><td>189</td><td>190</td><td>87</td><td>163</td></tr><tr><td>15</td><td>382</td><td>67</td><td>60</td><td>64</td><td>154</td><td>15</td><td>376</td><td>57</td><td>60</td><td>54</td><td>136</td></tr><tr><td rowspan="15">t2</td><td>1</td><td>0</td><td>0</td><td>27</td><td>0</td><td>45</td><td rowspan="15">t4</td><td>1</td><td>0</td><td>0</td><td>26</td><td>0</td><td>45</td></tr><tr><td>2</td><td>0</td><td>0</td><td>45</td><td>73</td><td>91</td><td>2</td><td>0</td><td>0</td><td>45</td><td>73</td><td>73</td></tr><tr><td>3</td><td>0</td><td>0</td><td>54</td><td>45</td><td>73</td><td>3</td><td>0</td><td>0</td><td>54</td><td>45</td><td>73</td></tr><tr><td>4</td><td>0</td><td>0</td><td>13</td><td>18</td><td>45</td><td>4</td><td>0</td><td>0</td><td>11</td><td>18</td><td>45</td></tr><tr><td>5</td><td>0</td><td>0</td><td>0</td><td>5</td><td>18</td><td>5</td><td>0</td><td>0</td><td>0</td><td>5</td><td>18</td></tr><tr><td>6</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>6</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>7</td><td>0</td><td>0</td><td>0</td><td>23</td><td>45</td><td>7</td><td>0</td><td>0</td><td>0</td><td>23</td><td>45</td></tr><tr><td>8</td><td>0</td><td>45</td><td>32</td><td>0</td><td>91</td><td>8</td><td>0</td><td>36</td><td>32</td><td>0</td><td>91</td></tr><tr><td>9</td><td>0</td><td>0</td><td>100</td><td>63</td><td>63</td><td>9</td><td>0</td><td>0</td><td>100</td><td>63</td><td>63</td></tr><tr><td>10</td><td>399</td><td>272</td><td>0</td><td>0</td><td>109</td><td>10</td><td>399</td><td>181</td><td>0</td><td>0</td><td>109</td></tr><tr><td>11</td><td>390</td><td>363</td><td>156</td><td>105</td><td>109</td><td>11</td><td>390</td><td>363</td><td>156</td><td>105</td><td>109</td></tr><tr><td>12</td><td>0</td><td>0</td><td>18</td><td>109</td><td>45</td><td>12</td><td>0</td><td>0</td><td>18</td><td>109</td><td>45</td></tr><tr><td>13</td><td>0</td><td>105</td><td>27</td><td>13</td><td>63</td><td>13</td><td>0</td><td>105</td><td>27</td><td>13</td><td>63</td></tr><tr><td>14</td><td>544</td><td>190</td><td>191</td><td>100</td><td>181</td><td>14</td><td>535</td><td>190</td><td>191</td><td>100</td><td>163</td></tr><tr><td>15</td><td>382</td><td>63</td><td>60</td><td>59</td><td>154</td><td>15</td><td>379</td><td>63</td><td>60</td><td>59</td><td>136</td></tr></tbody></table>