Scenario-Based Uncertainty Modeling for Power Management in Islanded Microgrid Using the Mixed-Integer Distributed Ant Colony Optimization

: Reliable droop-controlled islanded microgrids are necessary to expand coverage and maximize renewables potential. Nonetheless, due to uncertainties surrounding renewable generation and load forecast, substantial power mismatch is expected at off-peak hours. Existing energy management systems such as storage and demand response are not equipped to handle a large power mismatch. Hence, utilizing dump loads to consume excess power is a promising solution to keep frequency and voltage within permissible limits during low-load hours. Considering the uncertainty in wind generation and demand forecast during off-peak hours, the dump load allocation problem was modeled within a scenario-based stochastic framework. The multi-objective optimization with uncertainty was formulated to minimize total microgrid cost, maximum voltage error, frequency deviation, and total energy loss. The mixed-integer distributed ant colony optimization was utilized in a massive parallelization framework for the ﬁrst time in microgrids to solve the decomposed deterministic problem of the most probable scenarios. Moreover, a ﬂexible and robust load-ﬂow method called general backward/forward sweep was used to obtain the load-ﬂow solution. The optimization problem was applied to the IEEE 69-bus and 118-bus systems. Furthermore, a cost beneﬁt analysis was provided to highlight the proposed method’s advantage over battery-based power management solutions. Lastly, the obtained results further demonstrate the fundamental role of dump load as power management solution while minimizing costs and energy losses.


Introduction
Over the last two decades, a major shift from conventional forms of power generation toward renewable energy sources (RES) has occurred to minimize greenhouse gas emissions and facilitate the growth of greener and decentralized electricity supply industry [1].Likewise, developments in distributed generation (DG) technology, energy conversion efficiency, government incentive programs, and overall cost reductions have given rise to microgrids (MGs) as a viable solution for future smart grid projects [2].Nonetheless, there are economic and technical challenges that hinder the expansion of islanded MGs (IMGs) which have high RES penetration.The nature of these challenges is associated with the uncertainties surrounding variable renewable generation and demand forecast errors.Conversely, the inability to always match generation with demand, the insufficient transmission capacity to accommodate excess power generation, storage facilities' higher costs, and the need to minimize reliance on fossil fuels are considered the main obstacles that implicate the optimal operation and control of future IMGs [2].
Adequate control of IMGs is necessitated by international standards, such as IEEE Std.1547, to maintain voltage and frequency (V− f ) within acceptable margins [3].Moreover, a reliable control strategy is fundamental to ensure the autonomous participation of all DGs Energies 2023, 16, 4257 2 of 30 to cover any variations in demand.Therefore, droop control is often selected to enable the autonomous sharing of loads within an IMG due to its higher reliability and lower costs compared to other control schemes [4].As a result, IMGs that adopt the latter strategy are referred to as droop-controlled IMGs (DCIMGs).Conversely, with RES penetration levels exceeding 10% in emerging MGs, one major issue remains is how to relieve network congestion and grid integration difficulties.Hence, alleviating such issues arising from large power mismatches in DCIMGs shall maintain V− f safe limits [5].There are several energy management systems (EMS) that handle power deviation problems in DCIMGs.Such EMS include energy storage systems (ESSs) [6], demand response programs [7], and electric vehicle smart charging [8].Nevertheless, electric vehicle and demand response programs are often challenging to coordinate and execute [9], while battery ESSs (BESSs) suffer from high standing losses, environmental degradation, and expensive transportation costs [10].Therefore, EMS strategies that are based on electric vehicles and ESSs are suitable to manage small power mismatches during peak hours as a secondary solution, while different arrangements should be sought for at the event of prolonged wind generation exceeding two load cycles.
Dump loads (DLs) are considered as a good solution to absorb excess power to provide V− f support for synchronous and asynchronous generators [11,12].Moreover, the use of DL as dissipated heat with the help of an electronic load controller (ELC) has been utilized for hydro and wind self-excited induction generators [13,14].Conversely, DLs were put into heating and pumping applications to handle V− f deviations of the system [15,16] while DL application was distributed as an experimental multi-ELC project to provide hot water and heating services for households in a micro-hydro-driven MG [17].However, DLs are considered preliminary when it comes to power management solutions for DCIMGs; that is, studies seeking DL optimization within an IMG framework are scant.On the other hand, various attempts were made to find the optimal location and size of DG and ESS units serving various technical and economic objectives [18][19][20][21].Nevertheless, studies [18][19][20][21] did not account for V− f deviations at low load hours, nor did they consider generation and demand uncertainties.
To accommodate uncertainties in wind power and demand forecast errors within an EMS solution, the optimal droop settings of DGs were optimized in [22] to maximize loadability and minimize fuel costs.Likewise, the total MG cost was minimized in a real-time EMS as a stochastic solution considering demand response and BESS limits [23].Conversely, the authors of the stochastic EMS solution in [24] have sought emission reductions as a stand-alone objective alongside cost minimization while taking BESS as a backup to the IMG.Moreover, the work in [25] has taken similar approach as [24] by optimizing the droop setting of dispatchable DGs but with loadability maximization as an additional objective.In spite of the advancements in the uncertain EMS studies of [22][23][24][25], they did not provide a viable solution to the reliability issues of BESS as an off-peak primary management solution.Furthermore, the works in [22][23][24][25] did not account for the expected high generation/demand mismatch at low-load hours.Consequently, a novel attempt was made in [26] by allocating a DL to eliminate V− f deviations during off-peak hours.Thus, providing a legitimate solution to the efficiency problems of BESS-and demand-response-based EMS.Furthermore, the single DL allocation given in [26] was extended in [27] to allocate multi-DL units across a DCIMG to minimize the same objectives.Nonetheless, the DL problems provided by [26,27] have failed to address MG operational costs, emissions, and losses expected from DL allocation.Likewise, the work in [26][27][28] neither distinguished between dispatchable and non-dispatchable units nor considered uncertainties associated with wind and load powers.
Uncertainty in DCIMG operation roots back to the random parameters that define renewable energy and load diurnal state.The DL allocation into DCIMG is considered as a many-objective non-convex mixed-integer nonlinear programming (MINLP) problem.Moreover, by adding the uncertain non-decision variables into the optimization problem (i.e., wind and load uncertain powers), then the optimization problem is changed from a deterministic optimization problem with certain dimensions to become stochastic optimization problem with random boundaries.To that end, any chosen optimization technique must be able to accurately handle a significant batch of objective function evaluations within a reasonably small calculation time.This is considered fundamental to facilitate real-time EMS solutions within the smallest load cycles possible.Inversely, many of those stochastic EMS solutions [22][23][24][25], have neglected the enormous calculation time that will certainly render such solutions not fit for load cycles of less than 10 min.Ant colony optimization (ACO) is an acclaimed metaheuristic for a variety of engineering problems [29,30].However, the proposed optimization technique, mixed-integer distributed ACO (MIDACO), is based on the extension of ACO into the mixed-integer domains known as mixed-integer ACO (ACOmi) [31] collated with a robust penalty method known as the oracle penalty method (OPM) [32].The multi-objective handling by MIDACO is based on the utopia-nadir balance approach which aims to locate the best point on a Pareto front.Additionally, a massive parallelization strategy is utilized by fine-grained or co-evaluation parallelization to enable thousands of function evaluations in the shortest times possible [33,34].The main two significant advantages of MIDACO are its high-speed and accurate computational abilities.Thus, MIDACO is regarded as the state of the art in evolutionary and swarm intelligence computation for many-objective problems [33,35,36].
Due to the foregoing, it was deduced that ESS-and demand-response-based EMS strategies are not suitable to manage high power mismatch at off-peak hours.Moreover, neglecting uncertainty in wind power and load forecast has a negative impact on future DCIMG planning and optimal set-points.Previous DL studies as a power management solution during low demand hours [26][27][28], albeit scarce, have overlooked fundamental facets of the problem formulation in terms of adequate uncertainty modeling, MG running costs, DG emissions, and total energy losses.In this article, a novel methodology has been proposed to tackle shortfalls in previous studies by formulating a multi-objective stochastic optimization problem across the off-peak hours' horizon.Four objectives were minimized in total, viz., the expected total microgrid cost (TMC), the expected maximum voltage error (MVE), the expected frequency deviation, and the expected total energy loss (TEL).This was done using a state-of-the-art technique, MIDACO, combined with a newly developed and robust load-flow method.The stochastic optimization problem for optimal DL allocation and DG droop selection was modeled on IEEE 69-bus and 118-bus systems.So far and to the best of the authors' knowledge, no previous work has addressed the uncertainty in generation and demand during off-peaks hours within a DCIMG framework.Furthermore, all DCIMG optimization studies heretofore have neglected the impact of uncertainty during low-demand hours and their calculation burden.A summary of the main contributions of this work is as follows:

•
The DL allocation problem was tackled for the first time as a stochastic optimization problem within a scenario-based uncertainty framework.

•
A roulette wheel mechanism was utilized to generate 10,000 scenarios for generation and demand imbalance rather than a set of a few highly probable deterministic mismatch scenarios, as in [28].

•
The DL allocation problem was tackled considering technical, environmental, and economic objectives to provide comprehensive outlook on the optimal IMG operation during off-peak hours.

•
A robust, adaptive, and derivative-free method, the general backward/forward sweep (GBFS), was used to provide load-flow solution.GBFS offers more accurate representation of DG power updates compared to the load-flow method in [28].

•
The proposed optimization technique, MIDACO, has been utilized for first time in a stochastic optimization with uncertainty in microgrids.The proposed technique offers a massive parallelization framework to handle the calculation burden of uncertainties.Furthermore, it offers competitive selection criteria for the non-dominated solution compared to other acclaimed metaheuristics.

•
The advantage of DL-based EMS against battery-based EMS is further demonstrated via cost benefit analysis to provide yearly hot water demand for the microgrid.
This article is organized as follows: in Section 1, a brief introduction to the DL allocation problem in DCIMG with uncertainty in generation and demand.In Section 2, the methodology of droop control, uncertainty modeling, the proposed optimization method, and the proposed load-flow technique are described.In Section 3, the stochastic optimization problem is illustrated.Lastly, in Sections 4 and 5, a thorough discussion of the results and conclusions are presented, respectively.Additionally, to enhance the readability of this article, lists of all acronyms and symbols used herein are given in Tables A1 and A2, respectively, in Appendix A.

Methodology
In this section, the complete methodology for droop control, scenario-based uncertainty, optimization method, and load-flow technique as proposed in this article are elaborated in more detail.

Droop Control and Dump-Load-Based Energy Management System
The majority of DG units within IMG are interfaced as voltage source converters with bidirectional power flow where power electronics enable autonomous power sharing capability by droop control [4].To that end, active and reactive powers are shared by means of active power-frequency (P − f ) and reactive power-voltage (Q − V) droop relationships, respectively.For an inverter-based DG (IBDG) with inductive output impedance, the P − f and Q − V droop relations are given mathematically as in Equations ( 1) and (2) [28]: where |V 0 | and |V i | are the reference and operational voltage at bus i, respectively.Q Gi0 and Q Gi are the reference and generated reactive power at bus i, respectively.f 0 and f are the IMG reference and operational frequency, respectively.P Gi0 and P Gi are the reference and generated active power at bus i, respectively.m pi and n qi are the active and reactive droop coefficients at bus i, respectively.Moreover, the instantaneous droop control action of IBDG forms the foundation of an IMG hierarchal control strategy.The hierarchal control consists of three levels, viz., primary-level, secondary-level, and tertiary-level control.The primary level is achieved within a matter of seconds in response to the variation in load by the action of IBDG droop control.The secondary control, on the other hand, dictates the recovery of voltage and frequency to their nominal values.That is achieved by vertically shifting the curves of P − f and Q − V upwards or downwards depending on the under-or over-generation responses, respectively, by the primary control.The P − f and Q − V curves of an IBDG are depicted in Figure 1.The advantage of DL-based EMS against battery-based EMS is further demonstrated via cost benefit analysis to provide yearly hot water demand for the microgrid.
This article is organized as follows: in Section 1, a brief introduction to the DL allocation problem in DCIMG with uncertainty in generation and demand.In Section 2, the methodology of droop control, uncertainty modeling, the proposed optimization method, and the proposed load-flow technique are described.In Section 3, the stochastic optimization problem is illustrated.Lastly, in Sections 4 and 5, a thorough discussion of the results and conclusions are presented, respectively.Additionally, to enhance the readability of this article, lists of all acronyms and symbols used herein are given in Tables A1 and A2, respectively, in Appendix A.

Methodology
In this section, the complete methodology for droop control, scenario-based uncertainty, optimization method, and load-flow technique as proposed in this article are elaborated in more detail.

Droop Control and Dump-Load-Based Energy Management System
The majority of DG units within IMG are interfaced as voltage source converters with bidirectional power flow where power electronics enable autonomous power sharing capability by droop control [4].To that end, active and reactive powers are shared by means of active power-frequency (-) and reactive power-voltage (-) droop relationships, respectively.For an inverter-based DG (IBDG) with inductive output impedance, the - and - droop relations are given mathematically as in ( 1) and (2) [28]: where | | and | | are the reference and operational voltage at bus , respectively. and  are the reference and generated reactive power at bus , respectively. and  are the IMG reference and operational frequency, respectively. and  are the reference and generated active power at bus , respectively. and  are the active and reactive droop coefficients at bus , respectively.Moreover, the instantaneous droop control action of IBDG forms the foundation of an IMG hierarchal control strategy.The hierarchal control consists of three levels, viz., primarylevel, secondary-level, and tertiary-level control.The primary level is achieved within a matter of seconds in response to the variation in load by the action of IBDG droop control.The secondary control, on the other hand, dictates the recovery of voltage and frequency to their nominal values.That is achieved by vertically shifting the curves of - and - upwards or downwards depending on the under-or over-generation responses, respectively, by the primary control.The - and - curves of an IBDG are depicted in Figure 1.Lastly, tertiary control is dictated by the presences of an MG central controller (MGCC) [24].The MGCC is defined as the complete EMS strategy that is responsible for data collection, analysis, and optimization to determine the IMG optimal set-points for the controllable components.The MGCC is responsible for conducting a complete optimal load-flow of the system considering uncertainties during optimization cycles.The frequency of the optimization cycles is determined by the IMG needs and can be undertaken days, hours, or up to 15 min ahead.The shorter the optimization cycle is, the more reliable and accurate the optimal solutions provided by the MGCC are.The forecasted data presented to the MGCC are used to determine the optimal droop values and DL setting for the next IMG operational interval.This is ideally fed to DG and DL using a non-critical and low-bandwidth communication interface.This type of communication channel is only necessary at the end of each optimization cycle, thus increasing the reliability of the system.Likewise, according to the assumed notion of this study, large power mismatch at off-peak hours could not be sorted by relying only on BESS and demand response programs.
Consequently, IMG planning must be done first to optimize DL location, while a real-time DL-based EMS is utilized thereafter to deliver the optimal DG droop sets and DL size by virtue of MGCC.The idea is to balance the system's V − f fluctuations via dumping the extra power [28].Accordingly, this extra DL power can be applied in pumping and heating applications via electric boilers and water circulation systems.The role of BESS in this EMS strategy-or any other suitable auxiliary power management solution, for that matter-would be to simply manage any light deviation in generation and demand during peak hours.The proposed DL-based EMS of this study is depicted in Figure 2.
Energies 2023, 16, x FOR PEER REVIEW 5 of 31 Lastly, tertiary control is dictated by the presences of an MG central controller (MGCC) [24].The MGCC is defined as the complete EMS strategy that is responsible for data collection, analysis, and optimization to determine the IMG optimal set-points for the controllable components.The MGCC is responsible for conducting a complete optimal load-flow of the system considering uncertainties during optimization cycles.The frequency of the optimization cycles is determined by the IMG needs and can be undertaken days, hours, or up to 15 min ahead.The shorter the optimization cycle is, the more reliable and accurate the optimal solutions provided by the MGCC are.The forecasted data presented to the MGCC are used to determine the optimal droop values and DL setting for the next IMG operational interval.This is ideally fed to DG and DL using a non-critical and low-bandwidth communication interface.This type of communication channel is only necessary at the end of each optimization cycle, thus increasing the reliability of the system.Likewise, according to the assumed notion of this study, large power mismatch at off-peak hours could not be sorted by relying only on BESS and demand response programs.
Consequently, IMG planning must be done first to optimize DL location, while a realtime DL-based EMS is utilized thereafter to deliver the optimal DG droop sets and DL size by virtue of MGCC.The idea is to balance the system's - fluctuations via dumping the extra power [28].Accordingly, this extra DL power can be applied in pumping and heating applications via electric boilers and water circulation systems.The role of BESS in this EMS strategy-or any other suitable auxiliary power management solution, for that matter-would be to simply manage any light deviation in generation and demand during peak hours.The proposed DL-based EMS of this study is depicted in Figure 2. The load model chosen to represent the loads including the DL in the IMG of this study follows the static load model.According to the literature [37], the static load model is sufficient to describe the static components and approximate the dynamic components of the load.To that end, the constant power model was selected for this study by setting all load coefficients to zero [28].The load model chosen to represent the loads including the DL in the IMG of this study follows the static load model.According to the literature [37], the static load model is sufficient to describe the static components and approximate the dynamic components of the load.To that end, the constant power model was selected for this study by setting all load coefficients to zero [28].

