Multi-criteria decision making monarch butterfly optimization for optimal distributed energy resources mix in distribution networks

The optimal integration of distributed energy resources (DERs) is a multiobjective and complex combinatorial optimization problem that conventional optimization methods cannot solve efficiently. This paper reviews the existing DER integration models, optimization and multi-criteria decision-making approaches. Further to that, a recently developed monarch butterfly optimization method is introduced to solve the problem of DER mix in distribution systems. A new multiobjective DER integration problem is formulated to find the optimal sites, sizes and mix (dispatchable and non-dispatchable) for DERs considering multiple key performance objectives. Besides, a hybrid method that combines the monarch butterfly optimization and the technique for order of preference by similarity to ideal solution (TOPSIS) is proposed to solve the formulated large-scale multi-criteria decision-making problem. Whilst the meta-heuristic optimization method generates non-dominated solutions (creating Pareto-front), the TOPSIS approach selects that with the most promising outcome from a large number of alternatives. The effectiveness of this approach is verified by solving single and multiobjective dispatchable DER integration problems over the benchmark 33-bus distribution system and the performance is compared with the existing optimization methods. proposed model of DER mix the optimization technique significantly improve the system performance in terms of average annual energy loss reduction by 78.36%, mean node voltage deviation improvement by 9.59% and average branches loadability limits enhancement by 50%, and minimized the power fluctuation induced by 48.39% renewable penetration. The proposed optimization techniques outperform the existing methods with promising exploration and exploitation abilities to solve engineering optimization problems.


A B S T R A C T
The optimal integration of distributed energy resources (DERs) is a multiobjective and complex combinatorial optimization problem that conventional optimization methods cannot solve efficiently. This paper reviews the existing DER integration models, optimization and multi-criteria decision-making approaches. Further to that, a recently developed monarch butterfly optimization method is introduced to solve the problem of DER mix in distribution systems. A new multiobjective DER integration problem is formulated to find the optimal sites, sizes and mix (dispatchable and non-dispatchable) for DERs considering multiple key performance objectives. Besides, a hybrid method that combines the monarch butterfly optimization and the technique for order of preference by similarity to ideal solution (TOPSIS) is proposed to solve the formulated large-scale multi-criteria decision-making problem. Whilst the meta-heuristic optimization method generates non-dominated solutions (creating Pareto-front), the TOPSIS approach selects that with the most promising outcome from a large number of alternatives. The effectiveness of this approach is verified by solving single and multiobjective dispatchable DER integration problems over the benchmark 33-bus distribution system and the performance is compared with the existing optimization methods. The proposed model of DER mix and the optimization

Introduction
The global energy demand is rapidly increasing with the growing population and its high reliance on modern technology and infrastructure. According to the Department of Economic and Social Affairs United States [1], the global population was 7.713 billion by the year 2019 expected to reach 10.9 billion by the year 2100. Therefore, contemporary energy resources and systems should be strengthened to meet the forthcoming energy demand. World bank's sustainable energy for all study, including 264 countries across the globe [2] in 2017, has reported that about 14.67% of the global population is facing energy crises. Almost 3 billion people across the globe do not have access to modern and clean energy resources for cooking, heating and agriculture. These sectors still depend on conventional fuels like biomass/wood, kerosene and diesel [3] which also increases the global carbon footprint. Nowadays, the primary aim of ascendancy agencies across the world is to meet the globally increasing energy demand while minimizing carbon footprints. According to the United Nations' report on climate change published in Dec 2019 [4], the amount of CO 2 gas emission from 2015 to 2019 was 20% higher as compared to the previous five years. Therefore, the use of clean energy resources such as solar, wind, nuclear and other energies have been recommended [5]. The demand for fossil fuels is also expected to increase day-by-day and there will be hardly any dip in their prices if not replaced by alternative energy resources. Despite the unstable oil market, the global CO 2 emission has reached its all-time high in 2019 [6]. The share of renewables is accelerating to an extent but it remains lower than 2% of world primary energy-producing source combining all together [7]. Statistically, the countries having access to abundant fossil fuel resources are not much serious about the use of renewable energy, therefore, emitting the highest CO 2 across the globe [8]. Nonetheless, the world bank is continuously supporting many global clean energy development programs to increase the share of renewables. For example, it has committed five billion USD for energy access programs during the years 2014 to 2018, out of that 1.4 billion USD was bound in 2018 [9]. The world bank is also providing financial aids to support rising economy nations of south Asia for various energy efficiency programs. A market analysis and forecasting report published by International Energy Agency (from 2019 to 2024) [10] has reported that the world's renewable power generation capacity is estimated to grow by 50% between years 2019 to 2024 where solar power will lead by the expected growth of almost 60%. It has also been estimated that by the year 2050, nearly 80% of the world energy demand would be supplied by renewables only, provided that sustainable policies and infrastructure are timely adopted across the globe.

Rationale for distributed energy resources mix
The above discussion on world energy scenario at glance deduces that clean and sustainable energy resources penetration will grow in modern power systems to supply the increased global energy demand in the future. These resources are naturally dispersed thereby cost-effective to utilize at the place of their availability, especially renewables in distribution networks. In decentralized and deregulated power system, such energy resources are known as distributed energy resources (DER). Furthermore, a small-scale electric power generation connected to the consumer side of the system is defined as distributed generation (DG), connected to the consumer side of the system [11].
It could be the main source of power which may include both dispatchable and non-dispatchable DERs. Here, 'dispatchable' refers to the generations where power production or dispatch is controllable e.g. diesel generators (DE), micro-turbines (MT), fuel cells, various energy storage technologies such as flywheel energy storage (FES), pumped hydro energy storage (PHES), battery energy storage system (BESS), supercapacitors (SCap), supermagnetic energy storage (SMES), etc. [12]. Electric vehicles can also be considered as dispatchable once their availabilities are confirmed for a duration. Nonetheless, the power output of non-dispatchable DGs cannot be controlled such as wind turbines (WT) and solar photovoltaics (PV). Despite clean, free and abundant availability, the solar and wind energy resources involve many issues for the grid due to their fluctuation, producing indeterminacy and availability [13]. There could also be considerable financial and technical challenges [14]. One of the issues associated with solar power generation is its high non-availability throughout the year. Although it may be preferable by utilities and stakeholders due to its static operation, modular integration and high irradiation availability during peak load hours of the day.
On the other hand, the quick growth of DERs has completely changed the consumer only image of traditional distribution networks. These systems were designed considering uni-directional power flow topologies and protection schemes thus not able to accommodate largescale integration of non-dispatchable/renewables that introduces multidirectional power flow [15]. The hosting capacities of these networks are very limited due to over-voltages, over-currents, substation backfeeding and protection schemes. In order to maximize system abilities to host the large-scale integration of renewables and to manage increasing load demand, these networks need up-gradation that requires high investment cost. The distribution systems have experienced many reforms in the last few decades to improve power delivery performance and economics.
In this scenario, the outcome of a non-optimal DER integration or mix could be counter-productive [16]. The optimal synergy of dispatchable and non-dispatchable low-carbon technologies can overcome some of the limitations for such energy resources while generating an enormous amount of benefits for the system operators, DER stakeholders and consumers. Table 1 presents the individual and generalized benefits, challenges and issues associated with some of the promising DERs. It shows that many features of these resources are complementary since each energy resource or technology brings a different set of benefits, challenges and issues. For example, energy storage can ease the issues of intermittent renewable power generation and enhances the operational flexibility of active distribution networks by effective management of renewable generation and load demand. However, most of the promising storage technologies are still amateur and expensive with a limited life-cycles. The control and management strategies of storage systems are also very complex [17]. Similarly, the diesel/gasbased DGs have high running and emission costs with low investment costs. A brief comparison of some promising DERs is presented in the table. At the large-scale, the efficiency and economics of traditional power plants can also be enhanced by hybridization with renewables, as suggested in the literature. Some of the related research include steam power plants through solar re-powering [18], application of nanofluids to energy systems [19] and the use of gas turbines and heat recovery steam generators instead of the boilers [20].
The optimal DER mix and synergy of different DERs can overcome the limitations of high renewables penetration that leads to alleviation of emission [23], least investment cost [24], substations capacity release [25], less over-loadings of feeders and transformers [11] Power produced from solar PV system at time (kW) Total power loss of the system (kW) , Active, and reactive power injections of node at time (kW, kVAr) , , , Real, and reactive power generations at node (kW, kVAr) Total power intake from the main grid at time (kW) Mean of the power intake from the main grid over a time period of (kW) while improving reliability [26], stability [27], power quality [28] and many more, of distribution networks. The dispatchable DGs can also be used as backup power support for critical and non-critical loads depending on service providers' priority and available power reserve [29].

