Nature inspired methods for optimization
A Julia primer for process engineering

The problems considered in this book will all be formulated as the minimization of the desired objectives. There may be one or multiple objectives. In either case, the general formulation is

where \(d\) are the decision (or design) variables and \(\mathcal{D}\) is the search domain for these variables. \(\mathcal{D}\) might be, for instance, defined by a hyper-box in \(n\) dimensional space, \[ \mathcal{D} \equiv \left [ a, b \right ]^n \subset \mathbb{R}^n \] when only continuous real-valued variables are considered. More generally, the problem may include integer decision variables. Further, the constraints, \(g(d)\) and \(h(d)\), will constrain the feasible points within that domain.

The focus of this book is on solving problems that arise in process systems engineering (PSE). Optimization plays a crucial part in many PSE activities, including process design and operation. The problems that arise may have one or more of these challenges, in no particular order:

• nonlinear models
• multi-modal objective function
• nonsmooth and discontinuous
• distributed quantities
• differential equations
• small feasible regions
• multiple objectives

Each of these aspects, individually, can prove challenging for many optimization methods, especially those based on gradients to direct the search for good solutions. Some problems may have several of these characteristics.

Three chapters will present example design problems which exhibit some of the above aspects:

1. purification section for a chlorobenzene production process: discontinuous, non-smooth, small feasible region;
2. design of a plug flow reactor: continuous distributed design variable with behaviour described by differential equations;
3. heat exchanger network design: combinatorial, non-smooth, complex problem formulation.

These problems may be challenging for gradient based methods and therefore motivate the use of meta-heuristic methods. In this book, I will present meta-heuristic methods inspired by nature.

There are many optimization methods and different ways of classifying these methods. One classification used often is deterministic versus stochastic. The former include direct search methods [Kelley 1999] and what are often described as mathematical programming, such as the Simplex method [Dantzig 1982] for linear problems and many different methods for nonlinear problems [Floudas 1995]. The main advantage of deterministic methods is that they will obtain consistent results for the same starting conditions and may provide guarantees on the quality of the solution obtained and/or the performance of the method, such as speed of convergence to an optimal solution. If you are interested in using Julia for mathematical programming, the JuMP package is popular.

Stochastic methods, as the name implies, are based on random behaviour. The main implication is that the results obtained will vary from one attempt to another, subject to the random number generator used. However, the advantage of stochastic methods can be the ability to avoid getting stuck in local optima and potentially find the global optimum or optima for the problem.

Some popular stochastic methods include simulated annealing [Kirkpatrick et al. 1983] and genetic algorithms [Holland 1975] but many others exist [Lara-Montaño et al. 2022]. Simulated annealing, as the name implies, is based on the process of annealing, typically in the context of the controlled cooling of liquid metals to achieve specific properties (hardness, for instance). Genetic algorithms mimic Darwinian evolution, for instance, by using survival of the fittest to propagate good solutions and evolve better solutions. One common aspect of many stochastic methods is that they are inspired by nature.

In this book, we will consider three such nature inspired optimization methods: genetic algorithms, particle swarm optimization [Kennedy and Eberhart 1995], and a plant propagation algorithm [Salhi and Fraga 2011]. These are chosen somewhat arbitrarily as examples of stochastic methods. The selection is in no way intended to indicate that these are the best methods for any particular problem. The aim is to illustrate the potential of these methods, how they may be implemented, and how they may be used to solve problems that arise in process systems engineering.

This book includes all the code implementing the nature inspired optimization methods and most of the case studies, sufficient to allow readers to attempt the problems themselves.

Programming is at heart a practical art in which real things are built, and a real implementation thus has to exist. [Kay 1993]

The book has been written using literate programming [Knuth 1984]. The motivation for literate programming, to quote Donald Knuth, comes from:

[…] the time is ripe for significantly better documentation of programs, and that we can best achieve this by considering programs to be works of literature. [Knuth 1984]

Therefore, this book may be considered to be the documentation of the code used to implement the methods and to evaluate the problems defined in the case studies. The code presented in the book is automatically exported to the source code files, suitable for immediate invocation.

The technology supporting literate programming in this case is org mode (version 9.5.4). Org mode is a special mode in the Emacs editor (version 29.0.50.0). Org mode supports code blocks which may include programming code in a very wide range of languages. These code blocks are tangled into code files [Schulte and Davison 2011]. In the case of this book, all the code can be found in the author's github repository. As well as enabling tangling to create the source code files, org mode supports exporting documents to a variety of different targets, including PDF, epub, and HTML. The HTML version of the book is freely available.

Lastly, the literate programming support in org mode not only enables tangling, it also allows for code to be run directly within the Emacs editor while editing the document, with the results of the evaluation inserted into the document [Schulte and Davison 2011]. This book makes full use of this capability to process the results of the optimization problems, including extracting and generating statistical data with awk, grep, and similar tools, and plotting the outcomes with gnuplot.

Julia is a multi-purpose programming language with features ideally suited for writing generic optimization methods and numerical algorithms in general. To repeat the quote from the initial authors of the language,

Julia combines expertise from the diverse fields of computer science and computational science to create a new approach to numerical computing.[Bezanson et al. 2017]

For the purposes of this book, the key features of Julia are the following:

• dynamic typing for high level coding;

• multiple dispatch to create generic code that may be easily extended;

  

  [tweet]​

• functional programming elements for operating on data easily;

• integrated package management system to enable re-use and distribution; and,

• the REPL (read-evaluate-print loop) for exploring the language and testing code.

In this section, some of these aspects will be illustrated through code elements that will be used by all the nature inspired optimization methods presented in this book.

Julia version 1.7 has been used for all the codes in this book. I have tried to ensure that the style guide for writing Julia code has been followed throughout.

Different mathematical problem formulations of the same problem may have an impact on the solution process. Examples exist for the Simplex method [Hall and McKinnon 2004]. The choice of formulation may be a consideration for some problems. Further, for a given formulation, there may be alternative representations of the decision variables [Fraga et al. 2018] and these may affect the performance of individual optimization methods. Tailoring the formulation or representation to the method used may prove beneficial [Salhi and Vazquez-Rodriguez 2014].

For a given formulation, the representation should be chosen taking into account the following considerations:

1. The representation should be suitable for the optimization method or methods, enabling the implementation of the appropriate operations that are required by the methods. For instance, if the method is a genetic algorithm, the representation should be suitable for cross-over and mutation operators.
2. The representation should cover the complete space of feasible solutions to the problem or, at worst, ensure that the desirable potential solutions can be represented.
3. On the other hand, the representation should minimise the probability of the method's operators generating infeasible or undesirable solutions.

Multiple dispatch in Julia aids the process of considering different representations. For instance, different versions of the objective function implementation can be written that differ in the type of the first argument. This enables alternative representations without needing to change the implementation of the solution method. An example of this appeared recently [Fraga 2021c]. This latter example includes code in Julia using the Fresa implementation of a plant propagation algorithm [Fraga 2021b; Fraga 2021a] (see Section 2.5).

This book considers a selection of nature inspired optimization methods, specifically a selection of methods which are based on the evolution of a population of solutions. Three different solvers are described and fully implemented. These three have been selected as being representative of the class of evolutionary population based methods but the selection is not meant to be comprehensive, merely illustrative. The implementations presented in the book have been motivated by presentation and pedagogical advantages as opposed to aiming for the most efficient or feature-full implementations. References to production codes will be made where available. However, the codes presented are fully working implementations, as will be demonstrated through application to each of the case studies included in this book.

The three methods presented are

1. genetic algorithm, inspired by Darwinian evolution and survival of the fittest;
2. particle swarm optimization, inspired by flocks of birds and swarms of bees, looking for food or resting places; and,
3. plant propagation algorithm, inspired by the propagation method used by Strawberry plants.

All three methods fundamentally follow the same basic approach, illustrated in Algorithm 1.

All the methods presented in this chapter share some common functionality. This functionality is implemented in this section as a set of generic functions. These include fitness calculations for all the methods, including both single objective and multi-objective problems, selection procedures for those methods that use selection, and population management including function evaluation.

The codes in this section assume implicitly that the target optimization problem is one of minimization. If the problem considered were to be one of maximization, transformations of the objective function would be required. Many such transformations are possible including simply negating the value of the objective function.

Genetic algorithms are inspired by Darwinian evolution and, specifically, the survival of the fittest together genetic operators based on sexual reproduction and mutations [Holland 1975; Goldberg 1989; Goh et al. 2003; Katoch et al. 2021]. There are three key aspects that define a genetic algorithm:

1. selection, based on fitness;
2. crossover where new solutions are generated inspired by meiosis; and,
3. mutation where a new solution is generated by a small modification to an existing solution.

