Modified reptile search algorithm for optimal integration of renewable energy sources in distribution networks

This paper introduces a Modified Reptile Search Algorithm (MRSA) designed to optimize the operation of distribution networks (DNs) considering the growing integration of renewable energy sources (RESs). The integration of RESs‐based Distributed Generation (DG) systems, such as wind turbines (WTs) and photovoltaics (PVs), presents a complex challenge due to its significant impact on DN operations and planning, particularly considering uncertainties related to solar irradiance, temperature, wind speed, consumption, and energy prices. The primary objective is cost reduction, encompassing electricity acquisition, PV and WTs unit costs, and annual energy losses. The proposed MRSA incorporates two strategies: the fitness‐distance balance method and Levy flight motion, enhancing its searching capabilities beyond standard Reptile Search Algorithm and mitigating local optima issues. The uncertainties in load demand, energy prices, and renewable energy generation are represented through probability density functions and simulated using Monte Carlo methods. Evaluation involves typical bentchmark functions and a real 112‐bus Algerian DN, comparing MRSA's efficacy with other optimization techniques. Results indicate that the proposed DN optimization program with WTs and PVs integration reduces annual costs by 21.43%, from 6.2715E + 06 to 4.9270E + 06 USD, reduce voltage deviations by 21.67%, from 77.1022 to 60.4007 USD, and enhance system stability by 2.59%, from 2.3699E + 03 to 2.4314E + 03 USD, compared with the base case.

Keywords: optimal distribution networks planning; optimal operation; renewable energy sources; reptile search algorithm; uncertainty

The overall implementation diagram of the proposed Modified Reptile Search Algorithm for the optimal operation of distribution networks.

Renewable energy integration has become an increasingly important topic in recent years, as it offers many advantages in terms of reducing the environmental impact of electricity generation. One of the most promising approaches is the optimal integration of renewable energy resources (RERs)‐based Distributed Generations (DGs), such as wind turbine (WT), and photovoltaic (PV), into the distribution network (DN). This approach presents multiple advantages, such as reducing greenhouse gas emissions, minimizing energy losses, improving voltage stability, and reducing the dependence on fossil fuels.

In this context, the optimization of RERs integration in DNs has gained significant attention in the literature, as it represents a crucial step towards achieving a sustainable and dependable energy system. For this aim, the integration of RERs‐based DG units at different levels of the electrical networks has been recognized as a promising solution to address both the increasing demand for energy as well as to reduce the environmental impact. Additionally, the integration of such systems has shown multiple technoeconomic advantages reducing the total production cost of electricity, as well as enhancing system performance and voltage profiles.[1] However, the presence of DGs in DNs also introduces significant uncertainty, which greatly increases the intricacy of the optimal operation.[2] The main sources of uncertainty in DNs are the load demand, the prices of electricity, and the output power of renewable DGs, which are subject to variations depending on weather conditions (solar radiation, temperature, and wind speed). As a result, determining the optimal placement of DGs in DNs is a challenging and strenuous task.[3] The selection of appropriate combinations of PV and WT‐based DGs has the potential to increase the efficiency and reliability of the DNs by addressing the issues caused by their variable nature.[4] DGs can be strategically located and operated in the network to defer major system upgrades, improve voltage regulation, minimize distribution power losses, relieve heavily loaded feeders, and enhance equipment reliability.[5] However, inappropriate DGs location and sizing can lead to various negative impacts on the DN, such as voltage instability, increased power losses, harmonic distortion, and even equipment damage, which can affect the reliability and quality of power supply, potentially leading to economic losses and customer dissatisfaction.[6] Therefore, an optimal operation of the system is necessary to make full use of the benefits of DGs while mitigating their adverse effects.

A considerable amount of research has been conducted on optimal DG integration from various perspectives. In one study, Mansouri et al.[7] delved into numerous technological challenges that arise from the widespread penetration of PV into DNs, such as active power reduction, frequency regulation, energy storage, and reactive power injection. Another study[8] explored the technical obstacles that arise due to the extensive and intensive integration of PV into the DNs, such as frequency disturbances and voltage limit violations, as well as the resulting stability issues. The implications of integrating PV on a large scale into the DNs were also considered, along with potential solutions. The search‐based dragonfly algorithm was introduced as a proposed solution for the ideal allocation of DGs, taking into account uncertainties in load demand and DGs. In Prakash and Khatod,[9] various strategies for locating and sizing DGs in DNs were evaluated. In Zellagui et al.,[10] two popular methods, the particle swarm optimization (PSO) and the firefly method, were proposed as means of achieving operational, financial, and environmental benefits in DNs. In a recent study,[11] a combined approach utilizing the PSO and hybrid enhanced gray wolf algorithms was proposed for the optimum location for DGs. This approach aimed to minimize active power losses, system costs, and emissions while improving the voltage stability index (VSI) and reducing the voltage deviation index (VDI). A bioinspired algorithm based on Monte Carlo simulation was proposed by Hemeida et al.[12] In Rathore and Patidar,[13] the PSO coefficient was employed to minimize total energy losses by optimizing the sizing and deployment of DGs. In El Sehiemy et al.,[14] the Slap swarm optimizer was utilized to enhance the operational, financial, and environmental performance of power plants. The equilibrium optimizer method was introduced by Ahmed et al.[15] A highly effective algorithm was proposed for improving the size and placement of DGs within power networks and addressing the issues with microenergy management. In a study by Thokar et al.,[16] a bilayer approach to the placement of energy storage systems and PV in DNs was proposed. In Biswal and Shankar,[17] the Strength Pareto Evolutionary Algorithm 2 was proposed as a solution for the problem of capacitor and DGs placement with load uncertainty. In Hadidian‐Moghaddam et al.,[18] an innovative ant lion optimizer approach was proposed with the goal of lowering energy costs, reducing losses and voltage deviation (VD), and improving reliability. In Ahmed et al.,[19] the PSO approach was proposed, utilizing a probabilistic uncertainty modeling approach, for the siting and sizing of DGs. In Rao et al.,[20] Monte Carlo simulation was used for the best placement of DGs within power systems. In this research, Ullah et al.[21] create an energy optimization framework for intelligent microgrids with the aim of reducing operational expenses, minimizing emissions, and enhancing availability. In Ali et al.,[22] a demand side management strategy is proposed to optimize operational cost, pollution emission, and load coordination in smart grid using multiobjective wind‐driven optimization and modeling consumer behavior the paper uses probability density function (pdf) to forecast wind speed for integrating wind power. Hafeez et al.[23] propose an energy management strategy using price‐based demand response programs in an internet‐of‐things‐enabled smart grid to schedule smart home appliances. In Hafeez et al.,[24] a modular framework for efficient load scheduling in the smart grid, utilizing a hybrid algorithm and real‐time pricing data to reduce costs and peak demands, benefiting both residents and power companies.

