Next Article in Journal
Numerical Study on the Influence of Step Casing on Cavitating Flows and Instabilities in Inducers with Equal and Varying Pitches
Next Article in Special Issue
Computational Optimization of Porous Structures for Electrochemical Processes
Previous Article in Journal
Large Eddy Simulations of Reactive Mixing in Jet Reactors of Varied Geometry and Size
Previous Article in Special Issue
Highly-Efficient Caffeine Recovery from Green Coffee Beans under Ultrasound-Assisted SC–CO2 Extraction
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

MILP Formulation for Solving and Initializing MINLP Problems Applied to Retrofit and Synthesis of Hydrogen Networks

by
Patrícia R. da Silva
1,*,
Marcelo E. Aragão
2,
Jorge O. Trierweiler
1 and
Luciane F. Trierweiler
1
1
Group of Intensification, Modeling, Simulation, Control, and Optimization of Processes, Chemical Engineering Department, Federal University of Rio Grande do Sul (UFRGS), R. Eng. Luiz Englert, s/n, Campus Central, Porto Alegre 90040-060, RS, Brazil
2
Food and Chemistry School, Federal University of Rio Grande (FURG), R. Barão do Cahy, 125, Cidade Alta, Santo Antônio da Patrulha 95500-000, RS, Brazil
*
Author to whom correspondence should be addressed.
Submission received: 27 July 2020 / Revised: 27 August 2020 / Accepted: 2 September 2020 / Published: 4 September 2020
(This article belongs to the Special Issue Redesign Processes in the Age of the Fourth Industrial Revolution)

Abstract

:
The demand for hydrogen in refineries is growing due to its importance as a sulfur capture element. Therefore, hydrogen management is critical for fulfilling demands as efficiently as possible. Through mathematical modeling, hydrogen network management can be better performed. Cost-efficient Mixed-Integer Linear Programming (MILP) and Mixed-Integer Nonlinear Programming (MINLP) optimization models for (re)designing were proposed and implemented in GAMS with two case studies. Linear programming has the limitation of no stream mixing allowed; therefore, to overcome this limitation, an algorithm-based procedure called the Virtual Compressor Approach was proposed. Based on the MILP optimal solution obtained, the streams and compressors were merged. As a result, the number of compressors was reduced, along with the inherent investment costs. An operational cost reduction of more than 28% (example 1) and 26% (example 2) was obtained with a linear model. The optimal MILP solution after rearranging compressors was then provided as a good starting point to the MINLP. The operating costs were decreased by more than 31% (example 1) and 32% (example 2). Most of the cost reduction was obtained only with the usage of the MILP model. Besides, a higher level of cost reduction was only obtained when the linear model was used as the starting point.

Graphical Abstract

1. Introduction

Hydrogen has a prominent role in the refining industry, as both its production and its recovery are essential steps. Hydrogen consumption in oil refining increased from approximately 7 million tons in 1980 to 38 million tons in 2018 [1]. Its importance is sustained by three factors: (i) the increase in the processing of heavier oils with high levels of sulfur and nitrogen; (ii) the increase in environmental constraints; and (iii) the production of derivatives of higher added value [2,3]. Due to this trend, it is necessary to use more efficient hydrogen within the petroleum refining process.
A hydrogen network consists of hydrogen-producing units, hydrogen-consuming units, and purification units, capable of purifying hydrogen to achieve the required purity. The hydrogen generation units (HGU) have become increasingly present in refineries due to the importance of hydrotreatment units (HDT) because its function is to supply the hydrogen demand complementing those generated in the catalytic reform. The steam reform is the primary process used at the industrial scale to obtain hydrogen as a primary product. Catalytic reform and purge gas can be used as a secondary source of hydrogen. The main hydrogen-consuming units are hydrotreating, which uses hydrogen to improve the quality of naphtha, kerosene, solvents in general, diesel oil, heavy gas oils, paraffin, and lubricating oils [4]. The management of the hydrogen network in a refinery implies in the material balance at all these units.
The need for optimization of the hydrogen network in refineries was recognized in the 1990s, because, usually, the amount of hydrogen produced is higher than the amount consumed. This excess is usually incorporated into the fuel gas system or burned directly into the flare. Therefore, it is necessary to have greater control in the sources and consumers of hydrogen through network management as a whole, because it is not economically feasible to produce and burn the product with an excellent added value [5]. It is known that the cost of hydrogen is the second-highest cost in a refinery, behind only the cost of crude oil [6]. Therefore, savings in terms of the amount of hydrogen consumed and the operating cost of the network have great economic appeal.
Since then, many methodologies have emerged to accomplish it. In general, these methodologies can be divided into two categories: pinch methods and optimization methods (deterministic in this case) as mathematical programming approaches based on network design [7]. Graphical methods provide an essential insight into the integration of the refinery process and provide theoretical goals for minimum hydrogen use. As oil refining and the hydrogen network involve many restrictions, they must be considered during network modeling and optimization, such as pressure, impurities, and equipment capacity. However, in graphic methods, this is not possible, as only flow and purity restrictions are considered. Therefore, mathematical programming is the best alternative and the most used, providing more realistic results and networks [8].
In this work, a mathematical programming approach was used to develop a model to solve the problem of hydrogen network optimization based on operating costs and constraints. The idea is to apply the proposed model to existing networks. The optimization allows the possibility of including new equipment and finding better ways of connection between units. For linear optimization, a compressor rearrangement technique was proposed in this work to decrease the capital cost. It is called Virtual Compressor Approach (VCA). The methodology was proposed to make the linear model competitive and satisfactory for the retrofit of hydrogen networks, due to its advantages and characteristics. Besides, a nonlinear model was also developed for comparison, with an initialization strategy using the MILP solution. This proposal was developed to facilitate the resolution of nonlinear and obtaining more competitive hydrogen networks.

2. Literature Review

Previous works on hydrogen distribution management and analysis using a linear programming model, based on the graphical analysis of the pinch method, were found in the literature. Towler et al [9] proposed a linear model to optimize a hydrogen network, aiming to minimize the total hydrogen import as an external utility. Two procedures for problem relaxation were proposed. The disadvantages of this method are that pressure constraints are negligible, and the flow merging must be performed manually [9]. Fonseca et al. [10] employed the linear programming model to optimize the hydrogen network of a refinery taking account pressure considerations and achieved a 30% reduction in utility use with the objective function minimizing the total flow rate of fresh hydrogen from a hydrogen plant [10].
Considering nonlinear programming (NLP), Hallale and Liu [11], in addition to mentioning the graphical pinch method, developed a nonlinear mathematical model to reduce the hydrogen consumption of the network. The model took into account pressure constraints, existing compressors, and a strategy to install a purifier. The objective function was to minimize the total cost, including operating and capital costs [11]. Liu and Zhang [12] developed a systematic procedure for integrating purification in hydrogen network design. For this, an MINLP (Mixed-Integer Nonlinear Programming) model for purifier selection and integration was used, and with linear relaxation of bilinear forms MINLP model was solved first as MILP because of the advantages of using linear models for problem solution [12]. Kumar et al. [13] developed mathematical models (LP (linear programming), NLP, MILP (Mixed-Integer Linear Programming), and MINLP) to obtain the best optimization problem in two case studies. Comparing MINLP and NLP for case 1, MINLP showed a more significant reduction in operating costs and equal capital costs. For case 2, the formulations LP, NLP, and MILP were compared. The NLP model imports less hydrogen and features a more realistic network than the others. The conclusions were that mixed-integer linear and nonlinear programming models are considerably better than linear because it provides the less complicated and more realistic refinery system, and MINLP can include complexities as compressors, purity constraints, and pressure constraints [13].
Liao et al. [14] developed an MINLP model using an existing hydrogen network with a purifier. The objective function was the total annual cost, and the model was solved in GAMS (General Algebraic Modeling System) using DICOPT. The MINLP problem is decomposed into a series of NLP and MILP solvers. The total annual cost decrease by 22.6% and both the new compressor and PSA were incorporated [14]. Birjandi et al. [15] developed a methodology for the optimization of a hydrogen network based on a simultaneously resolved MINLP and NLP problem. Linearization techniques for nonlinear models were used to facilitate resolution by transforming nonlinear equality constraints into inequality constraints. Global optimization has reduced operating costs [15]. Matijasevic [16] presented a hydrogen network integration methodology for a case study of a local refinery. The minimum consumption of hydrogen was determined by pinch analysis. Then, the superstructure was modeled using a nonlinear mathematical model whose objective function was to minimize total operating costs. The problem was solved with the GAMS software [16].
Unlike what was found in the literature, this paper developed a cost-efficient MILP and MINLP optimization models for (re)designing of hydrogen networks or a new project. The main difference from the MILP model to the MINLP is that it is not possible to mix streams in the compressors as it generates nonlinearity. To reduce the cost of capital from the MILP, in this work, a compressor-retrofitting tool was proposed respecting the nominal capacities. Also, to facilitate the resolution of the nonlinear formulation, an initialization strategy was used using the linear solution as a feasible starting point.

3. Mathematical Programming Approaches