State-of-the-art on optimal DER integration and optimization
To exploit maximum benefits, many single and multiobjective DER integration problems are addressed and solved in the existing literature. In the first place, a specific single objective is optimized while allocating DER in distribution systems. Two or more objectives have been considered in the multiobjective DER integration problem formulations. Generally, the sites and sizes of the DERs are determined at the time of integration planning, just after the long-term economic assessment of a cumulative capacity of renewables and before the actual commission and operation. It may be observed that by considering the number, sites, sizes and types of energy resources, the optimal DER integration turns into a complex mixed-integer, non-linear and non-convex reallife multiobjective optimization problem. Some of the objectives are conflicting in nature therefore suitable optimization methods (OMs) and multiobjective approaches (MOAs) should be used to get the most compromising solution for all stakeholders. The detailed findings and shortcomings of OMs and MOAs identified in the existing literature are presented in the following sections.

Optimization methods
In literature, many optimization techniques and strategies are suggested to find optimal sites and sizes of DERs. Generally, the optimization techniques used for optimal deployment of DERs can be divided into the following categories: conventional, numerical, statistical and heuristic methods. The conventional techniques include 2/3 thumb rule [30], optimal power flow (OPF) [31], exhaustive search (ES) [32], analytical method (AM) [33], loss sensitivity factor (LSF) based methods [34], improved analytical (IA) and exhausted load flow (ELF) techniques [35], convex relaxation methods [36], efficient analytical (EA)-OPF method [37]. Similarly, some of the numerical techniques are mixed-integer program (MIP) [23], mixed-integer non-linear programming (MINP) [38], chance constrained-programming (CCP) [39], quadratic-programming, dynamic-programming, sequential quadraticprogramming (SQP) [40] and affine arithmetic (AA) [41]. Taguchi based methods can be considered as one of the statistical optimization techniques adopted to solve the DER integration problem of distribution systems [42]. These are found very effective for interaction analysis of variables [43]. It is found that analytical methods are computationally efficient but sometimes fail to find the global optimal solution [44], especially for multiobjective optimization problems. The possible reason could be the simplified assumptions made to ease the complexity of optimization problems. Similarly, the numerical optimizations are computationally fast and efficient but need very accurate problem formulation and modeling of the systems and problems [42]. Furthermore, most of these methods are not able to solve multiobjective optimization problems efficiently, the multiple objectives are generally converted into a single objective.
The artificial intelligence (AI) based meta-heuristic methods are very efficient to solve complex real-life engineering optimization problems. These techniques also have black box problem-solving abilities which are independent of the problem formulation and modeling. Some of the prominent nature-inspired and AI-techniques used to solve the DER integration problems include particle swarm optimization (PSO) [13], genetic algorithm (GA) [45], teaching learning-based optimization (TLBO) [46], ant colony optimization (ACO) [47], simulated annealing (SA) [48], invasive weed optimization (IWO) [49], -dominance-based evolutionary algorithm ( -DEA) [50], intelligent water drop (IWD) algorithm [51], backtracking search optimization algorithm (BSOA) [52], bacterial foraging optimization (BFO) [53], moth search optimization (MSO) [54], harmony search algorithm (HSA) [55], bat algorithm (BA) [56], artificial bee colony (ABC) optimization [57], ant lion optimization algorithm (ALOA) [58], firework algorithm (FWA) [59], etc. It has been analyzed that the standard variants of some of the popular metaheuristic methods show limitations when applied to solve complex real-life engineering optimization problems. Like GA and MSO show the tendency of premature conversion or local trapping. Similarly, SA and PSO methods are performing relatively poor to seek global optimal solutions [22,60]. Furthermore, the poor communication of TLBO in the second phase may lead to inadequate information sharing among the learners and causes local trapping [61]. To overcome the limitations of poorly performing standard variants of metaheuristic methods, various improvements and hybridization have also been suggested in the literature.
Recently, a computationally fast and very effective metaheuristic method named as 'monarch butterfly optimization (MBO)' was proposed P. Singh et al. by Wang et al. [72]. The method is inspired from the migration behavior of monarch butterflies and shows comparatively fast convergence to seek global optima for real-life optimization problems. To the authors' best knowledge, the promising abilities of the standard method have not been yet used to solve the complex multiobjective DER integration problems of distribution systems.