Scenario-Based Stochastic Uncertainty Modeling
There are a variety of issues that might affect how definitive a variable is.The majority of these issues are naturally occurring and related to inaccuracies in previous data collection, forecasting future data, and acquiring data specific to that variable [38].As a result, the presence of random variables alters the nature of IMG analysis and planning from a deterministic framework with fixed variables to one of a stochastic nature with uncertain variables.To that aim, planning issues for IMG must integrate a suitable uncertainty analysis tool.When it comes to a stochastic framework, probabilistic modeling is the favored option, where each uncertain variable is described as a probability density function (PDF) [2].These PDFs are used to determine the likelihood that a certain random variable will occur, wherein each PDF is divided into various unique probability levels [38].There are broadly three main probabilistic methods, viz., Monte Carlo simulation (MCS), analytical methods, and approximation methods [39,40].MCS methods are very efficient when problem solution search space is limited, while analytical methods are best when the number of random variables is low.Nonetheless, the optimization problem herein contains a huge solution search space with a large number of random variables.Thus, an approximation method is the best trade-off to attain an expected good solution.Out of the available approximation methods, scenario-based analysis is sought after when the number of random variables is large.To increase the accuracy of the method, the number of scenarios must increase, which will lead to higher computational burden.Nevertheless, this issue is eliminated via the proposed optimization technique; MIDACO can handle a large number of many-objective problems with a huge solution search space in a very fast calculation effort via massive parallelization framework [33,34,36].

Stochastic Load and Wind Modeling
Forecasted, historical, and measured load data are all determining factors toward the degree of load uncertainty.The study herein presumes that historical annual load data are available and that the system peak demand conforms with the IEEE reliability test system for load hourly shape [41,42].Additionally, an hour-by-hour prediction system is used to characterize load uncertainty as a normal-distribution PDF in accordance with standards within the literature [43,44].Typically, loads are distributed around the mean for any given hour, that is, by considering the load's percentage from the overall system demand in the form of an annual cumulative percentage.Additionally, for this study, a seven-hour window from 12 am to 7 am has been chosen to reflect the off-peak hours scenario with higher accuracy.Thus, normal distribution of each load was segmented into fifteen levels with lengths equal to half the standard deviation from the mean; that is, improved uncertainty modeling with a higher number of levels.Conversely, one of the most significant RESs in the world is regarded as wind energy.However, wind power is mostly dependent on wind speeds and, hence, is intrinsically intermittent in nature, much like many other unstable RESs.Therefore, it is essential to account for the influence of wind speed on wind power through suitable probabilistic uncertainty modeling.Wind speed, however, fits the Weibull distribution PDF rather than the normal distribution PDF of a load's probabilistic model [42,45].A wind speed (v) represented by a PDF, referred to herein as φ W , with a scale index c s and shape factor k s to conform with Weibull distribution is given by [45]: For a Weibull distribution, k s and c s may be calculated in a variety of ways.However, k s and c s may be roughly estimated using the following formulas [42, 45,46], provided that the average wind speed µ W and its standard deviation σ W of a certain site are known.
where the symbol Γ() refers to the gamma function.Detailed information regarding the wind speed PDF as well as k s ands c s derivations are found in [46,47].Nevertheless, as per the assumed notion of this study, the wind speed's mean and standard deviation are known in the form of historical annual data for a typical wind farm covering the previous three years [45].Additionally, according to the study's presumption, power is always available during peak hours, whereas wind speeds at higher elevations often tend to progressively rise during off-peak hours [48].Subsequently, a PDF for wind speed is also discretized into several levels or states, where each level reflects a range of projected wind speeds.
To achieve high accuracy for wind speed stochastic modeling, a 30 wind levels/states were constructed from zero wind speed increasing by a 1 m/s increment.Therefore, the probability of any given wind state (W st ) is equal to: where v u st and v l st are the upper and lower limits of wind speed for state W st .For a typical wind turbine (WT), the output power (P W (v st )) as a function of the mean wind speed (v st ) during the state (W st ) is given by [42,45]: where P Wr is the WT rated power.v co and v ci are the cut-off and cut-in wind speeds, respectively, with the former also known as the furling wind speed.v µ is the site's average wind speed.v r is the WT rated wind speed.

Scenario Generation, Reduction, and Aggregation
A sizable number of wind and load states is created as potential scenarios based on wind speed and load forecast PDFs.Moreover, by considering a roulette wheel mechanism (RWM) for this study, each unique PDF will be linked to a separate RWM.Subsequently, 10,000 scenarios are generated, where each (s) scenario is defined by the set (Ω s ) of uncertain variables.As a result, each Ω s consists of the values for each wind and load power random variable as well as their relative probability (Λ i s ).Likewise, the set Ω s is defined in an IMG for an N number of buses as: where P s Wi , P s Li , and Q s Li are the WT, load active, and load reactive powers at bus i during scenario s, respectively.NV is the number of uncertain variables in scenario s. lk and wk are the total number of loads and WT in the network, respectively.Each RWM is fragmented into many slots equal to the number of probability levels of the PDF, whereas an RWM slot length is defined by the random variable's corresponding normalized probability [24,38].Each uncertain variable is obtained from the RWM tool using a random number between 0 and 1.That is, an RWM slot is selected depending on the generated random number and its location on the RWM tool.Subsequently, the uncertain variable tagged with that slot along with its associated normalized probability are chosen.To attain the required number of scenarios, the slots selection process is repeated for all uncertain variables and for every s scenario.Nonetheless, many of the 10,000 created scenarios will have low probability and thus no bearing on the solution.As a result, eliminating repeated and low-probability scenarios must take place by selecting fewer and dissimilar scenarios with the highest probabilities.A reduced number of scenarios (NR) for this work was selected as 20 to conform with current norms in uncertainty modeling techniques [24,49].Lastly, each of the NR scenarios of the stochastic planning problem may then be utilized to construct a separate deterministic optimization problem.As a result, the probability of each deterministic problem's optimal solution within a scenario s will be equal to the convolved probability of all uncertain variables [50].It is important to note that this study makes the premise that the events linked to wind speed and load forecast are uncorrelated [42].In other words, neither the occurrence of a certain load state nor the occurrence of a particular wind speed state is affected by the other.This presumption is consistent with the notion herein that holds that generation will almost certainly be higher than demand during off-peak hours.Therefore, for a scenario s, the convolved normalized probability Λ N s is obtained by: Additionally, for each s scenario's deterministic problem, an aggregation technique based on weighted sum is used; that is, by accumulating the optimal solutions considering the importance of each scenario based on its higher probability effect.Thereby, the expected value for the objective function ( Fi (x)) of the stochastic problem taking into account all NR scenarios is provided by [38,49]:

The Proposed Optimization Algorithm
As mentioned in the introduction, the proposed DL stochastic optimization problem is a many-objective non-convex MINLP.This implies that the solution search space is challenging and vast.Such problems are known not to have a solution in polynomial time, making deterministic solvers inapplicable.Hence, to solve the stochastic many-objectives problem, a state-of-the-art metaheuristic algorithm, MIDACO, is proposed herein.This evolutionary algorithm is considered as an advanced hybrid optimization technique.It employs different heuristics for better exploration and a back-tracking local deterministic solver for enhanced exploitation [33].The main components of the MIDACO algorithm are the ACOmi for constructing iterates, or ants, using PDFs and the OPM for evaluating penalties from constraint violations [33].By employing the extended ACO for mixedinteger domains, the proposed method relies on a multi-kernel Gaussian PDF to construct solutions instead of a pheromone table as in the primordial ACO.For single-objective MINLP handling by MIDACO, three parameters are very influential, viz., ANTS, KERNEL, and ORACLE [33].That is, for constructing solutions, the method starts with a dynamic population of ants (N pop ) within a (k r ) number of kernels to generate and fine-tune ants as they get evaluated according to their oracle penalty value [31].A user-desired value for the objective function known as the oracle (Ω) estimates the penalty function value within OPM [32].Conversely, for multi-objective MINLPs, MIDACO uses the utopia-nadir balance concept to decompose the original problem into a series of single-objective sub-problems, each solved in the j-th dimension [51].This approach explores the Pareto front much faster and more efficiently leading to a better convergence of the multi-objective optimization.This differs from a non-dominated sorting approach that gives the entire Pareto front an equal importance making convergence much slower [33,51].Another advantage for the utopia-nadir balance approach is that, eventually, the best equally balanced solution is selected [33,51].Thus, an additional multi-criteria decision approach is not required in MIDACO.Accordingly, for any mixed-integer solution x that belongs to set of feasible solutions F, the utopia (U i ) indicates the fittest overall value for an objective function F i (x) Energies 2023, 16, 4257 9 of 30 whereas a nadir (N i ) is by far the least fit F i (x) value that relates to different utopia as follows [51]: Once the utopia-nadir information is known, then the solution x weighted d j i (x) and average D j (x) distances are calculated for an M-objectives MINLP.This is done for each objective i within the j-th dimension as follows [51]: Subsequently, by having d j i (x) and D j (x) for each sub-problem as well as the utopia- nadir information, a scalar function known as the balance function is created [51]: Each j-th sub-problem is solved as a single-objective problem using the target function T j (x) as defined hereafter [51]: Eventually, the original problem is re-constructed again using the set of target functions {T 1 , T 2 , . . . ,T M } evaluated in series or in parallel to attain the Pareto front.For CPU-time- intensive problems, such as the stochastic one proposed in this article, MIDACO offers a very efficient parallelization strategy.This strategy, known as fine-grained parallelization, handles the execution of objective and constraint functions in each individual ACOmi instance in parallel [33].To fine-tune multi-objective optimization, MIDACO utilizes the parameters, viz., BALANCE, EPSILON, and PARETOMAX.As for the former, it dictates the Pareto front search direction, thus expediting the convergence.The latter two parameters influence the precision and number of Pareto point collection, respectively.Other parameters such SEED and ACCURACY enhance the solution's global optimality and suitability, respectively [33].The proposed optimization technique flow chart is depicted in Figure 3.

The Proposed Load-Flow Method
There exist different load-flow techniques such as Gauss-Seidel, Newton-Raphson, and backward/forward sweep (BFS) [28].Nonetheless, many of them are not wellequipped to provide load-flow for IMG because of the low X/R ratio in distribution systems.Moreover, IMG has a variable system frequency and no constant slack bus voltage.Therefore, a special BFS (SBFS) method was proposed in [28] to account for droop Equations ( 1) and (2).However, the reactive power update in SBFS does not account for local DG voltage measurements.These local measurements are necessary for IMG with little or weak communication infrastructure.Conversely, GBFS was proposed in [52] which has two main stages, viz., BFS stage and update stage.These will be described briefly as follows.

The Proposed Load-Flow Method
There exist different load-flow techniques such as Gauss-Seidel, Newton-Raphson, and backward/forward sweep (BFS) [28].Nonetheless, many of them are not well-equipped to provide load-flow for IMG because of the low X/R ratio in distribution systems.Moreover, IMG has a variable system frequency and no constant slack bus voltage.Therefore, a special BFS (SBFS) method was proposed in [28] to account for droop Equations ( 1) and (2).However, the reactive power update in SBFS does not account for local DG voltage measurements.These local measurements are necessary for IMG with little or weak communication infrastructure.Conversely, GBFS was proposed in [52] which has two main stages, viz., BFS stage and update stage.These will be described briefly as follows.

Backward/Forward Sweep Stage
At this stage, all IMG variables are initialized, and all system quantities are calculated.The value 1 ∠ 0 • p.u. is given to voltages of all system buses including the virtual bus (VB).A VB mimics the behavior of a slack bus in grid-connected load flow; however, it has a variable voltage.Subsequently, inject currents as well as branch currents for all buses are calculated in forward sweep.Inversely, voltages on each system bus, except VB, are obtained in the backward sweep as follows [52]: where, for a radial system with n buses, and [B i ] are all column vectors of size n − 1 by 1 and correspond to at bus i as the apparent power inject, initial bus voltage, current inject, and branch current, respectively.
[BIBC] and [BCBV] are matrices of size n − 1 by n − 1 for branch inject-branch current as given in [53] and branch current-bus voltage as given in [28], respectively.[V 1 ] and [V in ] are column vectors of size n − 1 by 1 for VB voltage and new bus voltages, respectively.ζ 1 is the dynamic damping factor to suppress the oscillations in the system voltage error vector ∆V in as obtained stochastically by MIDACO in [52].V in is the new bus voltage vector as attained using ζ 1 and another sweep for V i .∆ f and ∆V 1 are the system frequency and VB voltage deviations, respectively.P G1 and Q G1 are VB-generated active and reactive powers, respectively.In the case that the VB has no DG unit connected to it, these values are set to zero.V 1 •B * 1 is the total apparent power exchanged at the VB.m pT and n qT are the system effective active and reactive droop coefficients at VB as defined in [28], and are given by:

The Update Stage
In this stage, the VB voltage and system frequency are updated along with the DGs' active and reactive powers as follows: where f c+1 and f c are, respectively, the system frequency at the c + 1 and c-th iterations.V 1c and V 1c+1 are, respectively, the VB voltage at the c + 1 and c-th iterations.ζ 2 is the dynamic damping factor to suppress the oscillations in the VB voltage error vector |∆V 1 | as obtained stochastically by MIDACO in [52].Calculation of ζ 1 and ζ 2 is presented in more detail in the following section.P Gi0 and Q Gi0 are the active and reactive power reference points at bus i, respectively.GK and N are sets for system generating buses and all system buses, respectively.Another distinct feature of GBFS compared to other BFS-based methods such as modified and nested BFS [52] is having one loop.This calculation loop with the counter c helps in minimizing computation burden by eliminating many internal loops.Moreover, according to GBFS implementation, the reactive power update of DG units considers the local voltage at each generating bus to mimic the nominal voltage |V 0 | recovery of DG units.However, such an approach will often cause divergence problems for ill-conditioned systems.Such issues are due to IMG with lower droop settings, high generation/demand mismatch, or higher line impedance.As a result, the correction vector for reactive power was proposed by GBFS in [52] and denoted here as γ i .The aim of γ i is to eliminate the reactive power error |∆Q Gi | in the IMG.
where Q c is the IMG correction factor based on average reactive power in the IMG, and is provided by [52]: Furthermore, to ensure adequate reactive power correction, β, a binary constant, was added to (32) to disable or enable the reactive power correction based on IMG needs [52]: Subsequently, the reactive power reference is adjusted to account for any power limit violations.Thus, a new reactive power value (Q Gi ) is used and given by [52]: If β = 0, this implies γ i is disabled and no change to reactive power is performed, i.e., Q Gi = Q Gi .In addition, due to frequency update, line impedance (Z i ) is updated.
where R i and X i are the resistance and reactance of the line connecting bus i and i + 1, respectively.Lastly, GBFS exits the loop when the condition for convergence is satisfied; in other words, by attaining an objective function (F(x 0 )) value below the voltage error threshold (ε Th ), as will be explained in the next section.That is, GBFS terminates when F(x 0 ) < ε Th , noting that ε Th = 10 −8 for all investigated cases in this article.A flow chart of GBFS is depicted in Figure 4.

Optimization Problem Formulation
In this section, the objective functions and constraints for the stochastic DL optimization considering uncertainty in wind and load power are described in detail.Energies 2023, 16, 4257 13 of 30

Optimization Problem Formulation
In this section, the objective functions and constraints for the stochastic DL optimization considering uncertainty in wind and load power are described in detail.

Objective Function for Dynamic Damping Factors in GBFS
The objective in the optimization problem of GBFS is to obtain a concurrent minimization of the two main voltage error vectors ∆V in and |∆V 1 |.This is achieved by stochastically adjusting the dynamic damping factors ζ 1 and ζ 2 .A weighted sum approach was elected to attain the dynamic adjustment of ζ 1 and ζ 2 without computational overhead and extra run time.This is attributed to weighted sum benefit in simplifying and transforming the multi-objective problem into a single objective [40,54].Furthermore, due to the existing knowledge of the objective function desired threshold, i.e., ε Th , the search efforts in MIDACO are greatly reduced.Therefore, the prosed optimization method is tuned accordingly using the parameters FOCUS and FSTOP.The former parameter can guide the search efforts of MIDACO towards a local region, which is known to have the desired threshold, whereas the latter parameter will terminate the algorithm once the objective function exceeds that threshold [33].Subsequently, as in Equation (37), the GBFS objective function to obtain ζ 1 and ζ 2 is provided as: where w 1 and w 2 are weights of ∆V in and |∆V 1 | and equal to 1 for enhanced conver- gence [52].x 0 is the decision variable for the objective function F(x 0 ) in the GBFS implementation.

Objective Function for Dump Load Stochastic Optimization
In this study, four objective functions are defined as the expected total microgrid cost (TMC), the expected maximum voltage error (MVE), the expected frequency deviation (|∆ f |), and the expected total energy loss (TEL).Those objectives are constructed consider- ing all reduced scenarios NR spanning across the chosen off-peak horizon (i.e., 12-7 am).Consequently, an objective function F i (x 1 ) with the decision variable x 1 is defined as: where P DL , Q DL , are, respectively, the dump load's active and reactive powers.N DL is the DL bus location.mn DL refers to the optimal droop setting considering the allocation of DL.The value of mn DL influences DG droops as follows: Having reactive droop slightly above the active droop value is common and known to further improve the load-flow convergence [19,55].Moreover, an aggregated sum for all NR deterministic problems within each hour that belong to the set (H) defining the off-peak horizon is given as total expected objective function Fi (x 1 ): where H is the total number of hours with a set H.An expected non-dominated solution using the Pareto front optimization technique, MIDACO, is used to find x 1 that satisfies all constraints and minimizes Equation (45).The stochastic objective function Fi (x 1 ) is solved by aggregating all deterministic optimization problems' objective function (F h,s i (x 1 )) in a solution matrix (SM) of size H × NR.The four objectives are detailed as follows:

Total Microgrid Cost
The four base cost functions making up the TMC objective described in this work which are also dispersed throughout the off-peak horizon for each reduced scenario are given as fuel cost (FC), maintenance cost (MC), emissions cost (EC), and technical costs (TC).Such a combination is the preferred idea when considering the typical technical, environmental, and economic goals for any IMG optimization task.Furthermore, we implicitly reduce emissions and losses by including active power generation from all dispatchable DGs as a component in the TMC objective.Therefore, it would be unnecessary to include a separate objective to address emissions due to the linear relation between active power and emissions.Nevertheless, network losses should be tackled as a separate objective since the relationship between apparent power and branch current components is not linear.Conversely, the TMC cost function considered herein is based on MG operational and running costs.Thus, it does not include capital or standing costs of the MG such as DL or storage installation costs.It is a common practice in the literature to neglect other form of costs apart from running costs, as they are mainly dependent on dispatchable unit generation [24,[56][57][58].Moreover, as per the notion of this study, DL application should at least offset the costs of storage installation for the same purpose.Nonetheless, a complete cost benefit analysis of DL and battery EMS is given separately in the discussion section.Given the foregoing, TMC is obtained as [56,57]: where for a scenario s during the off-peak hour h, the total microgrid, fuel, maintenance, emissions, and technical costs are given by TMC s h , FC s h , MC s h , EC s h , TC s h , respectively.The dispatchable DG's generated active power is given by P h,s Gi .In addition, the reactive and frequency costs are given by RC s h and FRC s h , respectively.η P is the efficiency in fuel consumption by a dispatchable DG.Ψ emis is the dispatchable DG's emissions rate.ψ f uel , ψ main , and ψ emis are, respectively, the fuel, maintenance, and emissions cost coefficients, noting that the technical cost is about MG reliability and stability as represented by RC s h and FRC s h costs.Reactive power is not connected with fuel, albeit leading to increased losses and penalties for generating reactive power [57].Conversely, when frequency deviations are out of the permissible range, the FRC s h costs are required to ensure MG supply quality.Although technical costs often refer to V− f deviations concurrently, in this article, only frequency penalty costs were elected.This is due to the fundamental role of frequency in IMG supply quality and the stringent permissible tolerance for frequency in the system (±0.4%)[59].Accordingly, RC s h and FRC s h are provided by [56,57]: where Ψ var is the dispatchable DG's reactive power coefficient [60].f h,s ss is the IMG steady state frequency at scenario s during the off-peak hour h.ψ f req is a penalty cost coefficient for frequency [56,59].Q h,s Gi is the reactive generated power by the dispatchable DG at scenario s during the off-peak hour h.Additionally, DCIMG quality of supply as well as safe operating constraints necessitate the provision of technical costs [56].It is known that a renewable DG does not cause emissions since it consumes zero fuel.Furthermore, each WT in this article was assumed in maximum power point tracking (MPPT) control, acting as an induction generator with a 0.9 leading power factor [24,61].

Maximum Voltage Error
In an IMG, a good stability and reliability indicator is the voltage.Therefore, many DG allocation studies have opted for stability indices such as total voltage variations and voltage stability index.Nonetheless, the allocation herein is for DL and not DGs.Thus, a flattened voltage profile is desired to ensure system voltages as close to nominal as possible.As a result, MVE was elected to attain voltage error minimization across the IMG, and is provided by [38]: where V h,s in is bus i voltage at scenario s during the off-peak hour h.

Frequency Deviation
On the other hand, the expected IMG frequency error is achieved as follows: where V h,s 1 , B h,s 1 , and P h,s G1 are, respectively, the voltage, branch current, and the active power at the VB within a scenario s during the off-peak hour h.

Total Energy Loss
As for the fourth objective, the expected total energy loss, it is calculated as follows: where, for each scenario s during the off-peak hour h, the active power loss and the branch current are given by P h,s loss and B h,s i , respectively.t h is the time duration at each s scenario which is equal to one hour.

Constraints for Dump Load Stochastic Optimization
In any IMG, there are different technical constraints that must be satisfied to ensure adequate and stable operation [3,62].Moreover, for each deterministic optimization problem, load flow must converge for any given scenario s and off-peak hour h.This very stable behavior signals power balance constraint adherence.Moreover, for each objective function evaluation the frequency, thermal, bus voltage, and dispatchable DG power limits must be adhered concurrently.These constraint functions (g h,s i (x 1 )) are handled during each scenario s at any given off-peak hour h as follows.
Bus i's voltage limits: Thermal limits: Steady-state frequency limits: Dispatchable DG power output limits: According to the presumption herein, non-dispatchable units operate with an MPPT algorithm (i.e., in PQ mode).Therefore, in any given scenario, the WT will always remain within power limits and thus were excluded from the constraint-handling function.Inversely, the limits for the decision variable x must be satisfied once at every function evaluation after the aggregated effect of all NR scenarios during low load hours (i.e., ∀h ∈ H) are taken into consideration.Thus, neither DL location, DL size, nor DG droop gains are affected by scenario or hour change.These are provided as follows [28]: DL size limits: 0.002 Droop coefficient limits: 10 −4 ≤ mn DL ≤ 0.05 (63) Noteworthy is that the per-unit system was used to write all numerical data using system nominal frequency as 50 Hz, power base 500 kVA, and 12.66 kV and 11 kV as voltage base for the 69-bus and 118-bus systems, respectively.

Results and Discussion
The case studies considered for the IMG in this work are the IEEE 69-bus and 118-bus systems which are shown in Figure 5a,b, respectively.
System line and load data were obtained from [63,64] for the 69-bus and 118-bus systems, respectively, while the 69-bus and 118-bus system generation bus locations were obtained from [26,65], respectively.All DGs and MG ratings are given in Table 1 for the dispatchable and non-dispatchable units which were obtained from [24,58], respectively.Furthermore, it was assumed that all non-dispatchable DGs were WT units, while dispatchable DGs, on the other hand, were installed as natural gas turbine (NGT) units.The DG arrangements for the no DL case (i.e., the base case) are given, respectively, in Tables 2 and 3 for the IEEE 69-bus and 118-bus systems.The pre-islanding reference generation for all NGT units were assumed as 2.545 + j1.909 p.u. at 0.8 lagging power factor, while WT units, on the other hand, were rated at 0.5 MW capacity with 0.9 leading power factor operation [24].Due to the expected large power mismatch as per the notion of this study, two and four identical DLs were considered in the stochastic optimization problem for the 69-bus and 118-bus systems, respectively.The DLs were rated at 500 kVA each and allocated simultaneously with the same size and to the same location for both test systems.
Droop coefficient limits: Noteworthy is that the per-unit system was used to write all numerical data us system nominal frequency as 50 Hz, power base 500 kVA, and 12.66 kV and 11 kV voltage base for the 69-bus and 118-bus systems, respectively.

Results and Discussion
The case studies considered for the IMG in this work are the IEEE 69-bus and 1 bus systems which are shown in Figure 5a,b, respectively.System line and load data were obtained from [63] and [64] for the 69-bus and 118bus systems, respectively, while the 69-bus and 118-bus system generation bus locations were obtained from [26] and [65], respectively.All DGs and MG ratings are given in Table 1 for the dispatchable and non-dispatchable units which were obtained from [58] and [24], respectively.Furthermore, it was assumed that all non-dispatchable DGs were WT units, while dispatchable DGs, on the other hand, were installed as natural gas turbine (NGT) units.The DG arrangements for the no DL case (i.e., the base case) are given, respectively, in Tables 2 and 3 for the IEEE 69-bus and 118-bus systems.The pre-islanding reference generation for all NGT units were assumed as 2.545 + 1.909 p.u. at 0.8 lagging power factor, while WT units, on the other hand, were rated at 0.5 MW capacity with 0.9 leading power factor operation [24].Due to the expected large power mismatch as per the notion of this study, two and four identical DLs were considered in the stochastic optimization problem for the 69-bus and 118-bus systems, respectively.The DLs were rated at 500 kVA each and allocated simultaneously with the same size and to the same location for both test systems.The simulation was modeled and executed in MATLAB ®® environment on hardware comprising an Intel core i7 2.6 GHz 8 GB RAM.Furthermore, the stochastic optimization problem was initialized for each deterministic problem with 0, 0, and 10 9 values for ANTS, KERNEL, and ORACLE, respectively.Meanwhile, PARETOMAX, BALANCE, and EPSILON were set to 1000, 0, and 0.01, respectively.The parameter selection is recommended at the respective values to allow for the best mix between exploration, exploitation, and speed [33].Lastly, to initialize the GBFS optimization problem [52,66], the parameters FSTOP and FOCUS were set to 10 −8 and 100, respectively.A hard-stopping criterion is chosen for MIDACO by fixing the maximum number for function evaluations, known as MAXEVAL, to 10,000.

Multi-Objective Optimization
As discussed previously, by concurrently considering all 20 reduced scenarios as per the desired accuracy rate, the formulated many-objective stochastic problem was created to minimize the expected TMC, MVE, |∆ f |, and TEL.Moreover, the chosen sampling rate for uncertainty modeling has resulted in NR deterministic optimization problems, each with constant load and wind power.It should be noted that the locations of DGs within each test IMG were selected to minimize losses while meeting network demand without relying on any external sources [65].Based on the forgoing, a non-dominated solution that considers all scenarios during the off-peak horizon for each expected objective function in the stochastic optimization problem is given in Table 4. Similarly, the utopia-nadir balance approach has resulted in finding the optimal non-dominated solution at the center of the Pareto front to satisfy all objectives.The non-dominated solution is highlighted in green as depicted in the Pareto front for the IEEE 69-bus and 118-bus systems in Figure 6a,b, respectively.According to results, the Energies 2023, 16, 4257 20 of 30 consideration of uncertainty in load forecast and wind speed in the DL allocation problem has resulted in costs savings, better voltage profile, and optimal supply quality if compared to the base case (i.e., without DL allocation while using droop sets of Table 2 or Table 3).Taking the TMC for the base case in particular, a significant increase in MG costs is observed against the DL allocation cases of both test systems.This is attributed to the higher fuel consumption and emissions by the dispatchable DGs in the IMG without DL.Furthermore, higher technical costs were incurred for operating the IMG without DL allocation, that is, to guarantee supply quality.Conversely, a smoother frequency profile is observed for both test systems utilizing the DL as a power management solution according to the improved expected |∆ f | during the off-peak hours as provided in Table 4. First step size only for |∆|.
Similarly, the utopia-nadir balance approach has resulted in finding the optimal nondominated solution at the center of the Pareto front to satisfy all objectives.The non-dominated solution is highlighted in green as depicted in the Pareto front for the IEEE 69-bus and 118-bus systems in Figure 6a,b, respectively.According to results, the consideration of uncertainty in load forecast and wind speed in the DL allocation problem has resulted in costs savings, better voltage profile, and optimal supply quality if compared to the base case (i.e., without DL allocation while using droop sets of Tables 2 or 3).Taking the TMC for the base case in particular, a significant increase in MG costs is observed against the DL allocation cases of both test systems.This is attributed to the higher fuel consumption and emissions by the dispatchable DGs in the IMG without DL.Furthermore, higher technical costs were incurred for operating the IMG without DL allocation, that is, to guarantee supply quality.Conversely, a smoother frequency profile is observed for both test systems utilizing the DL as a power management solution according to the improved expected |∆| during the off-peak hours as provided in Table 4.Moreover, flatter voltage profiles were achieved by DL's optimal sizing as well as DG's optimal droop setting as depicted in Figures 7 and 8 for the 69-bus and 118-bus systems, respectively.Meanwhile, the advantage of DL allocation in DCIMG has significantly reduced the expected averaged over-frequency by 0.1847 p.u. and 0.2114 p.u. for the 69- Moreover, flatter voltage profiles were achieved by DL's optimal sizing as well as DG's optimal droop setting as depicted in Figures 7 and 8 for the 69-bus and 118-bus systems, respectively.Meanwhile, the advantage of DL allocation in DCIMG has significantly reduced the expected averaged over-frequency by 0.1847 p.u. and 0.2114 p.u. for the 69-bus and 118-bus systems, respectively, against the base case.The advancements in frequency profiles are provided in Figures 9 and 10 for the 69-bus and 118-bus systems, respectively.As for the fourth objective (i.e., TEL), it is generally affected by the pre-existing reactive power mismatch in the IMG.In one hand, for an IMG with relatively small reactive power mismatch, the resultant TEL will be greater after the DL inclusion.On the other hand, when the pre-existing reactive power mismatch is relatively large, the resultant TEL will be lower after DL inclusion into the MG.This situation is explained after examining the reactive power compensation in predominantly capacitive networks and the role of DL allocation in IMGs.Notwithstanding the size of the existing reactive power mismatch before DL allocation into the network, considering TEL criterion within the many-objectives DL allocation will ensure loss reduction.Those refer to losses that would have been As for the fourth objective (i.e., TEL), it is generally affected by the pre-existing reactive power mismatch in the IMG.In one hand, for an IMG with relatively small reactive power mismatch, the resultant TEL will be greater after the DL inclusion.On the other hand, when the pre-existing reactive power mismatch is relatively large, the resultant TEL will be lower after DL inclusion into the MG.This situation is explained after examining the reactive power compensation in predominantly capacitive networks and the role of DL allocation in IMGs.Notwithstanding the size of the existing reactive power mismatch before DL allocation into the network, considering TEL criterion within the many-objectives DL allocation will ensure loss reduction.Those refer to losses that would have been otherwise incurred by the general DL utilization in distribution networks as seen from the single-and two-objectives results for the DL allocation presented in [28].

Comparison with Other Metaheuristics
The of the expected non-dominated solution offered by the proposed method herein is compared against other acclaimed evolutionary and swarm intelligence techniques, viz., multi-objective genetic algorithm (MOGA) [67], multi-objective particle swarm optimization (MOPSO) [68,69], and the non-dominated sorting genetic algorithm (NSGA-II) [70,71].The parameters for MOPSO, NSGA-II, and MOGA, were adopted from [28].Accordingly, 11 parameters were initialized for MOPSO, namely, population and repository size both 100, leader and deletion selection pressures both 2, and starting and ending inertia weights values 0.5 and 0.001, respectively.The cognitive and social learning coefficients were 0.1 and 0.2, respectively, grid per dimension 7, and mutation rate 0.1.As for NSGA-II, four parameters were set: population size 100, distribution index for mutation 20, distribution index for crossover 100, and mutation probability 0.25.Lastly, MOGA's population size was 100, crossover probability 0.8, and mutation probability 0.001.The results obtained by each metaheuristic technique are given in Table 5.According to results, an overall best for TMC and TEL was found by MIDACO for the 69-bus system, whereas MIDACO's TEL for the 118-bus system was substantially smaller if considered against the other methods.Inversely, taking the |∆ f | obtained by all methods as reference, the superiority of the proposed method is exemplified by the quality of the non-dominated solution achieved.Moreover, both MIDACO and NSGA-II have returned close values for TMC; nonetheless, the attained objectives of TEL and MVE were worst from NSGA-II compared to MIDACO's.Despite having much lower maximum function evaluation parameter, i.e., MAXEVAL, by MOPSO, MOGA, and NSGA-II, their recorded calculation times, albeit utilizing a parallelization strategy, were quite long, rendering them not practical.Noteworthy is that the DL locations attained by the four metaheuristics are not similar.This is attributed to the challenging and complex nature of many-objective optimization problems.Nonetheless, the handling of integer domains within MIDACO is enhanced to avoid premature convergence of the internal ACOmi instances [31].On the other hand, the competitive edge offered by MIDACO's fine-grained parallelization strategy was demonstrated by having much lower calculation times despite the huge number of function evaluations.The speed and accuracy advantage of MIDACO against the other metaheuristics proposed herein, shall indeed facilitate a competitive solution with decent and practical speeds.This is necessary for real-time optimization of challenging and stochastic MINLPs forming future IMGs.

Cost Benefit Analysis
As discussed in the optimization problem formulation, DL allocation is expected to reduce total running costs of the MG; that is, by considering fuel, maintenance, emissions, and technical cost reduction because of optimal dispatchable units' generation.Nonetheless, the assumption was that DL costs are offset by utilization of the power into useful operations rather than having it stored or dissipated.Therefore, to demonstrate the effectiveness of the proposed DL-based EMS (DLEMS), a cost benefit analysis (CBA) was conducted.The attained results of this analysis will signify the DL solution's advantage against a batterybased EMS (BEMS).By considering CBA for DLEMS and BEMS, the costs associated with DL and BESS installations will be analyzed, respectively.Thus, to put the CBA into practical perspective, the advantage of using DLEMS against BEMS considering power regulation in primary off-peak hours is investigated.This implies that the CBA has considered two unique EMS approaches to provide a system's hot water demand, namely, DLEMS and BEMS.In the former EMS approach, the overgeneration of excess power is absorbed by DL, while a portion of the system's demand for water-heating boilers and pumps is covered by the DL active and reactive powers, respectively [27].Furthermore, the remainder of hot water demand is supplied by electric-powered boilers from non-renewable sources.Conversely, for BEMS implementation, BESS is utilized to store the excess overgeneration power mismatch, while the system demand for hot water is supplied by natural-gas-fired boilers and pumps.
The motivation for the CBA herein is to provide a comprehensive economic dimension as a reflection from the DL solution to the V− f regulation in highly penetrated IMG.Therefore, a more realistic approach was adopted in this study contrary to that provided in [27].This was performed by considering all capital, fuel, maintenance, and running costs for hot water system installation and maintenance.Moreover, by considering uncertainties in wind and load powers, the amount of available power to be dumped was determined using stochastic scenario-based modeling.Additionally, in addition to assuming the levelized cost of energy (LCOE) for storage solutions in [27], the CBA herein has considered LCOE for gas and electric water-heating systems.Accordingly, when it comes to investment analysis for costs of producing energy, LCOE is widely used in the industry and academia to provide sufficient investment return overview for different energy technologies [72].Subsequently, the CBA is used to measure the cost-effectiveness of relieving high power congestion during off-peak hours with a DL solution rather than a storage solution.As a result, the analysis was conducted assuming two LCOEs for electric boilers, viz., from renewable and non-renewable sourced electricity.Furthermore, the provision for the total demand for hot water is scheduled during low-load hours by storing the hot water in dedicated cylinders to be used later via on-demand water circulation systems [27,73].The hot water volume for total daily demand (V tot hw ) in this article is assumed on average (considering winter and summer days) as 817.06 m 3 and 1293.67 m 3 for 69-bus and 118-bus systems, respectively [24,73].Similarly, the hot water volume (V hw ) equation is derived from the obtainable energy to produce heat as follows [73]: where P hw is the power utilized by the boiler for heating water during h hours.ρ w is the density of water.C w is the specific heat of water.η hw is the boiler efficiency for heating water and assumed as 0.99 and 0.80 for electric (η ele hw ) and gas (η gas hw ) boilers, respectively.
Energies 2023, 16,4257 ∆T is the difference between the set-point temperature T st , assumed at 60 • C, and the inlet temperature T in , assumed at 10 • C, and is given by [27,73]: Due to the foregoing, the utilization of electric or gas boilers for water heating will accumulate the following costs: where HC g hw and HC e hw are the BEMS and DLEMS costs for heating water, respectively.P gas hw and P ele hw are the power needed by gas and electric boilers to cover the total system demand for hot water, respectively.Note that P DL was deducted from P ele hw to represent renewable-based electric boiler costs.ψ gas LCOE is the gas boiler LCOE coefficient, which is valued at 57. 13 [72].Noteworthy is that the LCOE from [72] does not account for long distance thermal losses, nor the cost associated with long-distance thermal installation and maintenance.Nonetheless, the costs incurred from long-distance thermal operation are assumed the same for both EMSs in question and was omitted form the comparison.As for a BEMS strategy, P DL (i.e., the excess power) would require storage in an appropriate BESS.This will lead to extra costs which are called storage costs (SC BESS ).The value of SC BESS is dependent on the required storage power and the LCOE for the corresponding storage technology [74].Accordingly, two LCOE coefficients for storage technologies are considered herein to attain a better comparison against DLEMS.Thus, two BESS types are used in the BEMS implementation, viz., Nickel-Cadmium (Ni-Cd) and Lithium-ion (Li-ion).The cost coefficients for Ni-Cd (ψ Ni LCOE ) and Li-ion (ψ Li LCOE ) technologies are assumed as 691.06 USD/MWh (641.1 EUR/MWh) and 658.61USD/MWh (611 EUR/MWh), respectively [74].
All technical parameters and cost coefficients used in the CBA of this study are provided in Table 6.Lastly, the average daily hot water demand V tot hw is used to calculate the total cost for hot water per calendar year (HC tot hw ).This is done for both DLEMS and BEMS implementations as follows:  Provided in Table 7 are the results of the CBA for the 69-bus and 118-bus systems.According to the results, the costs associated with DLEMS to cover yearly hot water demand were much lower than they were for BEMS.This is true for both BESS technologies implemented, i.e., Li-ion and Ni-Cd.Likewise, the estimated costs for hot water demand per calendar year using the Ni-Cd-based BEMS were USD 2,535,058.14and USD 3,746,532.52for the 69-bus and 118-bus systems, respectively.Similarly, using the cheaper BESS option (i.e., Li-ion) in the BEMS implementation did not result in significant cost reductions, as they were USD 2,474,229.79and USD 3,662,771.75for the 69-bus and 118-bus systems, respectively.On the other hand, using DLEMS for the same purpose (i.e., heating systems hot water demand per calendar year) has resulted in significant and expected reductions in costs valued at USD 1,159,667.99and USD 1,850,319.83for the 69-bus and 118-bus systems, respectively.The major cost factors affecting this huge difference between DLEMS and BEMS are the storage and gas boiler costs.The latter is attributed to the higher power required for gas boilers to produce the same volume of hot water if compared with the more efficient electric boilers.The amount of power required by gas boilers were 7.4251 MW and 11.7563 MW for the 69-bus and 118-bus systems, respectively.Conversely, the power required by electric boilers were just at 6 MW and 9.5 MW for 69-bus and 118-bus systems, respectively.As for the former cost factor, BESS are still expensive considering different variables such as installation, end-of-life, and replacement of batteries [74].Moreover, despite that storage costs are dropping as technology advances, the overall cost elements of BESS are still not accounted for in medium-and long-term solutions [74].
Contrariwise, the notion of this study conforms with current norms that renewable energy penetration levels in DCIMG are rising in the foreseeable future.This will lead to higher periods of excess power mismatch that BEMS are not equipped to handle economically.Therefore, investing in more BESSs is not cost-effective as demonstrated by this CBA.Inversely, taking the DLEMS implementation for the 118-bus systems as an example, significant savings are achieved against the BEMS implementation.This is translated to huge saving of USD 1,812,452.0and USD 1,896,212.69for not using the Li-ion and Ni-Cd storage options, respectively.Furthermore, as mentioned previously, renewable energy is expected to increase drastically, leading to further reductions in electric boilers over all costs.Contrariwise, excess renewable power at off-peak hours will lead to more curtailment of renewable energy, higher costs of storage solutions, and loss of renewable transmission links.Therefore, an investment in DLEMS as a power management solution instead of BEMS is a more cost-effective venture.

Conclusions
In this article, the uncertainty surrounding demand forecast and renewable generation within the DL allocation problem to regulate V− f in highly penetrated DCIMG was investigated.A scenario-based stochastic modeling of uncertainty was used to model load forecast errors and wind generation.An RWM tool was used with higher accuracy to generate 10,000 scenarios by segmenting the load and wind PDFs into 15 and 30 probability levels, respectively.The optimal DL location and size as well as optimal DG droops were determined as a many-objectives problem to account for expected TMC, MVE, |∆ f |, and TEL.This was achieved using a state-of-the-art metaheuristic technique, MIDACO, combined with a robust and efficient load-flow method called GBFS.The optimization problem was applied on the IEEE 69-bus and 118-bus systems for validation.Moreover, a massive parallelization framework using the fine-grained parallelization strategy of MIDACO was utilized to handle the stochastic many-objectives problem.The proposed method has shown huge advancements in uncertainty handling via reduced calculation burden for probabilistic-based methods.Results have demonstrated the advantage of DL allocation as an off-peak power management solution for large power mismatch when considering uncertainty.Moreover, the efficacy of the proposed method was compared with different evolutionary and swarm intelligence techniques.Likewise, the obtained results by the proposed method showed better accuracy and significant speed advantage in calculation times.Moreover, the results provided by the CBA have provided a competitive economic advantage for DLEMS over BEMS as an off-peak power management solution.Future work could be extended to address a different mix of renewable and non-renewable generation, multiple dump load allocations at multiple different buses, hybrid power management systems using BESS and DL to handle power variation at peak and off-peak hours, respectively, and the modeling of other uncertainty factors such as loss of generation and inverter failure.

Figure 3 .
Figure 3. Flow chart of the proposed optimization method to solve the problem of DL allocation with uncertainty.

Figure 3 .
Figure 3. Flow chart of the proposed optimization method to solve the problem of DL allocation with uncertainty.

Figure 4 .
Figure 4. Flow chart of the proposed load-flow method.

Figure 6 .
Figure 6.Pareto front of the many-objectives stochastic problem.The optimal solution as derived with the utopia-nadir balance approach is highlighted by the green hexagon shape for (a) 69-bus; (b) 118-bus.

Figure 6 .
Figure 6.Pareto front of the many-objectives stochastic problem.The optimal solution as derived with the utopia-nadir balance approach is highlighted by the green hexagon shape for (a) 69-bus; (b) 118-bus.

Table 1 .
Technical characteristics of DG units and MGs under study.

Table 2 .
DG unit arrangements for 69-bus system.

Table 3 .
DG unit arrangements for 118-bus system.

Table 4 .
Multi-objective expected results considering stochastic scenarios.
First step size only for |∆ f |.

Table 6 .
Microgrid data for CBA of water heating system.

Table 7 .
Yearly hot water demand using different EMS strategies.

Table A2 .
List of Symbols and Nomenclatures.Wind speed probability density function U i , N i The objective function's utopia and nadir values k s , c s Weibull distribution shape factor and scale index w Solution x weighted and average distance µ W , σ W Average and standard deviation for wind speed B j (x) Total microgrid cost at scenario s during off-peak hour h v co , v ci Cut off and cut-in wind speeds for the wind turbine FC s h , MC s h , EC s h , TC s j i (x), D j (x) h Fuel, maintenance, emissions, and technical costs at scenario s during off-peak hour h, respectively loss Total microgrid's active and reactive power losses Fi (x 1 ) Dump load expected objective function ∆V 1 , ∆ f The virtual bus's voltage and frequency deviations g h,s i (x 1 ) Dump load constraint handling function at scenario s during the off-peak hour h

Table A2 .
Cont.Frequency at iterations c and c + 1, respectively T in , T st Inlet and set-point temperature for water heating boiler V 1c , V 1c+1 Voltage at iterations c and c + 1, respectively C w , ρ w Specific heat and density for water, respectively.