Sequential linear programming and particle swarm optimization for the optimization of energy districts

Abstract This article deals with the optimization of energy resource management of industrial districts, with the aim of minimizing customer energy expenses. A model of the district is employed, whose optimization gives rise to a nonlinear constrained optimization problem. Here the focus is on its numerical solution. Two different methods are considered: a sequential linear programming method and a particle swarm optimization method. Efficient implementations of both approaches are devised and the results of the tests performed on several energetic districts are reported, including a real case study.


Introduction
In this article the numerical solution of optimization problems arising in energy districts is investigated.
The electric system has been experiencing a dramatic evolution in recent years due to the growing use of unpredictable renewable energy sources and of the spread of distributed generation. On top of that, customers are increasing their attention towards a smarter approach to energy consumption, and paying more and more attention to convenient energy tariff schemes, choosing among different retailers. The problem of efficiently optimizing the resource management of a district is then considered, with the aim of providing the cheapest solution to the customers.
This problem may be addressed in different ways. One possibility is to formulate it as a mixed integer nonlinear programming (MINLP) problem, which can be addressed by adopting a proper linearization yielding a mixed integer linear programming (MILP) problem (Salgado and Pedrero 2008;Bischi et al. 2014). This in turn can be solved by one of the many commercial solvers available. In some cases the problem can even be directly formulated as an (MILP) (Schiefelbein et al. 2015). Another possibility is to formulate it as a multi-objective optimization problem (Ascione et al. 2016;Wang, Martinac, and Magny 2015), where the cost function is not the only objective that has to be minimized.
The research centre 'Enel Ingegneria e Ricerca' in Pisa, in collaboration with the Department of Civil and Industrial Engineering of the University of Pisa, addressed the problem of optimizing an energy district adopting an alternative approach. The problem has been formulated as a nonlinear programming problem (NLP). A software package has been developed (Pannocchia and Mancuso 2014), which takes information as an input such as energy market prices, specifications of CONTACT E. Riccietti elisa.riccietti@unifi.it *Present address: Institut de Recherche en Informatique de Toulouse (IRIT), Toulouse, France the machines, and current and forecast weather information, and builds a model of the energy district. It also computes the objective cost function (measuring the costs related to the exchange of energy with the network) and the functions modelling the physical constraints on the machines' variables. The aim of the optimization process is to find, one day ahead, the machines' parameters setting, providing the optimal dispatching of local energy resources for the following day in order to minimize the cost of energy at the customer site (Ferrara, Riccardi, and Sello 2014).
Here the focus is on the solution of the arising nonlinear minimization problem. In solving such a problem, a compromise between two needs should be found. On the one hand, one aims to find the machine asset corresponding to the lowest possible value of the objective function. On the other hand, the optimization process should be quick as it needs to be performed several times a day when the needs of the district or the weather conditions affecting the renewable resources change.
The objective is a nonconvex multimodal function with hundreds of variables. It does not have a simple analytical expression and it is expensive to evaluate as it is a sum over a large number of terms. Also, the function features some discontinuities. Thus, derivative-free methods are particularly attractive. As an alternative, methods using derivatives can be used if the objective function is approximated around discontinuities by a smooth surrogate function. While the function is still evaluated exactly, the derivatives of the surrogate function are employed. Only first-order methods can be used anyway as, due to such an approximation, one cannot rely on second-order derivatives. Also, the analytical form of only some components of the gradient is available and the remaining ones have to be approximated, for example, via finite differences. Then, approximating also the second derivatives would be computationally too heavy anyway.
With this in mind, it was decided to compare a Sequential Linear Programming (SLP) method and a Particle Swarm Optimization (PSO) method. SLP is a classical optimization method. It is a firstorder iterative procedure that builds approximations to the solution of the original nonlinear problem by considering a sequence of linear subproblems approximating it (Fletcher and de la Maza 1989;Byrd et al. 2003). The generated sequence usually converges to a local minimum close to the starting guess. PSO (Hu and Eberhart 2002;Lu and Chen 2008) is a derivative free method; part of the family of more recently developed metaheuristic algorithms, which aim at finding a global minimum of an optimization problem. This class comprises also Genetic Algorithms (Tessema and Yen 2009), Cultural Algorithms (Jin and Reynolds 1999) and Differential Evolution (Mezura-Montes, Miranda-Varela, and del Carmen Gómez-Ramón 2010). The amount of literature on PSO methods is huge. Many variants of the original approach have been proposed to successfully solve engineering problems. To cite just a few, in de-los-Cobos-Silva et al. (2018); Zhang and Xie (2003); Liu, Cai, and Wang (2010), the basic scheme has been improved by integrating it with Differential Evolution; in Mazhoud et al. (2013), a new constraint-handling mechanism (based on a closeness evaluation of the solutions to the feasible region) is proposed and in Cagnina, Esquivel, and Coello (2011), a double population and a special shake mechanism to avoid premature convergence are used. Many strategies have also been proposed that consist of using PSO as a global search method and coupling it with a fast local search algorithm to speed up its convergence (Fu and Tong 2010;Fan and Zahara 2007;I.F. and Vicente 2009;Martelli and Amaldi 2014).
The aim of this work is to investigate and compare the performance of SLP and PSO methods on these kinds of problems. Specifically, the study is intended to understand if it is worth using a method such as PSO, which aims at returning a solution close to the global minimum but is expected to converge more slowly than SLP, which is a first-order method.
Efficient and reliable implementations of the two procedures are also provided, designed especially for the problem at hand. The adopted SLP method is close to the one presented in (Robinson 1972), but is equipped with a trust-region approach to promote global convergence (Nocedal and Wright 2006). The PSO method is based on the constraint handling strategy introduced in Michalewicz and Attia (1994) and a new strategy to prevent stagnation in local minima, proposed here.
The two solvers have been inserted into the software package developed by the Enel research centre and the University of Pisa and tested on many different realistic examples of energy districts and on a real energy district provided by Enel. Preliminary analysis on this topic has been conducted by the authors in Riccietti, Bellavia, and Sello (2017), where the problem under study is introduced and preliminary numerical results are given. In this article the numerical solution of the arising optimization problem is analysed in depth.
The article is organized as follows. In Section 2 the district modelling is introduced and the machines that are part of the district are presented, focusing on the description of the optimization variables. In Section 3 the arising optimization problem is stated. In Sections 4 and 5 the SLP and PSO procedures implemented for the specific problem are described respectively. The convergence analysis for SLP is reported in Section 2 of the supplementary material. Finally, in Section 6 numerical results are presented. The two solvers are applied to examples of energetic districts and the results of the tests performed are shown. Further numerical results are reported in Section 1 of the supplementary material, where the two methods are applied to a well-known set of benchmark functions which are usually employed to test performance of metaheuristic methods (Liang et al. 2006). The aim is to provide a comparison of the implemented methods with the state-of-the-art metaheuristic methods.