In Sakr et al.,[25] a revised differential evolution algorithm was introduced for the best positioning of DGs. Additionally, the gorilla troops optimizer, which models gorilla social behavior and movement, has been applied to various engineering problems, including PV model extraction.[26]

The incorporation of the Levy flight distribution has been found to enhance the performance of optimization algorithms for engineering global optimization problems, leading to improved optimal solutions, as demonstrated in several studies.[27] Additionally, balancing fitness and distance has been shown to enhance the algorithm's capabilities for searching, as reported in various studies.[28] Overall, these findings highlight the potential of incorporating these techniques into optimization algorithms for engineering design problems.

Table 1 lists a comparison between the presented work and other related references for the Optimal Operation problem (OOP) of DNs with RERs.

1 Table Comparison between the presented work and other related references for the Optimal Operation of distribution networks.

1 Abbreviations: DG, Distributed Generation; DSTATCOM, distributed static compensator; PV, photovoltaic; WT, wind turbine.

This research paper proposes a modified version of the Modified Reptile Search Algorithm (MRSA) for the optimal operation of DNs. The proposed approach is compared against several established optimization techniques, including Sand Cat Swarm Optimization (SCSO),[29] PSO,[30] Dandelion Optimizer (DO),[31] Sine Cosine Algorithm (SCA),[32] improved harmony search (IHS),[33] and the conventional Reptile Search Algorithm (RSA).[28] The study considers the integration of PV and WT generators into DNs with the primary objective of cost minimization and the summation of VD and VSI over a 24‐h planning horizon. The model's performance is assessed in the presence of uncertainty in load, wind speed, solar irradiation, temperature, and energy prices for purchasing. Weather data, such as irradiance and temperature, are crucial input variables for simulating PV, whereas the output power of WT is influenced by several factors, such as turbine size and wind speed. Therefore, an accurate modeling approach is essential due to the sensitivity of both PV and WT. To this end, the study employs actual data on irradiance, temperature, and wind speed from the DN area, namely, from the real 112‐bus Algerian DN, to guarantee the precision and practicality of the findings.

The novelty and significant contributions of the developed algorithm compared with existing works are articulated as follows:

• 1. This study proposes an MRSA specifically designed to solve the optimal planning problem involving RERs within DNs.
• 2. Investigating the synergies between RERs (WTs and PVs), thereby enhancing grid stability and reducing reliance on expensive grid energy purchases.
• 3. The effectiveness and reliability of the MRSA algorithm are rigorously demonstrated through extensive statistical comparisons with established optimization techniques.
• 4. Rigorous uncertainty analysis to account for the stochastic nature of key parameters, including load demand, solar irradiation, wind speed, temperature, and energy pricing. This facilitates the robust optimization of the distribution grid operation.
• 5. The proposed MRSA algorithm is tested on standard benchmark functions and applied to address the optimal operation challenges of a real 112‐bus Algerian DN.

This paper is organized as follows: In Section 2, the optimal planning problem formulation. Section 3 outlines the modeling approach for uncertainty parameters. Sections 4 and 5 explain the mathematical formulation of the RSA and the proposed MRSA, respectively. Section 6 provides the simulation results. Lastly, the conclusions of this paper are depicted in Section 7.

This section describes the problem formulation that includes the objective functions and the corresponding constraints of the optimal operation.

In this paper, three objective functions are considered which include the following:

The objective function taken into consideration comprises the cost of electricity acquired from the network ( $${C}_{{Grid}}$$ ), the cost of PV units ( $${C}_{{PV}}$$ ), the cost of WT ( $${C}_{{Wind}}$$ ), and the yearly cost of energy loss ( $${C}_{{Loss}}$$ ), and it may be expressed as follows[34]: 1 $$$C=min({C}_{{Grid}}+{C}_{{PV}}+{C}_{{Wind}}+{C}_{{Loss}}).$$$

