A Novel Algorithm with Multiple Consumer Demand Response Priorities in Residential Unbalanced LV Electricity Distribution Networks

Demand Side Management (DSM) is becoming necessary in residential electricity distribution networks where local electricity trading is implemented. Amongst the DSM tools, Demand Response (DR) is used to engage the consumers in the market by voluntary disconnection of high consumption receptors at peak demand hours. As a part of the transition to Smart Grids, there is a high interest in DR applications for residential consumers connected in intelligent grids which allow remote controlling of receptors by electricity distribution system operators and Home Energy Management Systems (HEMS) at consumer homes. This paper proposes a novel algorithm for multi-objective DR optimization in low voltage distribution networks with unbalanced loads, that takes into account individual consumer comfort settings and several technical objectives for the network operator. Phase load balancing, two approaches for minimum comfort disturbance of consumers and two alternatives for network loss reduction are proposed as objectives for DR. An original and faster method of replacing load flow calculations in the evaluation of the feasible solutions is proposed. A case study demonstrates the capabilities of the algorithm.


Introduction
The shift from the old structure of the government-controlled electricity generation and tariff system existing in many parts of the world to the deregulated market model has aimed to incentivize innovation in the field, to lower the electricity cost, and to create new services for the consumers. In parallel, the climate change concerns have raised the need to push forward the generation of 'clean' electricity from renewable and non-polluting primary sources, which led to the necessity of creating the regulatory environment and tools required to enable wholesale and small-scale trading of electricity generated in large or small green electricity sources [1]. The proliferation of small-scale distributed generation sources in Medium Voltage (MV) and Low Voltage (LV) distribution networks changed the operation conditions and management requirements of these networks, that must now integrate new technologies and management procedures. Distribution System Operators (DSO) need to convert the existing supply infrastructure into smart grids with enhanced monitoring and communication capabilities. Consumers have to turn to HEMS in order to manage IoT devices and to sell the surplus of local generation to the market.
Load Control (DLC) methods. Here, multi-objective optimization problems are considered to solve DR for residential users, developing a Pareto Optimal Demand Response (PODR) algorithm [18], or using the Non-Dominated Sorting Firefly Algorithm (NDSFA) in [19]. The mixed problem of utility selection and DR supervision in a smart grid consisting of many electricity utilities and various clients is investigated in [20] with a reinforcement learning and Game-Theoretic Technique (GTT). In [21][22][23], particular approaches for the optimization of the PV-battery system of houses reducing both energy costs and DR activities based on Mixed-Integer Linear Programming (MILP) are presented. In other initial assumptions, the authors from [24,25] propose an efficient objective function to decrease the cost of microgrid operation considering large-scale plug-in electric vehicles and Renewable Energy Sources (RES). Here, the advantage of end-users is modelled by using incentives for including electric vehicle charging in DR programs.
This paper proposes a novel algorithm for DR optimization usable in radial three-phase, four-wire LV electricity distribution networks that supply one-phase residential consumers. Such algorithms were proposed before, in papers such as [9,16,17]. However, in these existing approaches the operation state of the network is assessed using load flow computations, which are time-costly, especially for networks with larger sizes [26][27][28], and they are limited to proposing two objective functions. Table 1 provides a comparison between the capabilities of the proposed approach with regard to the literature review discussed above. The notations from the last column, for the objective functions considered in the optimization, are: OF1-minimize the operation cost; OF2-minimize the consumer dissatisfaction; OF3-active power loss minimization; OF4-maximize the electricity consumption satisfaction; OF5-maximize the net revenue for consumers; OF6-air pollutants emission minimization; OF7-minimize the compensation costs of voltage management; OF8-maximize the load factor; OF9-minimize the curtailment weights. The type and number of consumers of the networks used in the case studies are given in the seventh column of the table, for comparison with the approach used in the paper. The cases with one consumer denote that the DR is performed at the household level.  [3,4] No Yes Yes No NSGA-II Test-1 cons. OF1+OF2 [5] Yes No Yes No PSO Real-151 cons. OF1+OF4 [6] No The proposed method offers more flexibility and ease of use, combining several features present in the literature review summarized in Table 1 and giving the choice of using more objective functions for optimization. The main characteristics of the proposed algorithm are:

•
The use of the Multi-Objective Particle Swarm Optimization (MOPSO) for identifying the optimal solutions • The solutions are provided as Pareto frontiers for optimizing two conflicting objectives

•
In the initial implementation, the users can choose between five possible objective functions, belonging to two distinct categories (technical and user comfort): Phase load balancing at the supply substation level The minimization of individual consumer disturbance caused by DR The minimization of overall consumer disturbance caused by DR An original, faster method for active power loss minimization A load-flow based method for active power loss minimization (for comparison purposes).
Except for the optional load flow method, the other objective functions require the knowledge of minimal data about the electrical network, that can be grouped in a matrix consisting of eight lines and a number of columns equal to the number of consumers existing in the network. The formulation of each DR objective has a simple mathematical model, which contributes to a fast execution time for the algorithm. The specific formulation of the input data allows the possibility of implementing individual consumer DR segmentation on independent levels according to comfort preferences. The optimization is based on the MOPSO algorithm because it allows an easy definition of other objectives for DR and has fast execution times. The results of the case study, performed on a real LV distribution system, show that the algorithm performs as expected for all the pairs of objective functions tested, and that the proposed loss minimization method combines increased computation speed and good accuracy when compared with the classic load flow approach.
The rest of the paper is structured as follows: Section 2, Materials and Methods, provides the detailed description of the proposed DR objectives, with their differences and practical usability, together with the mathematical model for each objective, followed by the general structure of the MOPSO algorithm in which they are used; Section 3 presents the electrical network used for testing the algorithm and the results of the case study; Section 4, Discussion, presents the advantages and limitations of the proposed approach, ending with the future research directions envisioned by the authors.

Materials and Methods
The proposed multi-objective DR optimization algorithm is aiming to provide to DSO and other entities a tool for enabling DR initiatives in distribution systems. As a working hypothesis, it is considered the case of low voltage four-wire distribution networks that supply mostly one-phase consumers connected through one-phase branchings to three-phase feeders. As it would be the case in real conditions, it is also considered that only a fraction of the total number of consumers connected in the network will participate in the DR program, due to technical (i.e., too low consumption) or personal (i.e., unwillingness) reasons. Out of the participating consumers, some can be equipped with HEMS systems which allow them to set custom DR levels to be accessed by the algorithm, that are modelling consumer preferences or comfort levels. Such systems are widely studied in the literature, one example being the HEMS proposed in [29].
When the utility initiates a DR request, given as a power amount P dec to be disconnected from the entire network, specified in kW or percent fraction from the total current network load, the algorithm has the task of distributing the power quantity requested between the consumers in the network such as a pair of two conflicting objectives to be fulfilled. The user (DSO) can choose from five available objective functions: • F1-Phase load balancing. The aim of this objective is to disconnect individual loads from the network so that phase load balancing is to be achieved on all three phases at the substation which supplies the entire LV network, or, in other words, to minimize the differences between the phase power flows at the supply point.
• F2-The minimization of individual consumer disturbance. This goal aims to distribute the load disconnection between the available consumers such as each individual consumer will have a minimal comfort reduction because of its participation to DR. • F3-The minimization of overall consumer disturbance, that will ensure that the DR request will affect the minimal number of consumers. • F4-Active power losses minimization in the network. Here, the aim is to disconnect loads ensuring that the operation of the network with the remaining load will result in minimal active power losses.
For this objective, an original approach is implemented, which prioritizes for disconnection the loads located near the farthest end of the distribution feeders and the loads where the power disconnected during the DR event is maximum. This approach aims to prevent the use of time-consuming load flow algorithms for computing the losses and to minimize the amount of input data used in calculations, as it does not require the knowledge of the network configuration or of the electrical parameters of its branches.
• F5-Active power losses minimization in the network. In this approach, the active power losses are computed using a load flow algorithm. This is the method most used in the literature, and is implemented for comparison with the approach proposed for F4.
All the fitness functions are taking into consideration the deviation between the total power disconnected in the DR event and the DR target set by the utility, P dec , which needs to be fulfilled as close as possible.
A more detailed description and the mathematical models for F1-F5 will follow in Section 2.3. The results are provided as Pareto frontiers, which are built using the Multi-Objective Particle Swarm Optimization (MOPSO) [30], described briefly in the following subsection.