The district model
The optimization problem to be solved arises form the district modelling described below, which was designed in Pannocchia and Mancuso (2014). The user has to specify which devices the district is a compound of, and the characteristic parameters of each of them to simulate their real behaviour. The variables that need to be optimized are the physical parameters of the machines that describe functionality, such as the electrical power produced by generators, the electrial power stored by batteries, and the thermal power provided by boilers. The decision variables chosen are all continuous (Pannocchia and Mancuso 2014), as described in the following subsections, so that the arising problem is a NLP. The decision variables are dimensionless variables in [0, 1] or [−1, 1], obtained by scaling the quantities under control with respect to their maximum value.
The aim of the optimization process is to individuate the optimal generation, load and energy storage profile for each device. A plan of the decision variables for each machine is built for the following day. The aim is to minimize the expenses related to the district management, taking into account real time information such as the price of the electricity (to sell or to buy), wind speed, solar radiation and ambient temperature. The time unit is set to be τ = 15 minutes, so that the day is divided into N = 24×60 15 = 96 quarters of an hour. The solver has to find the optimal decision variables for each device and for each time unit. Let i be the time index, i = 1, . . . , N = 96.
In the next subsections, the four different types of machines that can be included as part of an energy district are described. They are: electrical generators, accumulators, electrical loads and thermal configurations (i.e. groups of heat/cold generation units whose components are predefined) (Ferrara, Riccardi, and Sello 2014).