Multi-criteria decision making (MCDM) approaches
The real-life multiobjective DER planning is a large-scale optimization problem with many alternatives and objectives (criteria). It requires high memory and more computational time to implement such a big matrix. The performance of most of the multiobjective evolutionary algorithms deteriorates by the large size of the Pareto set, where the decision-makers have too many alternatives (Pareto front) to select from [73]. Most of the literature has formulated and solved the DER integration problem either as a single objective or multiple objectives transformed into to a single-objective optimization problem. A handful of literature only dealt with it as an actual multiobjective optimization problem. The popularly known MCDM or multiobjective approaches used to solve the optimal DER integration problems can include penalty function (PF) [22], -constrained ( -C) [26], Pareto front concepts (PFC) [45], weighted-sum (WS) [46,74], grey relational projection (GRP) [58], novel global harmony search II (NGHS-II) [75], goal programming (GP) [76], non-dominated sorting algorithm II (NSGA-II) [77] and fuzzy approaches (FA) [78].
It is identified that some of these MOAs also show limitations when applied to solve the real-life and large-scale multiobjective optimization problems. For example, two major drawbacks of WS approaches may include the non-uniform distribution of solutions and its inability to get the solution that lies in the non-convex region of the Pareto front [79]. The PF-based approaches are very popular to deal with constrained search-space by penalizing the infeasible solutions however need fine-tuning of the penalty factors. Therefore, the selection of optimal solutions is highly influenced by the prime objective(s) of higher value or a high degree of penalization [80]. Furthermore, a target needs to specify in GP-based approaches, where a -constrained approach requires a careful choice of master and slave objectives, which directly affect the final solution in both the techniques. The GRP approach combines the promising features of grey relational analysis and projection method to solve MCDM problems. One of the drawbacks of traditional projection methods includes the equal projections of multiple vectors on the ideal solution vector thus the alternatives cannot be compared [81]. These also lacking the characteristics of an affinity. Most of these approaches are only suitable for isochronous equidistant sequences therefore greatly limit their application [82]. Pareto-based optimization methods attempt to find non-dominated solutions, which cannot improve one objective without degrading some others and reflect different trade-offs between multiple objectives. The truncation strategy adopted in NSGA-II may destroy the solutions' diversity; lead to uneven distribution, as identified in [83], therefore some modified variants of these approaches are suggested in the literature.
The technique for order of preference by similarity to ideal solution (TOPSIS) is a robust and quite straightforward MCDM approach which is found to be very effective to select promising alternatives for various real-life multiobjective problems. The technique was developed by Hwang and Yoon in 1981 [84] however, it has not been much explored to solve MCDM problems of distribution systems.
After a thorough literature review in previous sections, the existing researches on DER planning can be categorized based on various modeling aspects such as DER and load types, load-levels, objective functions, problem formulations, optimization methods, etc. A comparative overview of the existing research works on DER integration in distribution systems is presented in Table 2 (see Refs. from the table). It briefly summarizes the contributions/findings of existing DER integration models and compares the problem formulations of DER planning based on objective functions, MOAs, OMs, load levels and DER types.
The comparison shows that most of DER integration models are solved either as a single objective optimization problem or by converting many goals into a single-objective optimization problem. Almost all MCDM problem formulations on mixed and multiple DER integration (dispatchable and non-dispatchable DER) are also formulated either by considering one goal or to maximize the overall monetary benefits of synergy by adding multiple financial objectives, therefore not completely modeled as a realistic MCDM problem with conflicting objectives. A handful of literature may be found that deals with the sophisticated real-life MCDM problems of the DER mix integration in the light of conflicting technical objectives and suitable MOAs.

Contributions
In this paper, a new MCDM problem is formulated and solved for the optimal integration of mixed DERs in distribution systems aiming to determine the optimal types, sites and sizes of multiple dispatchable and non-dispatchable energy resources.
The key contributions are listed here.
-The recently developed MBO method is introduced to solve the optimal DER integration problem of distribution systems. -A new deterministic MCDM problem is formulated to obtain the optimal DER mix by considering five vital objectives: the minimization of annual energy loss, node voltage deviation, grid substation back-feed and power fluctuations induced by renewables while maximizing the voltage stability index (VSI) of distribution networks. -A novel hybrid optimization method is proposed to solve the large-scale MCDM problem, i.e. MBO-TOPSIS. Whilst the metaheuristic technique generates non-dominated solutions (creating Pareto-front), the TOPSIS approach selects that with the most promising outcome from the enormous number of alternatives. -The effectiveness of these methods is verified by solving both single and multiobjective dispatchable DER integration problems over the benchmark 33-bus distribution system and comparing the performance with other existing optimization methods.
The proposed multiobjective DER mix problem is also applied to multiple scenarios. The results reveal that the optimal DER mix, under the proposed framework, produces more synergy benefits for network operators and consumers, in terms of annual energy loss and minimum power fluctuation caused by renewables, voltage profile and stability improvement, with minimized risks caused by back-feeding to substations.

Problem formulation
In this section, a deterministic MCDM problem is formulated for optimal mixed-integration of dispatchable and non-dispatchable DERs in distribution systems. Two fast-growing renewables have been considered for enhanced and mixed exploitation of non-dispatchable DER technologies in this work, i.e. WTs and PVs. The SCs are also considered to provide the reactive power support in some cases. The objective functions, constraints, TOPSIS based multiobjective problem formulation followed by individual deterministic modeling of each renewable power generation, shunt capacitor and load demand are presented in the following sections.

Objective functions
In modern distribution systems, the utilities aim to maximize their profits while meeting various techno-economic and social constraints simultaneously. The rapid growth of renewables along with artificial intelligence, advanced information and communication technologies offer many profit-making opportunities to DNOs, modern consumers and DER owners. In this situation, there are many key objectives that have to be optimized simultaneously while considering the need for each stakeholder in deregulated power system. To make the proposed model more realistic, some of the most practical, strategic, useful and operative objectives have been considered here. The detailed formulation of each objective function is presented below.

Minimization of annual energy loss
In the conventional power delivery system, the maximum power loss occurs in distribution networks directly affecting the annual revenue of utilities. A DNO could save millions every year by annual energy loss reduction [95]. It can be achieved by effective utilization and essential up-gradation of distribution system resources such as onload tap changers, voltage level standardization and feeder voltage regulators. Similarly, the DERs can also reduce the power loss when placed optimally however the main challenge of power loss minimization is conflicting with the goal of renewable penetration maximization. By considering the above-mentioned facts, the minimization of annual energy loss has been adopted as one of the objectives in the proposed formulation for optimal integration of renewable energy resources in distribution networks. The objective function of annual energy loss minimization is expressed as where, Eq.
(2) represents the total real power loss [16] of the system at time .
To determine the annual energy loss, the daily energy loss over the time period is multiplied with daily to annual conversion factor in Eq. (1). Here, and denotes the number of nodes in the network and resistance offered by the line connecting nodes and ; whereas, , , , and are representing the real and reactive power injections, voltage magnitude and angle of node at time respectively. P. Singh et al.

Table 2
A comparative overview of existing optimization methods, multiobjective approaches, energy resources and objectives used for optimal DER integration in distribution systems.

Minimization of node voltage deviation
The load demand of the network is changing all the time, so does the nodes voltage, which directly affects the equipment connected to the network. Some of the power consumers are very much attentive about the quality of power supplied to them because their machinery/equipment could be voltage-sensitive. Based on the regional requirements, regulators specify the supply voltage limits that have to be maintained by DNOs. The rapid growth of intermittent energy resources could affect the node voltages if not placed and sized adequately. Many studies have considered the nodes' voltage limits as the hard constraints, however, that directly affect the performance of optimization methods adopted. To mitigate some of such issues of renewables and optimization modeling, the nodes voltage deviation should be considered as one of the objectives when integrating DERs in distribution networks. The objective function of the node voltage deviation, measured from the nominal value, is expressed as In this equation, a maximum value of the square of each node voltage deviation ( 1 − ) 2 from the nominal value, i.e. unity, over a time period of is selected and then the sum of these voltage deviations for nodes is minimized. Thereby, the overall expression minimizes the node voltage deviation from nominal value whether it is overvoltage or under-voltage since square turn it into absolute value, as suggested in [83]. In addition to this objective, the minimum and maximum permissible node voltages are also ensured by the voltage limit constraint expressed in Eq. (16).