The Multi-Objective Particle Swarm Optimization Algorithm
The MOPSO algorithm, proposed in 2002 by C.A. Coello and M.S. Lechuga for solving Pareto-conditioned problems, is an extension for multi-objective optimization problems of the well-known Particle Swarm Optimization (PSO) algorithm [31]. Its basic flowchart for two objective functions is provided in Figure 1. Compared to the original PSO [32,33], the MOPSO evaluates the solutions, that are encoded as numeric vectors, for two conflicting objective functions, and accommodates the necessary modifications to keep the multiple solutions that make a Pareto frontier, using an external archive. The population members (particles) change their speed and position for the next iteration (p (it+1) ) according to the classic PSO equations: In Equations (1) and (2), v (it) and v j (it-1) are the speed of the particle in the previous (it-1) and current (it) iteration, rnd 1 and rnd 2 are vectors of random numbers, p best (it) and p (it) are the best personal and the current position of the particle, leader (it) is the position of the leader in iteration it, and w is an inertia coefficient decreasing over the iteration count to progressively narrow the search space around the current optimal solution. y modifications to keep the multiple solutions that make a Pareto frontier, using an e The population members (particles) change their speed and position for the next i cording to the classic PSO equations: igure 1. The flowchart of the multi-objective particle swarm optimization (MOPSO) algorithm After updating its position, each particle is subjected to a mutation procedure for increasing diversity inside the population. Moreover, since the optimal solution is not unique, the leader is selected for each particle from the external archive. Any mutation and selection procedure applied in the literature for evolutionary algorithms can be used for this purpose.
In the Pareto frontier of optimal solutions (external archive) found by the algorithm are allowed only non-dominated solutions. The conditions required for non-dominance when solving minimization problems are depicted in Figure 2. Usually, a fixed maximum number of solutions NS max is allowed in the Pareto frontier. For enforcing this constraint while encouraging a uniform distribution of solutions across the entire length of the frontier, the solutions from the most crowded sections are pruned with specific procedures, so that the number of solutions that are found at any time in the frontier, NS, is lower than NS max .
In Equations (1) and (2), v (it) and vj (it-1) are the speed of the particle in the previous (it-1) and current (it) iteration, rnd1 and rnd2 are vectors of random numbers, pbest (it) and p (it) are the best personal and the current position of the particle, leader (it) is the position of the leader in iteration it, and w is an inertia coefficient decreasing over the iteration count to progressively narrow the search space around the current optimal solution.
After updating its position, each particle is subjected to a mutation procedure for increasing diversity inside the population. Moreover, since the optimal solution is not unique, the leader is selected for each particle from the external archive. Any mutation and selection procedure applied in the literature for evolutionary algorithms can be used for this purpose.
In the Pareto frontier of optimal solutions (external archive) found by the algorithm are allowed only non-dominated solutions. The conditions required for non-dominance when solving minimization problems are depicted in Figure 2. Usually, a fixed maximum number of solutions NSmax is allowed in the Pareto frontier. For enforcing this constraint while encouraging a uniform distribution of solutions across the entire length of the frontier, the solutions from the most crowded sections are pruned with specific procedures, so that the number of solutions that are found at any time in the frontier, NS, is lower than NSmax. Several algorithms from the literature can be used for multi-criterial optimization, as the survey from Table 1 shows. The choice for the MOPSO considered several advantages. The algorithm was previously used with proven results in various engineering fields such as industrial manufacturing Several algorithms from the literature can be used for multi-criterial optimization, as the survey from Table 1 shows. The choice for the MOPSO considered several advantages. The algorithm was previously used with proven results in various engineering fields such as industrial manufacturing [34], photovoltaic panel design [35], water reservoir operation [36] or optimal dispatch of electricity generation [37]. Compared with other classic or metaheuristic optimization algorithms, MOPSO has a simpler mathematical model, which makes its more accessible and flexible for industry specialists looking to adopt new management tools with minimal training and costs. The simple mathematical model minimizes the computation time, improving the speed of the optimization and of the decision making process. Moreover, the PSO algorithm was previously used by the authors for solving single-objective DR optimization, in [38].