The basic algorithm is shown in Algorithm 2.

As presented, the algorithm does not include elitism. Elitism is a means of ensuring that a very good solution identified during the evolutionary procedure is not lost due to selection. In the implementation below, the best solution in any generation will automatically be included in the next population. This is known as elitism of size 1. For multi-objective optimization problems, the best solution will be one chosen arbitrarily from the set of non-dominated solutions, typically the solution that is the best with respect to one of the criteria. Other forms of elitism, especially for multi-objective problems, are possible and often desirable.

There are many variations on this algorithm, such as whether mutation should be applied potentially to the new solutions obtained by crossover and whether there should be a form of elitism where the best solutions in a population are guaranteed to be in the new population. Neither of those aspects are included in this algorithm. The aim of this book is not to explore the impact of the many variations possible but to illustrate the commonalities and differences between different nature inspired evolutionary population based methods.

The stopping criterion for this method, as for the others presented in the next chapters, will be the number of function evaluations. Again, many different stopping criteria could be considered including the number of cycles or generations to perform, the total time elapsed, and whether the population is stable in the sense that there is no longer any improvement in the best solution or solutions. As one of the case studies has significant computational requirements, the chlorobenzene process design, the number of function evaluations is an appropriate stopping criterion. Although it does not take into account the overhead of the actual method, the three methods considered in this book will have similar overhead.

From the algorithm, the outcome will be based on the implementation of the select, crossover, and mutate functions as well as the values of the parameters: the population size, \(n_p\), the crossover rate, \(r_c\), and the mutation rate, \(r_m\). These functions and parameters will typically require adaptation to the particular problem. The select method, however, will be generic and has already been defined above in Section 2.2.4. The crossover and mutate methods need to be based on the specific representation of solutions so may need tailored implementations.

The original method was developed by Kennedy & Eberhart [Kennedy and Eberhart 1995]. The basic idea is that each particle (a solution to the optimization problem) in a swarm has a velocity, allowing it to move through the search space. This velocity is adjusted each iteration by taking into account each particle's own history, specifically the best point it has encountered to date, and the swarm's overall history, the best point encountered by all particles in the swarm. See Algorithm 1.

The equation for updating the velocity of each particle, \(i\), is

where \(x_i\) is the current position of the particle, \(v_i\) its velocity, \(h_i\) the best place this particle has been (its historical best), and \(b\) the best point found by the whole population so far. The ω parameter provides inertia for non-zero values and the parameters \(c_1\) and \(c_2\) weigh the contributions from local history and global history towards the change in velocity.

Due to the calculations for updates to velocities and the use of the velocities to move the particles, this method is naturally directed at problems where the decision variables are real numbers. Although the method has been adapted to integer problems [Datta and Figueira 2011], the implementation here is for real-valued vectors as decision variables.

The plant propagation algorithm (PPA) [Salhi and Fraga 2011] is inspired by one of the methods used by strawberry plants to propagate: the generation of runners. In fertile soil, strawberry plants send of a number of short runners to exploit the current surroundings. In less fertile soil, fewer yet longer runners propagate to explore more broadly. The PPA therefore creates new solutions, based on existing solutions in a population, by generating a number of runners. The number of runners is proportional to the fitness and the distance these runners propagate is inversely proportional to the fitness. This is similar in motivation as an adaptive scheme for mutation in genetic algorithms [Marsili Libelli and Alba 2000].

The PPA is shown in Algorithm 1.

The plant propagation algorithm depends on three support functions: fitness, select, and neighbour. The first two have already been implemented as generic functions suitable for any method. The last, however, is specific to the plant propagation algorithm.

It is useful to have a control experiment to confirm that any or all of the above methods are indeed effective at search. Whether a fully random search could be considered a method inspired by nature or not, such a method is in the spirit of the other nature inspired methods in being stochastic. There are no parameters and the method relies on a function to generate new random solutions in the search space.

The actual search generates \(n_f\) random points, and to provide some statistics for comparison with the evolutionary methods, output is generated periodically with the best solution found so far.

 1: using NatureInspiredOptimization: Point, statusoutput
 2: function randomsearch(
 3:     # required arguments, in order
 4:     f;                    # objective function
 5:     # optional arguments, in any order
 6:     lower = nothing,      # bounds for search domain
 7:     upper = nothing,      # bounds for search domain
 8:     parameters = nothing, # for objective function 
 9:     nfmax = 10000,        # max number of evaluations
10:     output = true)        # output during evolution
11: 
12:     lastmag = 0           # for status output 
13:     lastnf = 0            # also for status output
14: 
15:     best = Point(randompoint(lower,upper), f,
16:                  parameters)
17:     nf = 1                # function evaluations
18:     while nf < nfmax
19:         p = Point(randompoint(lower,upper),
20:                   f, parameters)
21:         nf += 1
22:         if p ≻ best
23:             best = p
24:         end
25:         lastmag, lastnf = statusoutput(output, nf,
26:                                        best,
27:                                        lastmag,
28:                                        lastnf)
29:     end
30:     lastmag, lastnf = statusoutput(output, nf,
31:                                    best,
32:                                    lastmag, lastnf)
33:     best 
34: end

A simple benchmark problem is now present to illustrate how the nature inspired methods may be used. The benchmark problem is

where \(d\) is dimension of \(x\) and β typically 0.5. This is a continuous function defined on the full domain. It is smooth. As such, it is not the type of problem that requires black box optimization methods. However, this benchmark problem has been chosen to demonstrate the Julia code required to define a problem and solve it using the different solvers described previously.

The first step is to define a function that implements the objective function in Equation \eqref{orga2aeda6}:

1: function permdbeta(x, β = 0.5)
2:     d = length(x)
3:     z = sum( sum( (j^i + β)*((x[j]/j)^i - 1) for j ∈ 1:d )^2 for i ∈ 1:d )
4:     (z, 0)
5: end

This code illustrates some elegant programming elements of Julia for working with vectors. The summation signs in Eq. \eqref{orga2aeda6} translate directly into Julia: line 3. For instance, \(\sum_{i=1}^n x_i\) may be written as

sum( x[i] for i in 1:n )

in Julia.

The first line illustrates the use of default values for positional arguments, as opposed to optional arguments, when those arguments are not specified when the function is invoked. In this case, if only one argument is given, i.e. x, the value of β defaults to 0.5. This is the value typically used as a benchmark.

The code also highlights, line 4, that the objective function must return a tuple consisting of the actual value of the objective function, z, along with an indication of constraint satisfaction: 0 or less for a feasible point, a value greater than 0 for infeasible with the actual value ideally a measure of how far from feasible the point may be. For this example, as the problem is feasible over the full domain, the second element of the tuple is always zero.

For all the methods, the following initialization is common:

1: using NatureInspiredOptimization: Point
2: d = 4                           # dimension
3: x0 = zeros(d)                   # initial guess
4: a = -d * ones(d)                # domain
5: b =  d * ones(d)                # domain
6: p0 = [Point(x0,permdbeta)]      # initial population

where the dimension of the problem is set to 4, an initial point at the origin is defined, and the lower and upper bounds of the search domain defined. The initial population, p0, consists of the initial point evaluated with the objective function.

The GA and PSO methods require an initial population that is of the size they expect to evolve. The PPA can start with a population of size 1. However, for consistency, all methods will be started with the same size initial population:

1: using NatureInspiredOptimization: randompoint
2: n = 40                          # population size
3: while length(p0) < n
4:     push!(p0, Point(randompoint(a, b), permdbeta))
5: end

Given the random initial population, we can now apply each of the methods. First, the genetic algorithm:

1: using NatureInspiredOptimization.GA: ga
2: gabest, gapop = ga(p0, permdbeta;
3:                    np = n,
4:                    lower = a,
5:                    upper = b,
6:                    nfmax = 10_000_000)

Note that use of _ to denote separation of digits for large numbers, in line 2. 10_000_000 is much easier to read than 10000000 although, of course, we could have used exponential notation, 1e7.

The particle swarm optimization method is invoked similarly:

1: using NatureInspiredOptimization.PSO: pso
2: psobest = pso(p0, permdbeta;
3:               lower = a,
4:               upper = b,
5:               np = n,
6:               nfmax = 10_000_000)

Next, the plant propagation algorithm which differs in that the population size, n, is smaller than that used for the other methods. This is partly because the PPA generates multiple new solutions from each solution chosen for propagation, on average 2.5-3 for each, and also because of the results of the analysis we present in the next section of this book:

1: n = 10                          # population size
2: using NatureInspiredOptimization.PPA: ppa
3: ppabest, ppapop = ppa(p0, permdbeta;
4:                       np = n,
5:                       lower = a,
6:                       upper = b,
7:                       nfmax = 10_000_000)