Maximization of voltage stability index
It has been investigated that the quality of the node voltage profile is not a sufficient criterion to estimate the system's health over variable loading. To overcome the shortcomings of voltage profile, a voltage stability index was presented in [96]. The index helps to determine the voltage collapse points in radial distribution networks. It suggests that a branch with minimum VSI is more sensitive to voltage instability. To manage the node voltages within the threshold limits over a variable loading, VSI is also considered as one of the objective functions of DER integration. The objective function of VSI maximization is expressed as The mathematical expression for the VSI of a branch connecting nodes ' ' and ' ' is presented in Eq. (7).
In [96], the VSI was defined for the single system state that is now extended to states in (7). It represents the VSI of a line between nodes and at time , where denotes the line reactance. Here, the aim is to maximize the VSI value of a line that has the least VSI among all branches of the system over states, as expressed in Eq. (6).

Minimization of grid substation back-feeding
Governments across the globe are trying to increase the share of renewables in power generation that has risks to generate surplus power during light load hours. The excess power generation back-feeds the grid substation and traditional distribution systems were not designed to manage the reverse power flows. The protection schemes also need to be improved shortly to increase renewable penetrations and, to avoid system stability and security risks. In the present condition of systems, such back-feeds are rather restricted by system operators, therefore, it is also considered as one of the objective functions. In this paper, the back-feed power to the main grid, if any, is minimized, as expressed in Eq. (8). where, In Eq. (8), − represents the amount of back-feed power to main grid when renewable power is in surplus at time . As expressed in Eq. (9), − is considered zero when power intake of distribution system is positive at any time , i.e. ∑ =1 ( , − , ) + ≥ 0. The power generation and load demand of node at time are represented by , and , respectively.

Minimization of the fluctuations in power intake from the main grid
The load demand of distribution systems is highly fluctuating with the time therefore an accurate day-ahead demand scheduling is required. Moreover, the increasing penetration of renewables changes the magnitude of fluctuation in the power intake at grid substations that directly affects the grid frequency. Traditionally, warnings are issued to DNOs that overdraw power during a period when grid frequency falls below the nominal value and vice versa [97]. Sometimes, unscheduled interchange charges are also applied to energy purchases during the constraint violation period. It would be difficult to observe and mitigate such dynamic issues of frequency regulation in the planning stage, therefore the net load demand deviation of the distribution system has been considered as one of the performance measures in this paper. This measure helps to determine the optimal-mixed allocation of renewable energy resources that would restrain some amount of power fluctuation transferring to the upstream transmission networks. Therefore, a standard deviation of the power intake from the main grid, over states, is minimized in Eq. (11) when integrating the DERs in distribution networks.
The equation tries to maintain the power intake of the distribution system close to mean value thereby reduce the maximum variations in the scheduled power at grid substation. Here, and are representing the total power intake of distribution system from the main grid at time and its mean value over the time period of .
Eqs. (12), (13), (14), (15), (16) and (17) are expressing the constraints of active and reactive nodal power balances, DER hosting capacity of the network, DG active and reactive power hosting capacity limits of a node, node voltage limits, and thermal limits of the line connecting nodes and respectively. In these equations, , , , and max represent the Y-bus element, impedance angle, current at time , and maximum ampere limits of the line between nodes and . Similarly, the parameters , , ,̂, ,̂, max , min , , and denote peak real and reactive power loads of the system, real power generation at a node and its maximum limit, reactive power support at a node and its maximum limit, upper and lower permissible voltage limits of the nodes, binary decision variables of DG and SC deployment at node respectively. P. Singh et al.

MCDM problem formulation by using TOPSIS approach
The function of multiobjective or multi-criteria decision making problems can be declared as where, 2 is the number of objectives/criteria which must be greater than one. is representing a feasible set of decision vectors in the search domain, defined by constraint functions. An objective function with a decision vector is generally defined as ∶ → 2 ; = ( ), = 1, 2, … , 2 . There is no typical feasible solution for multiobjective optimization problems that minimize or maximize all objectives at the same time thus the concept of Pareto front or optimality is used. It states that the solution of any objective cannot be improved without deteriorating at least one criterion/objective from others. A feasible solution 1 ∈ is called to a Pareto or dominate to other solution 2 ∈ , only if: The solution 1 ∈ is said to be Pareto optimal if there exists no solution that dominates it. In the previous section, a multiobjective optimization problem is formulated for optimal allocation of different DERs. In order to solve such a complex multiobjective optimization problem, an MCDM approach known as TOPSIS has been adopted in this paper. The approach was developed by Hwang and Yoon in the early 80s [84] followed by some modifications by Yoon in 1989 [98] and Hwang and Liu in 1993 [99]. The basic idea of this technique is to select an alternative from the Pareto-front as a compromising solution that has the least geometric distance from positive ideal solution (PIS) and greatest geometric distance from negative ideal solution (NIS). The essential steps of the TOPSIS approach for solving an MCDM problem are given below.
Step-1: Create a decision matrix A 1 × 2 for 1 alternatives/solutions and 2 criteria or objectives, as expressed below.
Step-2: Normalize the decision matrix A developed in Eq. (19). The normalized decision matrix R is presented here. where, Step-3: In this step, the positive and negative ideal solutions are determined. PIS is a set of best solutions among 1 alternatives for each individual criterion . Similarly, NIS is a set of worst solutions among 1 alternatives for each individual criterion , as expressed below.
where, Step-4: Calculate the Euclidean distances of an alternative from PIS and NIS. The distance from PIS, i.e. + is determined as Similarly, the distance of th alternative from NIS, i.e. − is also calculated as where, ≥ 2. = 2 is the most commonly used value to determine Euclidean distances.
Step-6: Rank the alternatives according to their closeness indexes. An alternative with highest RCI value will be chosen as the most compromising solution.

Deterministic modeling of renewable power production, shunt capacitor and load demand
It is a challenging task to estimate wind and solar power generations accurately due to the intermittent behaviors of wind speed and solar irradiation. Similarly, changes in time make the load demand very uncertain. In practice, various forecasting techniques are being used to match the supply and load demand in daily system operations. The day-ahead hourly scheduling is very common in modern power system operations. However, the DER planning is done on a long-term basis therefore annual load and generation profiles are generally considered when optimally integrating the emerging technologies in the networks. Hourly deterministic power production from WTs, PVs and system load demand is presented in the following sections.