Electrical generators and thermal configurations
Thermal configurations are groups of heat/cold generation units such as cogeneration engines (or CHP, Combined Heat and Power), boilers and chillers. Examples of thermal configurations may be a CHP, a boiler and a heat pump connected in parallel, or an absorption refrigerator connected to a CHP and an electric chiller. The user can choose among 12 different predefined thermal configurations. Electrical generators comprise photovoltaic, wind turbine and fuel burning generators (Pannocchia and Mancuso 2014).
Let N gen denote the total number of generators and thermal configurations. For the kth device, α k ∈ R N denotes the decision variables vector, k = 1 . . . , N gen . Its component α k (i) is the load at ith time unit: namely the ratio between the actual power output and the maximum power output of the machine. The decision variables are then real and used to decide whether the unit is on or off through the computation of the unit status. The unit status is a discontinuous function θ k , defined as where α min,k is the minimum operating load for unit k. The cost functions for generators and thermal configurations depend on the unit status at each time instant i, which is evaluated as θ k (α k (i)). As a consequence, they feature some discontinuities. For fuel burning generators and thermal configurations with CHP a constraint on the number of ignitions NI(α k ) is provided for. This is an integer number counting how many times unit k is turned on within a day. The number of ignitions must be bounded above so that the following constraint arises: Note that values of NI(α k ) are integers, but NI(α k ) is not a variable of the problem. Subject to optimization, the decision variable is α k , which is real.

Electrical accumulators
Let β b ∈ R N be the decision variables vector for the bth accumulator, b = 1, . . . , N acc with N acc the total number of accumulators in the district. Its components β b (i) are the decision variables at the ith time unit, i = 1, . . . , N, where P b (i) is the power either supplied or absorbed by accumulator b at time i and PB b is the rated power output of the bth battery (Ferrara, Riccardi, and Sello 2014). Negative values of β b (i) means that the accumulator is providing power, and positive values of β b (i) means that the accumulator is storing power. For each accumulator and for each time unit there are also two physical restrictions on the state of charge SOC b , SOC b is a piecewise linear function of β b as it depends on the efficiency of the accumulator η b , which takes two different values depending on the sign of β b (charge and discharge mode), where c 1 is a constant depending on the state of charge of the previous time instant and the daily autodischarge of the accumulator, and c 2 is a constant depending on the nominal power and the capacity of the accumulator. SOC b is then non-differentiable when β b changes sign.

Electrical loads
Three different types of electrical loads are considered: L1 loads, which are mandatory electrical consumptions, L2 loads, which are electrical cycles that need to be completed one or more times at no specific time in the day and L3 loads, which are normally on and can be shut down for a limited amount of time without compromising the related process operation. L1 loads are given parameters for each time unit so these loads do not have decision variables associated and are not included in the count of total number of loads that is N loads = N L2 + N L3 where N L2 and N L3 are total numbers of L2 and L3 loads respectively. For the loads, the decision variables are associated with the indexes of starting times of cycles (for L2 loads) and with indexes of switch-off and switch-on times (for L3 loads). These should be integer variables. In order to avoid working with a MINLP, in the Enel package, scalar continuous variables in (0, 1] are considered as decision variables and they are then related to the integer quantities one actually wants to control (Pannocchia and Mancuso 2014). Let γ m denote the vector of decision variables for mth (L2 or L3) load. For L2 loads the decision variables are given by γ m ∈ R N m , where N m is the number of cycles that need to be completed by mth L2 load. γ m is used to compute the starting time i l of lth cycle as: For mth L3 load the vector of decision variables is γ m ∈ R 2NI m , with NI m maximum number of times that the electric load is activated. The odd components of γ m are related to switch-off times s l and the even ones to switch-on times a l , for l = 1, 2, . . . , NI m , On the decision variables of the loads, there are the bound constraints and also some physical nonlinear constraints. For L2 loads, one constraint ensures that the time between the starting points of two successive cycles is enough to complete the first of them. A second constraint ensures that the amount of time between the beginning of the last cycle and the end of the day is enough to complete the cycle. For L3 loads they guarantee three conditions. First, that the switch-off of the load precedes the switch-on. Second, that the load is not off for more than a given amount of time. Third, that if a load is off and is switched on, a suitable amount of time passes until it is switched off again (Ferrara, Riccardi, and Sello 2014). It is worth underlining that an approximation is introduced in the problem when real variables are considered rather than integer ones. However, taking into account that the time unit is rather small and that Enel is planning to significantly further reduce it in future work, it was preferred to allow such approximation, rather than that working with a MINLP.

