1 Introduction

Microgrid (MG) is the smallest part of the electrical system that consists of generators and loads. It is made up of a variety of loads, storage batteries, and a generating network. In a hybrid power system, electricity is created from renewable sources and stored in energy storage units to offer stabilization to the microgrid. A microgrid is a medium or low voltage power system with a distinct electric border that incorporates various distributed generators and demands [1]. Microgrid (MG) is an efficient low alternative for providing resilient and environmentally friendly power [2, 3]. In [4], the US Department of Energy defined a MG as "a group of connected distributed resources and loads." There are two main types of microgrids based on their connection to the main substation. The first type is the on-grid MG where the MG is connected to the electrical network. Therefore, the electrical network shares the connected distributed generators (DGs) to supply the required energy to the loads. The second type is the off-grid MG, where the MG is not connected to the electrical network. Therefore, the DGs should provide the required energy to the loads all the time. The energy management of these types of DGs is an important and challenging task for many reasons such as reducing the generation and operation costs, reducing the energy losses, enhancing the voltage profile, reliability, and the system stability. Continuous solar and wind energy is employed to develop the MG. Depending on the energy sources employed, the microgrid may run independently or integrated into an electrical system [5, 6], and [7]. The most important feature of a microgrid is to guarantee consistent electricity based on demand, which is challenging to predict, and energy output can change depending on the availability of renewable energy sources like solar and wind. Recently, it have been seen a fast increase in the expense of fossil fuels, as well as worries about climate change, prompting environmental effect studies and the usage of reliable and sustainable energy sources in electrical network [8].

Renewable energy has become the dominant source of generating in recent years, accounting for a significant portion of global generation. Renewable energy generation is gaining popularity owing to its little impact on the environment. It is utilized due to its low operating costs and extended life duration, although the construction cost is higher to save money on power in the long run. Renewable energy sources are advantageous since they are abundant in nature and environmentally beneficial. By 2030, the rate of increase in power consumption is predicted to be between 40% and 45 percent. RERs are appealing to environmentalists because they are pollution-free, readily available, and continuous in nature. The PV-based energy system is the most attractive RER commercially available due to its low cost, ease of installation, and mobility. It is coordinated with other RERs such as wind farms, fuel cells, diesel generator, electrolyzed, and set to construct a hybrid solar and wind energy system. In addition, it is having reliable electricity and productivity improvements. In most cases, this is paired with a battery storage unit. Thus, a wind and solar system is established when more than two renewable energy sources are utilized [9]. In residential microgrids, fluctuations in RES generation and the transition between grid-on and off-grid modes may raise concerns about system stability [10]. To meet the resources and load balancing, the MG technical and economic restrictions must be addressed. Many studies in the research, such as those in [11, 12], demonstrate that the correct planning of the hybrid MG needs an optimal component rating as per connected load.

Several efforts have been produced for the energy management of MGs. In [13], Bowen Hong et al. presented a useful general framework for designing the scheduling system's energy management core and its constituent parts. In [14], the author proposed a framework for the optimal configuration of a microgrid for minimizing its life cycle cost based on a genetic algorithm to satisfy cooling, heating, refrigeration, and auxiliary demands. The authors in [15] applied a binary integer linear programming for energy management hilly areas microgrid which includes renewable and conventional resources to reduce the cost and emission. A. Rezvani et al. in [16] implemented the lexicographic optimization and hybrid augmented weighted epsilon constraint for optimal and optimal scheduling of MG resources for cost and emission reduction considering market price variation. In [17], the normal boundary intersection (NBI) algorithm has been applied for optimal scheduling of a MG including the PV, FC, microturbine (MT), battery units, and WT for reducing the cost and emission. In this research, an MG is made up of an Energy Storage System (BESS), a PV, and a WT. Greenhouse gas emission reduction rules due to global warming are accelerating the transfer from fossil resources to RES [18].

Several optimization methods have been applied to design, develop, and optimize hybrid MGs. In [19], the overall cost and emissions of a microgrid are reduced by determining the best configuration of DGs while taking voltage stability and energy loss into account as local distribution system restrictions. In [20], an optimization technique is based on Mixed Integer Quadratic Programming (MIQP) for maximizing BESS rating in grid-connected commercial MG to maximize PV capacity and energy storage. In [21], performance the of the Grasshopper Optimization Algorithm (GOA) has been compared with the Cuckoo Search (CS) and Particle Swarm Optimization algorithms to optimize the ratings of various DGs in isolated MG. In [22], the author uses three meta-heuristic techniques to optimize the size and position of DERs in a distribution system, Imperialistic Competitive Algorithm (ICA), PSO, and Genetic Algorithm (GA) algorithms. The PSO is used in [23] to optimize DERs in a distribution system while taking into account different sorts of loads.

This paper presents a modified Capuchin Search Algorithm (MCapSA) for optimal energy management and planning of multi-microgrids (MMGs). The obtained results demonstrate the best performance of the proposed modified SCSO [24], PSO [25], WOA [26], SCA [27], HS [28], ALO [29], and the conventional CapSA [30]. It is important to note that PV and WT generators are incorporated as renewable energy sources, and 24-h planning is suggested with cost minimization as the primary objective as well as summation the VD and the VSI. This work examines how the planning model may be impacted by load, wind speed, and solar radiation uncertainty. As WT and PV are two popular renewable energy sources in modern electrical installations, weather data information like irradiance and temperature is necessary as input variables to simulate the PV. Because of the sensitivity of the PV output, an accurate model is required. As a result, many models are offered to demonstrate the PV's characteristics. Furthermore, the output power of WT is affected by several factors like as turbine size and wind speed.