Wind power production
As discussed earlier, wind power production is intermittent due to variable wind speed. WT is an effective technique to capture power from wind speed and transform into electric power. A fluctuating wind speed causes WT to produce variable electric power. This power generation depends on various parameters of a WT such as wind sweeping area, rotor diameter, pitch/blade angle, air density, wind speed, generator and gearbox efficiency etc., as shown in Fig. 1. The power production from a WT at time is expressed in Eq. (29), as devised in [100].
where, wind power production at any time is directly proportional to the third power of wind speed , i.e. ∝ ( ) 3 . In Eq. (29), , , and denote the air density, power coefficient and area swept by the rotor blades of WT respectively. For a given time and WT, remaining parameters and constants of Eq. (29) can be assumed constant. Without loss of generality and simplicity of the proposed model based on this assumption, the power production from a WT at time can be devised P. Singh et al.

as
here, , , _ , and _ denote the rated power capacity, rated, cut-in and cut-out wind speed of WT respectively.

Solar power production
In this section, the power production from solar photovoltaics is also modeled as done for WTs. Solar PV is one of the prominent renewable power generation technologies. Unlike wind energy, solar irradiation is only available in day time and also intermittent in nature. Ambient temperature and solar radiation are the dominating control variables for sizing the PV systems where solar panels are deployed in series and parallel combinations. The electric power production from a solar module at time can be expressed as where, power production of a solar panel is directly proportional to its efficiency , incident irradiation flux (W/m 2 ), and area of collector (m 2 ) at any time . If area of incident radiation flux and panel efficiency are assumed fixed at any time , the power generation from a panel/module can be determined by Eq. (32).
here, , and are the rated power generation and solar irradiation of module, and illumination at time respectively. A solar power system is comprised of many PV panels connected in series and parallel combination as shown in Fig. 2. Here, it is assumed that solar irradiation is uniformly available over all panels then the total power production from a solar power plant or × array at time can be calculated as P. Singh et al.

Load demand of the system
Like renewable power production, the load demand of the system also varies with time and depends on the type of consumers e.g., commercial, residential, industrial and mixed. Therefore, the hourly deterministic load demand of the system is also modeled here. An hourly load multiplying factor, is devised from the historical load data of a system. The real and reactive power loads ( , ) of a node at time are mathematically expressed as . . = 1.6 × 0 , and = 1.6 × 0 In this equation, 0 and 0 are representing the nominal real and reactive power loads of node . The peak load of a node is usually considered 1.6 times of the nominal loading, as suggested in [55].

Shunt capacitors
Apart from real power production or purchase, DNOs need to maintain the required amount of reactive power balance between supply and demand. Traditionally, fixed and variable shunt capacitors are employed to provide reactive power support in distribution systems. The DGs can also provide reactive power support but at a higher cost with respect to shunt capacitors. Moreover, the frequent adjustment of generation excitation control in order to supply or consume reactive power may damage the generator's field winding [75]. In this work, the reactive power support from a fixed SC at node is modeled as where, and denote the optimal number of capacitor banks and reactive power capacity of each bank to be deployed at node .

Monarch butterfly optimization for optimal mixed-integration of DERs
The optimization problem developed in the previous section is a large-scale, complex real-life, mixed-integer and non-linear optimization problem which is difficult to solve effectively by using conventional optimization methods. In this paper, a recently developed metaheuristic approach is adopted, i.e. monarch butterfly optimization. This technique is inspired by the migration behavior of the monarch butterfly of North America. It is developed by G.G. Wang et al. in 2015 [72]. This eastern North American species is notable for its annual southward migration from the north and central US and Canada to Florida and Mexico during the autumn or late summer. These butterflies change the location based on the weather conditions in these regions. Fig. 3 presents a migration behavior of monarch butterflies, when weather is not suitable in region-A they move to more favorable region-B. These flutters migrate thousands of miles with a multi-generational return during the spring. This migration behavior is modeled into some set of mathematical equations to introduce this new swarm-intelligence based optimization method, discussed here.

Monarch butterfly optimization
In this meta-heuristic technique, the butterfly flutter assumed to migrate from region-A to region-B in the month of April and further from region-B to region-A in September [72]. During this process of migration, they keep producing offspring which replace their parents. This complete process is mathematically bifurcated into two updating operators known as migrating operator and butterfly adjustment operator. During this migration, monarch butterflies follow the set of rules listed below.

Rule-1:
Whole flutter of monarch butterflies should stay in these two areas only so that the total number of butterflies in these two regions can form the complete population.

Rule-2:
The butterflies give birth to offspring only in one region by using the migration operator. Rule-3: The newly generated child replaces its genitor if it exhibits better fitness than its parent. This rule will keep the butterflies population constant.

Rule-4:
A butterfly individual with the best fitness value will automatically migrate to the next generation and will not be changed by any upgrading operators.

Butterfly migration operator (BMO)
This is the butterfly migration operator that helps to model the migration behavior of monarch butterflies. Suppose the population count of monarch flutter staying in region-A named as subpopulation-A is , which can be obtained as = ceil( × ). Similarly, the butterfly population for region-B called as subpopulation-B is = − . Here, and are representing the total number of monarch butterflies in both the regions and a ratio of butterflies in region-A respectively. In order to perform the migration operation, a decision variable, is generated as here, is a flutter migration period generally set to 1.2 by considering 12 months of a year. Whereas, is a random value generated from a uniform distribution function.
If ≤ then the new location of th butterfly is updated by using following migration operator, expressed as where , ( + 1) denotes the th element of that represents the location of th butterfly in generation + 1. Similarly, 1 , ( ) is the th element of 1 for the current generation . Furthermore, 1 is a randomly selected individual from subpopulation-A, i.e. . Suppose > then the th element of new location for th butterfly is calculated as here, 2 , represents the th element of location 2 for the current generation . Whereas, 2 is a randomly selected individual butterfly from subpopulation-B, i.e.
. It may be observed that the direction of butterfly migration can be controlled by adjusting the value of migration ratio . It will help to balance the migration in MBO method. For example, the high value of can force migration operator to choose more elements of butterflies from region-A and vice versa. Therefore, adequate value of this ratio is very important to provide a balance in migration process. Generally, the value of this ratio is considered to be = 5∕12 = 0.4166 as per the migration period. The pseudo-code of monarch butterfly migration process is presented in Algorithm 1.