Mathematical programming based on superstructure has advantages over pinch, in that it considers numerous limitations and variables when looking for solutions to the optimization problem. Limitations such as pressure, capacity, purity, operating costs, and investments in new equipment are some of the restrictions that may be included in the mathematical model formulation. The methodology to develop mathematical programming would be the development of the superstructure, including the sources, consumers, existing compressors, and purifiers. The formulation of the mathematical model also includes the objective function to be minimized or maximized subject to the set of constraints, the initialization strategy, and the resolution of the optimization problem. Typically, the objective function is the total annual cost of the hydrogen network [8].
Generally, the optimization problem can be formulated as a linear programming, mixed linear programming, nonlinear programming, or mixed-integer nonlinear programming problem. If linear combinations of variables can express the objective function and constraints, it is a linear optimization problem. Otherwise, the optimization problem is nonlinear. There are many optimization software used to solve optimization problems and already include algorithms called solvers [17].
Network management through mathematical modeling can be applied to an existing fixed topology, or to develop a new hydrogen network design. Thus, the approach of this article is based on the evaluation of the model developed for initial hydrogen network projects, through the validation with networks presented in articles already published. New equipment is considered, and the problem then becomes MILP or MINLP. Although the focus is operational, the problem addressed here is broader and has a significant industrial interest. The primary purpose of managing hydrogen networks is their production with minimum slack. Excess hydrogen production must be minimized, first because hydrogen is not easy to handle or store, and second, because it is not economically viable since the excess must be burned as fuel and furnaces and other processes.
The MINLP problems are more challenging to solve because they combine the NLP and MILP models and their characteristics. However, they result in more realistic networks and include several additional restrictions. According to the literature review, the use of MILP is not very recurrent, although when used, it presents significant results. Most articles found in the literature use nonlinear models for hydrogen network optimization. The advantages of using MILP is the linearity that facilitates the resolution of the optimization problem and the modeling of the logical constraints made in this article, which were not found in the literature. MILP problems are easier to converge to a global solution, since all the subproblems, for fixed binaries, are linear solved to global optimality [18,19].

3.1. Problem Statement

The problem to be addressed in this paper can be stated as follows: (i) a set of sources i ϵ hydrogen sources ( H S ), (ii) a set of consumers j ϵ hydrogen consumers ( H C ), and (iii) a set of purifiers k ϵ hydrogen purifiers ( H P = O H P N H P ), considering the existing purifiers, O H U , and the new purifiers, N H P . In the case of nonlinear formulation, there is still a set of compressors c ϵ hydrogen compressors ( H C P = O H C P N H C P , ) considering the existing compressors O H C P and new compressors N H C P . Figure 1 shows the two superstructures considered in this problem for the linear formulation (Figure 1a) and the nonlinear formulation (Figure 1b).
For each source, the maximum and minimum flowrate, as well as the hydrogen composition, and the outlet pressure are given. For each consumer, the inlet flowrate demand, pressure, and composition, the outlet purge flow, pressure, and composition are given. For each purifier, the maximum flow capacity, the composition of purified flowrate and purge flowrate, the pressure of purification, and the hydrogen recovery are given. It is also considered a fuel system in which waste streams can be burned and used as fuel to the process. For the existing networks, also given are the existing lines (unit connections), the distance between the units if informed, and the existing compressors (capacity and pressures) and purifiers.
The optimization problem is subject to the material balances and process operating constraints. For the retrofit case, process modifications are allowed to reduce the total operating costs (the objective function), despite the investment costs due to the installation of new pipelines, compressors, and possibly new purifiers.

3.2. Mathematical Model: MILP Formulation

Figure 1a shows the superstructure and all the possible connections among these four units between sources and consumers, sources and purifiers (existing and new ones), as well as flows between consumers and the purifying units for sources i and consumers j. The first step for the modeling development is to define which units are involved in the hydrogen network, for instance, which units provide hydrogen, which units consume hydrogen and the existing purifiers, and the potential purifiers that should be considered in the model.
The optimization problem of hydrogen network design in this work can be summarized as follows: the superstructure is formed by a set of sources of hydrogen i, a set of hydrogen consumers j and set of units of hydrogen purification k, account for the existing and new purifiers. The hydrogen sources have their minimum and maximum flow according to their capacity ( F H 2 I i , m i n   e   F H 2 I i , m a x ) as well as their hydrogen purity ( Y I i ). The hydrogen-rich stream can be sent to the consumers j ( F I J i , j ), to purification units k ( F I K i , k ) , or can be sent to the fuel system ( F I W i ) . The consumer’s units also have their known, and constant input required flows for the process ( F J j ), as well as its hydrogen purity ( Y J j ), in addition to the outflows ( F P j ) and hydrogen purity ( Y P j ), according to the hydrogen consumption of each specific process. The outlet flows from the consumers can be sent to purification ( F J K j , k ), can be used as a source for other consumers ( F J J j , j ) or can be sent to the fuel system ( F J W j ) to be used as the burning fuel. The purifying units have a known hydrogen recovery ratio ( r e c k ), as well as the maximum inlet flow capacity ( F P u r m a x , k ) and the constant purities of the hydrogen product pure streams ( Y K k ) and the composition for the stream of hydrogen not recovered stream ( Y K W k ) . The purified hydrogen stream from the purification can be used as a source for the consumers ( F K J k , j ) who need higher purity or can be referred to the fuel system ( F K W k ), if there is excess. The stream with the unrecovered hydrogen, F K W r e c k , has a small hydrogen composition, and it is sent directly to the fuel system. In this work, some considerations were made to simplify the model. The flowrates are considered only a binary mixture of hydrogen and methane. The partial pressure of the hydrogen and the flowrate are constant at the entrance and exit of the consuming units.

3.2.1. Sources

The overall material balance for each source is represented by Equation (1):
F H 2 I i = ( j H C F I J i , j + k H P F I K i , k + F I W i )   i H S
where F H 2 I i is the total flow from each source i, F I J i , j is the hydrogen flow from the source i to the consumer j, F I K i , k is the flow from the source i for the purification unit k, and F I W is the flow from source i sent to the fuel system. The available flow rate is limited by the capacity of the hydrogen generating units according to the following inequality constraints:
F H 2 I i , m i n F H 2 I i F H 2 I i , m a x i H S

3.2.2. Consumers

Equation (3) represents the overall material balance in the inlet of consumer units.
F J j = i H S F I J i , j + k H P F K J k , j + j H C F J J j , j   j H C
where F J j is the total flow directed to consumers, F J J j , j is the flow from one consumer j to another consumer j′ and F K J k , j is a flow rate from the purification unit k for the consumer units j. The index j’ which is used for cases where there is a connection between consumers. In this case, as it is not allowed between the same unit, j′ must be different from j. The hydrogen balance is then defined by Equation (4):
F J j Y J j = i H S F I J i , j × Y I i + k H P F K J k , j × Y K k + j H C F J J j , j × Y P j   j H C
where Y J j ,   Y I i ,   Y K k   and   Y P j are the volumetric fractions of hydrogen in the respective streams, consumer j, sources i, purifiers k, and purge of the consumer unit j. Besides, it is possible to calculate how much each consumer unit used hydrogen depending on the chemical process involved.
Equation (5) represents the overall material balance in the outlet of consumer units:
F P j = F J W j + k H P F J K j , k + j H C F J J j , j   j H C
where F P j is the total flow out of consumers, F J K j , k is the flow rate from the consumer unit j for the purification unit k, and F J W j is the surplus flow of consumers directed to the fuel system.

3.2.3. Purification Units

The purification unit is used, so that process streams are purified, providing hydrogen in a given purity, such as 99.99% in the case of PSA units. The overall material balance in these units is expressed as:
j H C F J K j , k + i H S F I K i , k = j H C F K J k , j + F K W k + F K W r e c , k   k H P
where F K W k the flow rate of the purifying unit k stream rich in hydrogen routed to burning and F K W r e c , k is the hydrogen flowrate not recovered by the purifying unit k sent to the burner. The hydrogen balance for each purifier is described as follows:
j H P F J K j , k × Y P j + i H S F I K i , k × Y I i = j H P F K J k , j × Y K k + F K W k × Y K k + F K W r e c , k × Y K W k   k H P
where Y K W is the fraction of hydrogen in the purge stream of purified k. The total flow entering the purifier is limited by the capacity of the purifying unit.
j H P F J K j , k + i H S F I K i , k k F P u r m a x , k   k H P
Given the hydrogen recovery of the purification unit, it is possible to calculate how much hydrogen is sent to the purge stream, i.e., the hydrogen not recovered.
( i H S F I K i , k × Y I i + j H P F J K j , k × Y P j ) × ( 1 r e c k ) = F K W r e c × Y K W k   k H P
The total flow through the PSA ( F K k ) can then be defined as:
j H P F J K j , k + i H S F I K i , k = F K k   k H P

3.2.4. Logical Constraints