In which, 2 $$${C}_{{Grid}}=365\times \sum _{h=1}^{24}{P}_{{Grid}(h)}\times {U}_{{Grid}(h)},$$$ where $${U}_{{Grid}}$$ represents the cost of purchasing electricity from the grid, and $${P}_{{Grid}(h)}$$ refers to the power withdrawn from the grid.[35]3 $$${C}_{{Loss}}=365\times {U}_{{Loss}}\times \sum _{h=1}^{24}{P}_{T{\rm{\_}}{Loss}(h)},$$$ where $${U}_{{Loss}}$$ is the cost of energy, and $${P}_{{T\_Loss}(h)}$$ refers to the total power losses. 4 $$${C}_{{PV}}={C}_{{PV}}^{{inst}.}+{C}_{{PV}}^{O{\rm{\&}}M},$$$ where $${C}_{{PV}}^{{inst}.}$$ is the cost of installing the PV, $${C}_{{PV}}^{O{\rm{\&}}M}$$ is the PV unit's operating and maintenance costs.[34]5 $$${C}_{{PV}}^{O{\rm{\&}}M}={U}_{{PV}}^{O{\rm{\&}}M}\times \sum _{h=1}^{24}{P}_{{PV}(h)},$$$ 6 $$${C}_{{PV}}^{{inst}.}={{CF}\times {U}_{{PV}}\times P}_{{rated}{\rm{\_}}{PV}},$$$ where $${CF}$$ is a factor affecting capital recovery, $${P}_{{rated\_PV}}$$ is the rated produced power of the PV.[34]7 $$${C}_{{Wind}}={C}_{{wind}}^{{inst}.}+{C}_{{wind}}^{O{\rm{\&}}M},$$$ where $${C}_{{Wind}}^{{inst}.}$$ is the cost of installing the WT, and $${C}_{{wind}}^{O{\rm{\&}}M}$$ is the wind's operating and maintenance costs.[34]8 $$${C}_{{Wind}}^{O{\rm{\&}}M}={U}_{{Wind}}^{O{\rm{\&}}M}\times \sum _{h=1}^{24}{P}_{{Wind}(h)},$$$ where $${U}_{{PV}}^{O{\rm{\&}}M}$$ and $${U}_{{Wind}}^{O{\rm{\&}}M}$$ represent the WTs and PVs operation and maintenance costs ($/kW).[34]9 $$${C}_{{Wind}}^{{inst}.}={{CF}\times {U}_{{WT}}\times P}_{{rated}{\rm{\_}}{wind}},$$$ where $${P}_{{rated\_wind}}$$ is the rated produced power of WT, $${U}_{{WT}}$$ and $${U}_{{PV}}$$ represent the WTs and PVs purchased cost ($/kW), and $${P}_{{Wind}(h)}$$ and $${P}_{{PV}(h)}$$ represent the yielded power from WTs and PVs systems. 10 $$${CF}=\frac{\beta \times {(1+\beta)}^{{NP}}}{{(1+\beta)}^{{NP}}-1},$$$ where $$\beta $$ and $${NP}$$ the interest rate and system lifetime of the winds turbine or PV unit.[36]

The generated power from the PV can be calculated as follows: 11 $$${T}_{c}={T}_{a}+\frac{{I}_{s}}{800}\cdot ({T}_{N}-20),$$$ 12 $$${P}_{{pv}}(t)={A}_{{pv}}\cdot {\eta }_{{pv}}(t)\cdot I(t).$$$

The surface area utilized by the set of PV, denoted as $${A}_{{pv}}$$ in (m2), is multiplied by a constant value representing the conversion efficiency of the PV panels, for the purpose of calculating the power generated by the PV, denoted as $${P}_{{pv}}$$ (kW), $${\boldsymbol{I}}$$ solar insolation in (kW/m2) and $${\eta }_{{pv}}$$ represent the instantaneous efficiency of PV panels. The instantaneous efficiency of PV panels is obtained using the following equation[37]: 13 $$$\begin{array}{c}\begin{array}{ccc}{\eta }_{pv}(t) & = & {\eta }_{r}\cdot {\eta }_{t}\times \left[1-\beta \cdot ({T}_{a}(t)-{T}_{r})-\beta \cdot I(t)\cdot \left(\frac{NOCT-20}{800}\right)\cdot (1-{\eta }_{r}\cdot {\eta }_{t})\right].\end{array}\end{array}$$$

The effectiveness of the maximum power point tracking equipment is denoted by $${\eta }_{t}$$ , and $${\eta }_{r}$$ represents the reference efficiency of the PV panels. The temperature coefficient of efficiency β, typically ranging from 0.004 to 0.006/°C for silicon cells, is also considered. The ambient temperature $${T}_{a}$$ (°C), PV cell reference temperature $${T}_{r}$$ (°C), and nominal operating cell temperature $${NOCT}$$ (°C) are also taken into account.