The following are the most important contribution of this paper:

  1. 1-

    Implementation of a new modified optimization technique for optimal energy management and planning of multi-microgrids (MMGs).

  2. 2-

    The energy management of the MMGs is solved with optimal integration (sizing and siting) of WTs and PV units at uncertain conditions.

  3. 3-

    Proposing an efficient modified Capuchin Search Algorithm (MCapSA) based on multiple strategies to boost the exploitation and the exploration process of the standard CapSA.

  4. 4-

    The superiority of the proposed algorithms is validated via the statistical results and comparisons with other algorithms on the standard benchmark functions, 33-bus and 69-bus multi-microgrids.

The rest of this paper is listed as follows: The energy management formulation is described in Sect. 2. Section 3 lists the uncertainty representation. The solution technique for solving energy management is provided in Sect. 4. The obtained results are covered in Sect. 5. Conclusions of this work are presented in Sect. 6.

1.1 Problem formulation

This section describes the problem formulation for the reduction in operating costs in the microgrid. We explain the total annual cost function, voltage stability, and constraints of the optimization operation management issue.

1.2 The objective function

In this work, three objective functions have been optimized which can be listed as follows:

1.2.1 The cost reduction

The overall cost reduction (\(C\)) is the first objective function that is taken into consideration, and it 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:

$$C=\mathit{min}({C}_{Grid}+{C}_{PV} +{C}_{Wind}+ {C}_{Loss})$$
(1)

where the following are the detailed elements of Eq. (1),

$${C}_{Grid}= 365\times {U}_{Grid}\times \sum_{h=1}^{24}{P}_{Grid(h)}$$
(2)

where \({U}_{Grid}\) is the purchasing electricity cost from the grid

$${C}_{Loss}=365\times {U}_{Loss}\times \sum_{h=1}^{24}{P}_{Loss(h)}$$
(3)

where \({U}_{Loss}\) is the cost of energy loss. \({P}_{Loss}\) are the active power loss.

$${C}_{PV}= {C}_{PV}^{inst.}+{C}_{PV}^{O\&M}$$
(4)

where \({C}_{PV}^{inst.}\) is the cost of installing the PV, \({C}_{PV}^{O\&M}\) is the PV unit's operating and maintenance costs,

$${C}_{PV}^{O\&M}={U}_{PV}^{O\&M}\times \sum_{h=1}^{24}{P}_{PV(h)}$$
(5)
$${C}_{PV}^{inst.}={CF\times {U}_{PV}\times P}_{\mathrm{rated}\_PV}$$
(6)

where \(CF\) is a factor affecting capital recovery, \({P}_{rated\_PV}\): is the rated produced power of the PV

$${C}_{Wind}= {C}_{wind}^{inst.}+{C}_{wind}^{O\&M}$$
(7)

where \({C}_{Wind}^{inst.}\) Is the cost of installing the wind turbine, and \({C}_{wind}^{O\&M}\) is the wind’s operating and maintenance costs.

$${C}_{Wind}^{O\&M}={U}_{Wind}^{O\&M}\times \sum_{h=1}^{24}{P}_{Wind(h)}$$
(8)

where \({U}_{PV}^{O\&M}\) and \({U}_{Wind}^{O\&M}\) represent the PV and Wind unit's operating and maintenance purchase cost ($/kW), respectively.

$${C}_{Wind}^{inst.}={CF\times {U}_{WT}\times P}_{rated\_wind}$$
(9)

where \({P}_{rated\_wind}\) is the rated produced power of wind turbine

$$CF=\frac{\beta \times {(1+\beta )}^{NP }}{{(1+\beta )}^{NP}-1}$$
(10)

where \(NP\) is the Wind turbine or PV unit lifetime

In addition, the generated power from the PV may be computed using the formula below [31]:

$$P_{{PV}} = \left\{ {\begin{array}{*{20}c} {P_{{rated\_PV}} \left( {\frac{{S_{i}^{2} }}{{S_{{std}} \times X_{s} }}} \right){\text{for}}\;0} & {if;S_{i} \le X_{s} } \\ {P_{{rated\_PV}} \left( {\frac{{S_{i} }}{{S_{{std}} }}} \right)} & {{\text{for}}\;X_{s} \le S_{i} \le S_{{std}} } \\ {P_{{rated\_PV}} } & {S_{{std}} \le S_{i} } \\ \end{array} } \right.$$
(11)

where (\({S}_{i}\)) stands for solar irradiance. The irradiation of the standard environment (\({S}_{std}\)) is 1000 W/m2 and (\({X}_{s}\)) for a specific irradiance point, which is equal to 120.

Also, the WT's produced output power may be computed as shown below:

$$P_{{WT}} \left( W \right) = \left\{ {\begin{array}{*{20}l} 0 \hfill & {for} \hfill & {W < W_{i} and > W} \hfill \\ {P_{{rated\_wind}} \left( {\frac{{W - W_{i} }}{{W_{r} - W_{i} }}} \right)} \hfill & {for} \hfill & {(W_{i} \le W \le W_{R} )} \hfill \\ {P_{{rated\_wind}} (} \hfill & {for} \hfill & {W_{r} < W \le W_{0} } \hfill \\ \end{array} } \right.$$
(12)

where \({W}_{r}\) is a rated wind speed; \({W}_{o}\) is a cut-out speed. \({W}_{i}\) indicates a cut-in speed.

1.2.2 Enhancement of voltage level

To improve the network performance, the voltage deviations should be kept to a minimum that can be defined as:

$$\sum VD=\sum_{h=1}^{24}\sum_{k=1}^{NB}\left|\left({V}_{k}-1\right)\right|$$
(13)

where \({V}_{k}\) denotes the voltage of the kth bus and \(\mathrm{NB}\) denotes the buses in the network

1.2.3 Improved system stability

$${VSI}_{k}={\left|{V}_{k}\right|}^{4}-4{\left({P}_{k}{X}_{km}-{Q}_{k}{R}_{km}\right)}^{2}-4\left({P}_{k}{X}_{km}+{Q}_{k}{R}_{km}\right){\left|{V}_{k}\right|}^{2}$$
(14)
$$\sum VSI=\sum_{h=1}^{24}\sum_{k=1}^{NB}{VSI}_{k}$$
(15)

where \({VSI}_{k}\) is the voltage stability index, \({R}_{km}\) is the resistance of the TL while the \({X}_{km}\) is its reactance. \({P}_{k}\) and \({Q}_{k}\) define the real and reactive powers, respectively. The following three objective functions are taken into consideration simultaneously:

$$f= {h}_{1}{f}_{1}+ {{h}_{2}f}_{2}+{h}_{3}{f}_{3}$$
(16)
$${f}_{1}=\frac{{C}_{RERs}}{{C}_{Base}}$$
(17)
$${f}_{2}=\frac{{\sum VD}_{RERs}}{{\sum VD}_{Base}}$$
(18)
$$f_{3} = \frac{1}{{VSI_{RERs} }}$$
(19)

when \({C}_{RERs}\), \({VD}_{RERs}\) and \({VSI}_{RERs}\) are the total cost, summation the voltage deviations and the voltage stability index with inclusion of the RERs while \({C}_{Base}\) and \({VD}_{Base}\) are the total cost and summation the voltage deviations at base case. The \({h}_{1}\), \({h}_{2}\), and \({h}_{3}\) are the penalty factors that are selected to be 0.5, 0.25, and 0.25, respectively, [32].

1.3 The constraints of the network

$${V}_{Min}\le {V}_{i}\le {V}_{Max}$$
(20)
$${P}_{rated\_ PV}+{P}_{rated\_wind}\le \sum_{i=1}^{NB}{P}_{Load,i}$$
(21)
$${PF}_{Min}\le PF\le {PF}_{Max}$$
(22)
$${I}_{n}\le {I}_{max,n} n =\mathrm{ 1,2},3\dots ,NT$$
(23)

where, \({V}_{min}\), and \({V}_{max}\) are the lower and maximum voltage limits. \({P}_{Load}\) and \({\mathrm{Q}}_{\mathrm{Load}}\) signified the real and reactive load, respectively; \({I}_{max,n}\) is the maximum allowable current limit of the n-th line; \({PF}_{min}\) and \({PF}_{max}\) are the minimum and the upper and the power factor limits of the WT, respectively; \(NT\) defines the number of the TL.

1.3.1 Restrictions on inequality

$${P}_{Grid}+{P}_{PV}+{P}_{WT}= \sum_{i=1}^{NT}{P}_{Loss,i}+\sum_{i=1}^{NB}{P}_{Load,i}$$
(24)
$${Q}_{Grid}+{Q}_{WT}= \sum_{i=1}^{NT}{Q}_{Loss,i}+\sum_{i=1}^{NB}{Q}_{Load,i}$$
(25)

where,\({P}_{Grid}\) and \({Q}_{Grid}\) refer to the active power and reactive power of the mean upstream station, respectively.

2 Modeling the uncertainties

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

2.1 The probabilistic model of solar irradiance

The Beta pdf has been employed to model the intermittent nature of the solar irradiance, as follows [37, 38]:

$${f}_{b}({g}_{s})=\left\{\begin{array}{ll}\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.$$
(26)

where the Beta pdf and the gamma function are denoted, respectively, by \(f_{b} \left( {g_{s} } \right){\text{ and }}\); \(\mu\) is the mean value of the irradiance while \(\sigma\) is its standard deviation of historical data which can be employed to calculate them as follows [39, 40]:

$$\beta =\left(1-\mu \right)\times \left(\frac{\mu \times \left(1+\mu \right)}{{\sigma }^{2}}-1\right)$$
(27)
$$\alpha =\frac{\mu \times \beta }{1-\mu }$$
(28)

The continuous pdfs are divided into several divisions. The probability of divisions occurrence can be given by:

$${prob}_{i}^{{g}_{s}}=\underset{{g}_{s,i}}{\overset{{g}_{s,i+1}}{\int }}{f}_{b}\left({g}_{s}\right)d{g}_{s,i}$$
(29)

where \({g}_{s,i}\) and \({g}_{s,i+1}\) are the limits of the interval i while \({prob}_{i}^{{g}_{s}}\) is its probability.

2.2 The probabilistic wind speed representation

For modeling the wind speed uncertainty, a Weibull pdf with a shape parameter of 2 is normally used [35, 41]. In this case, it is called Rayleigh pdf [38] which can described as follows:

$${f}_{r}(v)=\left(\frac{2v}{{c}^{2}}\right)\mathit{exp}\left[-{\left(\frac{v}{c}\right)}^{2}\right]$$
(30)

where v represents the wind speed while c denotes the scale parameters of Rayleigh pdf. By dividing the continuous pdf into several intervals, the probability of each interval is given by:

$${prob}_{i}^{v}={\int }_{{v}_{w1}}^{{v}_{w2}}{f}_{r}(v).d{v}_{i}$$
(31)

\({v}_{w1}\) and \({v}_{w2}\) are the limits of the wind speed of the interval i. \({prob}_{i}^{v}\) denotes probability of state i.

2.3 The probabilistic representation of load demand

Normal pdf is employed for representing the demanding uncertainty which is given using (32) [36].

$${f}_{n}\left(l\right)=\frac{1}{{\sigma }_{l}\sqrt{2\pi }}\times exp\left[-\left(\frac{l-{\mu }_{l}}{2{\sigma }_{l}^{2}}\right)\right]$$
(32)

where \({\mu }_{l}\) refers the mean load value while \({\sigma }_{l}\) is the standard deviation, respectively. The probability of a segment can be given by:

$${prob}_{i}^{l}=\underset{{l}_{i}}{\overset{{l}_{i+1}}{\int }}{f}_{n}\left(l\right)dl$$
(33)

where the \({l}_{i}\) and \({l}_{i+1}\) are the limits of the state i.\({prob}_{i}^{l}\) refers to the probability of interval i.

2.4 The combined probabilistic model

The combined probabilistic model (\({P}_{com,i}\)) is generated by combining the uncertainty models of irradiance, loading, and wind speed. For interval i, the combined model can be mathematically as below equation.

$${P}_{com,i}={{prob}_{i}^{v}\times prob}_{i}^{{g}_{s}}\times {prob}_{i}^{l}$$
(34)

3 Proposed technique

3.1 Capuchin search algorithm

CapSA is a new optimization technique that mimics the movement of the Capuchins in their foraging technique where the Capuchins climb, leap, and swing over the branches of the trees to obtain their food [30]. As the main aspect of the foraging motion, the Capuchins follow their leader and interact together in the same swarm which consists of females, males, and small apes. The motion of the Capuchins in the forest as a group to find the promising regions of food simulates the global searching phase of the optimization algorithms while in case the Capuchins follow their leaders, it mimics the exploitation phase of the algorithms, and they collect the other individuals in the group by barking and postures. The CapSA simulates three motions of the Capuchins which can be presented as follows:

3.1.1 Leaping motion

It is well known that the Capuchins jump between the trees at a large-scale distance like the projectiles which can be represented using the third law of motion according to (35).

$$x={x}_{0}+{v}_{0}t+\frac{1}{2}a{t}^{2}$$
(35)

Here, \(x\) reprints the new placement of the Capuchin while \({x}_{0}\) refers to its initial placement.\({v}_{0}\) denotes the initial speed. Equation (35) can be modified as follows:

$$x={x}_{0}+{v}_{0}^{2}\mathrm{sin}\left(2{\theta }_{0}\right)/g$$
(36)

where \(g\) represents the gravitational acceleration and \({\theta }_{0}\) refers to the leaping angle.

3.2 Swinging motion

The swing motion of the Capuchins over the branches of the trees simulates the local searching area of algorithms, and it can be modeled as follows:

$$x=L\mathrm{sin}\theta$$
(37)

where \(L\), \(\theta\) are the tail of length and the swing angle.

3.3 Climbing motion

The Capuchins climb to the trees to obtain the food; this motion is represented as follows:

$$x={x}_{0}+{v}_{0}t+\frac{1}{2}\left(v-{v}_{0}\right){t}^{2}$$
(38)

The initialization location of the Capuchins can be obtained as follows:

$${X}^{i}={X}^{\mathrm{Lb}}+\left({X}^{\mathrm{Ub}}-{X}^{\mathrm{Lb}}\right)\times rand$$
(39)

where \({X}^{Lb}\) denotes the upper boundary of the control variables and \({X}^{\mathrm{Ub}}\) denotes its upper boundary and \(0 \le rand\le 1\). The Capuchin swarm includes the leader or the Alpha leader and the follower where both of them look for the food resource location (F). The motion of the leader is represented using (40).

$$\begin{array}{c}{X}^{i}=F+\frac{{P}_{bf}{\left({v}^{i}\right)}^{2}\mathrm{sin}(2\theta )}{g}\\ i<n/2;0.1<\epsilon \le 0.20\end{array}$$
(40)

where \({P}_{bf}\) denotes the probability of the balance. \(F\) represents the position of the food. \(\epsilon\) denotes a random number. \(\theta\) denotes the Capuchin’s jumping angle, and it can be denoted as follows:

$$\theta =\frac{3}{2}r$$
(41)

where \(r\) refers to random value in range [0,1]. For balancing the transposition between the exploration and the exploitation of the algorithm and vector \(\tau\) is used for the pattern which can be represented as follows:

$$\tau ={\beta }_{0}{e}^{-{\beta }_{1}{\left(\frac{t}{T}\right)}^{{\beta }_{2}}}$$
(42)

where \({\beta }_{0}\), \({\beta }_{1}\), and \({\beta }_{2}\) are values number are set to be 2, 21, and 2, respectively, while \(T\) and \(t\) denote the maximum and current iterations, respectively. The velocity of the Capuchin is calculated using (43)

$${v}^{i}=\rho {X}^{i}+\tau {a}_{1}\left({X}_{best}^{i}-{X}^{i}\right){r}_{1}+\tau {a}_{2}\left(F-{X}^{i}\right){r}_{2}$$
(43)

where \({X}_{best}^{i}\) is the best position. \({a}_{1}\) and \({a}_{2}\) are two positive numbers. that control the effects of \(F\) and \({X}_{best}^{i}\) on the velocity. \({r}_{1}\) and \({r}_{2}\) are random umber with in [0,1]. \(\rho\) is an inertia coefficient equals to 0.7. The leader (Alpha) Capuchin moves in the ground in a large area to search the food location in case of the lake the food over the trees. The Capuchin motion on the ground is represented as follows:

$$\begin{array}{*{20}c} {X^{i} = F + \frac{{P_{ef} P_{bf} \left( {v^{i} } \right)^{2} {\text{sin}}\left( {2\theta } \right)}}{g}} \\ {i < n/2;0.2 < \in \le 0.30} \\ \end{array}$$
(44)

where \(P_{ef}\) denotes the probability motion of the Capuchins on the ground. The movement of the Alpha on the group is represented as follows:

$$\begin{array}{*{20}c} {X^{i} = X^{i} + v^{i} } \\ {i < \frac{n}{2} ; 0.3 < \le 0.50} \\ \end{array}$$
(45)

The swinging motion of the Capuchins on the branches of the trees to get the food in the local area is represented as follows:

$$\begin{array}{*{20}c} {X^{i} = F + \tau P_{bf} \times sin\left( {2\theta } \right)} \\ {i < n/2;0.5 < \le 0.75} \\ \end{array}$$
(46)

The Alpha and some Capuchins descend several times to search the food; this motion can be expressed as follows:

$$\begin{array}{*{20}c} {X_{j}^{i} = F_{j} + \tau P_{bf} } \\ {\left( {v_{j}^{i} - v_{j - 1}^{i} } \right)} \\ {i < n/2;0.75 < \le 1.0} \\ \end{array}$$
(47)

Capuchins may move randomly to find new regions of food resources. This motion is represented mathematically as follows:

$$\begin{array}{c}{X}^{i}=\tau \times \left[{X}^{\mathrm{Lb}}+\epsilon \times \left({X}^{\mathrm{Ub}}-{X}^{\mathrm{Lb}}\right)\right]\\ i<n/2;\epsilon \le Pr\end{array}$$
(48)

where \(Pr\) is a number selected to be 0.1. The followers flow the Alpha Capuchins; this motion is represented mathematically as follows:

$${X}_{j}^{i}=\frac{1}{2}\left({XL}_{j}^{i}+{X}_{j}^{i-1}\right)n/2\le i\le n$$
(49)

where \({XL}_{j}^{i}\) and \({X}_{j}^{i-1}\) are the location of the leader and the previous location of the followers, respectively.

3.4 The modified Capuchin search algorithm (MCapSA)

3.4.1 The quasi-oppositional-based learning (QOBL)

The oppositional-based learning (OBL) was applied for performance enhancement of several optimization techniques and boost the convergence rapidity. The opposition position has a superior probability to be nearer of the best solution than a random solution [39]. The opposite number means the mirror of this number which can be described as follows:

$${X}^{o}={x}^{min}+{x}^{max}-x$$
(50)

where \({x}^{min}\) is the lower margin of \({x}^{o}\) and \({x}^{max}\) is its upper limit. The opposite point of a multi-dimension point is represented as follows:

$${X}_{i,j}^{o}={x}_{j}^{min}+{x}_{j}^{max}-{X}_{i,j}, i=\mathrm{1,2},\dots ,n$$
(51)

The quasi-based learning is utilized to improve the oppositional-based method and represents the center point of the search area as follows:

$${C}_{i,j}=\left({x}_{j}^{min}+{x}_{j}^{max}\right)/2$$
(52)

The quasi-opposite point means the center and the opposite value. The QOBL combines quasi- and oppositional-based learning to improve which can be o find the global solution and represented as follows:

figure a

3.4.2 Levy flight random walk

The second modification is based on the Levy flight distribution (LFD) for allowing the populations to jump to new regions far away from the best so far solution. The primary feature of the LFD walk is the movement of the population in random and stochastic movement. Therefore, the application of the LFD to the motion of Capuchins can enhance the performance of conventional CapSA and avoid its stagnation. LFD can be described as follows:

$${X}_{i}(t+1)={X}_{i}(t)+\alpha \oplus \text{ L}{{\acute{\mathrm{e}}}}\text{vy }(\lambda )$$
(55)

where \(\alpha\) is to the step size. \(\text{L}{{\acute{\mathrm{e}}}}\text{vy }(\lambda )\) represents the levy distribution which can be described as follows [40]:

$${\text{Levy}} \, (\lambda )\sim \frac{\mu }{|v{|}^{1/\beta }}$$
(56)

in which

$$\mu \sim N\left(0,{\sigma }_{\mu }^{2}\right), v\sim N\left(0,{\sigma }_{\mu }^{2}\right)$$
(57)
$${\sigma }_{\mu }={\left\{\frac{\Gamma (1+\beta )\mathrm{sin}(\pi \beta /2)}{\Gamma [(1+\beta )/2]\beta {2}^{(\beta -1)/2}}\right\}}^{1/\beta },{\sigma }_{v}=1$$
(58)

here \(\Gamma\) refers to the Gamma function. \(\beta\) denotes a random value in the range [0,2]. The value of \(\beta\) is selected to be 1.5 [41].

3.4.3 Exploitation boosting based on DOA

The third modification is conceptualized from the exploitation phase of (DOA) where the prairie dogs have a special communication system [42]. Prairie dogs do a special sound in case of food availability or predator threats. The prairie dogs respond to the sound and go directly to the food source’s location while in case of the presence of predators or a hawk they make another warring sound and go into hiding. Application of the exploitation mechanism of the prairie dogs into the RBO is formulated using (60).

$$X_{{i,j}} = \left\{ {\begin{array}{*{20}c} {{\text{GBest}}_{{i,j}} - eC{\text{Best }}_{{i,j}} \times \varepsilon - CPD_{{i,j}} } & {if\;rand \ge 0.5} \\ {{\text{GBest}}_{{i,j}} \times PE} & {\;if\;rand < 0.5} \\ \end{array} } \right.$$
(59)

where \({\text{GBest}}_{i,j}\) is the best solution. \(\varepsilon\) is a small random value.

$$P{D}_{i+1,j+1}={ \, {\text{GBest}} \, }_{i,j}\times PE$$

where \({\text{GBest }}_{i,j}\) denotes the global best solution. \(eC{\text{Best }}_{i,j}\) denotes the effect of the current obtained solution and it can be calculated using (61) while \(CP{D}_{i,j}\) denotes the cumulative effect of all prairie dogs and it can be calculated using (62) and \(PE\) is the adaptive operator which can be calculated using (63). \(\Delta\) is a small number representing a small number.

$${\text{eCBest}}_{i,j}={ \, {\text{GBest}}}_{i,j}\times\Delta +\frac{{X}_{i,j}\times \mathrm{mean}\left({X}_{n,m}\right)}{{\text{ GBest }}_{i,j}\times \left(U{B}_{j}-L{B}_{j}\right)+\Delta }$$
(60)
$${\text{CPD }}_{i,j}=\frac{{\text{ GBest }}_{i,j}-rP{D}_{i,j}}{{\text{ GBest }}_{i,j}+\Delta }$$
(61)
$$PE = 1.5 \times \left( {1 - \frac{t}{{T_{{\max }} }}} \right)^{{\left( {2\frac{{{\text{ }}t{\text{ }}}}{{T_{{\max }} }}} \right)}}$$
(62)

The steps of applying the proposed algorithm for optimal integration of the RERs are depicted in Fig. 1.

Fig. 1
figure 1

Flow chart of the MCapSA for optimal energy management of MMGs

Full size image

4 The simulation results

The proposed MCapSA is applied to solving the energy management of MMGs of two standard networks IEEE 33-bus and 69-bus MMGs. Initially, the performance of the proposed method is tested on the standard benchmark functions and the simulation results were compared to other optimization techniques including SCSO [24], GWO [43], WOA [26], SCA [27], PF [44], DO [45], and the standard CapSA [30]. The program code of the proposed MCapSA and the other techniques for energy management was written in MATLAB software (MATLAB R2019b) and carried out on a PC with Intel i5, 2.5 GHz CPU, and 4 GB RAM. The studied cases are presented as follows:

4.1 Validation of the MCapSA technique on standard benchmark functions

Here, the proposed MCapSA is tested on 23 benchmark functions including the unimodal, multimodal, and fixed-dimension multimodal benchmark functions, which are listed in Appendix A [46,47,48]. For all cases, the parameters are set according to Table 1 and the obtained results are represented over 30 run times.

Table 1 Selected parameters of the optimizers
Full size table

4.1.1 Statistical results analysis

This section depicts the performance of the MCapSA compared with SCSO [24], PSO [25], WOA [26], SCA [27], HS [28], ALO [29], and the conventional CapSA [30]. Table 2 shows the statistical results including the mean, the worst, the average, the best, the standard deviation (SD) values, and the Wilcoxon P-value between the MCapSA and the other optimizers. According to Table 2, the proposed MCapSA optimizer is superior in terms of the mean, the worst, and the best values for the F1-F5, and F17-F23 while the obtained results for the reported algorithms are similar for F17 to F23. In addition to that for F12, the CapSA is better than MCapSA. The Wilcoxon P-values between the MCapSA and the other optimizers are presented in the 7th column of Table 2. When the P-value is less than 5%, it is a significant difference between optimizers [49]. On the opposite of that when it is more than 5%, there is 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 MCapSA has significant differences compared with SCSO, PSO, WOA, SCA, HS, ALO, and the standard CapSA for most of the studied benchmark functions.

Table 2 Statistical results of different optimizers for the standard functions
Full size table

4.1.2 The convergence curves analysis

The convergence curves of the MCapSA and other reported techniques including SCSO [24], PSO [25], WOA [26], SCA [27], HS [28], ALO [29], and the conventional CapSA [30] are shown in Fig. 2. According to the convergence curves, the proposed MCapSA has the best convergence speed for the fixed dimension, the unimodal, and the multimodal functions. In addition to that, the proposed MCapSA is superior and converged to the optimal solution faster than the conventional CapSA due to the proposed modifications which can boost the exploration and exploitation phases of the proposed optimizer.

Fig. 2
figure 2

Convergence curves of the test benchmark functions by different optimizers

Full size image

4.1.3 The boxplot analysis

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

Fig. 3
figure 3

Boxplot for test benchmark functions by different optimizers

Full size image

4.2 Solving the energy management of MMGs using the MCapSA

The proposed MCapSA is implemented for solving the energy management of MMGs with and without RERs. The proposed MCapSA method was tested on IEEE 33-bus and 69-bus MMGs. The system data of system are given in [50] while the specification of these networks is depicted in Table 3. It worth mentioning here that Table 3 lists the specifications of the studied systems (33-bus, 69-bus) at 100% loading. The studied distribution networks (DNs) including IEEE 33-bus and 69-bus grids have been divided into three microgrids as shown in Figs. 4 and 5, respectively. Each microgrid includes its RERs (WT and PV units). Therefore, the network has three WTs and three PV. The obtained results by MCapSA have been compared with those obtained by the classical CapSA, SCA, ALO, and PSO. For a fair comparison, the selected maximum number of maximum iterations and the populations are selected to be 80 and 18, respectively. The uncertainty of load, the wind speed, and the solar irradiance are considered as depicted in Sect. 3. The expected day-ahead load demand is shown in Fig. 6 while the market price of the purchasing energy is illustrated in Fig. 7 [51]. The expected wind speed and the irradiance are displayed in Figs. 8 and 9, respectively [52]. The cost coefficients for RERs, the purchasing power from the grid, and the constraints are presented in Table 4. The energy management studied cases are listed as follows:

Table 3 Specification of the 69-bus and 33-bus DNs
Full size table
Fig. 4
figure 4

Topologies of the first MMGs

Full size image
Fig. 5
figure 5

Topologies of the second MMGs

Full size image
Fig. 6
figure 6

Expected load profile

Full size image
Fig. 7
figure 7

Hourly market price

Full size image
Fig. 8
figure 8

Expected wind speeds

Full size image
Fig. 9
figure 9

Expected solar irradiance

Full size image
Table 4 Constraints and the cost coefficients
Full size table

4.2.1 The energy management of 33-bus MMGs system.

The key objective of the energy management of the MMGs is to reduce the total annual cost along with enhancing the system's performance. Initially, without the inclusion of RERs (base case), the total annual purchased from the grid is 3.0539E + 06 kWh, the energy losses are 1.4759E + 06 kW, the total purchase energy cost is 7.5653E + 06 USD, the energy loss cost is 8.855e + 04 USD, the summation of total annual cost is 7.6539E + 06 USD, the voltage deviations summation is 38.376 p.u., and the VSI is 6.291E + 02 p.u. Table 5 shows the simulation results with and without RERs which have been obtained by the MCapSA and other well-known algorithms. Table 5 summarizes the best obtained results including the total costs, summations of the VD, and the VSI. From Table 5, by application of the proposed optimization algorithm, the total annual cost is reduced from 7.65393E + 06 USD to 4.1019E + 06 USD. The summation of the voltage deviations is reduced from 38.376 p.u. to 14.374 p.u. Therefore, the voltage profile is enhanced considerably as illustrated in Fig. 10. Furthermore, the summation of the VSI has been enhanced from 629.02 p.u. to 724.323 p.u. The power losses of the system with and without the inclusion of the hybrid RERs are depicted in Fig. 11. It is obvious that the inclusion of the hybrid RERs in each microgrid optimal decreases the power losses considerably. The optimal locations of the hybrid units in the first, second, and third MGs are 6, 13, and 32, respectively. The rating of the WTs in kW of the first, second, and third MGs is 82, 1021, and 920, respectively, while the rating of the PVs is 151, 1481, and 54. It should be highlighted here that the total annual cost by application is better than the conventional CapSA. The statically results for the considered objective function that have been obtained by the application of different optimizers are recorded in Table 6. As depicted in Table 6, the obtained results by the proposed MCapSA are better than other optimizers in terms of the best, the mean, the worst, and the standard deviation (SD) values. The output powers of the PV units are indicated in Fig. 12. From Fig. 12, the yield powers of these units are varied continuously with the irradiance variations. Likewise, the generated powers of the WTs varied with the variation in the wind speed as shown in Fig. 13. Finally, the convergence curves by application of different optimization algorithms are shown in Fig. 14. Referring to this figure, the proposed algorithm has the best convergence characteristics compared with other optimizers.

Table 5 Energy management results of the first MMGs
Full size table
Fig. 10
figure 10

Voltage profile units for MMGs 1 a without hybrid system. b with hybrid system

Full size image
Fig. 11
figure 11

Power losses for the first MMG

Full size image
Table 6 Simulation results of objective function for applied algorithms
Full size table
Fig. 12
figure 12

Generated power of PV units of the first MMG

Full size image
Fig. 13
figure 13

WT units’ generated powers for first MMG

Full size image
Fig. 14
figure 14

Convergence curves of different optimizer for case 1

Full size image

4.2.2 The energy management of 69-bus MMGs system

The optimal sizing and ratings of the PV units and WTs have been assigned optimally for the three microgrids. In the base case, the total annual purchased from the grid is 3.1311E + 07 kWh, the energy losses 1.5716E + 06 kW, the total purchase energy cost is 7.7576E + 06 USD, the energy loss cost is 9.4297E + 04 USD, the summation of total annual cost is 7.8519 E + 06 USD, the voltage deviation is 39.0912 p.u. while summation of the VSI is 1.4866E + 03 p.u. Table 7 summarizes the simulation results with and without RERs which have been obtained by the MCapSA and the standard CapSA including the total costs, summations of the VD, and the VSI. According to Table 7, with implementations of the proposed algorithm, the total annual cost is reduced from 7.8519E + 06 USD to 4.3819E + 06 USD. The summation of the voltage deviations is reduced from 39.0912 p.u. to 23.5733 p.u. Therefore, the voltage profile is enhanced considerably as shown in Fig. 15. Furthermore, the VSI during the day ahead has been enhanced from 1.4866E + 03 p.u. to 1.5533E + 03 p.u. The system power losses with and without the insertion of the hybrid RERs are illustrated in Fig. 16. It is obvious that with the inclusion the hybrid RERs in each microgrids the power losses decrease considerably. The optimal sites of the hybrid units in the first, second, and third MGs are 12, 22, and 25, respectively. The WTs rating in kW of the first, second, and third MGs are 1119, 488, and 25, respectively, while the optimal sizes of the PV units are 1198, 1119, and 55. The generated power of the PV units WTs is shown in Figs. 17 and 18, respectively. It should be highlighted here that the total annual cost by application is better than the conventional CapSA. The statistical results that have been obtained by different optimizers for the considered objective function are tabulated in Table 8. It is clear that the proposed MCapSA is superior compared with other optimizers in terms of the best, the mean, the worst, and the standard deviation (SD) values. Finally, the convergence of the objective function by the application of different algorithms is shown in Fig. 19. Referring to this figure, the proposed algorithm MCapSA has the best convergence characteristics compared with other optimizers.

Table 7 Energy management results of the second MMG
Full size table
Fig. 15
figure 15

Voltage profile of units for the second MMG. a Without hybrid system. b With a hybrid system

Full size image
Fig. 16
figure 16

Power losses for the second MMG

Full size image
Fig. 17
figure 17

PV units’ generated powers

Full size image
Fig. 18
figure 18

WT units’ generated powers for the second MMG

Full size image
Table 8 Statical results for the second MMGs
Full size table
Fig. 19
figure 19

Objective function for ALO, CapSA, PSO, SCA, and McapSA for 69-bus system

Full size image

5 Outcomes and unique features

The summarized unique features can be explained as follows:

Proposing a novel Modified Capuchin Search Algorithm (MCapSA) for addressing the stagnation issue of the traditional Capuchin Search Algorithms (CapSA). The modification is based on: quasi-oppositional-based learning (QOBL), Levy flight random walk, and exploitation phase of DOA. The following outcomes can be listed as follows:

  1. 1.

    The purposed algorithm was applied for 33-bus and 69-bus multi-microgrid systems.

  2. 2.

    The suggested MCapSA has been evaluated on 23 different benchmark functions, including unimodal, multimodal, and fixed-dimension multimodal benchmark functions.

  3. 3.

    Obviously, the proposed algorithm is superior and converged to the optimal solution faster than the conventional, where the boxplots of the MCapSA are narrower and have significant difference compared with the other optimizers.

  4. 4.

    The proposed MCapSA algorithm was able to solve the energy management problem compared to the basic algorithm (CapSA) and other well-known algorithms.

  5. 5.

    The obtained results due to the objective functions compared to the basic case can be represented as follows:

    • For the IEEE 33-bus MMG

    • The total cost reduced from 7.65393E + 06 USD to 4.1019E + 06 USD which is the best compared with other techniques as illustrated in Table 5

    • The voltage deviation (VD) is reduced from 38.37576 p.u. to 14.3740 p.u.

    • The VSI is enhanced from 629.0769 p.u. to 724.3230 p.u.

    • For the IEEE 69-bus MMG

    • The total cost reduced from 7.851900E + 06 USD to 4.3819E + 06 USD.

    • The VD is reduced from 39.0912 p.u. to 23.5733p.u.

    • The VSI is enhanced from 1.4866E + 03 p.u. to 1.5533E + 03 p.u.

  6. 6.

    The optimal locations and sizes of the RERs in power system are one of the main features presented in this paper.

  7. 7.

    An optimal integration of the RERs can reduce the energy cost, enhance the voltage profile and the system performance.

  8. 8.

    A comparison due the voltage profile enhancement and the power loss has been applied between the technique and other well-known techniques to explain the superiority of purposed algorithm.

  9. 9.

    The total cost that has been obtained by application the proposed algorithm for the IEEE 33-bus MMG is 4.1019E + 06 USD while the total annual costs that have been obtained by application of CapSA, SCSO, PSO, WOA, SCA, HS, and ALO are 4.14E + 06 USD, 4.17E + 06 USD, 4.14E + 06 USD, 4.56E + 06 USD, 4.29E + 06 USD, 4.22E + 06 USD, and 4.16E + 06 USD, respectively. According to the previous results, the minimum cost can be obtained with the application of the proposed algorithm.

  10. 10.

    The total annual cost that has been obtained by application of the proposed algorithm for the IEEE 69-bus MMG is 4.1019E + 06 USD while the total annual costs that have been obtained by application of CapSA, SCSO, PSO, WOA, SCA, HS and ALO are 4.40E + 06 USD, 4.75E + 06 USD, 4.77E + 06 USD, 4.65E + 06 USD, 4.54E + 06 USD, 4.63E + 06 USD, and 4.69E + 06 USD, respectively. According to the previous results, the minimum cost can be obtained with the application of the proposed algorithm.

6 Conclusions

In this paper, the energy management issue for multi-microgrids (MMGs) was solved with optimal incorporation of WT and solar PV. A Modified Capuchin Search Algorithm (MCapSA) was presented for the energy management solution based on three mechanisms including QOBL, Levy flight motion, and the exploitation phase of the prairie dogs in the prairie dog optimization (PDO). The energy management was solved for cost reduction, and system performance enhancement at uncertainty of the solar irradiance, wind speed, and market price. The proposed algorithm was tested on standard benchmark functions and two distribution MMGs including the IEEE 33-bus and IEEE 69-bus MMGs, and the obtained results were compared with other well-known optimizers. The obtained result demonstrated that by the optimal energy management with inclusion the solar PV and WT using the proposed method, the total cost reduced from 7.65393E + 06 USD to 4.1019E + 06 USD or by 46.40%, the VD is reduced from 38.37576 p.u. to 14.3740 p.u. or by 62.54%, and the VSI is enhanced from 629.0769 p.u. to 724.3230 p.u. or by 15.1406% compared to the based case for the first MMG. Likewise, for the second MMG the total cost reduced from 7.851900E + 06 USD to 4.3819E + 06 USD or by 44.19%, the total VD is reduced from 39.0912 p.u. to 23.5733p.u. or by 39.70%, and the total VSI is enhanced from 1.4866E + 03 p.u. to 1.5533E + 03 p.u. or by 4.49%. In addition to that, the proposed MCapSA is superior compared to SCSO, PSO, WOA, SCA, HS, ALO, and the standard CapSA.