To consider the capital cost associated with new equipment, it is necessary to use constraint modeling, through logical propositions and disjunctions, so binary variables and logical inequality equations were included in the model with binary parameters. First, through the modeling of disjunctions, a binary variable z is associated with the existence of a particular flow F (e.g., F I J i , j , F K J k , j , F J K j , k , etc.). If the positive flowrate is greater than or equal to a small value ε , e.g., ε = 10 5 , the corresponding binary variable z assumes the value of 1. On the other hand, if the flowrate is lower than ε , the binary variable assumes the value of 0. F m a x are the flowrates between the units involved. These conditions are ensured by the following constraints:
{ F ε × z F ( min ( F m a x ) ) × z
A binary variable z c is associated with the installation of a compressor for the corresponding flow. For this case three events must hold simultaneously: (i) there is a non-zero flow, i.e., z = 1; (ii) there is no compressor previously installed identified by a binary parameter uc (1 if there is an existing compressor, 0 otherwise); and (iii) there is a pressure difference between the current unit and destination unit that requires a compressor identified by a binary parameter u d e l t a P (1 if the current pressure unit is lower than the destination pressure unit, 0 otherwise).
z c z + u d e l t a P + ( 1 u c ) 2
If any of these three events is false, then there is no need to install a compressor ( z c = 0), which is ensured by the set of constraints described in the set of Equation (13).
{ z z c 1 u c z c u d e l t a P z c
A similar procedure was used to consider the investment cost of piping. A binary variable z h is associated to the need of installing a new pipeline if two events hold: (i) there exists a non-zero flow in that connection, i.e., z = 1; (ii) there is no pipeline previously installed identified by a binary parameter u h (1 if there is a line, 0 otherwise).
z h z + ( 1 u h ) 1
If any of these two events do not hold, it must be ensured that no pipeline must be installed.
{ z h z z h 1 u h
There is also the possibility of installing new purification units. In this case, it is enough that there is any flow entering or leaving this unit. In this case, a binary variable z k n is associated with the installation of a new purifying unit and the logical constraints can be expressed by:
{ F K k ε × z k n F K k ( F P u r m a x , k ) × z k n k N H P
The same procedure for installing new compressors was also done (constraints in Equations (12) and (13)) if it is necessary to install new compressors on streams involving a new PSA.

3.2.5. Operating Costs

Operating costs include the production of hydrogen, the cost of electricity used in compressors, the operating cost of the purifying units, and the economic value corresponding to the burning gas in the fuel system. The cost of hydrogen production is assumed directly proportional to the flowrate, and it is defined as follows:
C H 2 I = i H S F H 2 I i × C i
where C i is the cost of producing hydrogen. The electricity cost of the compressor is directly proportional to the power ( W ):
W = F × w
where W is the power of the compressor with the flowrate being compressed F ,   w is the intensive power estimated from the stream properties ( C P , C V , z), the inlet and outlet pressure, and the compressor efficiency [11].
w = ( C P ¯ × T / η ) × ( ( P o u t P i n ) γ 1 γ 1 ) × ( ρ o / ρ )
where C P is the heat capacity, T is the stream temperature, η the efficiency of the compressor, P o u t and P i n are the outlet and inlet pressure, respectively, ρ o and ρ are the densities at design conditions and standard conditions, respectively, γ is the ratio of the heat capacity at constant pressure to that at constant volume. For a given connection, e.g., F I J i , j , the corresponding intensive power w i , j is previously calculated as a model parameter. For the complete model, the total electricity cost is calculated by the following Equation (20). The indices α and β represents the possible connections involved (i,j; j,k; k,j; j,j′; i,k; i-waste; j-waste; k-waste):
C H 2 C = ( α β F α , β × u d e l t a P α , β × w α , β ) × C e l e t r i c
where C e l e t r i c is the electricity cost. It is worth to note that each term is multiplied by the binary parameter u D e l t a P (1 if the pressure ratio is higher than one), for the cases in which the flowrate is not zero, but there is no need for compression. It does not matter if a new compressor is installed or an existing compressor is used, both consumes energy. Equation (20) will compute the energy cost correctly, and it takes into account the electricity used in existing and new compressors.
The cost of purifying unit is proportional to the feed flowrate:
C H 2 K = k H P F K k × C k
where C k is the cost of using the PSA purification units, new and existing ones.
The economy value corresponding to the burning of excess purge flows is corresponding to the cost of hydrogen and methane used as fuel and calculated as:
C H 2 F = C f u e l × F × ( y × Δ H ° H 2 + ( 1 y ) × Δ H ° C H 4 )
where C f u e l the cost per unit of energy, F is the gas flowrate, and y is the hydrogen composition. Assuming a binary mixture, 1 y represents the methane composition. The parameters Δ H ° H 2 and Δ H ° C H 4 are the standard heat of combustion of hydrogen and methane, respectively. For the complete model, taking into account the total contributions, the economic value corresponding to the total cost of fuel is calculated as follows:
C H 2 F T = C f u e l × α F α × [ y α × Δ H ° H 2 + ( 1 y α ) × Δ H ° C H 4 ]
The subscript α denotes all units sending streams to the fuel system (i, j, k). Since it corresponds to a saving cost, this value must be subtracted from the total operating cost. The operating cost parameters assumed in this work are presented in Table 1. The assumed values were the same used in example 1 [11], a case study of this work, also chosen based on the reviewed articles.

3.2.6. Investment Costs

The capital cost includes the cost of new compressors ( C n e w   c o m p r e s s o r ), new purification units ( C n e w   P S A ) and new pipelines ( C p i p i n g ) . Hallale and Liu [11] describe the cost for the inclusion of new compressors for a particular flowrate, with a fixed cost with a binary variable and a variable cost associated with the flow:
C n e w   c o m p r e s s o r = a × z c + b × W
W is calculated by Equation (18) and z c is the binary variable associated with the installation of a compressor for the corresponding flow and multiplied the fixed part of the new compressor cost, so it is considered only when the compressor is installed. The complete equation for accounting the new compressor cost is given by Equation (25). The indices α and β represents the possible connections involved (i,j; j,k; k,j; j,j′; i,k; i-waste; j-waste; k-waste).
C n e w   c o m p r e s s o r T = a × ( α β z c α , β ) + b × ( α β F α , β × u d e l t a P α , β × w α , β × ( 1 u c α , β ) ) × C e l e t r i c
The cost associated with the installation of new piping is described below, including a fixed part with a binary variable and a variable part dependent on flowrate. For these calculations, it is necessary to inform the distances between the already installed units of design.
C n e w   p i p i n g = ( c × z h + d × D 2 ) × L
With
D 2 = ( 4 × F / π × ϑ ) × ( ρ o / ρ ) = ( 4 × F / π × ϑ ) × ( T T 0 ) × ( P 0 P )
where L is the pipe length [m], c and d are constants, ϑ is the gas surface velocity (usually 15–30 m/s; assumed an average value of 22.5 m/s in this work), and D² is the equivalent square diameter [11]. The binary variable z h indicates the need to install the new pipeline. Equation (27) is replaced in Equation (26) in order to express the cost of piping as a function of the flowrate. The equation for the model (total cost of new piping) is represented by Equation (28). The indices α and β represents the possible connections involved (i,j; j,k; k,j; j,j; i,k; i-waste; j-waste; k-waste). Each term is multiplied by ( 1 u h α , β ) in order to consider only the cost of new piping.
C n e w   p i p i n g T = c × ( α β z h α , β × L α , β ) + d × ( α β F α , β × L α , β × w × ( 1 u h α , β ) × ( T T 0 ) × ( P 0 P ) )
There is also the possibility of installing new purification units. For this case, the cost of a PSA unit (purifier considered in this work) is a linear function of the unit flowrate (variable part) and include binary variable corresponding to the fixed installation cost:
C n e w   P S A = a P S A z k n + b P S A × F i n ,   P S A
where a P S A and b P S A are constants and F i n ,   P S A is the inlet flowrate of the PSA unit (MMscfd). The binary variable z k n is associated with the installation of a new purifying unit. The model equation is described as:
C n e w   P S A T = a P S A k N H P   z k n + b P S A × ( k N H P F K k )
This cost is only considered for new purifying units. The capital cost parameters used in this work are presented in Table 2. Different coefficients exist for calculating capital costs, including variations in temperature and materials involved. The most frequently used data in the reviewed papers were used, following Hallale and Liu [11]. The objective is to facilitate the comparison of the results obtained.

3.3. Formulation of the Optimization Problem

Based on all the costs involved in managing the hydrogen network described in the previous section, annual operating and annual capital costs are defined as:
C o p e r a t i n g = ( C H 2 I + C H 2 K + C H 2 C C H 2 F ) × t
C c a p i t a l = ( C n e w   P S A + C n e w   p i p i n g + C n e w   c o m p r e s s o r ) × A f
where A f is the annualizing factor, and t is the considered operating time of the plant in one year. The annualizing factor is defined by:
A f = f i × ( 1 + f i ) n / ( 1 + f i ) n 1
where n is the number of years of interest for the return on investment and f i is the interest rate. The Total Annual Cost (TAC) consists of the summation of the operating and investment cost:
T A C = C o p e r a t i n g + C c a p i t a l
For the retrofit case of existing networks, the economy saving used as economic criteria is calculated as:
E = C O P a c t u a l C O P n e w
where C O P a c t u a l   a n d   C O P n e w are the operating cost of the actual and new networks, respectively. The payback time is defined by the ratio of the total investment cost and the economy saving, and the following equation can estimate it.
p t = C c a p i t a l / A f E = ( C n e w   P S A + C p i p i n g + C n e w   c o m p r e s s o r ) C O P a c t u a l C O P n e w
The MILP model formulated in this work is described by the set of constraints (1, 2, 3–17, 20, 21, 23, 25, 28, and 30—HNS LM (Hydrogen Network Synthesis—Linear Model)). For process optimization, different objective functions can be chosen to be minimized. In this case, operating cost (31) for the retrofit case was chosen. The proposed model has the advantage of being a linear model, for which quite robust solvers can be used. However, the main drawback is that a compressor is associated with each possible connection individually in order to avoid nonlinear material balances to identify the composition of the stream being compressed. For this case, streams cannot be mixed to use the same compressor, and the resulting network may end up with more compressor units than an alternative nonlinear model, in which streams are allowed to mix.

3.4. Mathematical Model: MINLP Formulation

A nonlinear model was also developed. In this model, the compressors are considered as independent units that may be used to connect units that need compression (see Figure 1b). Different from the other units, the inlet and outlet pressure of each compressor are free variables. The maximum number of compressors to be considered is set in the superstructure modeling, and it is obtained in the model solution previously. In this model, streams are mixed to enter the compressor. Therefore, the hydrogen composition is unknown and must be treated as a variable. Besides, since no compressors are associated with each stream individually, the flowrates are only possible if the current origin pressure is higher than the destination pressure. For a particular flow F with upper bound F m a x , the constraints (37) ensure that flow is only possible for this case (higher pressure to lower pressure):
F F m a x × ( 1 u d e l t a P )
Despite the possibility of generating networks with fewer compressors, the nonlinearity comes up with a more difficult problem to be solved that is very dependent on the initial guess, as will be discussed later.
In the MINLP model, the superstructure is a bit different from the one presented, as illustrated in Figure 1b. In this case, the compressor is considered a unit of the network and, therefore, can have the same source (the compressor outlet) and consumer (the compressor inlet) functionality and must be present in the balance equations. The only nonlinearity in this model arises in the hydrogen balance in the inlet of the compressors because there is the merging of flows and, consequently, the product flow/purity. It is necessary to know the inlet composition because the outlet flow with this composition is sent to other units, and the hydrogen balances depend on this value.
The equations that describe the nonlinear model are described below. Equations (1), (3)–(9) of the linear model are replaced by the equations below, as compressors need to be considered in material balances. In sources, in addition to Equation (2), there is Equation (38), which describes the sum of flow rates from sources for consumers, purifiers, compressors ( F I C i , c ) and for burning. Hydrogen from the source can be sent to all these units.
F H 2 I i = ( j H C F I J i , j +   k H P F I K i , k + F I W i + c H C P F I C i , c )   i H S
For consumers, global and component material balances are made, where F C J c , j   is the flowrate from the compressor to the consumers and F J C j , c is the flow rate from consumers to compressors. The sum of the flowrate at the entrance of each consumer corresponds to the sum of the flowrate from the source, the purifier, another different consumer, and the compressor.
F J j = i H S F I J i , j + k H P F K J k , j + j H C F J J j , j + c H C P F C J c , j   j H C
The same is true for the hydrogen balance, where in addition to flowrates, purities are considered. Here there is the purity of the compressor ( Y C c ) .
F J j × Y J j = i   H S F I J i , j × Y I i + k H P F K J k , j × Y K k + j H C F J J j , j × Y P j + c H C P F C J c , j × Y C c   j H C
The sum of the outlet flowrate of each consumer corresponds to the sum of the flowrate that the consumer forwards to the burn (waste), to the purification unit, to another different consumer, and the compressor if necessary.
F P j = F J W j + k H P F J K j , k + j H C F J J j , j + c H C P F J C j , c j H C
The global material balance and for hydrogen is also applied for purifiers. The material balance corresponds to the sum of all flowrates at the entrance of the PSA, which include the flowrates from consumers, sources, and compressors. The purification unit, in turn, can send flow to consumers, compressors and can burn the excess (waste), which can be seen in Equation (42). Equation (43) corresponds to the hydrogen balance, considering the flows directed to the purifier and forwarded from the purifier. In addition to these equations, the purified flow rate must not exceed the PSA capacity (Equation (44)), and, through the recovery of the PSA, the flowrates that are sent for burning are obtained (Equation (45)).
j H C F J K j , k + i H S F I K i , k + c H C P F C K c . k = j H C F K J k , j + F K W k + F K W r e c , k + c H C P F K C k , c   k H P
j H P F J K j , k × Y P j + c H C P F C K c . k × Y C c + i H S F I K i , k × Y I i = j H P F K J k , j × Y K k + c H C P F K C k , c × Y K k + F K W k × Y K k + F K W r e c , k × Y K W k   k H P
j H P F J K j , k + i H S F I K i , k + c H C P F C K c , k F P u r m a x , k   k H P
( i H S F I K i , k × Y I i + j H P F J K j , k × Y P j + c H C P F C K c . k × Y C c ) × ( 1 r e c k ) = F K W r e c × Y K W k )   k H P
where F C K c . k is the flow rate from compressors to purifier, F K C k , c is the flow rate from the purifiers to the compressors. Also, as the compressors are like units in the hydrogen network, material balances are made. The sum of the flow that enters the compressors is called F C c , which consists of the sum of the flows from sources, consumers, and purifiers.
F C c = c     H C P F I C i , c + c H C P F J C j , c + c H C P F K C k , c   c H C P
Therefore, any flow that enters the compressor must be directed to the consumers and purifications units. If necessary, some part of the compressor flow that is not used can be sent directly for burning.
c H C P F I C i , c + c H C P F J C j , c + c H C P F K C k , c = c H C P F C J c , j + c H C P F C K c . k + F C W c   c H C P
It is also necessary to carry out the hydrogen balance in the flows that make up F C c .
F C c × Y C c = c H C P F I C i , c × Y I i + c H C P F J C j , c × Y P j c H C P + F K C k , c × Y K k   c H C P
In the same manner as in the MILP model, a binary variable z is associated with each possible flowrate, including the flowrates involving the compressor units, e.g., F I C i , c , F J C j , c , F K C k , c , F C J c , j , F K C k , c , and F C W c . The corresponding constraints are as described by Equation (10). Also, binary variables are associated with new pipelines (Equations (13) and (14)) and for new PSA (Equation (15)). The binary variable z c α , β are used to define if the compressor unit is installed assuming the value of 1, 0 otherwise. Differently from the MILP model, z c is not defined over a pair of streams; it depends only on the compressor unit.   F C c is associated with the flow of each compressor. Constraints (Equation (49)) is used to establish which compressors are used and their flow rates.
{ F C c ε × z c F C c F m a x × z c
As the pressures vary in the nonlinear model, pressure restrictions must be included, which guarantees the compressor inlet and outlet pressures. They are formulated in the same format as the logical flow restrictions. For a given compressor unit, the inlet pressure is set as lower than the minimum pressure among the pressure of the mixed streams entering the compressor (Equation (50)). The outlet pressure is set as higher than the maximum pressure among the pressure of the streams, leaving the compressor according to the pressure of the stream destination (Equation (51)). It is important to mention that, due to the minimization of the energy cost associated with the compressor in the objective function, which is proportional to the pressure ratio ( P C o u t , c / P C i n , c ), the inlet pressure is set as the minimum stream pressure entering the compressor c, and the outlet pressure as the maximum stream pressure leaving the compressor c.
{ P C i n , c P P j + ( P m a x P P j ) × ( 1 z j , c ) P C i n , c P I i + ( P m a x P I i ) × ( 1 z i , c ) P C i n , c P K k + ( P m a x P K k ) × ( 1 z k , c )
{ P C o u t , c P J j P m a x × ( 1 z c , j ) P C o u t , c P K k P m a x × ( 1 z k , j ) P C o u t , c P W P m a x × ( 1 z c , w )
where P C i n , c , and P C o u t , c are the compressor c inlet and outlet pressures, respectively, the binary variable z is associated with flowrates (i.e., z i , c ,   z j , c , z k , c ) and P m a x is the maximum pressure of the network used to make the constraints (50) and (51) redundant for the corresponding non-existent connection (the corresponding binary is set to zero due to the zero flowrate).
The operating and capital costs are calculated in the same way as in the linear problem, as well as the logical flow restrictions. The cost of hydrogen production is obtained by Equation (17), Equation (52) represents the electricity cost, Equation (53) represents the purification cost, and cost of fuel is represented in Equation (54).
C H 2 C = C e l e t r i c ×   c H C P F C c × ( C P T / η ) × ( ( P C o u t , c P C i n , c ) γ 1 γ 1 ) × ( ρ o / ρ )
C H 2 K = k H P ( j H C F K J k , j + F K W r e c + F K W k + c H C P F K C k , c ) × C k
C H 2 F T = C f u e l × α F α W α × [ y α × Δ H ° H 2 + ( 1 y α ) × Δ H ° C H 4 ]
The subscript α denotes all units sending streams to the fuel system (i, j, k, c). Equation (55) represents the cost of new compressors and Equation (56) the cost of new piping.
C n e w   c o m p r e s s o r = a × z c n e w c + b × F C c n e w c × ( C P T / η ) × ( ( P c o u t P c i n ) γ 1 γ 1 ) × ( ρ o / ρ )
C n e w   p i p i n g = c × ( α β   z h α , β × L α , β ) + d × ( α β F α , β × L α , β × w × ( 1 u h α , β ) × ( T T 0 ) × ( P 0 P ) )
The indices α and β represents the possible connections involved (i,j; j,k; k,j; j,j′; i,k; i-waste; j-waste; k-waste; i,c; j,c; k,c; c,j; c,k; c-waste). The MINLP model formulated in this work is described by the set of constraints (1, 2, 11, 14, 15, 16, 17, 37–56). The objective function is described in Equation (31). This MINLP model will be named to facilitate the description of the results by HNS NLM (Hydrogen Network Synthesis—Nonlinear Model).

3.5. Virtual Compressors

The main difference between the MILP model and the MINLP model is how the compressors are treated. In MILP, the compressors are associated with each particular flowrate. In this case, the streams are not mixed. However, in the MINLP, the compressors are treated as independent units, not associated with a flowrate. Then the stream can be mixed to enter the compressor and split leaving the unit. Besides the class of the resulting model (either linear or nonlinear), the linear model may result in a network with more compressors and pipelines than the nonlinear model. Both the linear and the nonlinear formulation are capable of representing the hydrogen network, so what differentiates them is the issue of allowed linearity (which can be improved through this proposed technique), the linear model is simpler to solve, and the global optimum solution is guaranteed.
To overcome a large number of compressor units and further investment cost reduction, a strategy to reduce the use of this equipment was carried out through an algorithm based on non-real streams or virtual compressors, i.e., it is possible to rearrange the streams and compressors if the compressor capacities were not reached. This developed technique is one of the contributions of this work. Through it, the linear model becomes competitive, compared to the nonlinear model, due to its advantages.
There are two cases where it is possible to perform this unit reduction: (Option 1) when there are streams with different composition being compressed and forwarded to the same unit or (Option 2) when streams coming from the same unit are compressed and forwarded to different units, as can be seen in Figure 2. In other words, it is possible to group streams and use the same compressor, thus decreasing the fixed part of the new compressor capital cost, since the variable part is flow dependent and does not change. It is worth nothing that the fixed cost of piping is also minimized due to the rearrangement of the streams.
For each option, the inlet pressure (in Option 1) and the outlet pressure (in Option 2) must be corrected according to the minimum and maximum pressure of the involved streams, respectively. In that case, the energy cost and the variable part of the investment cost must also be recalculated. It should be noted that using this procedure, the solution is not unique, and the best solution is that with the maximum total cost reduction. Despite eventually unfavorable pressure changes, the number of compressor units can be reduced. Therefore, when this procedure is performed, the investment cost is almost always reduced. In this work, since the number of possible rearrangements is small, this procedure was performed by enumeration.

3.6. Solution Strategy

In this work, the MILP and the MINLP model were used to the network (re)design. Compared to the linear models, nonlinear models are wholly dependent on the initialization, which has a more challenging convergence. Also, for MILP models, the global solution can be obtained without a high computational effort. The MILP solution can be rearranged to reduce the compressor units, with the virtual compressors approach. Besides, the solution obtained by the MILP model can be used as a good and feasible initial point for the MINLP model. It is crucial to the grassroots designs since, in the retrofit case, the existing network can be used as an initial point. All these possibilities were evaluated in this work, and further discussion is presented in the results section.
The initialization strategy used can be described as follows:
  • The flowrates are fixed according to the existing network for the retrofit cases, and an LP subproblem with Fobj = 0 subject to the material balances is solved to obtain a feasible solution.
  • The binary variables (z) are initialized according to the existing network, i.e., z = 1, where there is a non-zero flowrate, z = 0 otherwise. Also, the other binary variables (zc, zh, zkn) are fixed to zero, since they represent the installation of new compressors, piping, and purifying units.
  • The complete MILP model is solved. This result is defined as the existing network of each case study for later optimization (BASE CASE).
  • With all the variables values in the feasible solution defined by the existing network, the variables are set as free according to their lower and upper bounds. The complete MILP (HNS LM) is solved (objective function = minimize operating cost). The MINLP (HNS NLM) proposed model is also solved to compare with the item (6).
  • The optimized network obtained through the linear model is evaluated with the rearrangement of compressors. Here the values of operating cost and capital differ between them due to the decrease in the number of compressors and the possible increase in electricity.
  • This network design is used to initialize the MINLP nonlinear model (HNS NLM).
For all cases, it was possible to ensure that the starting point was a feasible point. Figure 3 summarizes the initialization techniques performed.

4. Results

The model described in the previous section was validated using two examples of hydrogen networks proposed in the literature. The mathematical programming model was implemented in the modeling system GAMS 22.2 on a 3.6 GHz Intel® Core™ I7 CPU ( GAMS Development Corporation, Washington, DC, USA).The solver used to solve HNS LM was CPLEX (CPLEX 10, GAMS Development Corporation, Washington, DC, USA, 2006), and for HNS NLM it was DICOPT (/DICOPT 2x-C, GAMS Development Corporation, Washington, DC, USA, 2006).
For the case studies, it was considered the retrofit design for existing hydrogen networks. Therefore, the existing structure was explored considering the installation of new pipelines, new compressors, and purifying units. The economy saving is obtained by the operating cost reduction compared to the original solution. However, there is also an investment cost associated with non-existing equipment and pipelines. The payback time, i.e., the investment cost divided by annual operating cost savings was also used as an economic indicator for comparing the model solution.
The original network was ensured as a feasible starting point for all optimization problems. It was accomplished by fixing all the values of stream flowrates according to the existing network, and the total operating costs were calculated according to the parameters listed in this work for each case study. For all cases, the original network was a feasible point. However, some authors have not presented the value of the parameters used to estimate the costs. Therefore, for a fair comparison, the costs were recalculated with the listed parameters in this work, and hence, despite the network configurations and flowrates are the same presented here, the costs are similar but not the same. Further discussion and considerations are given for each example.

4.1. Example 1

The first example is from Hallale and Liu [11]. The hydrogen network depicted in Figure 4 consists of a primary hydrogen production unit (H2plant) and a secondary source, which is catalytic cracking (CCR). In this process, there are six consumer units: HC (hydrocracker), JHT (kerosene hydrotreater), CNHT (cracked naphtha hydrotreater), DHT (diesel hydrotreater), NHT (naphtha hydrotreater), and IS4 (hydrodealkylation). Two previously installed compressors are used, and there are no purification units. Flowrates are expressed in MMscfd (million ft³/day, under standard conditions), stream purity, flowrates, and pressures are shown in Table 3.
The objective function chosen for the problem analysis was to minimize the operating cost of the hydrogen network, Equation (31), using the parameters listed in Table 3 and the network configuration depicted in Figure 4. A variation of ± 10% (vp) in the nominal flow of consumers was allowed, F J j and F P j   were allowed in the original article. For the installation of a new PSA, the purity of 99.99% with a maximum operating capacity of 50 MMscfd, a recovery rate of 90%, and purge purity of 40.2% was considered.
The annual operating costs for the original network were estimated at 39.819 $/year. This solution is referred here as Hydrogen Network -BASE CASE (HN0). The Hydrogen Network -BASE CASE corresponds to the existing basic topology.
The HNS LM model has about 763 single equations, 323 single variables, and 227 discrete variables. Through linear optimization, savings of $11 million per year were achieved with a total investment of $16 million. In this case, 9 new compressors and 16 new pipelines were installed, as well as a new PSA (HN1). Nearly a 28% reduction in operating cost was achieved. This network is shown in Figure 5a. The HN1 optimized network MILP model only imports 26.5 MMscfd, and the original network uses 44.9 MMscfd of hydrogen from H2 Plant, which represents a reduction of almost 41% in the amount of imported pure hydrogen.
In the HNS LM optimization, the merging of flows before the compressor units is not allowed. Therefore, the solution may result in a large number of installed compressors. However, the number of compressors can be reduced after the optimization, evaluating the obtained network, and, possibly, an even more significant cost reduction can be achieved. For the cases in which more than one stream leaving one unit is compressed and/or more than one stream is compressed to one unit, the streams can be rearranged to be compressed in a unique compressor unit saving the fixed cost associated to the compressor investment. According to the distance of the units, the cost of the pipeline must also be recalculated. As more than one alternative for the evolutionary network is possible, but they are only a few, this procedure can be executed manually by the designer. Therefore, we analyzed which compressors were already previously installed based on the units and purity involved and if their nominal capacity allowed them to receive more streams. If positive, the stream was directed to it, and the new associated compressor could be eliminated. The rearrangement technique using virtual compressors applied to the compressors of example 1 can be seen in Figure 5b.
According to the optimization result (HN1), 9 new compressors were installed, which can be rearranged, as explained in Figure 2. According to option 1, where different flow rates that go to the same unit are grouped, rearranging in only 2 new compressors and using the two existing ones. The total cost of these new compressors is $0.271 million ($0.230 million of the fixed cost and $0.041 million of the variable cost), and this represents an 86.4% reduction in the total investment in new compressors. The total cost of piping also reduces by 40% due to the rearrangement of the compressors. This impact on total investment is 12.4% less.
It should be noted that as the compressors are rearranged, the inlet pressure is the lowest pressure between the flows. Therefore, the cost of electricity is slightly changed due to this, so the cost of electricity increased by 14.5% (from $0.136 million to $0.156 million) and an increase of 0.06% in operating costs. The proposed new topology can be analyzed in Figure 5c, and HN1 will represent that network.
To compare the linear and nonlinear formulations, the original network was optimized through the nonlinear model HNS NLM, described in Section 3.4. The first initialization used here was the original network, in example 1 (Figure 4). The HNS NLM model has 944 single equations, 473 single variables, and 308 discrete variables. The operating cost obtained was $28.183 million per year and $7.846 million per year of capital cost (one PSA and 10 pipelines), called network HN2. The result obtained in the two proposed optimized networks is very similar; however, the nonlinear has fewer connections (Figure 6a). The most significant portion of the cost of capital corresponds to the quantity to be purified. The optimization of HN2 network is an integer solution (not an optimal as in HNS LM), which usually happens in nonlinear problems as it is not possible to guarantee optimum global optimization.
The second initialization made, which is the biggest contribution of this work, uses the result obtained from the HNS LM (HN1′- with compressors rearrangement) as the initialization of the nonlinear model HNS NLM, to facilitate the resolution of the nonlinear model. As already mentioned above, the HN1′ network with the rearrangement of the compressors has a significant reduction in the cost of new compressors. For this reason, it is an excellent point option for the nonlinear model. Besides, as can be seen in the results, since nonlinear optimization has great locations, this initialization helped to improve the result. The HN3 network (obtained using MILP as a feasible point in MINLP) resulted in the lowest operating cost, a reduction of 31.2% (Figure 6b). However, comparing the payback, which refers to both the economy and the necessary investment, the network with the lowest payback is HN1′. This shows that with the HNS LM model, good and significant results are achieved, but through nonlinear optimization, less complicated networks with lower operational costs are achieved. For this, it is important to evaluate the design of the proposed network through different initializations.
All the results obtained in the different optimizations are summarized in Table 4. It is observed that the most significant reduction in the operation cost was obtained in the HN3 network. However, taking into account the investment and the payback time, the HN1’ network proves to be an excellent alternative. Through the results obtained, it can be concluded that the two described models (linear and nonlinear) are efficient for the proposed optimization. The linear model is good enough and capable of providing considerably improved solutions. Besides, as an initial guess for the nonlinear model, it proved to be an even more competitive alternative. The compressor rearrangement technique provides a reduction in investments. When used to initiate the optimization of the nonlinear model, it provides designs with fewer lines and compressors.
As the original article of this case study does not present clear information about parameters and conditions used in the optimization [11], this work differs in values from the presented network. However, it is noteworthy that although the cost of the original network is different due to the explained, in this work, we considered the same calculation methodology for the original network (base case- HN0) and optimized network (HN1, HN2, HN3 …), with specific parameters and conditions chosen.
The result obtained from the optimization in Hallale and Liu [11] is a 26.6% reduction in operating cost and payback time of 1.6 years, whose objective function was to reduce operating costs, limiting the payback time to 2 years. The achieved results obtained here with the proposed methodology are satisfactory as HN1 (HNS LM) optimized network reduced by 28.1% the cost of operating with a payback of 18 months. The optimized HN2 (HNS NLM) network achieved a 29.3% reduction in the operating cost with a payback time of 16 months, while Hallale and Liu [11], reduced operating cost by 15%, with a 17 months payback. For this reason, the result obtained was better than that presented in the original article, as in percentage, a more significant reduction in operating cost and payback was achieved. With the proposal to use the linear solution as a feasible point, HN3 network, the reduction was even higher (31% in operating cost), which shows the efficiency of the proposed technique.

4.2. Example 2

The second example used is from Sardashti Birjandi et al. [15]. The network is made up of two hydrogen producing units, a catalytic cracking plant (CCR) and a hydrogen generating unit (H2plant), two purifying units (PSA), and 3 hydrogen consuming hydrotreating units (HDT I, HDT II, and HC), as illustrated in Figure 7. In addition to the information, some parameters described in Table 5 are required. This HNS LM model has 524 single equations, 224 single variables, and 158 discrete variables.
The annual operating costs for the original network were estimated at 44.017 $/year. This solution is referred here as Hydrogen Network -BASE CASE for example 2. This network corresponds to the existing basic topology (Figure 7).
Minimizing only the operating cost of the hydrogen network, savings around $12.4 million per year are achieved (HN4). For this design, the total investment of $22 million is paid off in 22 months. Six compressors, 10 new lines, and a new PSA were installed. The operating cost was reduced by 28.3%. To avoid the installation of a new PSA, the network has been further optimized (HN5), resulting in $11.7 million per year savings and with an even shorter payback time of approximately 2 months. Five new compressors and 6 new pipes were installed (HN5, Figure 8a). Almost a 26.5% reduction in operating cost was achieved.
In this case, when rearranging the compressors respecting the nominal capacity, there is a reduction from 5 new compressors to only 2 new ones and using the 3 already installed. The rearrangement technique using virtual compressors applied to the compressors of example 2 can be seen in Figure 8b. In terms of total compressor cost reduces from $1.194 million ($0.575 million fixed cost and $0.620 million variable costs) to $0.934 million ($0.230 million fixed cost and $0.704 million variable costs). It represents a 21.8% reduction in the total investment in new compressors. It is worth mentioning that the cost of electricity increased from $ 0.765 to 0.777 million per year due to the pressure drop in the rearrangement. The total investment cost reduces from $1.406 to $1.099 million. The proposed network design through the rearrangement of the compressors is represented by HN5’, as shown in the Figure 8c.
To make a more direct comparison with the retrofit results obtained in the original paper, the existing network was tested using the HNS NLM model described in Section 3.4. The HNS NLM has 787 single equations, 430 single variables, and 249 discrete variables.
The cost of operation in the nonlinear (HN6) problem is 4.7% lower than in the HNS LM problem (HN5). However, it is observed that the most significant difference is the amount sent to burning as fuel. The optimization of HN6 network is an integer solution, which usually happens in nonlinear problems as it is not possible to guarantee optimum global optimization. The design obtained in HN6 optimized network through an HNS NLM model is shown in Figure 9a. It is remarkable to highlight that the HN6 network has 3 new lines. However, the cost of piping in this problem is calculated as a percentage of the cost of capital, which in this case, is zero. Therefore, it is necessary to estimate an average value for the cost of piping, which can be obtained with the number of lines in previous examples. The average cost for 3 new lines is between $0.08 and $0.1 million per year, taking into account that the fixed part is the predominant value and does not vary much with the flow.
Using the same methodology as in example 1, the network optimized through the HNS LM (HN5′) was used as an initial value to solve the nonlinear problem. The idea of using the result obtained in the linear model to initialize the nonlinear model guarantees an even more significant reduction in operating cost, of 13.9%, with zero capital cost (despite 3 new lines). The initialization of the rearranged network generates better results, in addition to a network with fewer connections. The design network HN7 is shown in Figure 9b.
Table 6 summarizes the principal results obtained through linear and nonlinear models for example 2. The lowest operating cost is obtained with the initialization of the linear model in the HNS NLM resolution (HN7), in addition to presenting the advantage of easier convergence.
In their original article, Sardashti Birjandi et al. [15] proposed the hydrogen network optimization through an MINLP model and obtained a 12% reduction in TAC. Considering TAC, this proposed HNS NLM model was able to reduce TAC by 30% (HN7), and the proposed HNS LM model was able to reduce TAC by 25.3%, which is a promising result. It is important to note that this case study was adapted from the example taken from the literature and that as many parameters are not described, the results would not be the same.
This example also shows that optimization through the linear model achieves considerable savings. Besides, as an initial guess for the nonlinear model, it proved to be an even more competitive alternative, further reducing operating costs.

5. Conclusions

In this work, an HNS LM (Mixed-Integer Linear Model) and HNS NLM (Mixed-Integer Nonlinear) optimization model is proposed for designing hydrogen networks for efficient use of this resource with cost reduction and environmental benefits.
The mathematical model is based upon superstructures, and it accounts for hydrogen sources, consumer units, purifying units, a fuel system, pressure constraints, and existing equipment and pipelines. The model can be used for grassroots designs and the retrofitting case. In the former, all the structure must be installed with an investment cost. In the later, the existing infrastructure is explored to reduce costs allowing the installation of new compressors, purifying units, and pipelines with an inherent investment cost. For both cases, the operating costs and the investment costs are the standard objective function to be minimized. Economic issues such as economy savings, maximum investment available, the payback time can be considered while delivering the optimal network design.
The model is thoroughly described, with all constraints, including the logical modeling equations used to accomplish design decisions and a proper estimation of costs, and all the model parameters. Initialization strategies for new design and retrofit cases were developed, which showed satisfactory results and efficiency for this work, both for existing and new networks.
The model was implemented in the modeling system GAMS solved with the solver CPLEX, and DICOPT and case studies from the literature were used to validate and explore the model features. For all examples, the proposed model was able to represent the existing networks as a feasible point, as well as to optimize them. Significant economic savings have been achieved when compared to existing networks, which shows that it is possible to work towards minimum hydrogen production and with investments payable in short periods.
The main breakthrough is the assumptions made in the mathematical modeling resulted in a linear model, which always converges to a global optimum, and it is speedy and robust. On the other hand, the drawback is that the solution may end up with a large number of compressor units. This issue can be overcome with the proposed algorithm basis evolution strategy to reduce the number of compressor units and pipelines and, therefore, the investment costs. This strategy has presented an excellent performance for the examples considered in this work. Besides, this technique can be extended to other problems of mass integration, such as pumps in water reuse, where the structure could also be represented through a linear model to facilitate resolution.
For comparison purposes, an HNS NLM model was also developed, in which streams can be mixed to be compressed at the same compressor unit. In this case, the number of compressors units is reduced when compared to the HNS LM model. However, the solution is influenced by the initial value, and it does not always converge, leading to a poor local minimum. The HNS NLM model also satisfies the needs of this work for the retrofit case and presented good results. However, the nonlinearity increases significantly the time need to solve the optimization problem. It is noteworthy that the HNS NLM model uses a superstructure that is different from the HNS LM, as the compressors are seen as a unit. The results obtained through nonlinear optimization compared to the linear ones, it has more flexibility of operation, because of the possibility of merging flowrates and share compressors. Resource time is not one of the main advantages when comparing linear with nonlinear. However, in the future, this work will be applied for multi-scenario optimization combined with production scheduling, so faster and more efficient resolution will be a critical issue.
For each case, different networks were proposed with different constraints. In general, the results were better than the original works of the case studies. Even though it was explored, the model versatility design networks allowing different constraints generating alternative designs according to the process requirements.
Different comparisons were made between the optimized networks in this work. With that, it can be concluded that the HNS LM model is satisfactory to optimize the hydrogen networks, even more with the rearrangement of the compressors, capable of reducing the investment costs. A reduction of 28% (example 1) and 26% (example 2) was obtained in the operating cost. In terms of the nonlinear model, the best results were obtained with the initiation of the network obtained from linear optimization. As a result, the operating cost was reduced by 31.2% (example 1) and 32.8% (example 2). This initialization technique was not found in the literature and proved to be an excellent tool for the optimization of hydrogen networks.
In this work, the importance of optimizing hydrogen networks is evident, aiming to minimize the operational cost. In addition, it is known that networks actually operate not only under nominal conditions as considered here, but also operate under different scenarios and different uncertainties. Since several factors affect this process, it is essential that the network must be able to work in various conditions. Therefore, the importance of working with uncertainties and multi-scenario optimization is evident. The MILP formulation proposed here can be easily extended to a multi-scenario version. In our future works, the uncertainty level will be addressed.

Author Contributions

Conceptualization, P.R.d.S.; Methodology, P.R.d.S., M.E.A., J.O.T. and L.F.T.; Software, P.R.d.S. and M.E.A.; Validation and Formal Analysis, P.R.d.S. and M.E.A.; Writing-Original Draft Preparation, P.R.d.S.; Writing-Review & Editing, P.R.d.S., M.E.A., J.O.T. and L.F.T.; Funding Acquisition, J.O.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by ANP/Petrobras.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

List of Symbols

i , j , k , c Sets of sources, consumers, purifiers, and compressors
F H 2 I i Flowrate of hydrogen sources
F H 2 I i ,   m a x ,   F H 2 I i ,   m i n Maximum and minimum flow rate of hydrogen sources
F I J i , j Flowrate from source to consumer
F I K i , k Flowrate from source to purifier
F I W i Flowrate from source to waste (fuel system)
F J j Total consumer flowrate
F K J k , j Flowrate from purifier to consumer
F J J j , j Flowrate from consumer j to consumer j′
Y J j Consumer purity
Y I i Source purity
Y K k Purifier purity
Y P j Purge purity of consumer
F P j Total purge consumer flowrate
F J W j Flowrate from consumer to waste (fuel system)
F J K j , k Flowrate from consumer to purifier
F P u r m a x , k Maximum capacity of the purifier
F K k Total flowrate in the purifier
F K W k Flowrate from purifier to waste (fuel system)
F K W r e c ,   k Purge flowrate from purifier to waste (fuel system)
Y K W k Purity of purge flowrate from the purifier
F , F α , β Flowrate
F m a x Maximum flowrate
ΕParameter associated with the existence of flowrate
z Binary associated with flowrate
z c Binary of a new compressor
u d e l t a P Binary of the pressure difference between the units
u c Parameter associated with existence compressor
z h Binary variable from a new pipeline
u h Parameter associated with existence pipeline
z k n Binary variable from the new purifier
α, βRepresents the possible connections involved
r e c k Purifier recovery
C o p e r a t i n g Operating cost
C H 2 I , C i Total and hydrogen production cost
C H 2 K , C k Total and purification cost
C H 2 C , C e l e t r i c Total and electricity cost
WPower compressor
w Intensive power compressor
C p ¯ Heat capacity
T Temperature
ƞCompressor efficiency
γ Cp/Cv Ratio
ρ o Density in standard condition
ρ Density
P o u t Outlet pressure
P i n Inlet pressure
C H 2 F T ,   C H 2 F ,   C f u e l Cost of burning purge as fuel
y Hydrogen fraction in the purge flow
Δ H ° H 2 ,   Δ H ° C H 4 Combustion heat of hydrogen and methane
C n e w   P S A , C n e w   P S A T Cost of new purifier
a P S A ,   b P S A Parameters of new purifier cost
C n e w   p i p i n g ,   C n e w   p i p i n g T Cost of new pipelines
ϑ Superficial gas velocity
L Distance
c ,   d Parameters of piping cost
C n e w   c o m p r e s s o r , C n e w   c o m p r e s s o r T Cost of a new compressor
a ,   b Parameters of new compressor cost
t Annual operating time
A f Annualized factor
C c a p i t a l Capital cost
f i   Interest rate
T A C Total annual cost
EEconomy
C O P a c t u a l , C O P n e w Actual and new operating cost
p t Payback
F C c Total compressor flow
F I C i , c Flow from source to compressor
F C J c , j Flow from the compressor to consumer
Y C c Purity in compressor
F J C j , c Flow from consumer to compressor
F C K c . k Flow from compressor to purifier
F K C k , c Flow from purifier to compressor
P C o u t , c Outlet pressure in the compressor
P C i n , c Inlet pressure in the compressor
P m i n Minimum pressure
P m a x Maximum pressure
P I i Source pressure
P K k Purifier pressure
P W Waste pressure
P J j Inlet consumers pressure
P P j Outlet consumers pressure

References

  1. IEA. Global Demand for Pure Hydrogen, 1975–2018. 2019. Available online: https://www.iea.org/data-and-statistics/charts/global-demand-for-pure-hydrogen-1975-2018 (accessed on 1 December 2019).
  2. Ceric, E. Crude Oil, Processes and Products, 1st ed.; IBC: Saravejo, Bosnia and Herzegovina, 2012. [Google Scholar]
  3. Figueiredo, E.A.H. Aplicação do Diagrama de Fontes de Hidrogênio em Refinarias de Petróleo. Master’s Thesis, Universidade Federal do Rio de Janeiro, Rio de Janeiro, Brazil, 2013. [Google Scholar]
  4. Silva, R.; Marvulle, V.C. Arte da tecnologia do hidrogênio: Review. In Encontro Energia no Meio Rural; AGRENER GD: Campinas, Brasil, 2006. [Google Scholar]
  5. Borges, J.L. Diagrama de Fontes de Hidrogênio. Master’s Thesis, Universidade Federal do Rio de Janeiro, Rio de Janeiro, Brazil, 2009. [Google Scholar]
  6. Jiao, Y.; Su, H.; Hou, W.; Liao, Z. A multiperiod optimization model for hydrogen system scheduling in refinery. Ind. Eng. Chem. Res. 2012, 51, 6085–6098. [Google Scholar] [CrossRef]
  7. Jia, N. Refinary Hydrogen Network Optimization with Improved Hydroprocesso Modelling. Ph.D. Thesis, University of Manchester, Manchester, UK, 2010. [Google Scholar]
  8. Marques, J.P.; Matos, H.A.; Oliveira, N.M.C.; Nunes, C.P. State-of-the-art review of targeting and design methodologies for hydrogen network synthesis. Int. J. Hydrog. Energy 2017, 42, 376–404. [Google Scholar] [CrossRef]
  9. Towler, G.P.; Mann, R.; Serriere, A.J.L.; Gabaude, C.M.D. Refinery hydrogen management: Cost analysis of chemically-integrated facilities. Ind. Eng. Chem. Res. 1996, 35, 2378–2388. [Google Scholar] [CrossRef]
  10. Fonseca, A.; Sá, V.; Bento, H.; Tavares, M.L.C.; Pinto, G.; Gomes, L.A.C.N. Hydrogen distribution network optimization: A refinery case study. J. Clean. Prod. 2008, 16, 1755–1763. [Google Scholar] [CrossRef]
  11. Hallale, N.; Liu, F. Refinery hydrogen management for clean fuels production. Adv. Environ. Res. 2001, 6, 81–98. [Google Scholar] [CrossRef]
  12. Liu, F.; Zhang, N. Strategy of purifier selection and integration in hydrogen networks. Chem. Eng. Res. Des. 2004, 82, 1315–1330. [Google Scholar] [CrossRef]
  13. Kumar, A.; Gautami, G.; Khanam, S. Hydrogen distribution in the refinery using mathematical modeling. Energy 2010, 35, 3763–3772. [Google Scholar] [CrossRef]
  14. Liao, Z.; Wang, J.; Yang, Y.; Rong, G. Integrating purifiers in refinery hydrogen networks: A retrofit case study. J. Clean. Prod. 2010, 18, 233–241. [Google Scholar] [CrossRef]
  15. Birjandi, M.R.S.; Shahraki, F.; Birjandi, M.S.; Nobandegani, M.S. Application of global optimization strategies to refinery hydrogen network. Int. J. Hydrog. Energy 2014, 39, 14503–14511. [Google Scholar] [CrossRef]
  16. Matijašević, L.; Petric, M. Integration of hydrogen systems in petroleum refinery. Chem. Biochem. Eng. Q. J. 2016, 30, 291–304. [Google Scholar] [CrossRef]
  17. Petric, M. Integracija Sustava Vodika u Procesima Prerade Nafte. Diploma Thesis, Sveučilište u Zagrebu Fakultet, Zagreb, Croatia, 2014. [Google Scholar]
  18. Georgiadis, M.C.; Schilling, G.; Rotstein, G.E. A general mathematical programming approach for process plant layout. Comput. Chem. Eng. 1999, 23, 823–840. [Google Scholar] [CrossRef]
  19. Grossmann, I.E.; Guillén-gosálbez, G. Scope for the application of mathematical programming techniques in the synthesis and planning of sustainable processes. Comput. Chem. Eng. 2010, 34, 1365–1376. [Google Scholar] [CrossRef] [Green Version]
Figure 1. (a) Scheme of the Superstructure developed for the Mixed-Integer Linear Programming (MILP) problem. (b) Scheme of the Superstructure developed for the Mixed-Integer Nonlinear Programming (MINLP) problem.
Figure 1. (a) Scheme of the Superstructure developed for the Mixed-Integer Linear Programming (MILP) problem. (b) Scheme of the Superstructure developed for the Mixed-Integer Nonlinear Programming (MINLP) problem.
Processes 08 01102 g001aProcesses 08 01102 g001b
Figure 2. Virtual Compressor Approach–Possibilities of mixing streams in the compressors.
Figure 2. Virtual Compressor Approach–Possibilities of mixing streams in the compressors.
Processes 08 01102 g002
Figure 3. Summary of the methodology proposed in this article, through optimization via linear and nonlinear model.
Figure 3. Summary of the methodology proposed in this article, through optimization via linear and nonlinear model.
Processes 08 01102 g003
Figure 4. Existing hydrogen network for Example 1.
Figure 4. Existing hydrogen network for Example 1.
Processes 08 01102 g004
Figure 5. (a) Optimized network HN1 via HNS LM, for Example 1. (b) Virtual compressor approach applied to HN1 network. (c) Optimized network HN1’ with rearranged compressors.
Figure 5. (a) Optimized network HN1 via HNS LM, for Example 1. (b) Virtual compressor approach applied to HN1 network. (c) Optimized network HN1’ with rearranged compressors.
Processes 08 01102 g005aProcesses 08 01102 g005b
Figure 6. (a) Optimized network HN2 via HNS NLM for example 1. (b) Optimized network HN3 via HNS NLM with HNS LM as initialization, for example 1.
Figure 6. (a) Optimized network HN2 via HNS NLM for example 1. (b) Optimized network HN3 via HNS NLM with HNS LM as initialization, for example 1.
Processes 08 01102 g006aProcesses 08 01102 g006b
Figure 7. Existing hydrogen network for Example 2.
Figure 7. Existing hydrogen network for Example 2.
Processes 08 01102 g007
Figure 8. (a) Optimized network HN5 via HNS LM for Example 2. (b) Virtual compressors applied to HN5 network. (c) Optimized network HN5’ with rearranged compressors.
Figure 8. (a) Optimized network HN5 via HNS LM for Example 2. (b) Virtual compressors applied to HN5 network. (c) Optimized network HN5’ with rearranged compressors.
Processes 08 01102 g008aProcesses 08 01102 g008b
Figure 9. (a) Optimized network HN6 via HNS NLM for Example 2. (b) Optimized network HN7 via HNS NLM with HNS LM as initialization for example 2.
Figure 9. (a) Optimized network HN6 via HNS NLM for Example 2. (b) Optimized network HN7 via HNS NLM with HNS LM as initialization for example 2.
Processes 08 01102 g009
Table 1. Parameters used to calculate the operating cost.
Table 1. Parameters used to calculate the operating cost.
Hydrogen cost—H2 plant C i 0.07 $/Nm³
Hydrogen cost—CCR C i 0.08 $/Nm³
Electricity cost C e l e t r i c 0.03 $/kWh
Purification cost C k 0.0011 $/Nm³
Fuel cost C f u e l 2.5 $/MMBtu
Table 2. Parameters used to calculate the capital cost [11].
Table 2. Parameters used to calculate the capital cost [11].
Cost of new compressors (k$) 115 + 1.91 × W
W in (kW)
Cost of new piping ($) ( 3.2 + 11.42 × D 2 ) × L D2 (in2) and L (m)
Cost of new PSA (k$) 503.8 + 347.4 × F
F in (MMscfd)
Table 3. Flowrate, purity, and pressure information used in Example 1.
Table 3. Flowrate, purity, and pressure information used in Example 1.
Sources F H 2 I i
(MMscfd)
F H 2 I i , m a x
(MMscfd)
Y I i % P I i
(psia)
H2 plant45.0050.0092.50300
CCR23.5023.5075.00300
Consumers F J j
(MMscfd)
Y J j % P J j
(psia)
F P j
(MMscfd)
Y P j % P P j
(psia)
HC38.7892.00200011.2975.001200
JHT8.6575.005004.3265.00350
CNHT8.2186.535003.4775.00350
DHT11.3175.976008.6170.00400
NHT12.0871.443006.5560.00200
IS40.0475.00300
Table 4. Results obtained in the different optimizations models for example 1.
Table 4. Results obtained in the different optimizations models for example 1.
COST (×106)
HNS LMHNS LMHNS NLMHNS LM INITIALIZATION- HNS NLM
HN1HN1′HN2HN3
H2 production ($/year)38.65938.65940.43941.117
Electricity ($/year)0.1360.1560.2040.198
Fuel ($/year)10.57610.57612.93114.448
Purification ($/year)0.4290.4290.4700.568
Operating cost ($/year)28.64828.66728.18327.435
New compressor ($/year)0.9920.135-0.290
New piping ($/year)0.4150.4050.4190.341
New PSA ($/year)6.8016.8017.4268.937
Capital cost ($/year)8.2097.3427.8469.568
Total capital cost ($)16.41814.68415.69219.136
TAC ($/year)36.85736.00936.02937.003
Economy ($/year)11.21411.19511.67912.427
Payback (year)1.4641.3121.3441.540
Resource time (s)0.0400.0401.3375.427
Table 5. Flowrate, purity and pressure information used in Example 2.
Table 5. Flowrate, purity and pressure information used in Example 2.
Sources F H 2 I i
(Nm³/h)
F H 2 I i , m a x
(Nm³/h)
Y I i
(%)
P I i
(bar)
H2 Plant40,50090,00076.0022
CCR59,00065,00092.004.50
PSA F P u r m a x , k
(Nm³/h)
Y K k Y K W k Rec
PSA I80,00099.9038.000.85
PSA II50,00099.9967.80
Consumers F J j
(Nm³/h)
Y J j % P J j
(bar)
F P j
(Nm³/h)
Y P j % P P j
(bar)
HC54,30099.9919810,00075.0029.50
DHT750092.0055270024.007.50
NHT150092.0055128062.0010.00
Table 6. Results obtained in the different optimizations models for example 2.
Table 6. Results obtained in the different optimizations models for example 2.
COST (×106)
HNS LMHNS LMHNS NLMHNS LM INITIALIZATION- HNS NLM
HN5HN5′HN6HN7
H2 production ($/year)46.15346.15346.69047.714
Electricity ($/year)0.7650.7770.7900.826
Fuel ($/year)15.33915.33917.39719.761
Purification ($/year)0.7620.7520.7230.803
Operating cost ($/year)32.33132.34330.80629.583
New compressor ($/year)0.5970.467--
New piping ($/year)0.1050.082--
New PSA ($/year)----
Capital cost ($/year)0.7030.549--
Total capital cost ($)1.4061.099--
TAC ($/year)33.03432.89230.80629.583
Economy ($/year)11.68511.67413.21114.434
Payback (year)0.1200.094--
Resource time (s)0.0670.0674.4951.906

Share and Cite

MDPI and ACS Style

Silva, P.R.d.; Aragão, M.E.; Trierweiler, J.O.; Trierweiler, L.F. MILP Formulation for Solving and Initializing MINLP Problems Applied to Retrofit and Synthesis of Hydrogen Networks. Processes 2020, 8, 1102. https://0-doi-org.brum.beds.ac.uk/10.3390/pr8091102

AMA Style

Silva PRd, Aragão ME, Trierweiler JO, Trierweiler LF. MILP Formulation for Solving and Initializing MINLP Problems Applied to Retrofit and Synthesis of Hydrogen Networks. Processes. 2020; 8(9):1102. https://0-doi-org.brum.beds.ac.uk/10.3390/pr8091102

Chicago/Turabian Style

Silva, Patrícia R. da, Marcelo E. Aragão, Jorge O. Trierweiler, and Luciane F. Trierweiler. 2020. "MILP Formulation for Solving and Initializing MINLP Problems Applied to Retrofit and Synthesis of Hydrogen Networks" Processes 8, no. 9: 1102. https://0-doi-org.brum.beds.ac.uk/10.3390/pr8091102

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop