A parallelized hybrid genetic algorithm with differential evolution for heat exchanger network retrofit

Graphical abstract


Method details
Heat exchanger network synthesis (HENS) is an important tool to design energy efficient production plants in process industry. Furman and Sahinidis [1] showed, that HENS is N P -hard in the strong sense . The heat exchanger network retrofit formulation, is an extension of HENS with increased complexity. Complexity is further increased by the additional dimension of operating periods, the possibility to integrate bypasses or admixers, and practical constraints. It is unlikely that deterministic algorithms can provide a feasible solution for such problems on a large-scale. Stochastic algorithms are able to handle such problems. Therefore, a hybrid two-level evolutionary algorithm based on genetic algorithm (GA) and differential evolution (DE) is developed. The concept of GA was first introduced by Holland [2] and further extended by Goldberg [3] . DE was introduced by Storn [4] , Storn and Price [5] . In the developed algorithm, a GA handles the topology optimization in the top level and a DE in the sub level is used to optimize the operation parameter.

Topology optimization using genetic algorithm
In order to define the topology of a HEN, integer values are used to address the position of each heat exchanger (HEX). Genetic algorithms, which are based on the evolution process in nature, are suitable to handle such discrete variables. Thereby, individuals within a population are compared and only the fittest of them survive. During the process, new individuals are generated through the mating of the fittest individuals. Coping errors (mutations) may occur during mating. An overview of the algorithm is shown in Fig. 1 . The following sections describe the evolutionary operators in detail.

Genetic algorithm -initialization
In a first step, a population P t = { EAM 0 . . . EAM NCT } with NCT random individual topologies, hereafter called chromosomes, is initialized. Each of these chromosomes represents a HEN topology and is described using an exchanger address matrix (EAM). An EAM of an example topology is shown in Fig. 2   of the HEX. During the initialization, for heach HEX it is first randomly decided if the HEX exists. For all existing HEXs, a random hot stream, cold stream, and enthalpy stage is defined. To ensure only feasible solution, the preserving constraint handling technique is used (e.g., for the stream number a random value between zero and the number of hot streams NH is selected). Bypass and admixer are not initialized, thus their existence depends on the heat loads. Hence, the need of a mixer and its configuration is determined in the DE.

Genetic algorithm -evaluation
After the initialization, each chromosome needs to be evaluated. In order to define the required areas and thus, the resulting cost, the heat loads need to be defined first. The optimal heat loads are found by the DE. To reduce computation time, the penalizing constraint handling technique is used. Thereby, only feasible topologies are optimized using DE. Infeasible solutions are evaluated by a penalty function which is less computationally expensive. The fitness f gt GA ,ct of the topology chromosome ct in the topology generation gt is given by whereby C GA is the set of all GA constraints. A constant value , which is larger than the initial total annual cost of the existing process, is added to the sum of the squared distances to the feasible region.

Genetic algorithm -selection
To choose the parent chromosomes for mating, tournament selection is performed. This means that the fittest out of a given number of randomly selected chromosomes is chosen.

Genetic algorithm -crossover
With a crossover probability of CR GA , the parents mate. Thereby, a one-point crossover is performed. Both parent chromosomes are cut below a random selected HEX number. The HEXs above the cut are swapped between the two chromosomes. The two resulting children chromosomes replace the parent chromosomes and are evaluated.

Genetic algorithm -mutation
In nature, during the mating process in the crossover, coping errors (mutations) may occur. In the algorithm, with a mutation probability of MT GA , for each allele (a single scalar value within a gene), a mutation occurs. For the mutation in the existence of a HEX gene (ex), a random bit flip is performed.
In the genes (i), (j), and (k), the value is changed within the upper and lower boundary (e.g., for (i) the number of existent hot streams) of the gene using an uniform distribution. It is ensured by using the preserving strategy, none of these values exceed their boundaries.

Genetic algorithm -next generation
After the application of the evolutionary operators (selection, crossover, and mutation), the new generated topologies are evaluated and the population is updated by replacing the parent chromosome with the new children chromosomes for the next generation. In order to keep track of the best solutions, a Hall of Fame (HoF) list is created. This list contains the current best solutions. If in the evaluation step a fitter solution is found, the HoF list is updated. With this additional feature, the flexibility of the algorithm in use is increased as the engineer is now able to choose between the most promising solutions. The algorithm is terminated as soon the maximal number of topology generations NGT is reached.

Heat load optimization using differential evolution
In contrast to the topology optimization, the heat loads are continuous variables with upper and lower bounds. Differential evolution algorithms are best to deal with the continuous variables. Compared to deterministic algorithms such as gradient descent, DE does not require the model to be differentiable. DE algorithms use the same concepts of evolution as GAs. However, the order of the evolutionary operators is reversed. There are different options for the DE configuration. In this case, the standard configuration DE/rand/1/bin is used. This means that the individuals for mutation are selected randomly, only one difference for perturbation ( F P : perturbation factor) is considered, and a binomial crossover is performed. In Fig. 3 , an overview of the algorithm is shown. The following sections describe the algorithm in detail.

Differential evolution -initialization
For each feasible topology from the GA, a population P h = { X 0 . . . X NCH } with NCH random heat load chromosomes is initialized. Each of these chromosomes consists of all the heat loads of each heat exchanger in every OP , resulting in a two-dimensional array. Chromosome ch is given by whereby the heat load is an array of the size N E × N OP . The initialization of the heat loads of existing HEXs ˙ Q op ch,e is constraint by a minimal user defined value and a maximal value given by the enthalpy differences of the connected streams (case study dependent value).

Differential evolution -evaluation
For the evaluation of the population of heat loads, the fitness of each chromosome needs to be determined. Infeasible solutions are evaluated by a penalty function. In each generation gh the fitness f gh DE ,ch is given by whereby, the T AC is the total annual cost of the current topology and heat loads. TAC consists the yearly operating cost and the investment cost for the retrofit. The detailed model for calculating the TAC is formulated in the corresponding research paper [6] . In order to minimize the TAC, the reciprocal value represents the fitness of a chromosome. Each distance of a heat load constraint violation is given by n gh DE , v iol,c . For infeasible operation parameter, the fitness is defined by reciprocal value of the penalty function, given by whereby C DE is the set of all DE constraints. A constant value which is larger than the initial TAC is added to the sum of the squared distances to the feasible region.

Differential evolution -mutation
The first evolutionary operation in a DE is mutation. Thereby, a three non-equal chromosomes X gh r1 , X gh r2 , X gh r3 ( r 1 = r 2 = r 3 ) of the current generation gh are selected randomly. A new donor chromosome is generated by whereby one difference for perturbation is used and weighted with the perturbation factor F P ∈ [ 0 , 2 ] . After the mutation, for all infeasible heat loads a new random value within the boundaries is assigned.

Differential evolution -crossover
In the crossover operator, a trial chromosome U whereby, r gh ∼ U (0 , 1) has a uniform distribution. With a probability of CR DE a crossover is performed.
To ensure at least one crossover per operating parameter, a random index within the chromosome is chosen for which crossover is always performed.

Differential evolution -selection
In the selection operator, the new created trial chromosome is evaluated. To determine the new target chromosome for the next generation gh + 1 with a simple greedy selection is performed.

Differential evolution -next generation
For each generation, the evolutionary operators (mutation, crossover, and selection) are executed till one of the termination criteria is fulfilled. The first termination criterion is satisfied when the maximal number of heat load generations NGH is reached. The second termination criterion is reached when the number of consecutive generations without improvement of the fitness exceeds its limit.

Implementation and parallelization
The algorithm is implemented in Python 3.8.2 using the library DEAP -Distributed Evolutionary Algorithms in Python [7] . Evolutionary algorithms are predestined for parallel computing as they work with populations of chromosomes which can be evaluated separately by distributing them among multiple processors. Therefore, all feasible GA chromosomes are distributed to multiple processors on which the heat loads are optimized using DE. The source code of the algorithm is published under an open-source license [8] .