Butterfly adjustment operator (BAO)
The butterfly adjustment operator (BAO) is used to update the positions of butterflies in region-B. Similar to BMO, a random decision variable is generated in BAO as well, e.g.
. If ≤ then the position of th butterfly is updated as where, , ( + 1) denotes the th element of th butterfly position during the generation + 1. Furthermore, the , ( ) represents the th component of the fittest butterfly position , among both the regions A and B in current generation .
On the other hand, if > then the position of th butterfly will be updated as here, 3 , ( ) is representing the th element of 3 th butterfly position 3 ; ∀ 3 ∈ , randomly selected from region-B. In this condition, one more sub-condition is applied in which the random decision variable ' ' is compared with the algorithm parameter known as butterfly adjustment rate, i.e. . If > then following position updating rule is applied to further adjust the location of th butterfly.
where, is a weighting factor defined as = max ∕ 2 (42) here, max is the maximum walk step which the butterfly individual moves in one step and is the walk step size of th butterfly, calculated by using Levy flight operation.
In Eq. (41), the higher value of increments the significance of long step length thus will improve the exploration of MBO algorithm. On the other hand, the smaller value of reduces the influence of step length therefore improves the exploitation capability of the algorithm. The pseudo-code of butterfly adjustment operator is presented in Algorithm 2.

Optimal mixed-integration of DERs by using hybrid MBO-TOPSIS approach
In this section, the MBO algorithm is extended to solve the optimal DER integration problem of distribution systems. The decision variables of this optimization include sites and sizes of different DERs. The structure of an individual butterfly, i.e. , ; ∀ , ∈ containing the information of nodes and capacities of different DERs is presented in Fig. 4. Where, = size of optimization parameters or problem dimension. Therefore, the total number of decision variables in the proposed DER integration problem or dimension would be calculated as = 2 × ( + + ); where, , and are respectively representing the total number of wind turbines, solar power stations and shunt capacitors assumed to be installed in any distribution system. The total installation or rated capacity of WT, i.e.
; PV plant, i.e. = × × and SC bank, i.e. at node are optimally determined. For more DER mix, the structure of an individual shown in Fig. 4 can vary in the similar way.
The flowchart of proposed hybrid MBO method to solve DER integration problem is presented in Fig. 5.

Simulation results and discussion
In this section, the optimization problem developed in Section 2 for optimal DER mix in distribution networks is solved by using the proposed hybrid MBO-TOPSIS technique, i.e. Section 3.4. This real-life MCDM problem is investigated on a benchmark test distribution system of 33 buses [101]. It is a 12.66 kV radial network with nominal real and reactive power demands of 3.715 MW and 2.300 MVAr respectively. As stated, the purpose of this work is to introduce the promising MBO technique to solve the DER integration problem followed by the proposed multi-criteria decision-making problem for mixed integration of renewable-based DGs in distribution systems. The potential of MBO and proposed hybrid MBO-TOPSIS techniques is demonstrated before solving the developed multiobjective renewable DG integration problem. The performance of these techniques is compared with some of the existing optimization methods already available in the literature. These problems are formulated in the MATLAB environment on Intel(R) Core(TM) i5-6200 CPU@2.30 GHz processor with 8 GB RAM. A backward-forward load flow method is adopted to solve the power flow equations of active distribution networks. The validations of these methods are presented in the following sections, followed by the simulation results of the proposed multiobjective optimization problem for the DER mix.

Validation of standard MBO for single objective dispatchable DG integration problems
To establish the standard variant of MBO, a single objective dispatchable DG integration problem is formulated and solved for power loss minimization in 33-bus test distribution system. The optimal sites and sizes of MTs are determined with an aim to minimize the real power loss with nominal loading, expressed in Eq. (2). The optimal number of DGs for the 33-bus distribution system is three, as validated in [55]. Therefore the sites and sizes of three MTs, operated at unity power factor (UPF), are determined here. A comparison of simulation results obtained by MBO and some of the well-known optimization methods are presented in Table 3. The comparison shows that the MBO technique has great potential to find the optimal solutions for single objective DG integration problems of distribution systems. It provides maximum power loss reduction as compared to many of the existing optimization techniques, i.e. 64.78% reduction. This meta-heuristic technique outperforms analytical, numerical and heuristic methods such as LSF based approaches, FWA, IWO, analytical and improved analytical (IA) methods, ELF, TLBO and QOTLBO, hybrid HSA-PABC algorithm, efficient analytical method (EAM), SQP and HGWO. In addition, the node voltage profiles and branch voltage stability indices of the system, before and after DG integration, are presented in Figs. 6(a) and 6(b) respectively. The figures show that the proposed optimization model has also improved the voltage profile and stability index of the system despite solving a single objective DG integration (SODGI) problem. It has been observed that the minimization of power loss also improves the voltage profile and stability indices however up to a limited DG penetration.
Like other meta-heuristic methods, MBO also provides a different solution on each run. To further investigate the inherent characteristic and potential of MBO, some of its important performance parameters such as best, worst, mean and standard deviation of optimal solutions obtained in 50 independent trials are presented in Table 4. The population size and maximum generations considered in MBO for these trials are 50 and 100 respectively. An average CPU time of these independent runs is also presented in this table. It shows that the MBO method is efficiently seeking the optimal solution in each independent run since the mean fitness of these 50 autonomous runs is very close to its best fitness value. It has also been observed that the mean fitness of MBO itself is superior to the best fitness values of the existing methods presented in Table 4. The best and mean convergence characteristics of MBO, for power loss minimization, is also shown in Fig. 7. The convergence characteristic of the best fitness shows that the MBO technique is continuously seeking for the optimal solution and finds the global optima within 20 iterations. Similarly, the mean fitness characteristic in this figure demonstrates the extensive exploration capabilities of the method in the search space before exploiting the promising region, e.g. highly fluctuating mean fitness in the first 20 iterations due to some random variables of the method. The small spikes after 20th iterations show the local searching of this method in the promising area. This investigation on standard MBO and existing optimization P. Singh et al.  methods reveals that most of the loss sensitivity-based and standard meta-heuristic methods are unable to explore the global optima of DG integration problems. However, hybrid and improved optimization methods have shown some potential to explore the global region. Table 3 A comparison of optimal sites and sizes determined by MBO and some of the existing optimization methods for power loss minimization in 33-bus distribution system.    A comparison of optimal sites, sizes and objective functions determined by hybrid MBO-TOPSIS and some of the existing multiobjective optimization methods for 33-bus distribution system.