Finally, as the control for the experiment, the random search method:

1: using NatureInspiredOptimization: randomsearch
2: rsbest = randomsearch(permdbeta;
3:                       lower = a,
4:                       upper = b,
5:                       nfmax = 10_000_000)

The results from the above code are summarised in Figure 1. This figure shows the evolution of the best solution in the population. The stochastic nature of the methods means that every time one of the methods is applied, the results will differ. When we present a number of case studies below, this stochastic nature will be taken into account by solving each problem multiple times and performing some statistical analysis. For this illustrative example, however, only one attempt has been considered for each method, simply to illustrate how the methods are used and what results may look like. No conclusions about the efficacy of any of the methods should be drawn from this one experiment. The focus of this book is not on benchmark problems such as Equation \eqref{orga2aeda6} but on problems which arise in process systems engineering, as we shall see in the next chapter.

The case study is the optimisation of a purification section in the production of chlorobenzene. In the overall process, large quantities of benzene are used. Due to the partial conversion in the reaction section of the process, significant quantities of un-reacted benzene would be wasted if not recycled. To recycle the benzene, a stream consisting of primarily benzene needs to be further purified to ensure that the benzene sent back upstream in the process is pure enough to not affect the reaction. The case study we consider here is the design of this purification stage of the overall process.

The stage comprises three distillation units with a feed stream described in Table 1. Although the main aim is to purify the benzene recycle stream, the feed stream to this stage also contains some amount of the main product of the process, chlorobenzene. It is desirable to separate this product as well in sufficiently pure form and to avoid unnecessary loss. The process structure we consider is presented in Figure 1.

To assess the different alternatives, the economic criteria of capital and operating costs are used. Process models are required to determine the impact of design decisions on these criteria. For this problem, the process models, in Jacaranda, are based on the Fenske, Underwood, and Gilliland short-cut models and correlations often used for multi-component distillation column design [Rathore et al. 1974; Monroy-Loperena and Vacahern 2012].

These models are based on the concept of light (most volatile component in the separation) and heavy (least volatile of the components) key species. Essentially, if the species in a mixture are sorted according to their boiling points, the mixture can be separated between any two adjacent species. For instance, if there were a mixture with three species, A, B, and C, with increasing boiling points from A to C, there would be two possible separations for this mixture. The first would be between A and B, yielding two streams, one with almost all of the A but a little B and one with a little of the A, most of the B, and essentially all of the C. The second separation possible would be between B and C: one of the resulting streams from this separation would have all the A, most of the B and a little C and the other stream would have a little B and most of the C.

The Fenske equation determines the minimum number of stages, \(N_{\min}\), a distillation column will require to achieve a given separation identified by the recovery of the light and the heavy:

\[ N_{\min} = \frac {\log \frac {x_{D,l}x_{B,h}} {x_{B,l}x_{D,h}}} {\log \alpha_{l,h}} \]

where \(D\) refers to the distillate or tops product of the distillation unit and \(B\) to the bottoms product. \(x\) is the molar fraction of each species in the streams, \(\alpha_{l,h}\) is the relative volatility of the light key to the heavy key. Volatility is a measure of the vapour pressure, a property of each species, with higher values of volatility associated with the lighter species. The relative volatility of a pair of species is the ratio of their volatilities.

The compositions of the light and heavy keys in the product streams are a function of the recovery design variable. The recovery specifies the fraction of the key species that ends up in the desired output stream. The higher the recovery, the smaller the amount of each key that goes out with the other stream. Higher recovery is required for higher purity of the final product streams but there is direct relationship between the recovery and the size of the distillation column and hence its cost, both capital and operating.

The Underwood equation determines the minimum amount of reflux, the liquid sent back into the column from the top stream to aid in achieving the separation:

\[ R_{\min} = \sum_{i=1}^{n_{c}} \frac {\alpha_{i} x_{D,i}} {\alpha_{i}-\theta} \]

where \(\theta\) is the solution to the nonlinear equation

\[ \sum_{i=1}^{n_{c}} \frac {\alpha_{i}x_{F,i}} {\alpha_{i} - \theta} = 1 - q\]

where \(F\) indicates the feed stream, \(n_{c}\) is the number of species involved in the feed stream and \(q\) is the state of the feed stream. The state depends on the energy in the stream and that is a function of the temperature and the pressure of the stream.

The actual number of stages, above the minimum \(N_{\min}\), and the actual reflux, also above the minimum \(R_{\min}\), are correlated through the Gilliland equation. This is also a nonlinear equation and depends on another design variable, the reflux rate factor. The higher the factor, the larger the actual reflux rate used. The reflux rate affects the operating costs through the use of cooling and heating utilities: larger reflux means increased utility use. It also influences the capital cost: higher reflux leads to a wider column but one with less stages, usually leading to a lower capital cost overall.

In all of the above equations, the calculation of variables such as \(\alpha\) and \(q\) require the estimation of physical properties. Physical property models are typically nonlinear. For this problem, we have used the Antoine equation to determine vapour pressures, \(p^{*}\) used to calculate the relative volatilities, \(\alpha\), of the species:

\[ \log p^{*} = A - \frac {B} {T - C} \]

where A, B and C are the Antoine coefficients which can be found in reference books for many species. The \(log\) function used will depend on the source of the Antoine coefficients and may be base 10 or natural.

The relative volatilities of the species are the ratio of their vapour pressure to one reference species in the mixture. The temperatures will depend on the operating pressure, the third design variable for the distillation unit. Higher pressures will results in higher temperatures and smaller relative volatilities. The former will affect the choice of heating and cooling utilities; the latter will affect the number of stages required for the separation required.

The process design has several constraints beyond the physical models described above. Firstly, the benzene and chlorobenzene product streams must meet purity requirements. Secondly, there is a finite number of utilities available for meeting the heating and cooling requirements. For some combinations of operating pressure and stream compositions, utilities may not be available to meet the requirements imposed by the design of the distillation units. Finally, there is a limit on how much of the main product, chlorobenzene, can be lost through the bottom product stream of the third distillation unit. The constraints define a feasible region that is small, well under 1% of the box domain defined by the design variable bounds. The small feasible region together with the nonlinear models can pose difficulties to optimization methods.

A more detailed description of the full model used is given in [Zilinskas et al. 2006].

The Jacaranda system [Fraga et al. 2000] is an extensible object oriented framework for process design and optimization written in Java. It implements its own modelling suitable for the definition of complex process flowsheets with detailed processing step models. The process described in the previous section can be implemented in Jacaranda.

Jacaranda is available for download, with links and instructions for licensing detailed in the guide. Once downloaded, as a jar file, the CLASSPATH environment variable must be set to point to the jacaranda.jar file.

A package, Jacaranda.jl, is available to access Jacaranda and the objects created by that system. To add this package to your Julia configuration, simply invoke Julia in a command line or shell window:

$ julia
julia> ]
(@v1.7) pkg> add https://gitlab.com/ericsfraga/Jacaranda.jl
[...]
(@v1.7) pkg> BACKSPACE
julia> ^d
$

where the lines show the prompt you should be expecting at each point. The elided [...] output should show a number of packages being installed, those that Jacaranda.jl depends on. Once this has been done, you should be able to type the following in Julia:

1: using Jacaranda

In this section, we define the process model for the chlorobenzene recovery process. The model is written in the Jacaranda system modelling language. Each snippet of the model described here is tangled to an input file, clben-process-flowsheet.in, which will be used as input to the Jacaranda system. The structure of the process is described in section 3.1. It consists of a feed, two distillation columns, and three output streams. The Jacaranda model also defines an objective function based on specific defined criteria for optimization and the measure of constraint violations for infeasible solutions. This objective function provides the information necessary to be used by the nature inspired optimization methods developed in the previous chapter.

\sloppy The first part of the process model specification is administrative, asking for output from Jacaranda itself to be redirected to the specified file, clben-process-flowsheet.out and giving the project a name, clben-withconstraints. Two namespaces in the Jacaranda object framework are also specified. This allows us to write, for instance, Criteria instead of the fully qualified jacaranda.design.ps.Criteria. The first namespace is for process synthesis related objects and the second is for processing unit models.

output clben-process-flowsheet.out
project clben-withconstraints 
use jacaranda.design.ps
use jacaranda.design.units

The next part of the process model specification is the definition of constants used by the remainder of the model definition. Also, some variables used by Jacaranda's own optimization subsystem are defined. Jacaranda was originally designed to implement an implicit search procedure for finding optimal process designs, using discretization and dynamic programming [Fraga 1996; Fraga and McKinnon 1994]. Although these settings are not relevant for the use in this book, they may be useful for comparing results with other published results [Zilinskas et al. 2006].