The Input Data Required and the Encoding of a Solution
If load flow calculations are not used, the proposed algorithm requires minimal input data for performing the DR optimization. For an electrical network with NN buses or nodes, NC consumers and NB branches, one input matrix DataIN is required, with the structure presented in Table 2, where, for any consumer i, i = 1..NC, the following notations are used: index i -consumer index, phase i -the connection phase (1, 2, 3 or 0, denoting phases a, b, c, and three-phase consumer, respectively), dist i -the distance measured between the consumer and the supply source (MV/LV substation), PC i -the power demand of the consumer, in kW, in the interval of analysis and in the absence of DR load reduction, P DR,i,1 -P DR,i,4 -the power which the consumer is willing to disconnect in a DR event. The DR availability of any consumer, denoted as maxDR i, is described by a maximum of four levels of comfort that can be accessed in their decreasing order (indexes 4 to 1) by the utility in order to fulfill its DR objective. The consumers that do not participate to DR will not use any of these levels and will not be considered for optimization. The consumers with at least one active level are candidates for DR. They submit their maximum disconnection levels which they make available for load reduction, and the algorithm selects the actual number of levels that will be used for optimization. The consumers with one active level can also account for clients who do not have HEMS implemented, but wish to participate in the DR program. Thus, the real number of consumers from the network participating in the DR program is NC DR ≤ NC. A particle that describes a solution of the MOPSO algorithm is a vector of integer numbers of size NC DR : and has the structure described in Table 3, where level j is the highest available DR level remaining in use after optimization (i.e., the power available at levels level j +1 to maxDR j is disconnected during the DR event). An example of writing the matrix DataIN and interpreting a solution p j is given in Appendix A, at the end of the paper.

The Objective Functions Used for Optimization
The MOPSO algorithm can be used to build the Pareto frontiers for two out of five conflicting objective functions, previously defined as F1-F5. It should be noted that the objective function pool can be easily extended, with minimal changes to the algorithm, required only by possible new input data required by the new objective functions. The mathematical models of the objective functions are designed to be simple and intuitive, for minimizing the computation time and increasing the practical usability. The information contained in the input data is minimal and well known by network operators.

The Objective Function F1-Maximum Substation Phase Load Balancing (SPB)
The optimal operation conditions of three-phase electricity distribution systems require that the power flows on the three phases should be as close as possible to equal, on all the feeders that supply consumers. A substantial deviation from this constraint can lead to two undesired effects: • If the phase demand is unbalanced, power flows on the neutral wire increase, leading to supplementary power losses.

•
On the phase(s) with the highest load, supplementary power losses and higher voltage drops will occur, increasing with the length of the feeders and the load.
The normal operation conditions of these networks are usually with unbalanced phase load, because of the continuous variation of the demand and unoptimized phase connection of the consumers. Optimized DR can help alleviate the imbalance by finding the best candidates for DR, from the available pool, such as the demand remaining unaffected to be more balanced on the phases.
The algorithm used in the paper proposes that DR should be requested from the consumers in such a way that power is disconnected predominantly from the phases with the highest load, in order to equalize the power flows on the three phases. The balancing procedure used for DR optimization consists of several steps. F1.1. For each DR consumer j, read the DR request (the actual DR operation level which will remain active in the DR interval, as requested by the utility), according to the current particle p. This information allows the computation of the power disconnected from each consumer and phase.
where j = 1..NC DR. F1.2. Find the total active power which will be disconnected from each consumer, PD j , and subtract it from the total power consumption of the consumer, in order to determine the power requirement PCj remaining on each phase, after DR is applied.
Mathematics 2020, 8, 1220 9 of 24 F1.3. Using the consumer phase allocation phase i , i = 1..NC, compute the sum of demand P ph on each phase ph: F1.4. Compute the standard deviation of active power flows on the three phases of the network at the LV busbar. A low standard deviation of the three phase power flows (P a , P b , P c ) is a result of close values of these three quantities.
F1.5. Compute the difference between the total load disconnected in the DR event and the target set by the utility, because the phase balance must be achieved with the constraint that the amount of power disconnected via DR matches the target set by the network operator, P dec : F1.6. Compute the objective function F1 for particle p as a combination of the considered targets, balancing and amount of disconnected power, where a high deviation of the latter is seen as a penalizing factor for the first: The objective of the optimization is to find the solutions with minimal imbalance which adhere to the target quantity of power disconnected set by the network operator: One of the deterrents that would make the residential consumers avoid subscribing into DR initiatives proposed by utilities is the expected loss of comfort generated by the need to limit their electricity demand during the DR hours, which often will coincide to evening times, when the home power demand is at its peak. To alleviate this problem, the utility could consider to distribute its total DR requirement between the consumers such as the comfort level of each consumer would be minimally affected. This is similar to disconnecting power from each consumer on fewer DR levels as possible.
The computation of this objective function requires the following steps: F2.1. For each DR consumer j, read the DR request actualDR j , with Equation (4), according to the current particle p.
F2.2. Determine the number of affected DR levels, for each consumer j: F2.3. For affected levels greater than 1, simulate a penalizing factor. A high number of affected levels for a consumer means a greater loss of comfort. The higher the number of levels affected for an individual consumer, the higher the penalizing factor should be, increasing the value of a2. A low value of a2 means that, generally, consumers need to disconnect a minimal number of DR levels.
Mathematics 2020, 8, 1220 10 of 24 F2.4. Compute the difference between the total load disconnected in the DR event and the target set by the utility, with (9).
F2.5. Compute the objective function F2 for particle p as the average value of a2 for each DR consumer in the network, penalized by the deviation between the actual amount of disconnected power and the target set by the network operator: The two factors are averaged for the number of DR consumers and DR target for obtaining similar fitness values in networks with a different number of consumers and different DR targets set by the network operators.
The objective of the optimization is to find the solutions with the lowest value of F2: Another goal of an electricity utility when using DR in the network can be to minimize the number of affected consumers. Such an approach would ensure that the majority of consumers will not have to reduce their demand, and their comfort will not be affected. However, the comfort loss cost would be transferred to the consumers which will need to reduce their load, as they would be forced to minimize their demand. While this objective could be considered unfeasible to pursue on its own in real operation conditions, because of the high stress placed on the affected consumers, it becomes of interest in the context of Pareto optimality, when another objective such as power loss minimization is sought.
This objective is implemented in the MOPSO algorithm by searching for solutions that use the maximum number of DR levels for individual consumers. The following steps are used: F3.1. For each DR consumer j, compute the DR request actualDR j , with Equation (4). F3.2. Determine the number of affected DR levels, for each consumer j, with (12): F3.3. Count the number of zero elements computed in the previous step, that should be maximized, because this is the number of consumers which do not disconnect power in the DR event and will not suffer from loss of comfort: F3.4. Compute the difference between the total load disconnected in the DR event and the target set by the utility, with (10).
F3.5. Compute the objective function F3 for particle p, as the inverse of coefficient a3 averaged over the number of DR consumers, taking into account, as for F1 and F2, the compliance with the target of the network operator: The objective of the optimization is to find the solutions with a minimal number of affected DR levels for each individual consumer: min

The Objective Function F4-The Active Power Loss Minimization Using the Simplified Method Proposed in the Paper (SMP)
One of the main objectives of the DSO utilities is to operate their networks with minimal active power losses. At medium and high voltage levels, this goal can be achieved by optimally setting the operational configuration of the equipment or by optimal generation dispatch. At LV level, load management is one of the few measures available that do not necessitate supplementary investment costs. As it has been shown in the literature survey, this objective has been approached in the literature mainly using load flow calculations to determine the losses associated with the optimal DR solutions. While offering the best precision, this method is time consuming and requires the knowledge of the operating configuration of the LV network, together with the electrical parameters of all its elements, which are not always known in detail.
This paper proposes a much simpler approach to finding DR solutions which lead to network operation with minimal losses. Its advantages, with a marginally precision loss trade-off, are: • Significantly reduced computation time when compared to load flow methods, which makes the approach more attractive for implementation in real-time applications.

•
The knowledge of network configuration and electrical parameters is not needed.
The proposed method combines two principles widely used for loss minimization in radial electricity distribution networks, by encouraging DR load reduction for two categories of loads: • Consumers located far from the supply point • Consumers with a high demand and number of DR levels.
The associated objective function for this approach is using several steps: F4.1. For each DR consumer j, compute the DR request actualDRj, with Equation (4). F4.2. Find PD j and PC j , for each DR consumer, with Equations (5) and (6). F4.3. For each DR consumer j, find the network length separating it from the source, dist j , as given in the input matrix DataIN.
F4.4. Compute the total distance separating the consumers affected by the DR request from the source-the distance factor, the first component of the objective function: F4.5. Compute the total power disconnected from the consumers in the DR event-the maximum load disconnection factor, the second component of the objective function: F4.6. Compute the difference between the total load disconnected in the DR event and the target set by the utility, with Equation (9). F4.7. Compute function F4 for particle p combining the values for the three components (distance, high individual load, DR target): where P net is the power demand of the LV network before the DR event: The objective of the optimization is to find the solutions with the minimal value for F4:

Objective Function F5-Active Power Loss Minimization Using a Load Flow Method (LFM)
For the assessment of the performance of the load minimization method proposed in the previous subsection, a graph-theory based load flow method for radial distribution networks was implemented. Its steps follow.
F5. F5.2. Using the incidence matrix A written for the distribution network and the bus currents vector I b = I bus ∈ R 1x NC , find the branch currents vector I br , given in [A]: F5.3. If I branch are elements from I br , compute the branch losses for the entire network in [kW]: In Equation (28), K branch is a coefficient that accounts for the loss increase due to the phase load imbalance, computed with: where CUF branch is the current imbalance factor computed for each branch with: In (30) Each objective function is called a subroutine as functions F1 and F2, as seen in Figure 1. As an example, in Algorithm 1 the code used for the objective function F4, SMP is provided. Similar algorithms can be written for the other objective functions. In Algorithm 1, indexDR is a vector of size NC DR containing the indexes of the DR consumers in the network, determined using the first row from DataIN and the list of known DR clients from the network, pop is the population of the MOPSO algorithm, and NP is its size. Algorithm 1. The subroutine used to encode the calculation steps for the objective function F4-the simplified loss minimization method proposed in the paper.

Case Study
The algorithm aims to provide a flexible tool for DSOs and local microgrid operators who choose to implement DR programs for optimizing the operation of their network. For demonstrating its performances, a low voltage distribution network with real measured data and corresponding to the residential consumer profile mostly expected to adhere to DR (individual houses from a residential neighborhood) was chosen. The five objective functions were grouped in pairs and the optimization algorithm was run. The results show that the choice of objectives leads to different DR requirements from the consumers. The newly proposed objective function F4, aimed to replace time-consuming load flow computation, shows good adequacy.

The LV Network Used in the Case Study
The DR algorithm described in Section 2 was tested using a real LV, four-wire electricity distribution network, its one-line diagram being shown in Figure 3. The network supplies only one-phase residential consumers through Overhead Lines (OHL) feeders mounted on poles. In Figure 3, the black bullets represent pole indices, and the arrows point to the number of consumers supplied at each pole.
The DR algorithm described in Section 2 was tested using a real LV, four-wire electricity distribution network, its one-line diagram being shown in Figure 3. The network supplies only onephase residential consumers through Overhead Lines (OHL) feeders mounted on poles. In Figure 3, the black bullets represent pole indices, and the arrows point to the number of consumers supplied at each pole. The feeders are four-wire, while the consumer branchings are two-wire (active phase and neutral conductors). The supplied area consists of mainly one-and two-storied houses located at the outskirts of a small sized town. The main technical characteristics of this network are provided in Table 4. The consumers supplied by the network can be classified, based on their peak demand, as in Table 5. As described in Section 2, only the consumers with a peak consumption greater than a given limit (0.5 kW in this case) were included in the DR pool. In the time interval used for analysis, this approach resulted in 40.5 kW of power available for DR, from a maximum number of 84 consumers, The feeders are four-wire, while the consumer branchings are two-wire (active phase and neutral conductors). The supplied area consists of mainly one-and two-storied houses located at the outskirts of a small sized town. The main technical characteristics of this network are provided in Table 4. The consumers supplied by the network can be classified, based on their peak demand, as in Table 5. As described in Section 2, only the consumers with a peak consumption greater than a given limit (0.5 kW in this case) were included in the DR pool. In the time interval used for analysis, this approach resulted in 40.5 kW of power available for DR, from a maximum number of 84 consumers, at a total network power consumption of 116.6 kW. It should be noted that the LV substation transformer which supplies the network has a rated power of 180 kVA.

Algorithm Setup and Results
The algorithm was run for a maximum of 1000 generations, with 70 individuals in the population. The Pareto frontiers obtained for different pairs of objective functions are presented in Figures 4-8 and Tables 6-10.