Validation of the proposed MBO-TOPSIS for multiobjective dispatchable DG integration problems
After validation of the MBO method on a SODGI problem, the proposed hybrid MBO-TOPSIS is also validated for multiobjective DG integration (MODGI) problems of distribution systems. To compare with the literature, an additional DG allocation problem is formulated in this section by considering three conflicting objectives such as power loss and voltage deviation minimization while maximizing VSI of benchmark 33-bus distribution system, expressed in Eqs. (2), (5), and (6) respectively. The hybrid MBO-TOPSIS developed in Section 3.4 is used to solve this MODGI problem. The optimal sites and sizes of three dispatchable MTs are determined, as highly suggested in the existing literature. An optimal solution obtained by the proposed MBO-TOPSIS approach is compared with some of the well-known multiobjective optimization methods and presented in Table 5. It shows that the optimal solution of proposed hybrid MBO-TOPSIS is promising as it provides minimum node voltage deviation and maximum VSI for 33bus distribution systems as compared to existing methods like GA, PSO, improved NSGA-II and hybrid GA/PSO, TLBO and QOTLBO, CTLBO, QOCSOS and improved MOHSA. The value of system power loss is more than QOCSOS and IMOHSA however comparable with CTLBO. It seems that the QOCSOS and IMOHSA methods have had more focus on power loss reduction instead of voltage profile and stability improvement since voltage deviation and VSI values of these methods are very poor as compared to the proposed MBO-TOPSIS and other optimization techniques presented in the table. Furthermore, four DGs are deployed with IMOHSA [75] instead of three DGs in other case studies.
It is noticed that a higher power loss reduction has been obtained for the SODGI model in the previous section, this is because it aims only to minimize the power loss. To demonstrate the benefits of MODGI over SODGI models, the node voltage profiles and stability indices, obtained for the base case, SODGI and MODGI solutions are compared and presented in Figs. 8(a) and 8(b) respectively. The figures show that the solution obtained by MODGI provides a substantial improvement in the node voltage profiles and branch voltage stability indices as compared to the base case and SODGI model. The MODGI model and optimization method generates a more balanced and practical solution that reduces the power loss while simultaneously improving the nodes voltage profiles and branches voltage stability of the distribution system.
This comparison validates that the proposed hybrid MBO-TLBO method is reasonably efficient to solve MODGI problems of distribution systems. Now it can be adopted to solve the proposed multiobjective DER mix integration problem.

Simulation results of proposed multiobjective DER mix integration problem
The multiobjective optimization problem developed in Section 2, for multiple and mix DER integration in distribution system, is solved by the proposed hybrid MBO-TOPSIS method presented in Section 3.4. As Table 6 The values of parameters used in this case study. discussed, the proposed DER integration model is a complex combinatorial optimization problem therefore the most likely states of the generation and load demand are often used to minimize the computational burden [12,15]. Further, the consideration of most likely states or mean data would increase the probability of generating the most balanced solution and minimize the risk of non-optimal solutions, over highly variable states throughout the year. The proposed DER integration model has been designed in a deterministic environment to clearly show the technical impact of the DER mix, around the mean system states (low risk) without incorporating any type of system uncertainty. A concept of load multiplying factor, i.e. , is used to generate the hourly deterministic load demand of the system. It is obtained by dividing the hourly annual mean load with the peak load value and shown in Fig. 9(a). In this study, the considered load data is a mix of residential, commercial and industrial loads thereby the peak load demand appears between 11:00 to 16:00 h. Eq. (34) and load multiplying factor ( ) are used to estimate the hourly load of each node. For each solution that suggests the rated capacities of renewable DERs, the hourly mean of wind and solar power production is determined by using Eqs. (30) and (33), based on the historical hourly annual mean of wind speed and solar irradiation respectively. To demonstrate generation multiplying factors, the hourly power productions over the time are divided by the corresponding rated capacity of that DER. One sample of a solution is shown in Fig. 9(b) and (c). Fig. 9(b) shows that all WTs produce power at their corresponding rated capacities during 00:00 to 7:00 and 22:00 to 00:00 h since the average wind speed will remain between rated and cut-off speeds of the turbine. Due to low solar irradiation in the used data, the maximum average power generation of PV systems remains below the rated capacity. A maximum average power production of 80% is seen from Fig. 9(c). The hourly data of wind speed, solar irradiation and load demand are adopted from [12,15]. It is also assumed that all nodes have uniform availability of solar irradiation and wind speed since the distribution system is scattered in a small geographical area. The values of parameters used in the proposed model are summarized in Table 6.
To investigate the effect of different mixed renewable based DGs and SCs, following cases are framed and solved by using the proposed hybrid MBO-TOPSIS approach.   The simulation results of these cases are presented in Table 7. It includes the optimal values of decision variables such as nodes/sites, types and sizes of DERs and objective functions like DER penetration (DRP), annual energy losses (AEL), standard deviations of the power intake from the main grid (STDG), minimum of mean node voltages of the system (MMV) and a minimum of mean VSI values (MVSI) over time duration and percentage AEL reduction (AELR) of all cases. The best result of each objective among these cases are highlighted. From case I, it is observed that just with the integration of wind power generation operating with UPF, the annual energy loss has been significantly reduced compared to the base case, i.e. Case 0. The wind power generation has also improved the node voltage profile and branch stability indices. Despite these benefits, the standard deviation of the power intake from the main grid, over states, has increased. This case shows that wind power generation is highly volatile introducing fluctuations in power purchase from the main grid at grid-substation. This is one of the vital issues faced by utilities to balance the fluctuating power and demand mismatches.
To investigate the effect of mixed energy resources, Case I is extended by considering a PV power generation in Case II. All the DGs are assumed to be operated at UPF in this case. It can be analyzed that the mixed integration of two renewable energy generation technologies has further enhanced the annual energy loss reduction to 52.59%, higher than that of Case-I, i.e. 47.56%. Although PV generation is available in daytime only, its inclusion with WTs has improved the node voltage profile, VSI and STDG significantly. It basically supplies most of its power in peak load hours when wind power generation is low. Therefore, the mixed integration of solar PV and wind power generations compensates for the shortcomings of each other. It can be seen from the STDG value which is lower as compared to the wind only integration case, i.e. Case I. Nevertheless, the amount of compensation and mixed DER capacities are highly dependent on the wind speed and solar irradiation availability hours.
In cases I and II, the reactive power support from the DGs have not been considered. To investigate the effect of reactive power support from DGs, 0.85 lagging power factor (LPF) based WTs are considered in Case III when determining sites and sizes of WTs and one PV system. It may be observed that the optimal nodes are different than that of Case II. The total capacity of wind power generation is reduced P. Singh et al. Table 7 The optimal sites, sizes, and values of objective functions obtained for different DER mix integration in 33-bus distribution system.  but the size of solar PV increased since high power generation would be required in the daytime. The reactive power support from WTs has significantly improved the system performance in terms of annual energy loss, voltage profile, VSI and STDG. Especially, the AEL is significantly reduced as compared to Case II.
According to grid-codes in most of the countries, wind power producers are responsible to compensate the reactive requirements whether it could be supply or consumption. It may be observed that the excitation control to produce reactive power costs more than a capacitor. Moreover, it is not safe for generator field winding which can be damaged due to frequent and overuse of the excitation system to produce or consume the reactive power [75]. Considering the fact, a mix of PV, WTs and SCs are considered in Case IV to be simultaneously integrated in this system. The inclusion of SCs with DGs has shifted the WT from node-14 to node-8 with the same generation capacity obtained in Case III. On the other hand, the solar PV capacity is reduced without changing the node. Since this study has considered the commercially available fixed-set of WT sizes therefore the deployment capacity remains unaltered for a very small change of size in these cases. The commercially available WT sizes considered in this work are as follows: 500 kW, 850 kW, 1250 kW and 1500 kW [68]. A mixed DER integration has significantly reduced the annual energy loss while improving node voltage deviation and stability index. However, the presence of SCs has not reduced the STDG much. Finally, it can be concluded that the mixing of DERs, along with their operation at different power factors, has a significant effect on distribution networks performance. Fig. 10 presents the mean node voltage profiles of the system for all cases, over a time period of . It shows that the mean voltages of all nodes are significantly improved as compared to Case 0, i.e. base case. It is also observed that the WTs alone are not able to enhance the mean node voltages beyond a certain limit. For example, the mean voltages of the feeder connecting nodes from 25 to 33 are not uniformly improved as compared to the other nodes in this case, Case I. However, the inclusion of a PV system at node 29 uniformly improves the mean node voltage profile of the system in case II. The possible cause for this improvement could be the solar power production during the peak load demand in day time which is not true for WT in the previous case. With the same type of generation configuration, the reactive power support from WTs in case III further shift the mean node voltage profile towards the nominal value. Finally, the reactive power support from shunt capacitors highly improves the node voltages of the systems along with the same generation configuration of WTs and PV system. Similarly, Fig. 11 presents the net hourly power intake of active distribution system from the main grid for all cases. In case 0, the main grid supplies complete system demand therefore highest hourly load demand is observed in this case. The DER deployment in cases I to IV is significantly reduced the system load demand. The highest load demand reduction has been observed in case II, where a renewable mix of two WTs and one PV is operated at UPF. The ideal mix of PV and wind power generation has reasonably reduced the peak load demand. The net load demand of the system, in case III and IV, appears about the same due to the similar sizing of renewables with no effect of SCs and DER nodes on real power demand reduction in case IV. However, the SCs and different DER nodes in case IV are contributing to the real power loss reduction creating a very small difference in the net load demand of cases III and IV. Furthermore, no grid substation backfeeding or negative demand are observed in cases I, III & IV, as shown in Fig. 11. A very low back-feed or about zero load demand can be observed in case II in early morning hours.
Further, the voltage stability indices of the lines for all system states, after deploying various DERs in cases I to IV, are presented in Fig. 12. It shows that the stability indices of most branches in all cases are least during 10:00 to 20:00 h because these are the peak load demand duration, as observed from Fig. 9(a). However, this duration and number of branches in this region are decreasing with the growing P. Singh et al. mixing of DERs from cases I to IV. On the other hand, the maximum VSI of these branches are observed in 00:00 to 5:00 h when load demand is low.