variables
 const Base 0.001 "Base flow discretization"
 # flows of species in feed stream
 const Fbenzene  0.97 "Flow of Benzene"
 const Fc6h5cl   0.01 "Flow of C6H5Cl"
 const Fpc6h4cl2 0.01 "Flow of p-C6H4Cl2"
 const Fc6h3cl3  0.01 "Flow of C6H3Cl3"
 # set the cost base to 1993 based on the M&S index
 set msindex 964.2
 const nPs 4 "Number of discrete enthalpy levels"
end

For this process, we want to investigate the difference in operating versus capital costs so we define two criteria, each of which is the sum of the relevant criterion value generated by each unit in a flowsheet. The term in quotes for each criterion is the expression to evaluate for each unit. The sum means that the value of the criterion is the sum of expression evaluate for all the units in the resulting flowsheets.

Some criteria have been included which in fact are constraints on aspects of the process. The values of compositionspec (composition purity requirements) and flowspec (maximum flow of contaminants allowed) represent the maximum violation of the inequality constraints for product specifications. The value of UtilitySpec represents the temperature violation for heating and cooling utility requirements. All of these should be less than or equal to 0 for feasible solutions. Positive values give a measure of the infeasibility which can be used in assigning fitness to solutions in a population (see Section 2.2.3).

Criteria costs
 criterion CapitalCost sum "capcost"
 criterion OperatingCost sum "opercost"
 criterion CompositionSpec max "compositionspec"
 criterion FlowSpec max "flowspec"
 criterion UtilitySpec max "Tviolation"
end

Next we define the heat exchanger model for heating and cooling. There are two cost functions, one for the capital cost as a function of the area of the exchanger and the other the operating cost as a function of the duty (Q [kJ]), the cost of the utility used (Cu [$/kJ]), and the time (seconds in a year equal to h*hpy where hpy is the operating hours per year). To cater for the case where the temperature differences at both ends of the exchanger are the same, we use the Chen approximation [Chen 1987] to avoid numerical issues in the log function.

HeatExchanger hex
 model
  deltain = abs(T1in-T2out)
  deltaout = abs(T1out-T2in)
  eps = abs(deltain/deltaout-1)
  if eps < 4e-4
     lmtd = (1+eps/2-eps^2/12)*deltaout
  else
     lmtd = (deltain - deltaout)/log(deltain/deltaout)
  endif
  A = Q/U/lmtd
  capcost = 101964 + 92.9*A^0.65
  opercost = Cu*Q*h*hpy
 end
end

Utilities are specified which may be used to meet the heating and cooling requirements of the processing units. The parameters for costing the use of a utility are the inlet temperatures of both the process stream and the utility itself, the outlet temperatures, the overall heat transfer coefficient for the exchange, and the cost of the utility per unit of heat. Each utility definition therefore includes its inlet and outlet temperatures, the heat transfer coefficient, and the cost.

Along with the actual utilities, we define two special utilities solely to provide penalty function information when none of the normal utilities is sufficient to meet the requirements of the processing unit. These special utilities will give a positive value for Tviolation should none of the actual heating or cooling utilities be suitable. The Tviolation will be used as one of the criteria defined above and will therefore be available for the optimization system to assign fitness values to infeasible process designs.

DiscreteUtilities utils
  hot "HotPenalty"
    Tin 1000
    Tout 1000
    htc "5000 *W/m^2/K"
    cost "1.0246 /GJ"
    model
      Tviolation = T1out - 503.5
    end
  end
  hot  "Steam @28.23 atm"      503.5 503.5   "5000 *W/m^2/K" "1.0246 / GJ"
  hot  "Steam @11.22atm"       457.6 457.6   "5000 *W/m^2/K" "0.773824 / GJ"
  hot  "Steam @ 4.08atm"       417.0 417.0   "5000 *W/m^2/K" "0.573203 / GJ"
  hot  "Steam @ 1.70atm"       388.2 388.2   "5000 *W/m^2/K" "0.41796 / GJ"
  cold  "ColdWater @ 32.2degC" 305.2 305.2   " 500 *W/m^2/K" "0.0668737 / GJ"
  cold  "Ammonia @   1degC"    274.00 274.00 " 500 *W/m^2/K" "1.65035 / GJ"
  cold  "Ammonia @ -17.68degC" 255.32 255.32 " 500 *W/m^2/K" "2.96871 / GJ"
  cold  "Ammonia @ -21.67degC" 251.33 251.33 " 500 *W/m^2/K" "3.96226 / GJ"
  cold "ColdPenalty"
    Tin 0
    Tout 0
    htc "500 *W/m^2/K"
    cost "3.96226 /GJ"
    model
      Tviolation = 251.33-T1out
    end
  end
end

Note that the utilities have been defined in decreasing order of temperature.

All the data we need for this problem have been defined. Now set the global data variables that are used by the process evaluation procedure:

Data
 criteria costs
 hex hex
 utils utils
end

Given these data, we can now define the actual process flowsheet. This process consists of three distillation columns and four specific product tanks, as well as a feed tank (Figure 1).

The first unit is the feed tank. Associated with this tank is the actual feed stream. The feed stream consists of the output of the reactor and there are 4 different components (molecular species). All the components of the stream are defined using the KistaComponent class. This component class uses physical properties estimation methods [Ponton 1989; Boston and Britt 1978] based on the 1 atm boiling point of the component and the coefficients of the vapour heat capacity polynomial correlation.

vle.KistaComponent benzene
 bp 353.2
 cpv -33.917 4.743e-1 -3.017e-4 7.130e-8
 lhtc " 500 *W/m^2/K"
 vhtc "5000 *W/m^2/K"
end

vle.KistaComponent C6H5Cl
 bp 405.2
 cpv -33.888 5.631e-01 -4.521e-04 1.426e-07
 lhtc " 500 *W/m^2/K"
 vhtc "5000 *W/m^2/K"
end

vle.KistaComponent p-C6H4Cl2
 bp 447.3
 cpv  -14.344 5.534e-01 -4.559e-04 1.447e-07
 lhtc " 500 *W/m^2/K"
 vhtc "5000 *W/m^2/K"
end

vle.KistaComponent C6H3Cl3
 bp 486.2
 cpv  -14.361 6.087e-01 -5.622e-04 2.072e-07
 lhtc " 500 *W/m^2/K"
 vhtc "5000 *W/m^2/K"
end

Using these components, we can now define the feed stream. This stream consists of a single phase at specified temperature and pressure conditions. To define the stream, we first define the Phase object with the correct amounts of each component.

In defining the Phase object, we also specify the granularity of the component flows for each component. This is the basic discretization used by Jacaranda to ensure that the search graph generated implicitly is finite in size. The granularity should be appropriate for the particular problem and, in particular, should be consistent with discretizations performed in the separation units (for example) and the product specifications. In this problem, we have specified the same discretization factor (Base as defined above) for each component as a function of its flow in the feed stream. The individual flows of the components in the feed stream have been defined in the main input file as well. However, it should be noted that the discretization settings are not used in the case study presented below. They are here for completeness should an interested reader wish to experiment with Jacaranda itself directly.

vle.Phase feedPhase
 add benzene Fbenzene
 add C6H5Cl Fc6h5cl
 add p-C6H4Cl2 Fpc6h4cl2
 add C6H3Cl3 Fc6h3cl3

 base benzene   "Base * Fbenzene"
 base C6H5Cl    "Base * Fc6h5cl"
 base p-C6H4Cl2 "Base * Fpc6h4cl2"
 base C6H3Cl3   "Base * Fc6h3cl3"

 show
end

The stream definition includes the phase object defined above and the specification of the state of the stream. The state is defined by the initial temperature and pressure and the number of discrete pressure levels which can be represented. This number, nPs, is defined above.

vle.Stream feed
 add feedPhase
 nstates nPs
 T 313
 P "1 *atm"
end

Given the feed stream, we can now define the feed tank unit: this is a unit with no input streams and one output stream, the feed to the process.

vle.FeedTank feedtank-1
 Stream feed feed
end

The feed stream is sent to a distillation column, dist-2. This unit will separate the mixture in the feed stream, using benzene as the light key. The distillation model defines a number of design variables that will be used by the optimization procedures: rrf, the reflux rate factor, rec, the recovery of the light and heavy keys, and P, the operating pressure of the column:

vle.Distillation dist-2
  Component key benzene
  Real rrf 1.2 1.05 1.5 10
  Real trayefficiency 1.0 1.0 1.0 1
  Boolean eduljee false
  Boolean conceptual false
  Boolean sfa false
  Boolean totalcondenser true
  new tunable+specification "Recovery of light and heavy keys" Real rec 0.95 0.800 0.999 10
  LogReal P 5.039684199579494 1.0 32.0 16