The WT's produced output power ( $${P}_{{WT}}$$ ) is computed as follows: 14 $$${P}_{{WT}}(W){\rm{\hspace{0.33em}}}{\rm{\hspace{0.33em}}}=\left\{\begin{array}{ccc}0 & \mathrm{for}\unicode{x02007}W\lt {W}_{i}\unicode{x02007}\mathrm{and}\unicode{x02007}W\gt {W}_{o}, & \\ {P}_{{rated}{\rm{\_}}{wind}}\left(\frac{W-{W}_{i}}{{W}_{r}-{W}_{i}}\right) & \mathrm{for}\unicode{x02007}({W}_{i}\le W\le {W}_{r}), & \\ {P}_{{rated}{\rm{\_}}{wind}} & \mathrm{for}\unicode{x02007}({W}_{r}\lt W\le {W}_{o}), & \end{array}\right.$$$ the rated power of the used WT is 250 kW while its rated velocity $$({W}_{r})$$ is 15 m/s, the cut‐in speed ( $${W}_{i}$$ ) is 2.5 m/s, and cut‐out speed ( $${W}_{o}$$ ) the WT taken as 25 m/s.[38]

To improve the network performance, the voltage deviations should be kept closed to 1 p.u. value which can be defined as[[36], [39]]15 $$$\sum {VD}=\sum _{h=1}^{24}\sum _{k=1}^{{NB}}|({V}_{k}-1)|,$$$ where $${NB}$$ represents the number of buses in the network, and $${V}_{k}$$ represents the voltage of the kth bus.

The system stability index statement is as follows[[40]]: 16 $$${{VSI}}_{k}={|{V}_{k}|}^{4}-4{({P}_{k}{X}_{{km}}-{Q}_{k}{R}_{{km}})}^{2}-4({P}_{k}{X}_{{km}}+{Q}_{k}{R}_{{km}}){|{V}_{k}|}^{2},$$$ 17 $$$\sum {VSI}=\sum _{h=1}^{24}\sum _{k=1}^{{NB}}{{VSI}}_{k},$$$ where $${{VSI}}_{k}$$ is the voltage stability index, $${R}_{{km}}$$ represent the resistance of the transmission lines while the $${X}_{{km}}$$ is its reactance. $${P}_{k}$$ and $${Q}_{k}$$ define the real and reactive power at bus, respectively. The following three objective functions are taken into consideration simultaneously: 18 $$$F={\varepsilon }_{1}{F}_{1}+{{\varepsilon }_{2}F}_{2}+{\varepsilon }_{3}{F}_{3},$$$ 19 $$${F}_{1}=\frac{{C}_{{RERs}}}{{C}_{{Base}}},$$$ 20 $$${F}_{2}=\frac{{\sum {VD}}_{{RERs}}}{{\sum {VD}}_{{Base}}},$$$ 21 $$${F}_{3}=\frac{1}{{\sum {VSI}}_{{RERs}}},$$$ where $${RERs}$$ and $${Base}$$ are subscripts refers to with RERs and the base case, respectively. $${\varepsilon }_{1}$$ , $${\varepsilon }_{2}$$ , and $${\varepsilon }_{3}$$ are the weighted factors that were selected to be 0.5, 0.25, and 0.25, respectively.[42]

22 $$${V}_{{Min}}\le {V}_{i}\le {V}_{{Max}},$$$ 23 $$${P}_{{PV}{\rm{\_}}{rated}}+{P}_{{wind}{\rm{\_}}{rated}}\le \sum _{i=1}^{{NB}}{P}_{{Load},i},$$$ 24 $$${{PF}}_{{Min}}\le {PF}\le {{PF}}_{{Max}},$$$ 25 $$${I}_{n}\le {I}_{max,n},\,n=1,2,3,{\rm{\ldots }},{NT},$$$ where $${V}_{min}$$ and $${V}_{max}$$ are the lower and upper voltage limits. $${P}_{Load}$$ and $${Q}_{Load}$$ signified the real and reactive load, respectively; $${I}_{max,n}$$ is the maximum allowable current limit of the line; $${NT}$$ defines the number of lines; $${{PF}}_{max}$$ and $$P{F}_{min}$$ are the maximum and minimum of the WT power factor, respectively.

26 $$${P}_{S}+{P}_{{PV}}+{P}_{{Wind}}=\sum _{i=1}^{{NT}}{P}_{{Loss},i}+\sum _{i=1}^{{NB}}{P}_{{Load},i},$$$ 27 $$${Q}_{S}+{Q}_{{Wind}}=\sum _{i=1}^{{NT}}{Q}_{{Loss},i}+\sum _{i=1}^{{NB}}{Q}_{{Load},i},$$$ where $${Q}_{S}$$ and $${P}_{S}$$ are the reactive and real powers of the main network.

In this section, the considered uncertain parameters are represented as follows:

The Beta pdf has been employed to model the intermittent nature of the solar irradiance, as follows[[43]]: 28 $$$f({g}_{s})=\left\{\begin{array}{cc}\frac{\Gamma (\alpha +\beta)}{\Gamma (\alpha)\Gamma (\beta)}{s}^{(\alpha -1)}{(1-{g}_{s})}^{(\beta -1)}, & 0\le {g}_{s}\le 1,\alpha ,\beta \ge 0,\\ 0, & \text{otherwise},\end{array}\right.$$$ where $$\sigma $$ is the standard deviation while $$\mu $$ is the mean value which has been obtained from the historical data. $$\beta $$ and $$\alpha $$ can be obtained from the following equations[[45]]: 29 $$$\beta =(1-\mu)\times \left(\frac{\mu \times (1+\mu)}{{\sigma }^{2}}-1\right),$$$ 30 $$$\alpha =\frac{\mu \times \beta }{1-\mu }.$$$

For modeling the wind speed uncertainty, the Weibull pdf is utilized which can be described as follows[[47]]: 31 $$$f(W)=\left(\frac{k}{c}\right){\left(\frac{W}{c}\right)}^{k-1}\text{exp}\left[-{\left(\frac{W}{c}\right)}^{k}\right],$$$ where $$W$$ represents the wind speed. $$c$$ and $$k$$ are the scale and shape parameters of the Weibull pdf.

Normal pdf is employed to model the demanding uncertainty that can be mathematically described as follows[46]: 32 $$$f(L)=\frac{1}{{\sigma }_{l}\sqrt{2\pi }}\times \text{exp}\left[-\left(\frac{l-{\mu }_{l}}{2{\sigma }_{l}^{2}}\right)\right],$$$ where $${\sigma }_{l}$$ and $${\mu }_{l}$$ are the standard deviation and the mean value of the loading, respectively. The probability of the loading is divided into subsegments according to the following equation:

One of the most significant random characteristics in power systems is the price of electricity, which is an unreliable parameter that is acquired from the grid. To model the probability distribution of the electricity price, the Normal pdf can be used based on its mean value $${\mu }_{{EP}}$$ and standard deviation $${\sigma }_{{EP}}$$ , as shown in Equation (33)[49]: 33 $$$f(P)=\frac{1}{{\sigma }_{{EP}}\sqrt{2\pi }}\text{exp}\left[-\frac{{({EP}-{\mu }_{{EP}})}^{2}}{2{({\sigma }_{{EP}})}^{2}}\right].$$$

As the surrounding temperature varies continuously. Thus, the temperature is considered as uncertain parameter and it is assumed that the uncertainty of the temperature is represented using the normal probability distribution for modeling the uncertainty of the temperature as follows: 34 $$$f(T)=\frac{1}{\sqrt{2\pi {\sigma }_{T}^{t}}}\text{exp}\left[-\frac{(s{\rm{\mbox{--}}}{\mu }_{T}^{t})2}{2{\left({\sigma }_{T}^{t}\right)}^{2}}\right],$$$ where $${\sigma }_{T}$$ is the standard deviation while $${\mu }_{T}$$ is the mean value of temperature.

The RSA draws inspiration from the hunting and encircling behaviors of crocodiles in the wild. These behaviors involve crocodiles working together to surround and capture their prey.[28]

RSA uses mathematical modeling to mimic these behaviors and create an optimization process that is both gradient‐free and population‐based. This means that it can tackle optimization challenges of varying complexities, with or without specific constraints.

RSA begins the optimization by generating a set of possible solutions, represented as $$X$$ in Equation (35), through a stochastic approach. The algorithm then identifies the best solution obtained and considers it as an approximation of the optimal solution in each iteration. 35 $$$X=\left[\begin{array}{ccccc}{x}_{1,1} & \cdots & {x}_{1,j} & {x}_{1,n-1} & {x}_{1,n}\\ {x}_{2,1} & \cdots & {x}_{2,j} & \cdots & {x}_{2,n}\\ \cdots & \cdots & {x}_{i,j} & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {x}_{N-1,1} & \cdots & {x}_{N-1,j} & \cdots & {x}_{N-1,n}\\ {x}_{N,1} & \cdots & {x}_{N,j} & {x}_{N,n-1} & {x}_{N,n}\end{array}\right].$$$

The set of candidate solutions, $$X$$ , used in the RSA is generated randomly through the use of Equation (36), $${x}_{i,j}$$ represents the value of the solution at the $$j\mathrm{th}$$ position within the $$i\mathrm{th}$$ candidate solution. The number of candidate solutions generated is denoted as $$N$$ , while the dimension size of the problem is represented by $$n$$ . 36 $$${x}_{{ij}}=rand\times ({UB}-{LB})+{LB},\unicode{x02007}\unicode{x02007}j=1,2,{\rm{\ldots }},n.$$$

Two exploration search methods are utilized in the RSA. Each component is given a stochastic scaling factor to get more varied solutions and investigate different locations. The position‐updating equations used in the exploration phase are designed to imitate the surrounding behavior of crocodiles and are presented in Equation (37). It is noteworthy that the algorithm employs a simple rule to facilitate this behavior. 37 $$${x}_{(i,j)}(t+1)=\left\{\begin{array}{cc}Bes{t}_{j}(t)\times -{\eta }_{(i,j)}(t)\times \beta -{R}_{(i,j)}(t)\times rand, & t\le \frac{T}{4},\\ Best(t)\times {x}_{({r}_{1},J)}\times {ES}(t)\times rand, & t\le 2\frac{T}{4}\unicode{x02007}\text{and}\unicode{x02007}t\gt \frac{T}{4}.\end{array}\right.$$$

The values used in this calculation include $$Bes{t}_{j}(t)$$ , which denotes the position of the $$j\mathrm{th}$$ element in the optimal solution obtained until the current point. The integer $$rand$$ is created at random and ranges from 0 to 1, $${T}$$ , the upper limit of iterations and $$t$$ , the current iteration number. Equation (38) is used to calculate the hunting operator, denoted as ( $$i$$ ,), for the jth location in the ith solution. Iterations are conducted with a fixed sensitive parameter, β, set to 0.1, to regulate the exploration accuracy during the encircling phase. A reduction function, denoted as $${R}_{(i,j)}$$ , is employed to narrow down the search area and its calculation is determined by Equation (39). The random number $${r}_{1}$$ , generated between [1 $$N$$ ], is used to represent a stochastic location in the ith solution, represented by $${x}_{({r}_{1},J)}$$ . The value of $$N$$ represents the number of candidate solutions. The Evolutionary Sense $${ES}(t)$$ probability ratio varies randomly between 2 and −2 during each iteration according to Equation (40). 38 $$${\eta }_{(i,j)}=Bes{t}_{j}(t)\times {P}_{(i,j)},$$$ 39 $$${R}_{(i,j)}=\frac{Bes{t}_{j}(t)-{x}_{({r}_{2},j)}}{Bes{t}_{j}(t)+\epsilon },$$$ 40 $$${ES}(t)=2\times {r}_{3}\times \left(1-\frac{1}{T}\right).$$$

This equation describes some variables used in RSA, including $${r}_{2}$$ , which is a random number between [1 $$N$$ ], and $$\epsilon $$ , which is a small value The equation also involves the use of $${r}_{3}$$ , which represents a random value among −1 and 1. $${P}_{(i,j)}$$ is another variable used in the algorithm, which represents the proportional difference between of the $$j\mathrm{th}$$ position of the best‐obtained solution and the $$j\mathrm{th}$$ position of the current solution. This calculation is performed using Equation (41). 41 $$${P}_{(i,j)}=\alpha +\frac{{x}_{(i,j)}-M({x}_{i})}{Bes{t}_{j}(t)\times (U{B}_{(j)}-L{B}_{(j)})+\epsilon }.$$$

This excerpt explains the variables used in the hunting cooperation phase. The average position of the current solution is denoted by $$M({x}_{i})$$ , which is calculated using Equation (41). $$U{B}_{(j)}$$ and $$L{B}_{(j)}$$ represent the upper and lower boundaries of the $$j\mathrm{th}$$ position, respectively. The parameter $$\alpha $$ is fixed at 0.1 in this study. 42 $$$M({x}_{i})=\frac{1}{n}\sum _{j=1}^{n}{x}_{(i,j)}.$$$

The RSA algorithm uses two main search strategies, namely, hunting cooperation and hunting coordination, to investigate the search space and locate an ideal solution. These strategies are modeled in Equation (43) and are used for exploitation mechanisms. 43 $$${x}_{(i,j)}(t+1)=\left\{\begin{array}{cc}Bes{t}_{j}(t)\times {P}_{(i,j)}(t)\times rand, & t\le 3\frac{T}{4}\unicode{x02007}\text{and}\unicode{x02007}t\gt 2\frac{T}{4},\\ Bes{t}_{j}(t)-{\eta }_{(i,j)}(t)\times \epsilon -{R}_{(i,j)}(t)\times rand. & t\le T\unicode{x02007}\text{and}\unicode{x02007}t\gt 3\frac{T}{4},\end{array}\right.$$$

The $$j\mathrm{th}$$ location in the optimal solution obtained until now is represented by the term $$Bes{t}_{j}(t)$$ . The hunting operator applied to $$j\mathrm{th}$$ location in the $$i\mathrm{th}$$ solution, denoted as $${\eta }_{(i,j)}$$ , is calculated using Equation (38). The disparity percentage of the $$j\mathrm{th}$$ element in the current solution compared with the $$j\mathrm{th}$$ element in the best‐obtained solution is denoted by $${P}_{(i,j)}$$ , and is calculated using Equation (41). $${R}_{(i,j)}$$ is a parameter utilized to shrink the search space, and its value is determined using Equation (39).

The suggested MRSA is based on two improvement methods. The first modification aims to improve the exploration phase of the RSA using the fitness‐distance balance (FDB) approach. The FDB is an efficient selection method that can guide the algorithm to the global solution.[[50], [52]] The FDB selection method is based on score values of the solution candidates. Then, the locations of the crocodiles are modified based on the fitness values and the distance between the candidates and the best solution. The distance values are calculated as follows: 44 $$${{DS}}_{i}=\sqrt{{\left(({x}_{i}^{1}-Bes{t}^{1})\right)}^{2}+{\left(({x}_{i}^{2}-Bes{t}^{2})\right)}^{2}+\cdots +{\left(({x}_{i}^{d}-Bes{t}^{d})\right)}^{2}}.$$$

Then, construct the vectors of the fitness and the distance values as follows: 45 $$${DS}=[D{S}_{1},D{S}_{2},{\rm{\ldots }},D{S}_{n}],$$$ 46 $$$F=[{F}_{1},{F}_{2},{\rm{\ldots }},{F}_{n}].$$$

Then the distance and the fitness value can be normalized as follows: 47 $$$normD{S}_{i}=\frac{D{S}_{i}-min({DS})}{max({DS})-min({DS})}m,$$$ 48 $$$norm{F}_{i}=\frac{{F}_{i}-min(F)}{max(F)-min(F)},$$$ where $$min$$ and $$max$$ are the minimum and maximum in the distance and the fitness vectors. The following formula is used to calculate the FDB score: 49 $$$FDB\,scor{e}_{i}=\alpha \,\ast \,(1-norm{F}_{i})+(1-\alpha)\,\ast \,normD{S}_{i}.$$$

In which 50 $$$\alpha =0.5\,\ast \,\left(1+\frac{t}{{T}_{max}}\right).$$$

The second modification is based on boosting the exploitation phase of the RSA technique by updating the locations of the crocodiles around the best solution using the Levy flight distribution according to (51).[[54], [56]]51 $$${x}_{(i,j)}(t+1)={r}_{3}Best(t)-{r}_{4}{x}_{(i,j)}(t)+{C}_{1}\cdot {L}_{F}\cdot ({x}_{(r,j)}(t)-{x}_{(i,j)}(t)),$$$ where $${r}_{3}$$ and $${r}_{4}$$ refer to a random value in the range (0, 1). $${x}_{(i,j)}(t)$$ and $${x}_{(r,j)}(t)$$ represent the current location of the crocodiles and the best location, respectively. $${C}_{1}=2{r}_{4}(1-T/{T}_{max})$$ represents an operator that measures the intensity of the Levy flight. $${L}_{F}$$ is the Levy flight function which can be measured as follows: 52 $$${L}_{F}=0.05\times \frac{u\times \sigma }{|v{|}^{1/\beta }}.$$$

In which 53 $$$\sigma ={\left(\frac{{\rm{\Gamma }}(1+\beta)\times \sin (\pi \beta /2)}{{\rm{\Gamma }}((1+\beta)/2)\times \beta \times {2}^{(\beta -1)/2}}\right)}^{1/\beta },$$$ where $$u$$ and $$v$$ denote random values that can be obtained from the normal distribution. β refers to a constant that is 1.5.

Figure 1 shows the MRSA for the optimal operation solution. The methodology for solving the stochastic optimal operation of a distribution grid is depicted in Figure 2.

The suggested MRSA is applied to solving the real 112‐bus Algerian DN. Initially, The proposed approach's performance was assessed through a comparison with other established optimization methods, such as SCSO,[29] PSO,[30] DO,[31] SCA,[32] IHS,[33] and the conventional RSA.[28] All the algorithms, including the proposed MRSA, were programmed in MATLAB software (MATLAB R2019b) and executed on a computer with an Intel i7, 2.5 GHz CPU, and 6 GB RAM. The studied cases are presented as follows:

Here, the proposed MRSA is utilized for 23 classic functions, including the multimodal, unimodal, and fixed‐dimension multimodal functions which have been listed in Appendix A.[[58]] For all cases, the parameters are set according to Table 2 and the obtained results are represented over 30 run times.

2 Table Selected parameters of the optimizers.

2 Abbreviations: DO, Dandelion Optimizer; HMCR, harmony memory considering rate; IHS, improved harmony search; MRSA, Modified Reptile Search Algorithm; PAR, pitch adjusting rate; PSO, particle swarm optimization; RSA, Reptile Search Algorithm; SCA, Sine Cosine Algorithm; SCSO, Sand Cat Swarm Optimization.

This section depicts the performance of the proposed MRSA compared SCSO,[29] PSO,[30] DO,[31] SCA,[32] IHS,[33] and the conventional RSA.[28] Table 3 shows the statistical results, including the mean, the worst, the average, the best, and the Wilcoxon p value between the proposed MRSA and the other algorithms. According to Table 3, the proposed MRSA optimizer is superior in terms of the mean, best, and worst values for the F1–F4, and F17–F23 while the obtained results for the reported algorithms are similar for F17–F23. In addition to that for F15 and F5, the RSA is better than MRSA in the value of the best score but the mean value of the MRSA is better. The p values of the Wilcoxon test compared with MRSA and the other optimizers are included in the 7th column of Table 3. When the p value is less than 5%, it is a significant difference between algorithms.[60] On the opposite of that when is more than 5%, no significant difference between optimizers. In addition to that if the results of different algorithms are identical, the p value will be N/A. From the p value, it is clear that the proposed MRSA has significant difference compared with SCSO, PSO, DO, SCA, IHS, and the standard RSA for most of the studied benchmark functions, the population size and maximum iteration for all algorithms are 30 and 300, respectively.

3 Table Statistical results of different optimizers for the standard functions.

• 3 Note : Bold values indicate the best obtained solutions.
• 4 Abbreviations: DO, Dandelion Optimizer; IHS, improved harmony search; MRSA, Modified Reptile Search Algorithm; PSO, particle swarm optimization; RSA, Reptile Search Algorithm; SCA, Sine Cosine Algorithm; SCSO, Sand Cat Swarm Optimization.

The convergence carves of the MRSA and other reported techniques including SCSO,[29] PSO,[30] DO,[31] SCA,[32] IHS,[33] and the conventional RSA[28] can be realized from Figure 3. According to the convergence carves, the proposed MRSA has the best convergence speed for the fixed‐dimension functions, the unimodal functions, and the multimodal functions. In addition to that, the suggested MRSA is superior and converged to the optimal solution faster than the conventional RSA due to the proposed modifications which can boost the exploration and exploitation phases of the proposed algorithm.

Boxplots are ideal for displaying data distributions in quartiles that can be used to show the characteristics of data distribution. Figure 4 depicts the Boxplots of the studied MRSA and the other reported optimization algorithms. It is obvious that the boxplots of the MRSA are narrower compared with the other optimizers.

In this section, the proposed MRSA is applied to solve the optimal operation and determine the best location for solar PV unit ratings and the WTs on 112‐bus Algerian DN. The system, shown in Figure 5, consists of 112 buses and 111 branches with a cumulative load of 3367.60 kW and 3725.70 kVAR. The meteorological data for the wind speed, temperature, and solar radiation are collected over a period of 3 years from Ehsan and Yang.[61] Also, the load demand data are collected for 3 years which given in Kaur et al.[62] To validate the effectiveness of the proposed MRSA algorithm, the obtained results have been with those by other meta‐heuristic optimization algorithms, including RSA, SCSO, DE, PSO, SCA, and IHS. To ensure a fair and valid comparison, the maximum number of iterations and populations selected for the proposed algorithms were set at 60 and 25, respectively. The presence of uncertainties, including load variations, wind speed, energy price, temperature, and solar irradiance, were considered in accordance with the methodology outlined in Section 3. Figure 6 illustrates the expected day‐ahead load demand, while Figure 7 displays the market price for purchasing energy. In addition, Figures 8–10 present the expected wind speed, irradiance, and temperature, respectively. Table 4 details the cost coefficients for RERs as well as the operational limitations.

4 Table Constraints and the cost coefficients.

5 Abbreviations: PF, power flow; PV, photovoltaic; WT, wind turbine.

The main objective of this paper is to optimize the total annual cost while improving system performance. Initially, the base case was considered without any RERs integrated into the DN. The results showed that the total annual purchased energy from the grid was 7.2657E + 04 kWh, with energy losses of 2.4901E + 03 kW. The total purchase energy cost was 6.2170E + 06 USD, and the energy loss cost was 5.4533E + 04 USD, resulting in a total annual cost of 6.2715E + 06 USD. The summation of the VD was 77.1022 p.u. and the VSI was 2.3699E + 03 p.u. Table 5 shows the simulation results which have been obtained by different algorithms. According to Table 5, with the application of the proposed MRSA with the ideal allocation of the RERs, the total annual cost has been from 6.2715E + 06 to 4.9270E + 06 USD and the voltage deviations have been reduced from 77.1022 to 60.4007 p.u., while VSI has been improved from 2.3699E + 03 to 2.4314E + 03 p.u. The system voltage profiles with and without RERs are provided in Figures 11 and 12. Judging from Figures 11 and 12, the voltage profile was enhanced considerably with the inclusion of the PV units and the WTs. Additionally, as seen in Figure 13, the power losses have been significantly decreased. The optimal locations of PV units in the DN are identified as 81, 82, and 102, and the optimal locations of WTs are identified as 5, 95, and 112, respectively as shown in Figure 14. The WTs ratings are 1000, 250, and 250 kW, respectively, while the corresponding area of the solar modules are 4310, 4000, and 5000 m2. The statistical outcomes for the objective function achieved through various optimization techniques are tabulated in Table 5. As depicted in Table 6, the results obtained by the proposed MRSA algorithm outperform other optimization techniques in terms of the mean, best, and worst values. The power output of the PV units is illustrated in Figure 15, indicating that their yields fluctuate continuously with irradiance variations. Similarly, the generated powers of the WT vary with changes in wind speed, as demonstrated in Figure 16.

5 Table Energy management results of Algerian distribution network.

• 6 Note : Bold values indicate the best obtained solutions.
• 7 Abbreviations: MRSA, Modified Reptile Search Algorithm; PF, power flow; PV, photovoltaic; RESs, Renewable Energy Sources; RSA, Reptile Search Algorithm; WT, wind turbine.

6 Table Simulation results of objective function for applied algorithms.

• 8 Note : Bold values indicate the best obtained solutions.
• 9 Abbreviations: DO, Dandelion Optimizer; IHS, improved harmony search; MRSA, Modified Reptile Search Algorithm; PSO, particle swarm optimization; RSA, Reptile Search Algorithm; SCA, Sine Cosine Algorithm; SCSO, Sand Cat Swarm Optimization.

This paper solved the OOP of a real 112‐bus Algerian DN with ideal integration of the RERs, including the PV units and the WTs using a new MRSA. The proposed MRSA is based on the FDB method and the Levy flight motion selection strategies. The optimal operation was solved by considering the uncertainties of the price, the load, the wind speed, temperature, and the solar irradiance. The proposed MRSA has been applied and tested on standard benchmark functions and the obtained results were compared with other optimization algorithms, including the SCSO, DE, PSO, SCA, IHS, and the standard RSA. The numerical results of optimal inclusion of the RERs using the proposed MSA show that the total cost has been reduced from 6.2715E + 06 to 4.9270E + 06 USD, the VD has been reduced from 77.1022 to 60.4007 p.u., and an enhancement in VSI has been enhanced from 2.3699E + 03 to 2.4314E + 03 p.u. in comparison to the base case. In addition, the findings demonstrate that the proposed MRSA algorithm for solving the optimal operation surpasses other optimization techniques, including SCSO, DE, PSO, SCA, IHS, and the standard RSA, in terms of performance and effectiveness. However, this study has certain limitations, such as the absence of considerations for energy storage systems and electric vehicles (EVs). Future work will expand to encompass energy management in distribution systems, including various energy storage systems, like, batteries, compressed air, and Superconducting Magnetic Energy Storage, as well as incorporating EV charging stations.

See Tables A1–A3.

A1 Table Unimodal functions.

A2 Table Multimodal functions.

A3 Table Fixed‐dimension multimodal benchmark functions.