Arising optimization problem
The objective cost function f represents the overall daily district expense for exchanging electricity with the network (Pannocchia and Mancuso 2014). Iff i denotes the cost at each time instant i, the objective function is Eachf i can be expressed asf i = c i f i , where c i is the cost of electricity at instant i and f i is the net electrical power exchanged between the district and the network at instant i, which is positive if the district is selling electricity to the network or negative if the district is buying it. The exchange of power with the network can be expressed as the sum of all the partial exchanges of district devices for the ith time unit where f i,k (x) is the power generated by the kth generator, f i,m (x) is the power absorbed by the mth electrical load and f i,b (x) is the power released by the bth accumulator, which is negative when the accumulator is charged (Pannocchia and Mancuso 2014).
If n is the problem dimension, with then f : R n → R is a nonlinear, nonconvex multimodal function that features some discontinuities. As discussed in the previous section, a number of devices have physical constraints on their decision variables (see equations (3), (5)). The constraints are nonlinear and relaxable (minor violations are allowed). They can be stacked together obtaining a vector of constraints g : R n → R p , with p the total number of constraints, where g j (x) for j = 1, . . . , N gen + N acc + N loads is a vector containing the constraints on the jth device.
The resulting optimization problem is the following nonlinear constrained problem: where x ∈ R n is the stacked vector of all devices decision variables and x min ∈ R n and x max ∈ R n denote the bound constraints from (1), (4), (7), i.e.

Sequential linear programming
Sequential Linear Programming (SLP) is an iterative method to find local minima of nonlinear constrained optimization problems by solving a sequence of linear programming problems, (Robinson 1972;Byrd et al. 2003). The adopted SLP method is close to the one presented in (Robinson 1972), but is equipped with a trust-region approach to promote the global convergence. Moreover a boundconstraints handling strategy is adopted, which is similar to the one used in (Herty et al. 2007) where an SLP approach without trust-region is presented for equality constrained problems. The proposed approach is also different from those presented in (Fletcher and de la Maza 1989;Byrd et al. 2003) as second-order information is not used. At each iteration k, f and g are approximated in a neighbourhood of the current solution approximation x k , with first-order Taylor series where d = x − x k . It is worth remembering that functions depending on the status function θ k and the constraints functions for accumulators feature some discontinuities. Then, while they are evaluated exactly, for the computation of gradients they are approximated in a neighbourhood of those points by a surrogate function obtained replacing the discontinuous status and efficiency functions with continuous and differentiable sigmoid functions (Pannocchia and Mancuso 2014).
To obtain a globally convergent method, a trust-region strategy is employed (Conn, Gould, and Toint 2000). Then, a new constraint is added that consists of a bound on the step-length of the following form: d ∞ ≤ k , where k is called the trust-region radius. The new constraint is added to the bound constraints so that at each iteration the following problem is solved: At each iteration the computed solution d k of (11) is used as a step to define the new solution approximation: x k+1 = x k + d k . It is not possible to work directly with problem (11) as the linearized constraints may be inconsistent, i.e. it might be impossible to find a d for which (11b) holds, or the solution found d k might be such that the new approximation x k + d k does not satisfy the nonlinear constraints To deal with this, the constraints are usually added in the objective function as a penalty parameter, i.e. a new term is added to function f that is positive if the constraints are not satisfied, and is zero otherwise. The resulting new objective function is called a penalty function. This term can be chosen in many different ways so that different penalty functions are obtained. Following (Fletcher and de la Maza 1989;Byrd et al. 2003) the l 1 penalty function was chosen and the following penalized objective function is obtained: where ν > 0 is the penalty parameter. If ν is sufficiently large, i.e.
where λ * i , i = 1, . . . , p are the Lagrange multipliers of the inequality constraints g(x) ≤ 0, it is possible to show that l 1 is an exact penalty function (Nocedal and Wright 2006).