end

The top product of the first distillation will go into a product tank, intended to receive pure benzene. The benzene is eventually destined to be recycled back to the reaction section in the upstream process. The recycle is not considered in this problem. The specification is that the benzene be at least 98 mole% pure (compositionspec) and for there to be no more than 0.005 kmol s^-1 of chlorobenzene (flowspec).

vle.ProductTank pureBenzene-3
  Expression spec "1"
  model
    compositionspec = 0.98-x[benzene]
    flowspec = F*x[C6H5Cl] - 0.005
  end
end

The second distillation column, which processes the bottom product of the first distillation unit, also has benzene as the light key and the same three design variables for optimization:

vle.Distillation dist-4
  Component key benzene
  Real rrf 1.2 1.05 1.5 10
  Real trayefficiency 1.0 1.0 1.0 1
  Boolean eduljee false
  Boolean conceptual false
  Boolean sfa false
  Boolean totalcondenser true
  new tunable+specification "Recovery of light and heavy keys" Real rec 0.999 0.800 0.999 10
  LogReal P 32.0 1.0 32.0 16
end

The top product stream, intended to be pure benzene, is sent to another product tank that has the same composition and flow specifications: 98 mole% benzene with no more than 0.005 kmol s^-1 of chlorobenzene:

vle.ProductTank pureBenzene-5
  Expression spec "1"   # always true
  model
    compositionspec = 0.98-x[benzene]
    flowspec = F*x[C6H5Cl] - 0.005
  end
end

The third and final distillation unit separates chlorobenzene from the waste products.

vle.Distillation dist-6
  Component key C6H5Cl
  Real rrf 1.2 1.05 1.5 1
  Real trayefficiency 1.0 1.0 1.0 1
  Boolean eduljee false
  Boolean conceptual false
  Boolean sfa false
  Boolean totalcondenser true
  new tunable+specification "Recovery of light and heavy keys" Real rec 0.95 0.800 0.999 10
  LogReal P 1.5874010519681996 1.0 32.0 16
end

The top product from this distillation unit is sent to a product tank with the specifications for the chlorobenzene, specifically a purity of at least 90 mole% (compositionspec).

vle.ProductTank pureC6H5Cl-7
  Expression spec "1"  # always true
  model
    compositionspec = 0.90-x[C6H5Cl]
  end
end

The bottom product from the third distillation unit is sent to another product tank. This tank has a constraint that there be no more than 10 mole% benzene or 10 mole% chlorobenzene.

vle.ProductTank otherProducts-8
  Expression spec "1"  # always true
  model
    compositionspec = max({x[benzene]-0.10, x[C6H5Cl]-0.10})
  end
end

Having defined all the processing steps in the process, including feed and product tanks and the actual processing steps, we define the actual Flowsheet object which consists of the units and their connections:

Flowsheet clben
  build
    feedtank-1 dist-2
    dist-2 pureBenzene-3 dist-4
    dist-4 pureBenzene-5 dist-6
    dist-6 pureC6H5Cl-7 otherProducts-8
  end
  continuous
  evaluationmethod breadthfirst
end

The continuous specification means that the discretization parameters defined in this model are not taken into account. Jacaranda's Flowsheet object can model the process in both continuous and discrete spaces. For the purposes of this book, only the continuous modelling is considered.

This completes the definition of the clben process flowsheet in Jacaranda. In the next section, the model is accessed through the Jacaranda-Julia interface, retrieving the clben object to be able to evaluate the process with different values for the design variables.

The first step in defining the objective function for the chlorobenzene process design problem is to use Jacaranda to load the process flowsheet object, defined in the previous section:

1: function process()
2:     Jacaranda.init("clben-process-flowsheet.in", "clben")
3: end

The Jacaranda.init function loads the input file given as the first argument. The second argument is the name of the Jacaranda variable that implements the ObjectiveFunction interface (in the jacaranda.util.opt namespace). The Flowsheet object implements that interface. The Jacaranda.init method returns the following 5-tuple:

(flowsheet, x0, lower, upper, nf) 

where flowsheet is a handle on the Jacaranda process flowsheet object. The remaining entries in the tuple relate to the flowsheet itself. x0 is an initial guess for the design variables. These initial values are those specified in the definition of the flowsheet and its processing units in the Jacaranda input file. lower is the vector of lower bounds on the design variables, upper the vector of upper bounds, and nf the number of objective function values that will be returned when the process flowsheet is evaluated.

The design variables for this process are the design variables for each distillation unit: the operating pressure, P \(\in [1, 32]\) atm, the reflux rate factor, rrf \(\in [1.05, 1.5]\), and the recovery of the light and heavy keys, rec \(\in [0.8, 0.999]\). As there are three distillation units, there are 9 design variables.

The objective function values, as returned by the evaluation of the process flowsheet, include not just the criteria values used to assess the particular design alternative but also the violations of any constraints, e.g. product purity. These are defined in the Criteria object. For the chlorobenzene process, there are 5 criteria in total. Two of these, the capital cost and the operating cost, are actual objective values for optimization and the other three are feasibility constraints. Therefore, when we evaluate the process flowsheet, we split the values returned into the objective function values, as a vector with two values, and a single feasibility indication, as the maximum of the constraint violations:

1: function evaluate(v, flowsheet)
2:     f = Jacaranda.evaluate(flowsheet, v)
3:     # Jacaranda returns, for this process, two
4:     # objective function values and three
5:     # constraint violations so put these in the
6:     # form that an optimization system can use
7:     (f[1:2], maximum(f[3:5]))
8: end

The following code illustrates the use of the evaluate function. The process flowsheet is evaluated at four different points in the search space: the initial guess defined in the Jacaranda input file, x0, the lower and the upper bounds, and at a mid-point between these bounds.

1: using NatureInspiredOptimization.Chlorobenzene: process, evaluate
2: flowsheet, x0, lower, upper, nf = process()
3: [evaluate(x0, flowsheet)
4:  evaluate(lower, flowsheet)
5:  evaluate(upper, flowsheet)
6:  evaluate(lower + (upper-lower)/2.0, flowsheet)]

The results, below, show the results for each of the four evaluations. Each line consists of the two objective function values, within square brackets, and the measure of constraint evaluation.

4-element Vector{Tuple{Vector{Float64}, Float64}}:
 ([2.1074435175185073e6, 3.512848134293994e6], 149.3968323935136)
 ([1.989911840550072e6, 1.3787081622759001e6], 0.7885017421602787)
 ([2.6923280544667e6, 6.012089599489303e6], 463.7476422252121)
 ([1.7892859471243918e6, 4.264066021597366e6], 236.37180766744075)

The first line tells us that the value returned is a vector with four elements. Each element, as expected, is a tuple (with the objective function values) and a real number (the constraint violation). None of the four points evaluated is feasible. One of the challenges of this problem is finding feasible points and any optimization search procedure will need to be able to start with infeasible solutions and hopefully eventually identify feasible solutions on its route to finding a good design.

Having defined the Jacaranda process flowsheet object, and tested it, we can now consider applying each of the optimization methods to find the best design. In the cases that follow, all codes have the same preamble for this problem, which we define here:

1: using Jacaranda
2: using NatureInspiredOptimization.Chlorobenzene: process, evaluate
3: 
4: # define the process flowsheet
5: flowsheet, x0, lower, upper, nz = process()

This preamble, as in the example above, creates the Jacaranda process object and defines variables for the optimization methods, including lower and upper bounds for the search domain. An initial point in the search, x0, has been defined and will be used although it has no special significance. As we have seen above, this point does not yield a feasible design. Further, it has objective function values (in so much as these mean anything for an infeasible design) that are not necessarily good when compared with the values obtained for the lower bounds, upper bounds, or mid-point values of the design variables.

The design of a process using simulation software for the evaluation of alternative designs poses a challenge for gradient based optimization methods. Nature inspired optimization methods, however, are well suited to this problem. The application of the three nature inspired optimization methods to this problem has illustrated the ability of these methods to handle infeasible designs, evolving to feasible solutions quickly. The GA and PPA, in particular, appear to obtain good results with a reasonable computational effort. The PSO does appear to converge quickly to a good solution but then fails to evolve further.

Some process designs require identifying a distributed quantity, i.e. one which varies in either space or in time. An example of this is a plug flow reactor (PFR) in a fixed bed with catalyst distributed along the bed. The catalyst improves or, in some cases, enables the desired reactions. An example [exercise 2, Vassiliadis et al. 2021:311] of the design of a PFR with a catalyst is used to illustrate problems with distributed quantities.