Discussion
In the previous sections, the proposed MBO and MBO-TOPSIS techniques solved the single and multiobjective optimization problems aiming to find optimal sites and sizes of dispatchable and non-dispatchable DERs in distribution networks. The MBO validation for solving the single goal DER integration problem in Section 4.1 reveals that this newly developed method has strong abilities to solve complex real-life engineering optimization problems quickly and efficiently as compared to existing techniques. This technique outperforms some well-known analytical (LSF, IA, ELF, EAM), numerical (SQP), and meta-heuristic methods, used to solve the complex single objective optimization problem of dispatchable DER in distribution systems. It also shows faster and better solution searching abilities as compared to many techniques in this category, e.g. FWA, IWO, HSA, TLBO, HGWO. Further, the proposed hybrid MBO-TOPSIS method outperforms many of the existing multiobjective techniques when applied to solve a multiobjective DER integration problem of the distribution system in Section 4.2. It provides more accurate results in comparison to existing approaches such as GA, PSO, INSGA-II, hybrid GA/PSO, QOTLBO, CTLBO, QOCSOS and IMOHSA. The MBO and hybrid MBO-TOPSIS techniques have shown great potential to solve complex single and multiobjective optimization problems with well-balanced exploration and exploitation abilities. These methods significantly optimize the network performance in terms of power loss reduction, voltage profile and stability improvement. In the future, these techniques can be explored to solve various operation management problems of distribution networks, which require computationally fast optimization algorithms.
The proposed multiobjective DER mix integration problem developed in Section 2 has demonstrated the high potential of synergy benefits in Section 4.3. The distribution system performance is significantly improved in terms of annual energy loss reduction, voltage profile and stability improvement, and demand deviation reduction. A maximum annual energy loss saving of 78.36% has been achieved when an optimal mix of WT, PVs and SCs deployed simultaneously. The mean node voltage profile of the network has improved by 9.59%. At the same time, the voltage stability margin (VSM) of distribution lines are also enhanced. It is defined as the difference between the current and maximum loading limits of the network. Fig. 13 demonstrates the increased VSM for all cases, varying from a minimum of 30% in case I to a maximum of 50% in case IV. The additional VSM has increased the network ability to connect 50% extra load in the future, without affecting the voltage stability limits of the system. It directly generates enormous benefits by deferring the network reinforcement cost. This margin can be increased further by deploying more flexible DERs in the networks, e.g. energy storage, electric vehicles. Furthermore, the reduced fluctuations in the power fed from the main grid would be helpful to manage the frequency-related services of the network since high real power demand mismatch directly affects the system frequency.

Conclusions
This paper reviews DER integration models, optimization techniques and MCDM approaches applied to determine the optimal sites, sizes and mix of DERs in distribution networks. Monarch butterfly optimization, a recently developed optimization method, is introduced to solve single and multiobjective problems of power distribution systems. A new deterministic MCDM problem is also formulated to find the optimal DER mix (dispatchable and non-dispatchable) considering five key objectives: minimization of annual energy loss, node voltage deviation, grid substation back-feed and the fluctuation induced by renewables while maximizing the voltage stability index of distribution networks. A novel hybridization of MBO and TOPSIS techniques is also proposed to solve the real-life large-scale MCDM problems. The effectiveness of the proposed methods is verified by solving both single and multiobjective dispatchable DER integration problems of a benchmark 33-bus distribution system. And then the simulation results are compared with some of the well-known optimization methods.
The key findings of this work include: • The introduced MBO method outperforms some of the wellknown analytical, numerical and meta-heuristic methods used to solve the complex single objective optimization problem of dispatchable DER in distribution systems. The method shows faster and better solution searching abilities as compared to many existing techniques in this category. • The proposed hybrid MBO-TOPSIS method has also found promising when its performance compared with some of the existing multiobjective optimization methods. It provides more accurate results than many of the existing optimization methods. • The MBO and the proposed hybrid MBO-TOPSIS techniques have shown great potential to solve complex single and multiobjective optimization problems with well-balanced exploration and exploitation abilities. • Under the deterministic MCDM optimization framework of mixed DER integration, the proposed method has significantly enhanced the performance of the 33-bus distribution system over different scenarios. It has minimized the annual energy loss by 78.36%, improved node voltage profile by 9.59% and increased the lines voltage stability margin by 50% while minimizing the fluctuations in power supplied by the main grid and renewables back-feeding to grid substation.
The single and MCDM MBO can be adopted to solve other large-scale, real-life engineering optimization problems in the future. Furthermore, the proposed MCDM problem, for the DER mix, can be extended to other flexible energy resources such as energy storage, electric vehicles.
Future work can also address the uncertainty modeling of renewable generation, load demand, EV growth and energy markets.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.