Definition 4.1 (Exact Penalty Function):
A penalty function (x, ν) is said to be exact if there exists a positive scalar ν * such that for any ν ≥ ν * , any local solution of the nonlinear programming problem (8) is a local minimizer of (x, ν).
Since it is necessary to take into account the bound constraints, the following problem needs to be considered: If (13) holds, local minimizers of (14) are equivalent to local solutions of (8) to a large extent (see Fletcher (1987); Han and Mangasarian (1979); Janesch and Santos (1997) for details). The general scheme of the algorithm is the same as before, with the difference that instead of linearizing both f and the constraints, function is linearized and at each iteration a problem with just bound constraints is solved, that surely has a solution. As choosing the correct value of ν is difficult, a sequence of linearized penalized problems is solved, adjusting the penalty parameter during the course of the computation. Then, at iteration k given the current iterate x k and the current penalty parameter ν k , the following linear programming problem has to be solved: After the solution d k is found, it is necessary to decide whether to accept the step or not. The step acceptance is based on the agreement between the model function l k and the objective function , which is measured by the ratio between the actual reduction and the predicted reduction (Conn, Gould, and Toint 2000) If ρ k > ρ bad , where ρ bad is a tolerance to be fixed, typically ρ bad ∈ 0, 1 4 , the step is accepted. In this case the trust-region radius is left unchanged or it is possibly enlarged. Otherwise the step is rejected, the trust-region is shrunk and (15) is solved again.
In Algorithm 1, the kth iteration of the SLP method is sketched.
(3) Solve the linear programming problem (15) obtaining a candidate step d k .
Otherwise reject the step: It is possible to prove the global convergence of the sequence {x k } generated by Algorithm 1 to a stationary point of problem (14), as stated in the following theorem. See Section 2 of the supplementary material for the proof.

Theorem 4.2 (Global convergence of Algorithm 1):
Let f and g be C 1 functions and let {x k } be the sequence generated by Algorithm 1. Then, either there exists an iteration indexk such that xk is a stationary point for problem (14), or there exists a subsequence S of indexes such that {x k } k∈S has an accumulation point x * which is a stationary point for problem (14).
Theorem 4.2 states the existence of an accumulation point x * of the sequence {x k } generated by Algorithm 1 that is a stationary point of problem (14), regardless of the starting point. In Section 4.1 the chosen penalty parameter update strategy is described, which ensures that the penalty parameter ν is large enough to let x * be a solution of the original problem.

Implementation issues
In this section, three important features of the algorithm are discussed: the solution of the subproblem (15), the stopping criterion and the updating strategy for the penalty parameter. Function l k in (15a) is non-differentiable, but problem (15) can be written as the following equivalent smooth linear programming problem, introducing the vector of slack variables t (Byrd et al. 2003): Then, the solution d k to the above problem is sought, along with the Lagrange multipliers vectors λ k , π k andλ k of constraints (17b), (17c) and (17d) respectively. These multipliers are employed to implement a reliable stopping criterion for Algorithm 1. Algorithm 1 is stopped whenever a pair (x k , λ k ) satisfies the following conditions: with a tolerance to be fixed. This stopping criterion provides a measure of the closeness of the computed solution to a point satisfying first-order optimality conditions for problem (8). See Theorem 2.5 in Section 2 of the supplementary material for a theoretical support for the employment of such criterion. As far as the penalty parameter is concerned, it is desirable to have penalty function (12) to be exact, according to equation (13) and Definition 4.1. Then, at Step 5 of Algorithm 1, in case of successful iteration, the following updating strategy is adopted: is set. Both in (18) and (20) current Lagrange multiplier estimates of the LP subproblems (17), provided by the function used to solve the LPs, are used as an approximation to those of the original problem.