The aim is to maximise the yield of the third species, the product of reactions involving 2 other species, with the feed consisting of the first species alone. The amount of catalyst along the bed of the reactor is \(u(x)\) and the reactor is 2 units of length long. The design problem is to identify \(u(x)\) for \(x \in [0.0,2.0]\) which maximises this yield:

where the right hand side is the molar composition of the third species at the end of the reactor. The value of \(y_3\) is obtained by solving differential and algebraic equations. These equations describe the behaviour of the three species, \(y_i, \, i=1,2,3\) as a function of the distance, \(x\), along the reactor

with initial conditions

and constraints and bounds

As discussed in Section 1.6, the choices of both formulation and representation for the design variables may have an impact on the performance of the optimization method and the resulting outcome. The first decision for this particular case study is how to formulate the continuous distribution of the catalyst over the length of the reactor.

The most straightforward formulation may be to discretise the spatial domain into a number of intervals and specify a step function with a specified constant catalyst amount for each interval in the spatial domain. Another straightforward alternative would be a linear function for the catalyst amount in each interval with the amounts matching at the end-points. The latter is a piece-wise linear formulation. A further alternative would be to define more complex functions, either for each of the intervals or over the full domain, functions which could have desirable continuity and smoothness properties.

A more complex function representation was used in [Rodman et al. 2018] for describing a time varying quantity (the operating temperature). The use of a continuous formulation with smoothness conditions was chosesn as being appealing to process operators, those charged with implementing the operating conditions identified by the solution of the operating design problem. A combination of spline functions for the operating temperature of the fermentation reactor was used to describe the continuous design variable. The same motivation may apply to the design of the PFR as well.

For the purpose of illustration, I will present first a step function to describe amount of catalyst as a function of spatial position, \(x\). Define \(n_x\) intervals over the full the spatial domain, \(x \in [0, 2]\): \[I_i = [x_{i}, x_{i+1}], \qquad i=1,\ldots,n_x \] with \(x_1=0\) and \(x_{n_x} = 2\). For each interval, a constant amount of catalyst, \(u_i\), will be present. Therefore, the decision variables for this optimization problem are the catalyst amounts in the \(n_x\) intervals, \(u_i\), and the intermediate points defining the intervals, \(x_i\), for \(i=2,\ldots,n_x-1\).

The first element to represent is the set of interior points defining the intervals. The simplest would be a uniform discretisation. A uniform discretisation leads to the need for a large number of intervals to enable sufficiently fine granularity to capture changes in the distributed quantity. A large number of intervals implies a large number of values for the distributed quantity to be represented. This can affect the computational efficiency but, more importantly, it will reduce the effectiveness of the optimization methods by increasing the dimension of the search space. Therefore, it may be advantageous to consider an adaptive placement of the points; see [Rodman et al. 2018] for an example of an adaptive representation.

Again, for clarity of exposition, I will present a uniform discretization first, specified by \(n_x\) implicitly: \( x_i = (i-1) \delta x\) for \(i \in 1,\ldots,n_x\) and \(\delta x = \frac{2}{n_x}\). For this discretisation, what remains is to find values of \(u_i \in [0,1]\); these will be the decision variables for the optimization problem.

The results using the discrete step representation, see Figure 1, are designs that an engineer may consider uninteresting or even not realistic for actually building the reactor. Such almost erratic changes in the amount of catalyst to place along the reactor may not be appealing. A similar comment was made, by Rodman et al. [Rodman et al. 2018], when discussing results obtained with a similar discrete representation for the operating temperature of a fermentation unit in brewing. The suggestion in that work was to use a representation for the optimization procedure that is guaranteed to lead to more pleasing outcomes in terms of the actual design.

In this section, I present a similar approach, one which illustrates the use of multiple dispatch in Julia. The challenge is to find a representation that describes designs with desirable attributes. The main desirable attribute would be more gradual changes in the amount of catalyst on the bed of the reactor. Therefore, a representation that is based on a continuous and smooth function would satisfy this desire.

The results obtained by the various methods for larger numbers of intervals tend to have a similar shape: a high value at the start, a plateau along the length of most of the reactor, and then no catalyst for the last part of the reactor would seem to be the aim. However, it is useful to look at the results obtained, with the PPA, for \(n_x=5\), where all the attempts return designs such as that shown in Figure 1. This shows a slight increase in the amount of the catalyst in the plateau region.

Consider the function illustrated in Figure 1:

for \(x \in [0:1]\). Picture two curves like this, each scaled in both dimensions, connected by a straight line. The result could be used to cover the space of designs where the shape of the resulting curve is similar to that observed in Figure 1. Other functions, such as Bezier curves, could be used to represent the possible designs.

The advantage of using a function to represent potential solutions to the design problem is that these solutions will exhibit not only a particular shape but will be continuous over the whole spatial domain and also smooth everywhere except at the two connection points.

To stretch and place the curves, we parameterise the \(a(x)\) function,

by replacing the 10 in the version used for illustration with α. This parameter changes the steepness of the curve. Then we can scale and shift this function by

We want \(f(x) \in [0,1]\) for all values of α, β, and γ. β does not change the range of \(a\); if we restrict α ∈ [0,1], the product is ∈ (0,1). The amount of shifting will depend on α: \(\gamma \in [0,1-\alpha]\). This combination of bounds on these parameters will ensure that the result is suitable for a representation of \(u(x)\) for the reactor design model.

To create the composite function which includes a first component based on \(a(x)\), a connecting line, and then the second case of \(a(x)\), we have to shrink and position the curves appropriately. We define two interior points \(x_1\) and \(x_2\) and then defined the three components as follows, to give a value of \(u \in [0,1]\) for all \(x \in [0,2]\):

where \(f_1\) and \(f_2\) are the same as \(f(x)\) with their own values of α, β, and γ.

This representation of \(u(x)\) has 8 parameters in total, 3 for \(f_1(x)\), 3 for \(f_2(x)\), \(x_1\), and \(x_2\), all with bounds \([0,1]\). The advantage of this representation is that it provides a continuous function throughout the interval with a small number of parameters.

Implementing this alternative representation requires the definition of a function to retrieve the value of \(u(x)\) for any \(x \in [0,2]\) and to be able to access this function from the function that implements the right hand side of equation \eqref{orge7c096d}. The code below uses Julia's multiple dispatch capabilities to implement this.

Process integration, where excess heat in one part of a process is used to meet heating demands in another part, has a long history due to the potential to save significant costs by reducing the amount of external utilities used for heating and for cooling. Heat is transferred between processing steps in a process through a heat exchanger network. A good review of this topic can be found in Chapter 3 of [Towler and Sinnott 2013].

Identifying the best such network is a challenging optimization problem. This problem has been the focus of research in optimization methods and modelling for a long time [Fraga et al. 2001b; Zamora and Grossmann 1998; Morton 2002; Floudas et al. 1986; Linnhoff and Hindmarsh 1983; Linhoff and Flower 1978]. The problem is typically formulated as a mixed integer nonlinear programme (MINLP) representing a superstructure of all the possible structures or networks.

In this chapter, I will present a graph theoretic approach to modelling heat exchanger networks. This approach will use only real valued decision variables and will therefore be suitable for solution using the three nature inspired optimization methods as presented and implemented in this book. This problem is significantly more complex than the previous two and involves the following:

1. definition of data structures for the various aspects of heat exchange network designs;
2. the representation of the superstructure, a directed graph which represents the set of potential structures to be considered by the search procedures; and,
3. the evaluation of a specific network design, requiring traversal of a complex graph and evaluation of the individual units.

Associated with the above is a simple visualization function for presenting solutions to the heat exchanger network design problem.

The solution to a heat exchanger network design problem is a directed graph that represents the flow of process streams and heat. An example network is shown in Figure 1. A graph which includes all possible exchanges between all the process streams and external utilities is known as a superstructure. The aim of the design method is to select, from this superstructure, the best actual network with exchangers designed to exchange heat as required. Defining a superstructure is challenging and so it is worth being able to represent such superstructures easily.

Given the parsed information for each stream, we wish to represent the network in a data structure that allows us to traverse it easily and evaluate the network, given values for amounts to exchange and the split fractions (or binary choices).

The actual network is a directed graph, implementing the state-task network, with states represented by streams (as defined above) as the contents of edges and tasks represented by steps as the content of nodes [Emmerich et al. 2001].

Edges

connect one node to another and carry process streams:

1: mutable struct Edge
2:     from       # node
3:     to         # node
4:     stream     # state
5:     id         # unique id
6:     traversed  # for evaluation
7: end

Nodes

have any number of in and out edges:

 1: mutable struct Node
 2:     in         # edges
 3:     out        # edges
 4:     step       # task
 5:     type       # :hot or :cold
 6:     id         # unique id
 7:     evaluated  # for costing
 8:     z          # objective
 9:     g          # infeasibility