Algorithm Setup and Results
The algorithm was run for a maximum of 1000 generations, with 70 individuals in the population. The Pareto frontiers obtained for different pairs of objective functions are presented in Figures 4-8 and Tables 6-10.

Algorithm Setup and Results
The algorithm was run for a maximum of 1000 generations, with 70 individuals in the population. The Pareto frontiers obtained for different pairs of objective functions are presented in Figures 4-8 and Tables 6-10.      . Figure 8. The Pareto frontier for the pair F1-F5.    . Figure 8. The Pareto frontier for the pair F1-F5.    In Appendix B, Table A1, the solutions obtained by the algorithm corresponding to the best values obtained for each objective function are given, compared to the initial solution (no DR applied). For comparison purposes, the synthetic data (functions F1-F5, network losses and unbalance, number of affected consumers) are provided in Tables 11 and 12. Several aspects can be observed from the data presented above. First of all, the number of solutions in the Pareto frontier is small for any pair of objective functions chosen for optimization. This can be attributed to the fact that, as an analysis on the matrix DataIN shows, there are multiple solutions that can involve different consumers chosen for DR, resulting in the same value of the fitness function. As an example, disconnecting 0.2 kW alt DR level 2 from consumers 146 or 167 would have the same effect on the value of objective function F1 (MPB, Maximum Phase Balancing at substation level) or F3 (OCD, Overall Consumer Disturbance). Due to its size, the matrix DataIN is provided only as a Supplementary File with the paper.
The objective function F1 (MPB) has the same effect on power losses as F4 (SMP) and F5 (LFM), because phase balancing leads to lower current flows on the neutral wire and lower maximum phase current flows, which both result in lower active power losses. This is why it can be observed that the pair F1-F5 has the lowest number of solutions, because objectives F1 and F5 are not always conflicting (Table 10, Figure 8). The same can be observed for pair F5-F4 (LFM, SMP), where the extremities of the Pareto frontier are very narrow (between 6.23 and 6.3 for F5 and 0.49-0.52 for F4). If the load flow method is used as reference, the error in loss calculation is of 3.34% for the proposed method, while the loss reduction achieved by DR is of 33% (best LFM).
The values from Figure 7 and Table 9 show that the proposed method can be a good tradeoff between speed and accuracy. While the results regarding the active power loss in the network are marginally worse than those obtained using the load flow method (F5) (5.99-6.06%, as opposed to 5.79%, and compared to the reference no DR case of 8.64%), they are better than the active power losses obtained after balancing with the objective function F1-6.15%, as seen in Table 12. The proposed method, SMP (F4), has also the advantage of better speed. In MATLAB, on a Ryzen 5 1600X 6-core machine with 16 GB of RAM, the time cost of running an iteration of the algorithm with the function pair F3-F4 is 0.214 s, while running the function pair F3-F5 takes 0.654 s, which makes the proposed method three times faster than the classic graph theory-based load flow. Running the algorithm with the settings specified in Section 3.2 (1000 iterations, population size 70, Ryzen 5 PC with MATLAB), the average running times for combinations of F1 to F4 varied between 37 s (pair F1-F4) and 33 s (pair F2-F3). If the time-costly load flow objective is used, the time can amount up to 235 s (pair F1-F5). These times can be further improved by using distributed, parallel computation amongst the available computing cores and also depend on the number of DR consumers in the network. In the paper, parallel computation was disabled and the number of DR consumers (thus the length of the MOPSO particles processed) is high, at 84, for simulating the use of hardware and network situations that are common in real operation conditions.
Regarding the robustness characteristics of the proposed method, the computation times obtained above, the literature comparison from Table 1, and the network data from Table 4 show that the proposed methods are suitable for large size networks, with a high number of consumers. Using as reference the data from [39], the network chosen for the case study is representative for the type of residential consumption currently existing in Europe, and comparing it with the information from Table 1, it is the amongst the largest used, at a length of 1.5 km and a number of DR consumers of 84 out of 176 possible.
Using any of the objective functions F1-F5 leads to reducing the active power losses in the network, because in all the cases DR reduces the network load, but the most significant reductions are obtained using the objective functions F1, F4, or F5, which are mostly technical in nature. The other two objective functions, F2 and F3 are comfort based, aiming to reduce the disturbance for the individual consumers affected by DR (F2) or the number of consumers affected (F3). As Table 12 shows, the number of affected consumers is the largest for F2 and the smallest for F3, a behavior consistent with the sought objective (in order to achieve the same disconnected power amount, a larger number of consumers will be affected less and a smaller number of consumers will be affected more). From Table A3, it can be seen that when optimizing with function F3, the algorithm chooses for DR preferentially the consumers with the highest number of DR levels available.
The economic aspects of implementing the DR initiatives in LV residential networks can be viewed from two perspectives. First, there is the effect on the network. Reducing power losses leads to savings in primary energy resources. For the best cases found by the algorithm, losses are reduced via DR with 4.31 kW when using the objective function F4 and 4.51 kW when using the objective function F5. If this value is extrapolated for 365 days and 3 h a day (the peak evening hours), the saving amount to 4.72 and 4.94 MWh of electricity. Their cost can vary according to the primary source. At a reference price of 0.2 EUR/kWh for households, as seen in the EUROSTAT statistics, the yearly savings for this network and DR scenario would be up to 988 EUR.
On the other hand, the consumers agree to participate in DR in exchange for other services or tariff reductions offered by the network operator. These advantages are difficult to quantify, because they can imply both economic and comfort components.
Moreover, for all the fitness functions used, the algorithm complies with the DR limit set by the network operator (20.25 kW) with an error of 0.05 kW determined by the one-decimal load representation on the DR levels. The only exception is the optimal solution found for F3, which is a 0.15 kW error.