Particle swarm optimization
Particle Swarm Optimization (PSO) is a stochastic evolutionary method designed to converge to a global minimum of a function f, (Kennedy 2011). It is inspired by the behaviour of bird swarms. Following the natural metaphor, PSO evolves a population of individuals, referred to as particles, within the search space, which behave according to simple rules and interact to produce a collective behaviour to pursue a common aim; in this case the localization of a global minimum. The swarm is composed of s particles, each of them representing an approximation of the global minimum of the optimization problem, and it is represented by a vector x ∈ R n . Each particle is associated with a velocity vector v. The method is an iterative procedure where at each iteration k vectors x k and v k are updated as follows: where x i k and v i k are the position and velocity vector of the ith particle, i = 1, . . . , s, at the kth iteration, p i best,k and p g best,k are respectively the best position reached by the ith particle so far and the best position reached by the whole swarm (the best position is the one that corresponds to the lowest value of the objective function), c 1 , c 2 and w are positive weights, r 1 and r 2 are random variables with uniform distribution in [0, 1], r 1 , r 2 ∼ U(0, 1). At each iteration, vectors p i best,k and p g best,k are updated too, ). The solution approximation provided by the procedure is p g best,k * where k * is the last iteration index. Bound constraints are considered by bringing a particle x back to the nearest boundary that has left the search space. It is also necessary to change the particle velocity, otherwise at the next iteration it is likely to have a new violation: k and r is a random variable uniformly distributed in (0, 1). Originally, PSO methods were developed to deal with problems with only bound constraints, and later they were employed to also solve constrained problems (Aziz et al. 2011;Parsopoulos, Vrahatis et al. 2002). See also (Mezura-Montes and Coello 2011; Jordehi 2015) for a more recent review on constraints handling strategies for PSO methods. Penalty function approaches are widely employed to make the method suitable for the solution of that kind of problem. Following (Michalewicz and Attia 1994;Michalewicz and Schoenauer 1996), the following quadratic penalty function is employed: with τ being the penalty parameter that is decreased at each iteration, so as to increase the weight of the penalty term in the objective function, to penalize with increasing severity constraints violations. Then, given a penalty parameter τ k , function (x; τ k ) is used in (23) in place of the objective function f (x). The kth iteration of the PSO procedure is sketched in Algorithm 2.

Algorithm 2: The kth iteration of the PSO algorithm
(1) Given c 1 , c 2 , w min , w max , k max , τ k , x max , x min , x i k , v i k , p i best,k , p g best,k for i = 1, . . . , s, perform the following steps.
In the implemented PSO method a new strategy to help the swarm escape from local minima is proposed. When the swarm appears to be stuck, i.e. when after a certain number of iterations the value of the objective function is not decreased, the velocity updating scheme (21) is modified adding a new term in the equation, which is proportional to the distance of the global best position from the particle best position, where c 3 is a weighting parameter to be set, and r 3 ∼ U(0, 1). When coefficient c 3 = 0 the weight of the term (p i best,k − x i k ) is increased while that of (p g best,k − x i k ) is decreased. In this way, the particles are less attracted by the global best position, fostering the swarm to escape from a stalemate. When a new improvement in the objective function is obtained, coefficient c 3 is set back to zero and the standard update is employed. Numerical experiments have shown that in some tests the addition of this term indeed helps the swarm to find a better solution approximation. In the following, this PSO variant will be addressed as PSOc 3 .
The main advantages of PSO methods are that they do not require regularity assumptions on the objective function or computation of the first derivatives, and therefore they are suitable when little information on the objective function is available. Clearly the fact that little information on f is used leads to a slow method that requires many iterations to converge. However, the method could be efficiently implemented on a parallel architecture.
These methods are heuristic and standard convergence results. Similarly to those proved for exact optimization methods, they are not usually provided. However, a different kind on analysis of the algorithm can be performed. Using results from dynamic system theory, it is possible to provide an understanding about how the swarm searches the problem space through the analysis of a single particle trajectory or of the swarm seen as a stochastic system (Trelea 2003;Clerc and Kennedy 2002). The analysis provides useful guidelines for the choice of the free parameters to control the system's convergence tendencies.