10: end

Network

consists of sets of nodes and edges

 1: mutable struct Network
 2:     inlets :: Vector{Node}  # root nodes
 3:     outlets :: Vector{Node} # leaf nodes
 4:     edges :: Vector{Edge}   # all edges in graph
 5:     nodes :: Vector{Node}   # all nodes in graph
 6:     ne :: Int               # next edge ID
 7:     nn :: Int               # next node ID
 8:     # empty network constructor
 9:     Network() = new(Node[], Node[],
10:                     Edge[], Node[], 1, 1)
11: end

where some of the nodes have been identified as root or inlet nodes and leaves or outlet nodes.

Each of the above types has specific constructors and some ancillary functions:

 1: function Edge(n :: Network)
 2:     e = Edge(nothing, nothing, nothing, n.ne, false);
 3:     add(n, e);
 4:     e
 5: end
 6: function Edge(n :: Network, s :: Stream)
 7:     e = Edge(nothing, nothing, s, n.ne, false);
 8:     add(n, e);
 9:     e
10: end
11: function Node(n :: Network, type)
12:     node = Node(nothing, nothing, nothing,
13:                 type, n.nn,
14:                 false,
15:                 zeros(3),
16:                 0.0);
17:     add(n, node);
18:     node
19: end
20: function Node(n :: Network, s :: Step, type)
21:     node = Node(nothing, nothing, s,
22:                 type, n.nn,
23:                 false,
24:                 zeros(3),
25:                 0.0);
26:     add(n, node);
27:     node
28: end
29: iscold(n :: Node) = n.type == :cold
30: ishot(n :: Node) = n.type == :hot
31: function add(network :: Network, edge :: Edge) 
32:     push!(network.edges, edge)
33:     network.ne += 1
34: end
35: function add(network :: Network, node :: Node) 
36:     push!(network.nodes, node)
37:     network.nn += 1
38: end

Of special note is the initialization of nodes to have three values for the objective function (lines 15 and 24). For the HEN design problem, we will consider capital cost, operating cost, and exchange area potentially as separate criteria for optimization.

In Appendix 8, some functions for printing the structure of a HEN network are presented. I have moved these out of the way as they are long and messy and do not add much to the description of the design of such networks.

The next code block takes the parsed segment structures and creates the network graph:

 1: function segment2graph(network :: Network,
 2:                        edge :: Edge,
 3:                        segment :: Segment,
 4:                        type)
 5:     s = 1                   # count through steps
 6:     ns = length(segment.steps)
 7:     first = true
 8:     let root = nothing, leaf = nothing
 9:         while s ≤ ns
10:             step = segment.steps[s]
11:             s += 1
12:             # ignore Noop steps
13:             if ! (step isa Noop)
14:                 node = Node(network, step, type)
15:                 edge.to = node
16:                 if first
17:                     first = false
18:                     root = node
19:                 end
20:                 node.in = [edge]
21:                 if step isa Exchange
22:                     edge = Edge(network)
23:                     node.out = [edge]
24:                     edge.from = node
25:                 elseif step isa Split
26:                     # first branch
27:                     e1 = Edge(network)
28:                     e1.from = node
29:                     # next step should be a segment
30:                     step = segment.steps[s]
31:                     s += 1
32:                     # create a sub-graph from segment
33:                     (root, leaf1) = segment2graph(network, e1, step, type)
34:                     e1.to = root
35:                     # second branch
36:                     e2 = Edge(network)
37:                     e2.from = node
38:                     # connect branches to split node
39:                     node.out = [e1, e2]
40:                     step = segment.steps[s]
41:                     s += 1
42:                     # create a sub-graph from segment
43:                     (root, leaf2) = segment2graph(network, e2, step, type)
44:                     e2.to = root
45:                     # now join branches up: Mix
46:                     step = segment.steps[s]
47:                     s += 1
48:                     @assert step isa Mix "Expected a Mix after segments but got $step" 
49:                     node = Node(network, step, type)
50:                     # out edges of each branch leaf
51:                     # node are in-edges of the Mix
52:                     # node
53:                     node.in = vcat(leaf1.out,
54:                                    leaf2.out)
55:                     leaf1.out[1].to = node
56:                     leaf2.out[1].to = node
57:                     # prepare to connect to what
58:                     # follows
59:                     oldedge = edge
60:                     edge = Edge(network)
61:                     node.out = [edge]
62:                     edge.from = node
63:                 elseif step isa Utility 
64:                     edge = Edge(network)
65:                     node.out = [edge]
66:                     edge.from = node
67:                 else
68:                     error("We should not get here")
69:                 end
70:                 # the last node will be a leaf of
71:                 # the directed graph
72:                 leaf = node
73:             else
74:             end
75:         end
76:         # return root and leaf of segment
77:         return (root, leaf)
78:     end
79: end

The actual heat exchanger network superstructure, as a directed graph, can now be created using the network structures defined above along with the generation of segments from the representation of the network:

 1: function createnetwork(streams, segments)
 2:     @assert length(streams) == length(segments) "Need a segment for each stream and vice versa."
 3:     network = Network()
 4:     for i ∈ 1:length(streams)
 5:         # create root (inlet) node for stream
 6:         node = Node(network,
 7:                     Inlet(streams[i]),
 8:                     streams[i].type)
 9:         push!(network.inlets, node)
10:         # first edge from root node is the stream
11:         edge = Edge(network, streams[i])
12:         edge.from = node
13:         node.out = [edge]
14:         # and this edge is connected to the graph
15:         # corresponding to the segment for this
16:         # stream
17:         (root, leaf) = segment2graph(network, edge, segments[i], streams[i].type)
18:         # the leaf is also remembered to allow for
19:         # backward traversal
20:         push!(network.outlets, leaf)
21:     end
22:     # link up all process stream exchanges
23:     matchexchanges(network)
24: end

The network structure created by the function above is not fully connected. Missing are the actual exchange links between the hot and cold streams. Each end of the links is present in the hot and cold streams but they are not matched up yet. This is done on line 23, invoking the function defined below.

In this function, the sub-graph associated with each hot stream is traversed from root to leaf. If an exchange node is encountered, the cold streams are traversed in reverse order to find the appropriate matching exchange node. The cold streams are traversed in reverse order because we assume counter-current operation. The links we wish to create correspond to the exchanges illustrated by the top dashed line between the hot and cold streams in Figure 1.

 1: function matchexchanges(network)
 2:     for node ∈ network.nodes
 3:         if ishot(node) && node.step isa Exchange
 4:             # for each hot exchange node, look for
 5:             # corresponding cold node, in reverse
 6:             # traversal of cold stream
 7:             for i ∈ length(network.nodes):-1:1
 8:                 n = network.nodes[i]
 9:                 if iscold(n) && n.step isa Exchange
10:                     # for cold side, link up the two
11:                     # nodes if this cold node is not
12:                     # already linked
13:                     if n.step.hot == node.step.hot && n.step.cold == node.step.cold && nothing == n.step.match
14:                         # link up
15:                         n.step.match = node
16:                         node.step.match = n
17:                         break
18:                     end
19:                 end
20:             end
21:             @assert node.step.match != nothing "Did not find a match"
22:         end
23:     end
24:     network
25: end

A heat exchanger network design problem is defined by the hot and cold streams, as illustrated in Table 1, the utilities available for meeting any unmet demands through external heating and cooling, the superstructure to consider for possible designs (section 5.2), and the cost model described in section 5.1 and implemented by the code in section 5.1.2.4.

To encapsulate the problem, we define an appropriate structure:

 1: mutable struct HENS
 2:     streams :: Vector{Stream}
 3:     representation :: Vector{String}
 4:     network :: Network
 5:     hexmodel # f(A)
 6:     tolerance :: Real
 7:     utilities :: Vector{ExternalUtility}
 8:     function HENS(s, r, h, t, u)
 9:         network = createnetwork(s,parse(r))
10:         new(s, r, network, h, t, u)
11:     end
12: end

The evaluation of a network is a function of the design variables. These include the actual amounts to exchange between process streams, and implicitly with external utilities for any unmet demands, and the split fractions for branches in the network.

The actual amount to exchange between two streams will be defined by a fraction, \(q\), of the heat available in a hot stream for exchange to the cold stream, subject to the cold stream having sufficient need for this amount of heat. Therefore, for some values of \(q\), an instance of the network defined by the values of the design variables may be infeasible. Note that the cold streams will not have \(q\) design variables associated with them; only the hot streams will have these.

The split fraction, \(\eta\), will is the other design variable and every split operation will have one associated with it.