Discussion
The algorithm presented in the paper can serve DSOs and other entities such as microgrid operators as a tool to test and implement DR services which take into consideration multiple optimization goals, combining technical and user comfort criteria. The higher number of fitness functions that are available in its current implementation give it an advantage over other alternatives from the literature. Furthermore, there is the possibility to include with minimum modifications more types of objectives, as each objective is being implemented as a subroutine called independently.
In real world implementation, the algorithm is, however, dependent on the existence of intelligent grids, capable of two-way communication with the network operator, and intelligent Home Energy Management systems are required to enable the multi-level DR approach. This can be currently considered a limitation in its applicability, until Smart Grids will become more common.
The results of the case study have shown that there is a dependence between the chosen optimization and the results. For instance, the DR requests that give lowest power loss values are not detected by combinations of fitness functions that include user comfort-oriented objectives, which, in turn, lead to solutions with higher than average loss values.
The results from Tables 6-10 suggest also that the algorithm is more likely to find compromise solutions, rather than favoring the cases that belong to the extremities of the Pareto frontier, weighted towards one of the objective functions. Some reasons for these results could be found in the significant length of the particles needing to be optimized by the MOPSO (84 in this case) and in the typical behavior of metaheuristics, which have simpler optimization models, with the tradeoff of finding near-optimal solutions.
The simplified approach to load minimization proposed in the paper, while does not reach the global optimal loss value, shows the advantage of speed coupled with near-optimal results, which can constitute an acceptable compromise for real time or short-term optimization. This objective function can have an advantage over the classical load flow approach in scenarios when network operators prioritize speed in their decision-making processes.
Future research directions envisioned for the algorithm are the testing of other optimization methods such as the NSGA-II for generating the Pareto frontiers, further improvement of the mathematical model of the DR objectives, the extension of the optimization model to include several consecutive DR intervals and the DR rebound, and further expansion of the number of concurrent optimization criteria, by using a many-objective approach such as the NSGA-III algorithm [40]. Funding: This research received no external funding.

Acknowledgments:
The paper is a result of an international collaboration started in the framework of the Erasmus Mundus gLINK-Sustainable Green Economies through Learning, Innovation, Networking, and Knowledge Exchange program.

Conflicts of Interest:
The authors declare no conflict of interest.

∆P
The The current position of the particle P a , P b , P c The three phase power flows p best The best personal position of the particle PC The power demand of the consumer PD The total active power which will be disconnected P dec Power disconnected P DR The power which the consumer is willing to disconnect in a DR event ph The phase phase The The bus rated voltage v The speed of the particle w The inertia coefficient

Appendix A
For the LV distribution network from Figure A1, with NN = 9 buses, NB = 9 branches, and NC = 12 consumers, out of which only NCDR = 7 consumers are participating in the DR program, the matrix DataIN can be written as in Table A1.