Implementation issues
The choices that have been made regarding the method of implementation take several factors into account. The problems under consideration have hundreds of variables and the objective function is a sum over a large number of terms. Also, the algorithm must return results quickly to be useful in practice.
Regarding the number of particles, it is impossible to use wide swarms because each particle of the swarm requires two function evaluations; the objective and the constraints functions, at each iteration. On the other hand, using few particles means low exploration capability. After several numerical tests, it was found that a swarm of 20 particles represents a good compromise between solution quality and execution time. For parameter w in (21), instead of using a fixed value, a linearly decreasing scheme is adopted (Shi and Eberhart 1998), where k max is the maximum number of allowed iterations. With this choice, convergence velocity is slowed down ensuring a more accurate exploration of the search space. The process is stopped when a maximum number of iterations k max is performed or when there are no improvements over a fixed number of iterations κ. An improvement is measured in terms of a decrease in the objective function, and it is judged not to be sufficient when the following condition is satisfied for κ consecutive iterations and for a tolerance to be fixed:

Numerical tests
Both solvers were implemented in MATLAB and numerical tests were performed to study their practical performance. First, the proposed methods have been benchmarked against some of the state-of-the-art metaheuristic methods on a set of functions commonly used to test the performance of such approaches. It emerged that the proposed methods can compare with other metaheuristic approaches recently proposed in the literature. The results of these tests are reported in Section 1 of the supplementary material. Then, 13 different problems arising from the described industrial application and corresponding to models of synthetic energy districts were taken into account, plus one arising from the modelling of an existing district in Pisa, provided by the Enel research centre. To solve these problems, the two solvers were inserted in the software tool previously developed by the Enel research centre and University of Pisa.
Here, the values of the free coefficients used in the procedures are specified. For PSOc 3 in (21) c 1 = 1.3 and c 2 = 2.8, in (25) c 3 = 1, in (26) w max = 0.6, w min = 0.1, the tolerance for the stopping criterion (27) is = 10 −3 , κ = 20 and k max = 1700. The swarm size is 20. The choice of this rather small value is motivated by the fact that in energy districts test cases the objective function is expensive to evaluate and the use of wider swarms would lead to prohibitive computational costs. In the PSO constraints handling strategy (24) it was set τ k+1 = (1 − 0.01)τ k and τ 0 = 0.1. For SLP in (12) it was set ν 0 = 1. In Algorithm 1 ρ bad = 0.10 and ρ good = 0.75, and the linear subproblems were solved using the Matlab function linprog(Interior Point Algorithm) with default choice parameters. Notice that the structure of the subproblems is not taken into consideration and a special purpose code could lead to a more efficient solution to the subproblems, but this is outside of the scope of this article.
Results shown in the following tables are the average of those obtained over 100 runs, varying the starting point for SLP solver. In the tables for each one of the test cases, the following statistics are reported:f andk arithmetic mean of the function values and the number of iterations, σ f the standard deviation on function values, min f and max f minimum and maximum function values obtained, time(·) total time in seconds or minutes.

Examples of synthetic energy districts
In this section the results gained by the two methods when applied to 13 synthetic test examples of energy districts are presented. The arising optimization problems have different dimensions, the number of variables is comprised between 294 and 494, the number of bound constraints between 588 and 796 and the number of process constraints is 10 for the first test and 213 for the others. The results of the tests performed are reported in Table 1. At first glance, from these results it is possible to deduce the following remarks. The convergence rate of the SLP method is much higher then that of the PSO method. This is a consequence of the fact that the SLP solver is a gradient-based method. Indeed, as the PSO method does not employ first-order information, it requires a high number of iterations to allow the swarm to carefully explore the search space. In term of computational time, an iteration of the PSO method is cheaper than one of SLP. Each SLP iteration, indeed, requires the solution of a Linear Programming problem and the computation of some derivatives by finite differences, which is computationally expensive. Despite this, the total time required by SLP is much lower than that required by the PSO method, as SLP performs far fewer iterations. As expected, SLP is more suitable for real time optimization.
Regarding function values, it can be observed that in many tests PSO presents higher means compared to SLP, but generally the values on different runs are close to the mean value as the standard deviation is really low. On the other hand, SLP provides lower means but high standard deviation, which means both higher worst values and lower best values. Then, the probability of finding a worst result on a single run using SLP is high, while for PSO it is more likely to have a result close to the mean value. This is highlighted in Figure 1(a), where for each test case the mean values together with the standard deviation for the two solvers are compared. The left bars (light grey) refer to PSOc 3 method and the right ones (dark grey) to SLP method, the standard deviation is highlighted by black vertical lines. Note that overall PSO provides good results, even if a small number of particles are employed. Its performance could be improved increasing the swarm size to better explore the search space. The difference between PSO and SLP is indeed especially evident in Test 5, which is the one in which the search space has the highest dimension. However, from the numerical experience it emerged that increasing the number of particles would be prohibitive from a computational point of view.
Then, a deeper statistical analysis of the results was performed. First it was checked if the results obtained in the tests satisfy the conditions for the use of a parametric test, cf. (García et al. 2009;Derrac et al. 2011). They were tested with the Kolmogorov-Smirnov test, via the Matlab function kstest. It emerged that the data are not normally distributed, meaning that a non-parametric analysis would be more meaningful. Therefore, the data were analysed by the Wilcoxon signed ranks test, cf. (Derrac et al. 2011), based on the mean values found, through the Matlab function signrank. The Wilcoxon test is a nonparametric test that is used to determine whether two independent samples were selected from populations having the same distribution, therefore it can be employed to detect significant differences between the behaviour of two algorithms. The null hypothesis is that the difference between two sample means is zero. From the test it emerges that the data produce evidence that is sufficient to reject the null hypothesis, according to significance level α = 0.05. If R + denotes the sum of ranks for the problems in which SLP outperforms PSOc 3 , and R − the sum of ranks for the opposite, one obtains R + = 84, R − = 21. Then, the statistical tests confirm that SLP shows an improvement over PSO on these problems if the mean values are considered.
In the next section, results of the tests performed on a real energy district are shown.