The following helper functions retrieve design variables, for those nodes that have such, and enable setting the values of those design variables.

 1: getVariable(step :: Exchange) = step.Q
 2: getVariable(step :: Split) = step.η
 3: getVariables(problem :: HENS) = getVariables(problem.network)
 4: function getVariables(network :: Network)
 5:     [getVariable(n.step)
 6:      for n ∈ network.nodes
 7:          if ((ishot(n) && n.step isa Exchange)
 8:              || n.step isa Split)]
 9: end
10: 
11: setVariable(step :: Exchange, Q) = step.Q = Q
12: setVariable(step :: Split, η) = step.η = η
13: setVariables(problem :: HENS, vars) = setVariables(problem.network, vars)
14: function setVariables(network :: Network, vars)
15:     map((n,v) -> setVariable(n.step,v),
16:         [n for n ∈ network.nodes
17:              if ((ishot(n) && n.step isa Exchange)
18:                  || n.step isa Split)],
19:         vars)
20: end

Note the use of filtering for the getVariables function and both filtering and mapping for setVariables. Both are functional programming features of the Julia language. Functional programming makes it easy to traverse and work with iterable collections of data, such as the vector of nodes in the network. Filtering is used to extract specific elements from a collection. In lines 5 through 8, the values of the design variables from those nodes in the network that have design variables are extracted. Mapping applies a function to elements from one or more collections. On lines 15 through 19, the map function is used to assign each variable in vars to nodes in the network which are of the right type. These nodes are identified using filtering.

The filtering relies on introspection where the code can interrogate the runtime system as to the types of entities, such as nodes in the network. Also note that use of multiple dispatch for the setting and getting design variable values for the different steps, making the network setting and getting functions more straightforward, not requiring any check of the type of the step in each node.

Once design variables can be set, the evaluation is a matter of traversing the network and evaluating the individual nodes. The hot streams are traversed first as these will define the specified amount of head to exchange between hot and cold streams. Then the cold streams are traversed and the actual exchangers designed. If the amount of heat to exchange is more than the cold streams require, the design will be considered to be infeasible with the excess amount a measure of the infeasibility.

Three different measures of the performance of the network are defined: the total annualised cost, the total area for exchangers, both integrated and utility, and the total amount of heat removed or added using external utilities. This provides the capability to design to minimise any of these quantities or to consider the problem in a multi-objective sense.

Note that the evaluation code updates data in place for the network (the edges in particular) which means that this evaluation cannot be used with multi-processing. This can be addressed by making a deep copy of the network for the evaluation. I have not done so for simplicity.

 1: function evaluate(vars, problem)
 2:     network = problem.network
 3:     setVariables(network, vars)
 4:     # reset network: flags, costs, designs
 5:     map(e -> (e.traversed = false), network.edges)
 6:     for n in network.nodes
 7:         n.evaluated = false
 8:         n.z = [0.0, 0.0, 0.0]
 9:         n.g = 0.0
10:         if n.step isa Exchange || n.step isa Utility
11:             n.step.hex = nothing
12:         end
13:     end
14:     # evaluate whole network with breadth first
15:     # traversal: hot streams first
16:     while reduce(|,
17:                  [(ishot(n) && !n.evaluated)
18:                   for n ∈ network.nodes])
19:         for n ∈ network.inlets
20:             if ishot(n) && !n.evaluated
21:                 evaluate(n, problem)     # recursive
22:             end
23:         end
24:     end
25:     # now cold streams
26:     while reduce(|,
27:                  [(!ishot(n) && !n.evaluated)
28:                   for n ∈ network.nodes])
29:         for n ∈ network.inlets
30:             if !ishot(n) && !n.evaluated
31:                 evaluate(n, problem)     # recursive
32:             end
33:         end
34:     end
35:     # collect and add up objective function and
36:     # feasibility values
37:     (reduce(+, [n.z for n ∈ network.nodes]),
38:      reduce(+, [n.g for n ∈ network.nodes]))
39: end

The network evaluation code also makes use of functional programming aspects of Julia. Specifically, in lines 5, 16-18, 26-28, and 37-38 the map and reduce operations are used. The map operation in line 5 shows the use of an anonymous function to initialise the traversal status of all the edges in the network. The other lines listed use reduction. Reduction applies the given operator, such as | for logical or and + for addition, to the elements of the collection given as second argument. In lines 16-18, the result of the reduce function is true if any element of the filtered network is true. Likewise in lines 26-28. Line 37 uses reduction to sum up the individual values of the objectives; line 38, similarly, adds up the constraint violations.

The evaluation of the network requires evaluation of all nodes (lines 21 and 31). Evaluation of nodes involved designing heat exchangers. A node can only be evaluated if all incoming edges have been traversed (breadth first traversal). This means that the stream associated with each incoming edge is up to date with the information required (temperatures, heat duty available or required). Once the node is evaluated, the procedure continues recursively to evaluate all nodes for the particular stream. The evaluation of a hot side of an exchange simply updates the temperatures and the actual heat duty to be exchanged with the cold stream. The evaluation of a cold stream exchange node designs the actual exchanger for the duty, as specified by the hot stream side.

 1: function evaluate(node :: Node, problem :: HENS)
 2:     # check if node is ready for evaluation
 3:     if node.in == nothing ||
 4:         reduce(&, [edge.traversed
 5:                    for edge ∈ node.in])
 6:         # node can be evaluated which means designing
 7:         # the processing step associated with it.
 8:         design(node.step, node, problem)
 9:         for i ∈ 1:length(node.out)
10:             node.out[i].traversed = true
11:             # recursively traverse graph
12:             if nothing != node.out[i].to
13:                 evaluate(node.out[i].to, problem)
14:             end
15:         end
16:         node.evaluated = true
17:     end
18: end

All nodes are traversed. Traversing a node determines the state of the outputs given the inputs by designing the process step associated with that node.

Using multiple dispatch, the appropriate design method is invoked by the type of processing step (line 8). The Inlet step copies the stream definition to the first Edge for the stream:

1: function design(step :: Inlet, node, problem)
2:     # copy stream to the first (& only) edge
3:     node.out[1].stream = Stream(step.inlet)
4: end

Nodes associated with two specific steps, however, require an actual design: Exchange and Utility steps. For an integrated heat exchange, involving two process streams, the hot stream side of the match specifies the amount of heat to exchange. The cold stream side takes care of designing the actual exchanger, assuming that the match is feasible.

 1: function design(step :: Exchange, node, problem)
 2:     @assert length(node.in) == 1 "Exchange $type node expects one inlet"
 3:     in = node.in[1].stream      # inlet stream
 4: 
 5:     # actual exchange duty: if hot, it is defined by
 6:     # the design variable (step.Q) as a fraction of
 7:     # the available duty (in.Q); if a cold stream,
 8:     # the value should have been set already when
 9:     # the hot side was processed.
10:     if ishot(node)
11:         Q = step.Q * in.Q
12:         # tell cold side the amount to (potentially)
13:         # exchange
14:         step.match.step.Q = Q
15:     else
16:         Q = step.Q
17:     end
18:     # define the stream for the edge leaving this
19:     # node as being the same as the stream coming
20:     # in.  This stream will be modified by the
21:     # exchange if an exchange actually takes place.
22:     # Very small values of a desired exchange are
23:     # ignored.
24:     node.out[1].stream = Stream(in)
25: 
26:     # only do any work if we have an actual amount
27:     # to exchange; note that a negative value
28:     # indicates that the hot stream wished to give
29:     # away more heat than the cold stream can
30:     # accept.  This will be caught in the design of
31:     # the exchanger.
32:     if abs(Q) > problem.tolerance
33:         if Q ≤ in.Q
34:             T = in.Tin + (Q/in.Q)*(in.Tout-in.Tin)
35:             # duty for stream is adjusted
36:             node.out[1].stream.Q = in.Q-Q
37:             # as is the inlet temperature
38:             node.out[1].stream.Tin = T
39:             # if cold stream, design actual exchanger
40:             if ! ishot(node)
41:                 try
42:                     node.step.hex = Exchanger(
43:                         step.match.in[1].stream,
44:                         in,
45:                         Q,
46:                         problem.hexmodel)
47:                     # three criteria
48:                     node.z = [node.step.hex.cost,
49:                               node.step.hex.A,
50:                               0.0]
51:                 catch e
52:                     # the design may fail (Q<0, ΔT<0)
53:                     if e isa Infeasibility
54:                         node.g = e.g
55:                     else
56:                         # something else went wrong
57:                         # so propagate error upwards
58:                         throw(e)
59:                     end
60:                 end
61:             end
62:         else
63:             # infeasible exchange: violation is
64:             # difference between heat required and
65:             # heat available
66:             node.g = abs(Q-in.Q)
67:         end
68:     end
69: end