Pisa district
This is a real district in Pisa provided by the research centre 'Enel Ingegneria e Ricerca' that comprises: a photovoltaic generator with rated power 14 kWe, a wind farm with rated power 3 kWe and 2 L1 loads related to lighting and heating consumptions. Moreover, a thermal group is also part of the district, which is composed of a CHP characterized by rated power 25 kWe and rated thermal power 75 kWt, a gas boiler with rated thermal power 35 kWt and a tank for the storage of hot water with capacity 9400 kJ/deg. This was modelled choosing one among the predefined thermal configurations provided by the software tool. The decision variables for each of these machines are described in Section 2. As for synthetic energetic districts, the time horizon is one day, which is divided into N = 96 time steps. The arising optimization problem has 288 variables, 576 bound constraints and 1 physical constraint. The detailed results of the optimization process are reported in the last row of Table 1.
This test is much less constrained than those presented in the previous section; it has just one mild nonlinear constraint and the dimensionality of the search space is also lower. In this case PSO algorithm performs better than SLP algorithm, providing a lower mean value of the objective function and also a much smaller standard deviation. Notice that in this case also for PSOc 3 algorithm the execution time is quite reasonable.
This test case is particularly interesting because data referring to the cost of unoptimized management of local resources (i.e. the management that is actually running the district) is available. Specifically, the daily cost function of the district for energy consumption is f 0 = 89.3, as it is depicted in Figure 1(b), left bar. By means of the proposed procedures, it is possible to predict values of the machines' decision variables which minimize the cost function. The district management gained by setting the machine parameters to these optimal values will be referred to as optimal management. In Figure 1(b), right bars, the optimal values of the objective function found by SLP and PSO respectively are reported. They represent the cost associated to the optimal management and can be compared to f 0 (the value of the unoptimized cost function) to evaluate savings arising from the employment of the optimized management. The comparison shows that the optimized approach enabled the district to save 18% of the overall daily expense to buy energy from the network.

Conclusions
In this article a model of energy districts is considered that gives rise to a NLP. The main focus of the article is the numerical solution to such a problem. To this aim, a particle swarm optimization solver and a sequential linear programming solver were implemented and compared. The solvers were tested first in constrained tests taken from the literature, and it emerged that they can compare with other recently proposed methods. Then, the solvers were tested on many different examples of synthetic energy districts, including a real case study. From this study it was highlighted that PSO provides values less distributed around the mean, but lower means are generally provided by the SLP method. The results suggest that a hybrid PSO/SLP method may be a promising alternative solver, which will be the objective of future work. Noticeably, from the test performed on the real case study it arises that the optimized management of resources gained by the optimization package provides considerable savings in the energy bill.