In this equation, *t* represents the time, *Pch* and *Pdch* denote the charge and discharge rates, and *SoC* stands for the state of charge of the storage unit that has a maximum storage capacity limit and a minimum storage capacity boundary owing to the long-run operation and maintenance facets. Each storage unit can be charged or discharged by a limited amount within an hour. The last equation also ensures that the storage unit cannot operate at charge and discharging modes simultaneously. In this equation, *I*(*t*) is a binary variable.

#### 2.3.8. The Uncertainty Modeling for PV Sources

The nature of solar power and photovoltaic irradiance can be represented by various types of distribution functions. One of the most efficient types for this purpose is the lognormal probability distribution function (PDF), which effectively and precisely characterizes the intensity of irradiance on a typical day. The PDF of solar irradiance (ζs) trailing the lognormal distribution function with a mean value of *μs* and standard deviation of σs that is presented in Equation (20).

$$f\left(\zeta\_s \left| \mu\_s, \sigma\_s^2 \right.\right) = \frac{1}{\sigma\_s^2 \sqrt{2\pi}} \exp\left\{-\left(\frac{\left(\log\left(\zeta\_s\right) - \mu\_s\right)^2}{2\sigma\_s^2}\right)\right\} \qquad \forall \zeta\_s > 0 \tag{20}$$

The solar power generation (*PPV*) depends on the solar radiation intensity, also known as solar irradiance, which is shown by *ζ*s. The generation curve can be different based on different items such as the location and type of installation, the technology of solar panels, as well as the ambient temperature. A typical photovoltaic irradiance-power curve is given in Equation (21).

$$P\_{PV}(\mathcal{J}\_s) = A\_{PV} \eta\_{pv} \mathcal{J}\_s \tag{21}$$

$$
\eta\_{pv} = \eta\_0 [1 - 0.0042 \left( \frac{\mathcal{I}\_s}{18} + T\_a - 20 \right)] \,\eta\_{\text{inv}} \tag{22}
$$

The amount of power generation by a solar panel depends on different factors, which can be estimated by Equations (21) and (22). Above, *η* stands for efficiency, and *Ta* denotes the ambient temperature.

In order to model the intermittencies, a dataset of historical records associated with the solar radiation. Such irradiance datasets can effectively be matched with a lognormal or Weibull distribution function. These distribution models usually have high mathematical compatibility with various natural phenomena. The PDF helps to generate randomly derived samples within the occurrence range while specifying a confidence range to eliminate low-probable cases. For each scenario, a random value is derived from the lognormal model. If the number of scenarios is high, various scenario reduction techniques can be employed that merge similar cases into one class. Based on the model objectives, the operation state can be adjusted within the confidence range with respect to the level of risk chosen by the operator. The operator can arrange the operation point in a risk-averse zone to ensure a reliable and secure operation, although it incurs higher costs. On the other hand, the operation of the grid, with respect to the whole uncertainties, can be stated in a high-risk zone by ignoring less-probable incidences to boost the profitability and the economy of the model.

#### 2.3.9. Load Forecast Uncertainty

Load forecasting is a critical item in many planning applications in power systems. Load forecasting can also be performed in different intervals for various purposes. Load forecasts are dependent on weather conditions (temperature, humidity, air brightness, wind speed). In addition, each day of the week has its own load curve. Load consumption curves on holidays and non-holidays are also different from each other. In different seasons of the year, according to the factors specific to each season, such as day length, the load consumption pattern will be different. Short-term load forecasting is a pivotal data needed for running a day-ahead market. The results of STLF are employed by GENCOs to discover prices with regard to the way other market participants may act. Price forecast signals are also applied to bidding strategy techniques to maximize profitability. The system operator will clear the market for the next-day delivery according to the forecasted values along with the bids from GENCOs. ISO is responsible for balancing the market between generation level and the forecasted demand at each interval. However, with regard to the inaccuracy of forecasts, there is always a violation that must be redressed by storage units or other controlling actions such as load shedding.

There are multiple approaches to deal with load forecasting. The most common methods are the time series method, regression methods, and intelligent methods. Intelligent methods are also classified into several categories such as artificial neural networks, fuzzy logic-based methods, also ANFIS models, wavelet transformation methods, as well as

support vector machines. In order to deal with uncertainties, probabilistic models can be integrated into the employed load forecasting paradigm. In this study, the Gaussian process quantile regression (GPQR) is employed to deal with uncertainties [38]. Due to the stochastic nature of load patterns as well as various external factors such as weather conditions, calendar effects, and seasonal factors, the power demand signals exhibit intermittent and volatile characteristics. Hence, a prediction scheme is needed to provide the most probable distribution of power demand patterns rather than a crisp value. A Gaussian process model is a set of random variables, in which any member corresponds with a probability distribution that can be represented by the following function:

$$f(l) \sim GP(\mu(l), COV(l, l'))\tag{23}$$

In this equation, *μ*(*l*) and COV(*l*,*l* ) denote the mean and covariance functions that can be calculated as follows:

$$
\mu(l) = E\{f(l)\}\tag{24}
$$

$$\text{COV}(l, l') = E\left\{ \left[ f(l) - \mu(l) \right] \times \left[ f\left(l'\right) - \mu\left(l'\right) \right]^T \right\} \tag{25}$$

The covariance function illustrates the similarity between data points. One of the most widely-used functions for describing covariance is the squared exponential (SE) model that can be represented as follows:

$$COV\_{SE}\left(l, l'\right) = \theta\_f^2 \exp\left(\frac{||l - l'||}{\theta\_{l \text{length}}^2}\right) \tag{26}$$

The parameters *θf* and *θl* control the scale and length. With regard to the differentiability of this covariance, which implies that the Gaussian process is very smooth, a simplified alternative, known as Matern covariance, can be replaced because, in practice, no phenomena do not reflect such a strong smoothness.

$$COV\_{Mat} \left( l - l' \right) = \sigma^2 \frac{2^{1-v}}{\Gamma(v)} \left( \sqrt{2v} \frac{l - l'}{i} \right)^v B\_v \left( \sqrt{2v} \frac{l - l'}{i} \right) \tag{27}$$

In the above, *Bv* represents the modified Bessel function, in which *v* and *i* are both positive hyperparameters. To simplify this covariance function, the value of *v* is supposed to be *v* = *p* + 1/2. In this parameter, *p* is a non-negative integer. The value of *v* is set as 5/2 or 3/2 in most previously conducted research studies, which subsequently and, respectively, are named Mat5 and Mat3. The probabilistic prediction model also has another covariance function named the period covariance, which is useful to model periodic phenomena and can be described as follows:

$$COV\_{Mat}\left(l, l'\right) = \theta\_f^2 \exp\left(-\frac{2}{\theta\_{length}^2} \sin^2\left(\pi \frac{l - l'}{p}\right)\right) \tag{28}$$

Typically, the load forecasts are highly influenced by a variety of features that can be defined by Equation (29). According to this equation, load at time *t* depends on the similar hours (*t*∈{0, 24}) throughout the historical records, the day of the year (*d*∈{1, ... , 365}), the value of the load at similar intervals, the value of weather variables such as temperature, and the price at similar intervals.

$$
\hat{y} = f(t, d, \upsilon\_l, \upsilon\_{w\_l}, \text{price}) \tag{29}
$$

GPQR method seeks for the relationships and correlation between inputs and output based on a probabilistic framework. Quantile regression (QR) delineates a type of regression analysis that detects and exploits the relationships between quantiles of the conditional distribution of a response variable and input variables. The least absolute deviations

regression integrates the median of the conditional distribution, which is called the norm regression case of quantile regression shown by *L1*. Unlike least-squares regression and norm regression, the quantile regression encompasses minimization of the summation of asymmetrically weighted absolute residuals. Hence, QR can find more sophisticated relationships between inputs and outputs with better precision. QR has more flexibility and compatibility to deal with large datasets, such as market analysis or econometrics. With regard to the accuracy of predictions, the loss function can be described as follows:

$$L\_{\pi}(\varepsilon\_{i}) = \begin{cases} \tau \varepsilon\_{i} & \text{if } & \varepsilon\_{i} \ge 0 \\ (\tau - 1)\varepsilon\_{i} & \text{if } & \varepsilon\_{i} \ge 0 \end{cases} \tag{30}$$

The required quantile is defined by *τ*∈[0, 1] and *εi* = *yi* − *y*ˆ*i* so that *yi* is the actual model and *y*ˆ*i* denotes the predicted quantile model. So far, different linear programming methods are employed to achieve the desired quantiles through direct loss function minimization, which leads to the maximization of a likelihood. To solve this drawback, the Gaussian process is incorporated into the QR model. The density function of loss based on this model can be represented as follows:

$$L(t|\mu, \sigma, \tau) = \frac{\tau(1-\tau)}{\sigma} \exp\left[-\frac{t-\mu}{\sigma}(t - I(t \le \mu))\right] \tag{31}$$

where *τ*∈[0, 1] is responsible for shaping and controlling the skewness of the distribution curve, *μ* stands for the mean value, and *σ* denotes the standard deviation, which should always be positive. The binary variable *I* takes the value of 1 when the condition is true; otherwise, it takes 0.

$$\mathcal{U}\_{\mathsf{T}}(y,q) = Z \exp\left[-\sum\_{i=1}^{N} L\_{\mathsf{T}}(y\_i, q\_i)\right] \tag{32}$$

In Equation (32), *q* denotes the predicted value of the *τ* quantile, Z is the normalization constant. Afterward, a Gaussian process is placed on the QR function:

$$p(q) = GP(q|0, K) \tag{33}$$

The GPQR model training can be conducted by integral maximization. The expectation propagation algorithm can be employed to locally approximate this integral.

$$\underset{q}{\text{argmax}} \int\_{q} \mathcal{U}\_{\text{r}}(y, q) p(q) \, d(q) \tag{34}$$

Suppose a dataset of historical load records as *l* = {*x1*, *x2*, ... , *xN*}, which are independently distributed samples. The estimated shape of this distribution function can be obtained through a Gaussian kernel density estimator model.

$$\mathcal{K}(\mathbf{x}\_1, \mathbf{x}\_2, \sigma) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{\left(\mathbf{x}\_1 - \mathbf{x}\_2\right)^2}{2r^2}} \tag{35}$$

Mean absolute percentage errors (MAPE) and root mean squared error (RMSE) are two employed evaluation metrics that can exhibit the performance of the forecasting model.

$$RMSE(y\_i, \mathcal{Y}\_i) = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} (y\_{i\prime} \mathcal{Y}\_i)^2} \tag{36}$$

$$MAPE(y\_{i\prime}, \hat{y}\_{i}) = \frac{1}{N} \sum\_{i=1}^{N} \left| \frac{y\_{i\prime} \hat{y}\_{i}}{y\_{i}} \right| \tag{37}$$

#### *2.4. Energy Management Paradigm by Microgrid's Controller*

The controlling framework of a microgrid is a vital element in this planning scheme. The microgrid controller, also known as the microgrid operator, is responsible for collecting raw data and analyzing them to find the optimum operation point with respect to all operational constraints of both sides. In other words, it has to provide an accurate estimation of generation sources for all intervals and take a sensible range of risk for being preserved from the detrimental zone. Then it can dedicate the level of power exchange with the upstream grid for various hours. To implement such an intelligent autonomous mechanism, a smart environment is needed that is only feasible through IoT infrastructures [39]. The following flowchart concisely describes the proposed scheme. The paradigm of the proposed scheme is depicted in Figure 2.

**Figure 2.** The paradigm of the proposed scheme.

#### **3. Lightning Attachment Procedure Optimization (LAPO) Method**

This algorithm was inspired by the procedure of lightning the attachment process to the ground or any type of earthed object. In order to simulate this optimization algorithm, at first, some test points are considered on a cloud or the ground as the initial population. By increasing the electric field at the points on the cloud, the streamer channels begin to break down the air and move toward the ground, which is called a downward leader. As these downward leaders move toward the ground, the opposite charges will be enhanced in the ground and produce upward leaders. The striking point at the final jump step is a point where the upward leader collides with the downward leader [40–42]. The algorithm is composed of some sections that will be presented as follows:

#### *3.1. Initialization*

At first, the test points on the surface of a cloud as well as on the ground are defined. These dedicated charges are called the initial population.

#### *3.2. The Next Jumping Node*

Each branch of the lightning has some probable points in front of itself, which can move toward one of them. The choice of the next jumping point of the lightning is highly correlated with the intensity of the electric field between the point and the probable points. Of course, the next jump will not necessarily be toward the point with the maximum

field, and part of the process of selecting the next point will occur randomly. For the mathematical modeling of this section, a random node is selected for each node in the search space. If the electric field in this point (merit value of the point or the fitness value at the intended point) is greater than the background field (which is supposed to be the mean field), the lightning moves toward this point; otherwise, the path of lightning will be in another direction. It is noticeable that the term opposite direction does not imply an upward movement exclusively rather than an angled motion toward this point. Hence, for determining the next jump, the following equation can be employed:

$$P = \text{sign}(F\_{av} - F(r))\tag{38}$$

$$X\_{\rm new}(i,j) = X(i,j) - P \times rand \times \left( Xav(1,j) - X(r,j) \right) \tag{39}$$

#### *3.3. Streaming Forward and Elimination of Branches*

When a new lightning branch is formed, this branch can stream downward until the branch's charge cuts down and becomes lower than a specific value (critical value). The critical value is assigned to be 1 μC because the air cannot be broken down for the lower levels of charges, and the branch will be faded.

#### *3.4. Upward Leader Propagation*

When the branches start to move downward, the upward leaders start to move upward. The motion speed of the upward leader is correlated with the amount of charge aggregated at this point as well as the total charge of the downward branch. The distribution of charge in the downward branch is so that there are fewer charges near the cloud and more charge at the tip of the leader. This is why the tip of the leader is brighter. This incremental charge aggregation over the branch can be considered as a linear or exponential function. Here, the branch's charge is considered exponentially, as pointed out by the following two equations:

$$S = 1 - \left(t/t\_{\text{max}}\right) \times \exp\left(-t/t\_{\text{max}}\right) \tag{40}$$

$$\mathcal{C}\mathcal{x} = \mathcal{S} \times \left( X\_{\text{min}}(i\_{\prime}f) - X\_{\text{max}}(i\_{\prime}f) \right) \tag{41}$$

Above, *CC* denotes the branch's charge, *X*min is the charge of the branch in the tip of the leader, *X*max indicates the branch's charge at the beginning point. The upward leaders will be propagated according to the following equation:

$$X\_{new}(i,j) = X(i,j) + rand \times \mathbb{C}c \tag{42}$$

#### *3.5. Convergence*

When an upward leader reaches a downward leader, the collision point of the lightning will be determined. In other words, when the optimum point is specified, the algorithm will be terminated. The abovementioned procedure will be iteratively executed until the optimum point is found.

#### **4. Simulation and Results**

In this section, the optimal placement for the microgrid is carried out on the targeted test system. The targeted system is a radial 69-bus IEEE test system, which is depicted in Figure 3 in the form of a single-line diagram. In addition, in order to have a better appraisal of the effectiveness of the proposed method, the proposed method is also tested on a radial 33-bus IEEE test system, which is depicted in Figure 4. However, the input data is the same for both grids. In this paper, all the methods are implemented in MATLAB 2018a in a Core i5 PC with 3 GHz processing frequency of CPU and 8 GB of RAM. The convergence behavior is another aspect with which the methods are compared to each other.

**Figure 3.** The single-line diagram of the 69-bus test system.

**Figure 4.** The single-line diagram of the 33-bus test system.

The characteristics corresponded with the lines and the buses of this system refer to [43,44], respectively. The normalized demand of the microgrid is shown in Figure 5. The demand for this microgrid at the peak is considered to be 250 kW. The peak of the load curve of the grid is supposed to be the mid-peak demand obtained from [43,44]. Moreover, the hourly solar insulation profile is normalized and illustrated in Figure 6.

**Figure 5.** The normalized demand of distribution network for both systems.

**Figure 6.** The normalized daily solar irradiance for both distribution systems.

The maximum generation capability of the PV panel per month is shown in Table 1. The placement mechanism is executed in the worst condition. Thus, the calculations are conducted for the month of August. This month indicates the lowest level of mean solar radiation.

**Table 1.** The mean solar radiation for each month for both radiation systems. [45].

The economic data for different sources were extracted from reference [44]. The investment cost of the PV panel is considered to be 7.44 €/kW, and the O&M cost is supposed to be 40 €/year. The cost of each battery is determined as 0.283 €/Ah, and each battery incurs 50 €/year for O&M costs. The operation cost and maintenance cost of the diesel generator is 0.55 €/W and 0.2 €/h, respectively.

The demand in the targeted microgrid is supposed to be composed of residential and industrial loads. Table 2 delineates the price of electricity corresponded with the types of loads and the demand level.

**Table 2.** The electricity price for different load types and various demand levels for both distribution systems.

#### *4.1. The IEEE Test Systems (33-Bus and 69-Bus)*

The optimization for the 33-bus test network is conducted with regard to the acquired economic data. The results are demonstrated in Table 3. In addition, the convergence behavior of two targeted methods in reaching the optimal solution is illustrated in Figure 7. Accordingly, the LAPO method has reached the optimal solution with fewer iteration. In terms of accuracy, the results of LAPO are better than the ABC algorithm. As can be seen, at the early hours of the time horizon (at midnight), when there is no solar power generation, a battery energy storage unit can supply the demand of internal loads of the microgrid as it is clear that the battery power is negative during these hours, which indicates that it is discharged. The battery starts charging after 8:00 and reaches its maximum charge level at 11:00. Therefore, no power is stored in it until it starts discharging again at 18:00. It is recharged at 19:00 and fully discharged at 20:00 due to the high volume of consumption. The diesel generator commits at 8:00 and stops power exchange at 9:00. These conditions continue until 20:00, and after that, it continues to serve the loads until midnight. Due to the significant amount of solar production, a large volume of this production is sold to the distribution network, and the obtained profit improves the objective function. As it is evident, the network is designed in such a way that no load is interrupted during the day and all the loads are completely satisfied.

**Table 3.** The design and size determination of microgrid for the distribution test systems (33-bus and 69-bus).

**Figure 7.** The comparison of the convergence behavior of different methods in seeking for the best solution for the 33-bus distribution test system.

With regard to the economic data provided in the previous section, the optimization problem is simulated on the 69-bus test system. The optimized results of the simulation are presented in Table 3. In addition, the convergence behavior of both optimization techniques employed in this study is depicted in Figure 8. As can be seen, it is evident that the LAPO approach not only has reached the optimal solution in less iteration, but also it has found a better solution in comparison with the ABC method.

**Figure 8.** The comparison of the convergence behavior of different methods in seeking for the best solution for the 69-bus distribution test system.

The behavior of the microgrid in August is illustrated in Figure 9. As can be perceived, in the early hours at midnight, while the sunrise has not occurred, the battery is responsible for maintaining power for the microgrid. Thus, the power battery is negative in these hours, which implies that the battery is discharging. The battery starts to be charged at 8:00 in the morning, and it will be fully charged at 11:00. Hence, the batteries will not be charged anymore and remain idle until 18:00. At 18:00, the batteries are to some extent discharged. However, they are charged again at 19:00. With regard to the high level of consumption at 20:00, the batteries deliver all of their stored energy to the microgrid. The diesel generator is connected to the microgrid at 8:00, and it is disconnected at 9:00. It continues until 19:00. Again, at 19:00, the diesel generator starts to supply the load until the end of the time horizon of this study (24:00). With respect to the high level of solar generation, a considerable share of this production is sold to the distribution grid. The obtained economic benefit has improved the objective function. As it is clear, the microgrid is so designed that the loads are supplied entirely and load shedding has not occurred. The microgrid behavior in August is depicted in Figure 9.

**Figure 9.** The behavior of the microgrid in August.

The power purchased from the main grid is shown in Figure 7. As it shows, during the day hours when there is solar energy, the demanded power from the upstream grid is drastically reduced. It should also be noted that the network losses and the corresponding costs were 1061 kW and 55.86 € per day, which reached the values of 956.79 kW and 51.78 € per day, after the installation of microgrids in the network, respectively. Figure 8 shows the amount of load shedding in the distribution network at different times of the day. The results figure out that during the hours when there is solar energy, solar unit injects power into the distribution network, the amount of shed load has reached zero, which is very beneficial in terms of cost and efficiency for the distribution network as well as reliability improvement from microgrid point of view. The model is also tested on the 69-bus test system. The results imply that the bought power from the main grid by the 69-bus distribution network is illustrated in Figure 10. As it is evident, during the hours that the PV panel can generate power, the absorbed power from the sub-transmission grid is considerably diminished. The power losses of the distribution grid and the subsequent cost of them were about 1350 kW equivalent with 65.76 €/day before the implementation of the microgrid scheme, and it is alleviated to 989.32 and 56.28 €/day after implementation of the microgrid. In Figure 11, the hourly shed load in the distribution network is shown. It is obvious that whenever the PV generation exists, all of the loads in the distribution grid are satisfied, and no load-shedding measure in the distribution network is executed. This matter significantly improves the performance of the distribution network.

**Figure 10.** The bought power from the main power system by the distribution network.

#### *4.2. Incorporation of Uncertainty in the Model*

The demand for a grid and the daily radiation intensity are not deterministic parameters, and they usually have uncertainties. Hence, these parameters must be forecasted, and the scheduling must be conducted based on the forecasts. Therefore, the existence of forecasting error is inevitable, and the error must be included in the optimization model in order to mitigate the risk. A common way to model the uncertainties corresponded with solar radiation and demand of the grid is to employ a normal probability distribution. In

other words, the mean value and the standard deviation of these parameters are calculated, and the design and the scheduling must be carried out based on these forecasts.

**Figure 11.** The amount of shed load in the distribution network in the presence and absence of microgrid.

The most famous and the most accurate approach for dealing with uncertainties and probabilistic problems is the Monte Carlo algorithm. In this approach, a large number of probabilistic samples are defined within a specified range, and the scheduling is performed for all of these samples. Then, the probability density function of targeted parameters is extracted. In this study, an upper and a lower boundary are dedicated for the demand of microgrid and distribution network along with the daily solar radiation. It is supposed that the radiation and demand will be materialized within the defined range. These boundaries are depicted in Figures 12–14.

**Figure 12.** The upper bound, the lower bound, and the mean value of demand of microgrid in the probabilistic study.

**Figure 13.** The upper bound, the lower bound, and the mean value of demand of distribution network in the probabilistic study.

**Figure 14.** The upper bound, the lower bound, and the mean value of daily solar irradiance in the probabilistic study.

In order to solve the stochastic problem, the Monte Carlo approach is employed. In this respect, for each uncertain parameter, 1000 stochastic samples are generated randomly. To solve the stochastic problem, the following step-by-step procedure is considered:

Table 4 outlines the results of the design of the microgrid with the incorporation of uncertainties.

**Table 4.** The results of optimal sizing of test grids incorporating uncertainties of demand and solar radiation and the consideration of the maximum cost of stochastic scenarios as the objective function.

In the case of the 33-bus test system, the cost of power generation is greatly increased. This increase is due to the conditions in which the load may be maximum and production may be minimum. Therefore, these conditions must also be taken into account in the problem. However, as can be seen, the cost has increased by nearly 49% compared to the case where uncertainty is not considered. Thus, it should be considered that in the presence of sources of intermittency, the cost of power supply from the network may increase remarkably, which conveys the importance of stochastic solution.

The increase in the system cost pertaining to the uncertainty will be 27%, which is nearly 20% less than the previous case. Therefore, the uncertainties are applied to the model using the average intermittencies in the system cost in all possible scenarios as an objective function. As can be seen, in the stochastic model, the cost of the system has increased by 26% compared to the deterministic case.

As the results of the 69-bus test system figure out, the cost of supply in this condition is significantly increased. This increase conveys a condition that the demand is at a high level and the solar radiation is at a low level. Therefore, the worst conditions must be contemplated in the design of the microgrid. However, as it is evident, the cost is drastically increased by 50% in comparison with the case of the deterministic model (regardless of uncertainties). Thus, it should be noted that the overall cost of the microgrid scheme may be enhanced extensively for uncertain problems.

It is important to pay attention to a critical question. The cost of supply has increased by 50% in order to include uncertainties. It conveys a high level of consumption and a low level of generation. The occurrence of such generation and consumption levels for this system at the targeted time is not a definitive and certain event. The question is whether it is sensible to increase the cost of the system by 50% for a forecast that has an occurrence probability of 10%.

In the procedure of seeking the optimum solution of uncertain parameters, for each suggested answer provided by the optimization program, all possible scenarios are taken into account, and the highest cost incurred to the system, among all scenarios, is regarded as the cost of the suggested answer. If the suggested answer does not satisfy even one or more constraints, that answer should be dismissed. Such a procedure profoundly mounts the system cost.

In order to tackle this problem, no answer should be dismissed, and the system cost must be obtained for all answers in as scenarios, even when they were unable to meet some constraints. Finally, the average cost of all scenarios is dedicated to the answer. Hence, if the suggested answer does not meet the constraints in some scenarios, the answer will not be dismissed. It is obvious that, for the problems with a few scenarios, the inclusion of scenarios, which used to be disregarded, does not have a remarkable impact on the solution of the problem. With consideration of the proposed condition, the optimal results of the simulation can be expressed in Table 5. This assumption has been led to an increase of 26% in the system cost for the incorporation of uncertainties in the model. It is just over 20% fewer than the normal approach. Thus, the uncertainties can be included in the design by incorporation of the average cost of the system. As can be seen, the system cost has increased by 26% in comparison with the deterministic approach.

**Table 5.** The results of optimal sizing of test grids incorporating uncertainties of demand and solar radiation and the consideration of the average cost of stochastic scenarios as the objective function.

#### **5. Conclusions**

In this study, a microgrid scheme including loads, PV panels, energy storage system, and the backup supply system of a diesel generator was designed. The object of the design was to minimize the cost of supply. In this microgrid, it was intended to not absorb any power from the distribution grid. However, if there is an excess of generation by intrinsic power resources of the microgrid, the excess power would be sold to the distribution network. The costs pertaining to PV panels, batteries, and diesel generator consists of purchasing and installation cost, operation and maintenance cost, and replacement cost. If the proposed scheme cannot supply the loads (generation deficiency) at some hours, the load-shedding measure must be imposed for the surplus demand. Just the same, the excess generation must be injected into the upstream network. In order to boost the performance of the design, the placement of the microgrid is so determined that the network power losses are minimized. Thus, the optimal design of the microgrid was performed at the same time as the placement of the microgrid. The problem solution was targeted to be stochastic. Hence, two upper and lower boundaries were dedicated to the amount of demand and solar irradiance. In order to solve the problem, the Monte Carlo approach with 6222 random samples is employed. The results of the 33-bus test system imply that the incorporation of uncertainties in the model has drastically increased the supply cost by about 49%. In order to avoid unnecessary risks, the averaging method is deployed, which

boosts the risk-averse scheduling. In this case, the system cost has increased by 27% in comparison with the deterministic approach. The results of the 69-bus test system imply that the incorporation of uncertainties in the model has drastically increased the supply cost by about 50%. In order to avoid unnecessary risks, the averaging method is deployed, which boosts the risk-averse scheduling. In this case, the system cost has increased by 26% in comparison with the deterministic approach. The optimization is solved using the LAPO method, and the results are compared with the ABC algorithm. The results indicate the proper performance of the suggested algorithm in terms of convergence speed and accuracy.

**Author Contributions:** Conceptualization, A.F.N., H.S., H.N., B.V., Y.A. and M.B.; Methodology, A.F.N., H.S., H.N., B.V., Y.A. and M.B.; Supervision, A.F.N., H.N., B.V.; Validation, A.F.N., H.N., B.V.; Visualization, A.F.N., H.N., B.V.; Writing—original draft, A.F.N., H.S., H.N., B.V., Y.A. and M.B.; Writing—review & editing, A.F.N., H.S., H.N., B.V., Y.A. and M.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

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

#### **References**

#### *Article*

## **Development and Implementation of a Novel Optimization Algorithm for Reliable and Economic Grid-Independent Hybrid Power System**

#### **Mohammed Kharrich 1, Omar Hazem Mohammed 2, Salah Kamel 3, Ali Selim 3,4, Hamdy M. Sultan 5, Mohammed Akherraz 1 and Francisco Jurado 4,\***

Received: 17 August 2020; Accepted: 16 September 2020; Published: 21 September 2020

**Abstract:** Recently, fast uptake of renewable energy sources (RES) in the world has introduced new difficulties and challenges; one of the most important challenges is providing economic energy with high efficiency and good quality. To reach this goal, many traditional and smart algorithms have been proposed and demonstrated their feasibility in obtaining the optimal solution. Therefore, this paper introduces an improved version of Bonobo Optimizer (BO) based on a quasi-oppositional method to solve the problem of designing a hybrid microgrid system including RES (photovoltaic (PV) panels, wind turbines (WT), and batteries) with diesel generators. A comparison between traditional BO, the Quasi-Oppositional BO (QOBO), and other optimization techniques called Harris Hawks Optimization (HHO), Artificial Electric Field Algorithm (AEFA) and Invasive Weed Optimization (IWO) is carried out to check the efficiency of the proposed QOBO. The QOBO is applied to a stand-alone hybrid microgrid system located in Aswan, Egypt. The results show the effectiveness of the QOBO algorithm to solve the optimal economic design problem for hybrid microgrid power systems.

**Keywords:** economic energy; Bonobo Optimizer; hybrid renewable energy system; microgrid; PV panels; wind turbine; energy storage

#### **1. Introduction**

Despite the steady increase in electric power production, it is still below the required level, due to the increase in load demand caused by the population increase as well as the increased use of technology in the residential, industrial and agricultural fields. According to the International Energy Agency (IEA), the global electricity demand will grow at an annual rate of 2.1% until 2040. This increases electricity's share in the total energy consumption to 24% in 2040. It is expected that renewable energy sources (RES) will face a significant increase in global investment in the coming years, to cover more than half of the energy consumption in the world by 2040. These energies will make up for the shortfall in electrical energy production and contribute to a reduction in carbon dioxide emissions in the atmosphere, thereby reducing pollution significantly [1–3].

In order to invest in RES to optimize electrical energy production and raise the efficiency of the systems, many studies in the world recommend combining different technologies to form hybrid renewable energy systems (HRES) [4,5]. Consequently, these sources complement each other, support the national grid, and reduce the use of traditional power plants depending on fossil fuels that release greenhouse gases and pollute the environment [6]. However, the design of these hybrid systems needs sophisticated programs and smart algorithms capable of reaching the optimal solution taking into consideration all the conditions and constraints such as reliability aspects, economic cost, sensitivity factors, availability of RES, etc. [1,2,7,8].

Several studies have been conducted on the technical and economic feasibility of hybrid systems in past years to determine their viability. Many of these studies have used different modeling of HRES, and they have applied different algorithms and various software tools to achieve their goals. According to the literature, these challenges still exist and are the focus of a lot of research, especially on finding the best algorithms and modern techniques in reaching the optimal solutions of the optimization problem of finding the optimal sizing of the installed capacities of the components of HRES [9–15].

In [16], the pre-feasibility analysis of a stand-alone energy system using HRES including renewable and conventional energy sources was applied using HOMER software in Newfoundland, Canada. In one of the earlier studies [17], the authors conducted a feasibility study of generating electricity using RES for a hybrid system in a stand-alone village in Chhattisgarh, India. In [18], the authors introduce a realistic solution for energy demand from a hybrid power system consists of wind turbines (WT), photovoltaics (PV), and battery energy storage systems (BESS). Through a real measurement of meteorological data in 2017, concerning especially the wind speed, solar radiation and temperature, the output power of the proposed hybrid system is calculated. Load satisfaction is considered to evaluate the feasibility of the system. The optimum solution is found using the linear TORSCHE optimization technique, while a comparative study between PV/WT/battery and PV/WT has been accomplished and an economic analysis was presented. As a result, the hybrid PV/WT/battery is proved more economical than using each system individually.

Xiao Xu et al. [19] designed and investigated a hybrid PV/WT/hydropower/pump storage as a case of study. The optimal configuration of the HRES is found using a techno-economic index that respects the maximum Loss of Power Supply Probability (LPSP) and minimum investment cost. The Multi-Objective Particle Swarm Optimization (MOPSO) is used to trade off analysis between two objectives. Besides, the curtailment rate (CR) of the WT and PV are taken into consideration due to policy requirements. The authors in [20] proposed an optimized design of an energy system featuring the highest penetration of renewable energy. This system is composed of WT, PV, geothermal, diesel, and BESS; otherwise, the system is obtained respecting the technological and financial feasibility constraints. The model developed is based on weather and electric demand data measured to reach the optimal sizing of the hybrid system. Three objective functions are conflicting, which are the Net Present Cost (NPC), renewable energy fraction and the energy index of reliability.

In [21], the authors implemented and compared three algorithms to find the optimal design of a hybrid WT/PV/Biomass/BESS energy system. Based on the obtained results, the Harmony Search Algorithm (HSA) was faster and efficient in the convergence, compared to Jaya and PSO optimization algorithms. The techno-economic study has been implemented to have the optimal unit sizing of the HRES, which guaranteed a cost-effective, efficient, and reliable power supply for the customers of electric energy. The constraints are chosen to enhance the reliability and efficiency of the hybrid system, using the LPSP and the energy fraction factors.

In this paper, a new smart algorithm named Bonobo Optimizer [22] was employed and improved using a quasi-oppositional method, and the modified Quasi Oppositional BO (QOBO) was utilized for optimal economic designing of a stand-alone microgrid hybrid system in Aswan, Egypt, where the hybrid system consists of RES (PV panels, WT and BESS) with diesel generators. Then, the results were compared between the traditional and improved BO. This proved the ability of the QOBO algorithm to reach the optimal solution in a shorter time and with better efficiency compared to the traditional BO algorithm. Other algorithms, namely Harris Hawks Optimization, Artificial Electric Field Algorithm and Invasive Weed Optimization are applied, and the results are compared where the efficiency of the QOBO algorithm has been proved. Additionally, a sensitivity analysis of the proposed systems scenarios was performed to obtain the optimal solution.

#### **2. Mathematical Description of the Proposed Hybrid System Components**

The schematic diagram of the suggested HRES is shown in Figure 1. Four scenarios are applied, which include the PV power plant, WT power plant, diesel generator, Biomass, BESS and inverter.

**Figure 1.** Configuration of the proposed microgrid hybrid energy system.

Two strategies are adopted in this paper; the first is the biomass/PV as shown in Figure 2 and the second uses the PV or WT or both as in Figure 3. The main strategy steps for the operation of the proposed system can be explained as follows:

#### *2.1. PV System*

The PV system is considered as a number of cells connected in series. The output power of the PV system is presented based on many parameters as introduced in Equation (1) [23]:

$$P\_{pv} = I(t) \times \eta\_{pv}(t) \times A\_{pv} \tag{1}$$

where *I* represents the solar irradiation, *Apv* represents the area covered with PV modules and η*pv* is the efficiency of the PV system that can be calculated as follows:

$$\eta\_{\overline{\eta}\overline{\eta}}(t) = \eta\_{\overline{\tau}} \times \eta\_{\overline{t}} \times \left[1 - \beta \times (T\_4(t) - T\_{\overline{\tau}}) - \beta \times I(t) \times \left(\frac{\text{NOCT} - 20}{800}\right) \times (1 - \eta\_{\overline{\tau}} \times \eta\_{\overline{t}})\right] \tag{2}$$

where *NOCT* is the nominal operating cell temperature (◦C), η*r* is the reference efficiency, η*t* is the efficiency of the maximum power point tracking (MPPT) equipment, β is the temperature coefficient, *Ta* is the ambient temperature (◦C), *Tr* is the solar cell reference temperature (◦C).

**Figure 2.** Power management of the PV/Biomass hybrid renewable energy sources (RES).

**Figure 3.** Power management of the PV/WT/diesel/battery energy storage system (BESS), PV/diesel/BESS and WT/diesel/BESS hybrid RES.

#### *2.2. Wind Energy System*

Based on the basics of aerodynamics, wind power can be presented as [24]:

$$P\_{\rm wind} = \begin{cases} 0, & V(t) \le V\_{ci\prime} \; V(t) \ge V\_{co} \\ a \times V(t)^3 - b \times P\_{r\prime} & V\_{ci} < V(t) < V\_r \\ & P\_{r\prime} & V\_r \le V(t) < V\_{co} \end{cases} \tag{3}$$

where *V* represents wind speed, *Pr* is the rated power of wind, *Vci*, *Vco* and *Vr* are the cut-in, cut-out, and rated wind speed, respectively. *a* and *b* are two constants, which can be expressed as:

$$\begin{cases} \quad a = P\_r / \left(V\_r^{\;3} - V\_{c\bar{l}}^{\;3}\right) \\ \quad b = V\_{c\bar{l}}^{\;3} / \left(V\_r^{\;3} - V\_{c\bar{l}}^{\;3}\right) \end{cases} \tag{4}$$

The rated power of wind is calculated as given in the following equation:

*Pr* = 1 2 × ρ × *Awind* × *Cp* × *Vr* 3 (5)

where ρ represents the air density, *Awind* is the swept area of the wind turbine, *Cp* is the maximum power coefficient ranging from 0.25% to 0.45%.

#### *2.3. Biomass System*

The biomass system is a renewable energy system, which produces power as given in Equation (6) [23].

$$P\_{\rm BM} = \frac{Total\_{\rm bio} \times 1000 \times CV\_{\rm bio} \times \eta\_{\rm bio}}{8760 \times O\_{\rm tim}} \tag{6}$$

where *Totalbio* is the total organic material of biomass, *CVbio* is the calorific value of the organic material (20 MJ/kg), η*bio* is the biomass efficiency, which is taken as 24% and *Otime* presents the operating hours each day.

#### *2.4. Diesel System*

The diesel generator is used as a back-up, working just in case there is a need, is connected directly with the load, and starts when the battery is fully discharged and the load is more than 30% of its rated capacity. The model of the diesel generator regarding its output power is presented by the following Equation [25]:

$$P\_{d\text{g}} = \frac{F\_{d\text{g}}(t) - A\_{\text{g}} \times P\_{d\text{g,out}}}{B\_{\text{g}}} \tag{7}$$

where *Fdg* is fuel consumption, *Pdg*,*out* is the output power of diesel generator, *Ag* and *Bg* are the constants of the linear consumption of the fuel.

#### *2.5. BESS System*

The battery energy storage system (BESS) is a mandatory element for the isolated hybrid systems. BESS is charged in the periods of power excess and discharged when the load increases. The capacity of the BESS is expressed as follows [25]:

$$\mathcal{C}\_{\text{bat}} = \frac{E\_l \times AD}{DOD \times \eta\_{inv} \times \eta\_{b}} \tag{8}$$

where *El* is the load demand, *AD* represents the autonomy daily of the battery, *DOD* is the depth of discharge of the battery system, η*inv* and η*b* are the battery and inverter efficiency, respectively.

#### **3. Formulation of the Optimization Problem**

#### *3.1. Net Present Cost*

The objective function in the optimization model is the minimization for the Net Present Cost (NPC) which is the pillar factor considered for any project design; it is counted as a sum of all components costs including the capital (*C*), operation and maintenance (*OM*) and replacement costs (*R*), considering also the fuel cost of the diesel *FCdg* , taking into account the interest rate (*ir*), inflation rate (δ), and escalation rate (μ) and the predefined project lifetime (*N*). The NPC modeling is expressed as follows [23,24]:

$$NPC = \mathbb{C} + OM + R + FC\_{d\mathfrak{g}} \tag{9}$$

#### 3.1.1. PV and WT Costs

The costs of PV and WT are presented in a similar concept, their capital cost is expressed based on its initial cost (λ*PV*,*WT*) and its area (*APV*,*WT*), the capital cost is as follows [26]:

$$C\_{PV,WT} = \lambda\_{PV,WT} \times A\_{PV,WT} \tag{10}$$

The operation and maintenance costs are expressed as [26]:

$$OM\_{PV,WT} = \theta\_{PV,WT} \times A\_{PV,WT} \times \sum\_{i=1}^{N} \left(\frac{1+\mu}{1+i\_r}\right)^i \tag{11}$$

where θ*PV*,*WT* is the annual operation and maintenance costs for any components. The replacement costs are considered null because the project lifetime and the *PV* or *WT* lifetime are the same.

#### 3.1.2. Diesel Generator Costs

The costs of the diesel generator are presented as follows [27]:

$$C\_{d\S} = \lambda\_{d\S} \times P\_{d\S} \tag{12}$$

$$\text{OM}\_{d\S} = \theta\_{d\S} \times N\_{rm} \times \sum\_{i=1}^{N} \left( \frac{1+\mu}{1+i\_r} \right)^{i} \tag{13}$$

$$R\_{\text{diesel}} = R\_{d\text{g}} \times P\_{d\text{g}} \times \sum\_{i=7,14\dots} \left(\frac{1+\delta}{1+i\_I}\right)^i \tag{14}$$

$$C\_f(t) = p\_f \times F\_{d\lg}(t)\tag{15}$$

$$F\mathbb{C}\_{d\S} = \sum\_{t=1}^{8760} \mathbb{C}\_f(t) \times \sum\_{i=1}^{N} \left(\frac{1+\delta}{1+i\_I}\right)^i \tag{16}$$

where *Cdg* is the capital cost, λ*dg* is the initial cost of the diesel generator for each KW, *OMdg* represent the actual O&M cost, θ*dg* is the annual O&M cost of diesel, *Nrun* is the number of operating hours of diesel generator per year, *Rdiesel* is the diesel generator replacement cost, *Rdg* represents the annual replacement cost of diesel generator, *pf* is the fuel cost, *Fdg* is the annual consumption of fuel and *FCdg* is the total fuel cost.

#### 3.1.3. BESS Costs

The capital and O&M (containing the replacement) costs of the BESS are expressed as follows [26]:

$$
\mathbb{C}\_{BESS} = \lambda\_{\text{but}} \times \mathbb{C}\_{\text{but}} \tag{17}
$$

$$OM\_{BESS} = \theta\_{\text{hit}} \times \mathbb{C}\_{\text{bat}} \times \sum\_{i=1}^{T\_B} \left(\frac{1+\mu}{1+\delta}\right)^{(i\_-1)N\_{\text{bat}}} \tag{18}$$

where λ*bat* is the BESS initial cost and θ*bat* is the annual O&M cost of BESS.

#### 3.1.4. Biomass Costs

The biomass costs are presented as follows [28]:

$$C\_{\rm bg} = \lambda\_{\rm bg} \times P\_{\rm bg} \tag{19}$$

$$OM\_{b\S} = \theta\_1 \times P\_{b\S} \times \sum\_{i=1}^{N} \left(\frac{1+\mu}{1+i\_r}\right)^i + \theta\_2 \times P\_w \times \sum\_{i=1}^{N} \left(\frac{1+\mu}{1+i\_r}\right)^i \tag{20}$$

where λ*bg* is the biomass initial cost, θ1 is the annual fixed O&M cost and θ2 is the variable O&M cost of the biomass system, and *Pw* is the annual energy generated by the Biomass system (kWh/Year).

#### 3.1.5. Inverter Costs

The inverter capital and O&M costs are presented as follows [27]:

$$\mathbf{C}\_{inv} = \lambda\_{inv} \times P\_{inv} \tag{21}$$

$$OM\_{Inv} = \theta\_{Inv} \times \sum\_{i=1}^{N} \left(\frac{1+\mu}{1+i\_r}\right)^i \tag{22}$$

where λ*inv* is the inverter initial cost and θ*Inv* is the annual O&M cost of the inverter.

#### *3.2. Levelized Cost of Energy*

The Levelized Cost of Energy (*LCOE*) is a critical factor. The consumers do not care about project cost or its lifetime, but their interest is to know how much to pay for each kilowatt-hour of consumption. Therefore, the *LCOE* is a measure of the average *NPC* over its lifetime, its equation is expressed as follows [25]:

$$LCOE = \frac{NPC \times CRF}{\sum\_{t=1}^{8760} P\_{load}(t)} \tag{23}$$

where *Pload* is the load demand; *CRF* is the capital recovery factor used to convert the initial cost to an annual capital cost, and is expressed as follow:

$$\text{CRF}(ir, R) = \frac{i\_r \times \left(1 + i\_r\right)^R}{\left(1 + i\_r\right)^R - 1} \tag{24}$$

where *R* denotes the lifetime of the hybrid system.

#### *3.3. Loss of Power Supply Probability*

The loss of power supply probability (LPSP) is a technical factor used to express the reliability of the system. The LPSP is expressed as follows [25]:

$$LPSP = \frac{\sum\_{t=1}^{8760} \left( P\_{load}(t) - P\_{pv}(t) - P\_{wind}(t) + P\_{dg,out}(t) + P\_{bmin} \right)}{\sum\_{t=1}^{8760} P\_{load}(t)} \tag{25}$$

#### *3.4. Renewable Energy Fraction*

The transfer from classical electricity production to renewable energy projects was not easy. The majority introduced RES partially, while the objective is to use all projects with 100% renewable energy. Therefore, the renewable energy factor is dedicated to calculating the percentage of the renewable energy used. The renewable energy fraction (*RF*) is expressed as follows [25]:

$$RF = \left(1 - \frac{\sum\_{t=1}^{8760} P\_{d\lg out}(t)}{\sum\_{t=1}^{8760} P\_{r\varepsilon}(t)}\right) \times 100\tag{26}$$

where *Pre* represents the total power from RES.

#### *3.5. Availability Index*

The availability index (*A*) is calculated to predict customer satisfaction. The availability index measures the energy converted to the load while confirming the ability of the designing system of the project. The availability index is calculated as follows [23]:

$$A = 1 - \frac{DMN}{\sum\_{t=1}^{8760} P\_{\text{load}}(t)}\tag{27}$$

$$DMN = P\_{b\text{min}}(t) - P\_b(t) - \left(P\_{pv}(t) + P\_{w\text{mid}}(t) + P\_{dg\text{,out}}(t) - P\_{load}(t)\right) \times u(t) \tag{28}$$

while, *u* will be equal to 1 when the load is not satisfied, and 0 when the load is satisfied.

#### *3.6. Constraints*

The constraints are presented to achieve the desired system design. In this microgrid system, the constraints are shown as follows:

$$\begin{array}{l} 0 \le A\_{\textit{pv}} \le A\_{\textit{pv}}^{\text{max}}, \\ 0 \le A\_{\textit{wind}} \le A\_{\textit{wind}}^{\text{max}} \\ 0 \le P\_{d\textit{gl}} \le P\_{d\textit{gl}}^{\text{max}}, \\ 0 \le P\_{\textit{Cap},\text{flat}} \le P\_{\textit{Cap\\_bat}}^{\text{max}} \\ \underline{LPSP} \le \underline{LPSP}^{\text{max}}, \\ \underline{RY}^{\text{min}} \le \underline{RY}, \\ \underline{A}^{\text{min}} \le \underline{A} \\ \underline{AD}^{\text{min}} \le \underline{AD} \end{array} \tag{29}$$

#### **4. Algorithms**

In this section, the conventional BO and proposed QOBO are illustrated. In addition, both algorithms are compared with well-known optimization techniques (HHO, AEFA and IWO) which are briefly described in Appendix A.

#### *4.1. Bonobo Optimizer*

Bonobo optimizer is a new optimization algorithm that was proposed in [22]. In BO, the social reproductive behavior of the bonobo is modeled based on four mating strategies: promiscuous, restrictive, consortship, and extra-group mating. These mating strategies are subjected to the living condition of the bonobo, hence two terms named positive phase (PP) and negative phase (NP) have been used to present the situations of this life. In this framework, PP describes the peaceful living in which the mating can be done. On the contrary, NP expresses a hard life. In the BO, each solution is called *XB* and the best solution is *X*α *B*. The mathematical modeling of the BO algorithm is presented in the following subsections.

#### 4.1.1. Bonobo Selection Using Fission–Fusion Strategy

The solutions update of the BO algorithm depends on the mating strategies subjected to the current phase. However, a bonobo should be selected before each mating based on the fission–fusion social group strategy. As noted, the bonobo community lives in small groups with different sizes (random and unpredictable) for a few days and the communities rejoined again to the main community. Hence, based on this behavior, a bonobo for mating can be selected. The mathematical formulation for the maximum number of these temporary subgroups *Nsub* can be expressed as follows:

$$N\_{\rm sub} = \max\left(2, \left(\varepsilon\_{\rm sub} \times N\right)\right) \tag{30}$$

where *N* is the total number of the population and ε*sub* denotes the sub-group size factor. To find the selected bonobo *XP B* for mating with *Xi B* to create a new bonobo *Xnew B* , if the best bonobo in the subgroup in terms of the fitness function is better than the *Xi B*, then it is selected as *XP B*, else a random one should be selected form the subgroup.

#### 4.1.2. Creation of New Bonobo

After achieving the selected bonobo *XP B*, four mating strategies are used in the BO algorithm to create a new bonobo *Xnew B* based on the current phase (PP or NP). For the PP case, promiscuous and restrictive mating have higher probabilities ρ*ph* for occurrence. On the contrary, in NP, the probabilities

ρ*ph* of consortship mating and extra-group mating are higher.

Promiscuous and Restrictive Mating

In this mating strategy, the new bonobo can be created using the following equation:

$$X\_B^{\text{new}} = \left. X\_B^i + r\_1 \times \mathcal{S}\_{\text{ceef}}^a \times \left( X\_B^a - X\_B^i \right) + (1 - r\_1) \times \mathcal{S}\_{\text{ceef}}^p \times \mathcal{C}\_{f \text{lig}} \times \left( X\_B^i - X\_B^p \right) \tag{31}$$

where *r*1 is a random number between [0, 1]. *S*α *coe f* and *SP coe f* are the sharing coefficients for the alpha bonobo *X*α *B* and the selected bonobo *XP B*, respectively. *Cflag* is a flag value that equals −1 or 1 for restrictive and promiscuous mating, respectively. A controlling parameter in terms of the phase probability ρ*ph* is used to adopt the mating strategy. Initially, ρ*ph* is set to 0.5. Hence, if a random number *r* is found to be less than or equal to ρ*ph*, a new bonobo is created based on promiscuous and restrictive mating, otherwise, consortship mating and extra-group mating can be used.

Consortship and Extra-Group Mating

If *r* is greater than ρ*ph*, consortship and extra-group mating can occur. However, a new random number*r*2 between [0, 1] is used with a probability of extra-group mating ρ*xg* to represent the occurrence of extra-group mating when *r*2 is less than or equal to ρ*xg* as follows [22,29]:

$$X\_{B}^{\text{new}} = \begin{cases} X\_{B}^{i} + \beta\_{1} \times \left(X\_{\text{max}}^{i} - X\_{B}^{i}\right), & X\_{B}^{\alpha} \ge X\_{B}^{i}, and \, r\_{4} \le \rho\_{d} \\\ X\_{B}^{i} - \beta\_{2} \times \left(X\_{B}^{i} - X\_{\text{min}}^{i}\right), & X\_{B}^{\alpha} \ge X\_{B}^{i}, and \, r\_{4} > \rho\_{d} \\\ X\_{B}^{i} - \beta\_{1} \times \left(X\_{B}^{i} - X\_{\text{min}}^{i}\right), & X\_{B}^{\alpha} < X\_{B}^{i}, and \, r\_{4} \le \rho\_{d} \\\ X\_{B}^{i} + \beta\_{2} \times \left(X\_{\text{max}}^{i} - X\_{B}^{i}\right), & X\_{B}^{\alpha} < X\_{B'}^{i}, and \, r\_{4} > \rho\_{d} \end{cases} \tag{32}$$

$$\beta\_{1} = e^{\left(r\_{4}^{2} + r\_{4} - \frac{2}{r\_{4}}\right)} \tag{33}$$

$$\beta\_{2} = e^{\left(-r\_{4}^{2} + 2r\_{4} - \frac{2}{r\_{4}}\right)}$$

where *r*3 and *r*4 are random numbers between [0, 1] and *r*4 - 0. ρ*d* is a directional probability with initial value which equals 0.5. β1 *and* β2 are intermediate parameters between [0, 1]. *Xi min* and *Xi min* are the values of the upper and lower boundary.

If *r*2 is greater than ρ*xg*, a new bonobo can be created using the consortship mating strategy as follows:

$$X\_B^{new} = \begin{cases} X\_B^i + \mathbb{C}\_{f\text{lag}} \times e^{-r\_5} \times \left(X\_B^i - X\_B^P\right), & \mathbb{C}\_{f\text{lag}} = 1 \text{ or } r\_6 \le \rho\_d\\ & X\_{B'}^P \quad \text{Otherwise} \end{cases} \tag{34}$$

where *r*5 and *r*6 are two random numbers.

#### 4.1.3. Parameter Updating

The BO's parameters are updated during the iterative process based on the best solution *X*α *B* at each iteration, where if there is an improvement in the final solution compared to the previous iteration, the BO's parameters can be updated in the following way.

The negative phase count is set to zero (*NPcont* = 0) and the positive phase count grows by increments of one (*PPcont* = *PPcont* + 1). In addition, ρ*xg* = ρ*xg*\_*initial* and ρ*ph* = 0.5 + *Cp* where *Cp* is the amount of the change in the phase, and can be calculated as *Cp* = min(0.5, *PPcont* × *rcp*) where *rcp* is the rate of the change in the phase. Moreover ρ*d* = ρ*ph* and

$$
\varepsilon\_{\text{sub}} = \min \left( \varepsilon\_{\text{sub\\_max}}, \left( \varepsilon\_{\text{sub}\_{initial}} + \text{PP}\_{\text{cont}} \times \text{rcp}^2 \right) \right) \tag{35}
$$

where ε*subinitial* = 0.5 ∗ ε*sub*\_*max*.

On the other hand, if there is no improvement, the BO's parameters are updated as follows:

$$\begin{aligned} NP\_{\text{cont}} &= NP\_{\text{cont}} + 1 \text{ and } PP\_{\text{cont}} = 0, \\ \mathbb{C}p &= \min(0.5, \ NP\_{\text{cont}} \times rcp), \\ \rho\_{\text{x}\emptyset} &= \rho\_{\text{x}\emptyset, \text{initial}} \min(0.5, \ \rho\_{\text{x}\emptyset, \text{initial}} + NP\_{\text{cont}} \times rcp^2), \\ &\quad \text{and} \\ \varepsilon\_{sub} &= \min\{\varepsilon\_{sub\bot\text{sup}\bot} \left(\varepsilon\_{sub\bot\text{inf}} - NP\_{\text{cont}} \times rcp^2\right)\}. \end{aligned}$$

The overall steps of the BO algorithm are presented in Algorithm 1.

#### **Algorithm 1**: BO

Initialize a set of random search bonobo *Xi B* = *X*1 *B*, *X*2 *B*, ... , *XN B* within the limits *Xi min* ≤ *Xi B* ≤ *Xi max*. Initialize the BO's parameters Evaluate the objective function for all bonobos Identify the alpha bonobo *X*α *B* While (k < *Kmax*) Determine the actual size of the temporary sub-group Choose a bonobo using fission-fusion society strategy Create a new bonobo *Xnew B* as follows: if *r* ≤ ρ*ph* Create new bonobo using promiscuous or restrictive mating strategy else *r* > ρ*ph* Create new bonobo using consortship or extra-group mating strategy end if Calculate the objective function Update alpha bonobo *X*α *B* and the BO's parameters. *K* = *K* + 1 end while Return the final best solution *X*α *B*

#### *4.2. Improved Quasi-Oppositional BO (QOBO) Algorithm*

As with any population-based algorithm, BO has some problems such as falling in the local optima. However, in this work, an improved BO based on three leaders' selection and quasi-opposition-based learning is developed.

#### 4.2.1. Three Leaders

In this method, instead of using the best solution (alpha bonobo *X*α *B*) for updating the new bonobo *Xnew B* and ignoring the other best solutions, three leaders can be used to increase the diversity of the solutions as follows

$$X\_B^{\alpha} = w\_1 \times X\_{\text{best}\_1} + w\_2 \times X\_{\text{best}\_2} + w\_3 \times X\_{\text{best}\_3} \tag{36}$$

where

$$w\_1 = \frac{r\_7}{r\_7 + r\_8 + r\_9}, \; w\_2 = \frac{r\_8}{r\_7 + r\_8 + r\_9}, \; \text{and } w\_1 = \frac{r\_9}{r\_7 + r\_8 + r\_9}$$

*r*7, *r*8, *and r*9 are random values between [0, 1].

#### 4.2.2. Quasi-Oppositional

Opposition-based learning (OBL) [30] has been widely used to improve many optimization techniques such as quasi-oppositional teaching-learning (QOTLBO) [31,32], Quasi-oppositional swine influenza model-based optimization with quarantine (QOSIMBO-Q) [33] and Oppositional Jaya Algorithm [34]. In the OBL, improvements can be achieved by using the candidate solution and its opposite at the same time. Hence, in this work, the opposite solution of the BO algorithm *Xi B* can be expressed as presented in [35]:

$$X\_B^{\text{quocw}} = \mathbb{C} + r\_{10} \left( \mathbb{C} - \overline{X\_B^{\text{new}}} \right) \tag{37}$$

where *r*10 is a random number between [0, 1], and *C* is a middle point between *Xi min* and *Xi max* which can be calculated as follows:

$$C = \frac{X\_{\text{min}}^i + X\_{\text{max}}^i}{2} \tag{38}$$

Additionally, *Xnew B* is the opposite solution which can be calculated as

$$
\overline{X\_B^{ncw}} = X\_{min}^i + X\_{max}^i - X\_B^{ncw} \tag{39}
$$

The overall steps of the improved BO based on three leaders and the quasi-oppositional method are presented in Algorithm 2.

#### **Algorithm 2**: QOBO

Initialize a set of random search bonobo *Xi B* = *X*1 *B*, *X*2 *B*, ... , *XN B* within the limits *Xi min* ≤ *Xi B* ≤ *Xi max*. Initialize the BO's parameters Evaluate the objective function for all bonobos Determine the alpha bonobo *X*α *B* using three-leader method While (k < *Kmax*) Determine the actual size of the temporary sub-group Choose a bonobo using fission-fusion society strategy Create a new bonobo *Xnew B* as follows: if *r* ≤ ρ*ph* Create new bonobo using promiscuous or restrictive mating strategy else *r* > ρ*ph* Create new bonobo using consortship or extra-group mating strategy end if Calculate the objective function for all new bonobos *Xnew B* Find quasi-oppositional model for all new bonobos *Xqnew B* Calculate the objective function for all new bonobos *Xqnew B* if *f Xqnew B* ≤ *f Xinew B* , *Xnew B* = *Xqnew B* Else *Xnew B* = *Xnew B* end if Update alpha bonobo *X*α *B* using three leader method and the BO's parameters. *K* = *K* + 1 end while Return the final best solution *X*α *B*

#### **5. Case Study**

To validate the robustness of the QOBO algorithm, it has been applied for addressing the studied problem of optimal configuration of the proposed multiple scenarios HRES, i.e., the PV/WT/diesel generator/BESS, PV/biomass, PV/diesel generator/BESS and WT/diesel generator/BESS. The proposed hybrid systems have been introduced in the isolated mode for satisfying the load requirements in the proposed site.

The project is applied in Aswan, Egypt as shown in Figure 4. The annual load curve over a time interval of one hour is shown in Figure 5. Figures 6–9 present solar irradiation, temperature, wind speed and atmospheric pressure in the studied region. Four standalone scenarios of the hybrid system will be evaluated for covering the load demand in that site. These configurations are: (1) PV/WT/diesel/BESS, (2) PV/biomass, (3) PV/diesel/BESS and (4) WT/diesel/BESS. The proposed QOBO is validated on optimal sizing of these four hybrid systems and the optimization results are comprehensively compared with the corresponding ones obtained from BO, HHO, AEFA and IWO algorithms.

**Figure 4.** Location of the case study (Aswan) on the world map.

**Figure 5.** Annual load curve over a time interval of one hour with a peak demand of 70 kW.

**Figure 6.** Solar irradiation over the studied region.

**Figure 7.** Temperature variation in Aswan.

**Figure 8.** Annual variation of wind speed over the year in Aswan.

**Figure 9.** Atmospheric pressure variation in Aswan.

#### **6. Results**

The main object of this research paper is to find the optimal design of the proposed hybrid system and to validate the accuracy of the proposed QOBO optimization method. The optimal sizing is based on the objective functions introduced in (9) and the parameters of optimization are: (i) the area of PV system, (ii) the area swept by the WT, (iii) the rated power of diesel generator, (iv) the nominal capacity of the battery, (v) the consumption of the biomass fuel. To confirm the suitability of the QOBO in addressing such optimization problem, QOBO, BO, HHO, AEFA and IWO were launched 100 times for each configuration and statistical study was conducted based on the best minimum value of the fitness function. For a deep analysis of the obtained results and to ensure the sensitivity analysis, four indices were chosen, namely, NPC, LCOE, LPSP and the availability index. In the next subsections, the optimization results are provided for the standalone system with multiple scenarios. Modelling and simulation of the optimization problem were accomplished using MATLAB 2015a program, while the adjusting parameters for the three algorithms are the same, i.e., the number of maximum iterations is taken as 100 iterations and the search agents' number is 30 agents. The input technical and economic data for the system components are presented in Table 1. The results of the statistical measurements for the modified QOBO and the conventional BO with HHO, AEFA and IWO algorithms are listed in Tables 2 and 3. From the previously mentioned tables, the reader can conclude that the QOBO technique generates the best minimum value of the fitness function in all cases. The convergence curves of the 100 iterations implemented for all the studied configurations using QOBO, BO, HHO, AEFA and IWO are presented in Figure 10a–d.

#### **Table 1.** Units for magnetic properties.

#### *6.1. Validation of QOBO Algorithm*

The results of the statistical measurements for the modified QOBO, the conventional BO, HHO, AEFA and IWO algorithms are listed in Tables 2 and 3. From Table 2, the reader can find the results of the optimal sizing for the multiple scenarios studied, as well as the convergence time of each simulation, and conclude that the QOBO algorithm finds the best results with a short time compared with the other algorithms. From Table 3, the reader can compare between algorithms and the different scenarios of the proposed hybrid system using multiple factors. Briefly, it is noticed that the hybrid PV/biomass system is highly competitive, mainly using the developed QOBO algorithm, the optimized system is calculated with \$110,807, which means an LCOE of 0.1053 \$/kWh, the constraints are satisfied and the project is 100% supplied by renewable energy sources. In this scenario, the performances of the QOBO and the BO are almost equal, while in other scenarios, the difference is clearly noticed.

**Table 2.** Sizing results of different scenarios obtained from different optimization methods.

#### **Table 3.** Factor results for all scenarios.

#### *6.2. Combinations of the Studied System Components*

In this section, the results obtained in the convergence simulation of the NPC as a fitness function using the QOBO are presented. The optimized parameter results (i.e., *Apv*, *Awind*, *Pdgn*, *PCap*\_*bat*, *PBM*) for all suggested combinations are listed in Table 3 with the rating of the inverter that takes the value of the peak load demand. From Figure 10, the reader can notice that using QOBO, BO, HHO, AEFA and IWO algorithms, the best minimum values of fitness function (NPC) is obtained for the second configuration, i.e., hybrid PV/biomass energy system. From the table, it is obvious that QOBO generates the minimum value of LCOE in all cases.

The reliability of the proposed scenarios of the proposed HRES are respected and the availability of power is highly assured, the penetration RES is considered in this paper, while different results are

obtained. The minimum penetration of 70% is obtained for the WT/Diesel/Battery scenario while the maximum penetration of 99.75% is obtained for the PV/WT/Diesel/BESS scenario. The daily battery autonomy is also influenced by the configuration of the HRES, the best autonomy is achieved for the WT/Diesel/BESS scenario taking nearly 4 days, while the minimum autonomy is obtained in PV/WT/Diesel/BESS case with only 6 h. The last system is composed of the different energy resource which explains the independence for a specific resource. Table 4 presents a detailed overview of all costs needed, for all scenarios presented and for all proposed algorithms.

**Figure 10.** *Cont.*

**Figure 10.** Convergence of the objective function of all algorithms for different scenarios; (**a**) PV/WT/Diesel/BESS, (**b**) PV/Biomass, (**c**) PV/Diesel/BESS, (**d**) WT/Diesel/BESS.

#### *Appl. Sci.* **2020** , *10*, 6604

#### *6.3. Sensitivity Analysis*

RES is intermittent which can be affected by any variation of sizing, meteorological or economic data. The sensitivity analysis is a method that helps to select and/or to expect the optimal configuration of the hybrid system. The sensitivity analysis in this paper is implemented on the best scenario of the proposed, i.e., the PV/Biomass in the Aswan region. The selection of the sensitivity variables is based on the sizing of components in order to analyze the effect of sizing variation on four factors which are NPC, LCOE, LPSP and the Availability index.

Figure 11 shows the effect of variation in the sizing of PV and biomass units on the NPC. The PV sizing is highly impacted by the total cost of the hybrid PV/Biomass system, which means that in the case of reducing the area of PV units the NPC is reduced too. On the other hand, if the area covered by PV modules is increased, the NPC increases too. The variation in the sizing of biomass unit is increased throughout the interval −20 to 20 slowly and it has no noticeable impact on the NPC anyway. Figure 12 shows the effect of variation of PV and biomass sizing on the LCOE. The NPC and the LCOE are linked with a linear equation which means that they have the same shape. The LCOE reached 0.08 \$/kWh when the area of the PV system is reduced by 20%. Figure 13 shows the impact of variation in the sizing of PV and biomass systems on the LPSP. The impact of PV size is very important for the LPSP, because when the size of the PV system is increased the LPSP is enhanced, mainly in the −20% to 0% interval. When the PV size is changed in the interval of 0% to +20%, the LPSP is increased to 2% while when the PV size is changed to −20%, the change in LPSP equals 16.4% which is a very bad sign for system building. The Biomass system does not affect the value of the LPSP and the transition between −4% to 0% is explained as the obtained sizing of the system is optimum. Figure 14 shows the impact of the variation of PV and Biomass sizing on the availability index. The availability index enhanced exponentially with the increase in the PV sizing. In the interval between −20% and 0, availability progresses quickly, while after zero, the availability begins to be stabilized and it is clearly shown in the interval between +12% and +20%.

**Figure 11.** Sensitivity analysis application for net present cost (NPC).

**Figure 12.** Sensitivity analysis application for Levelized Cost of Energy (LCOE).

**Figure 13.** Sensitivity analysis application for Loss of Power Supply Probability (LPSP).

**Figure 14.** Sensitivity analysis application for the Availability index.

The PV system through this analysis is demonstrated as a very important element for having a high hybrid systems criterion. However, the Biomass system helps the PV to satisfy the constraints and its variation does not have a serious impact on the performance of the hybrid system.

#### **7. Conclusions**

With the increased penetration level of RES into electrical energy production in the microgrid systems, new challenges have emerged on the international scene. These challenges are represented in finding ways to optimize the design of the hybrid system by using smart algorithms and software. Among these dilemmas, the economic cost and feasibility of installing systems in different locations in the world is considered the most important challenge. Therefore, this research proposes a developed algorithm called Quasi-Oppositional Bonobo Optimizer (QOBO) for the optimal economic design of a stand-alone hybrid microgrid system in Aswan, Egypt. Four configurations of the hybrid system have been implemented, which consist of RES (PV panels, WT and biomass) with diesel generators and battery storage systems. The obtained results showed that the PV/Biomass scenario is the most cost-effective system with an NPC of \$110,807 and LCOE of 0.1053 \$/kWh; otherwise, the best configuration of the microgrid system contained 293.971 m2 of PV and 1020.18 ton/year consumed by the biomass system; the PV/Diesel/BESS scenario is also cost-effective with NPC of \$153,401 and LCOE of 0.1457 \$/kWh. On the other side, the LPSP and availability index are satisfied and without the need for traditional resources. Additionally, the results showed the ability of the QOBO algorithm to reach the optimal solution in a shorter time and with better efficiency compared to the traditional BO, HHO, AEFA and IWO algorithms in all cases studies. Furthermore, a sensitivity analysis of the proposed systems scenarios was performed to obtain the impact of unit size on the performance of the hybrid system, where it has been emphasized that PV system sizing is very important and has a great impact on the overall performance of the system. The obtained results from this study would be useful material for decision makers working on the development of the renewable energy sector in Egypt. In future studies, it is suggested to apply the proposed QOBO in other engineering problems.

**Author Contributions:** Conceptualization, S.K. and M.K.; Data curation, O.H.M., H.M.S. and A.S.; Formal analysis, F.J., M.A. and H.M.S.; Methodology, M.K., A.S. and S.K.; Resources, O.H.M., S.K. and H.M.S.; Software, A.S., M.K.; Supervision, F.J.; Validation, H.M.S. and O.H.M.; Visualization, S.K., A.S. and M.A.; Writing—original draft, M.K., O.H.M., A.S. and H.M.S.; Writing—review & editing, S.K., M.A. and F.J. All authors together organized and refined the manuscript in the present form. All authors have approved the final version of the submitted paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the NSFC, China-ASRT, Egypt, Joint Research Fund, under Grant<PHONE_NUMBER>6.

**Acknowledgments:** The authors gratefully acknowledge the contribution of the NSFC (China)-ASRT (Egypt) Joint Research Fund, Project No.<PHONE_NUMBER>6 for providing partial research funding to the work reported in this research.

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

#### **Nomenclature**

#### **Appendix A. Algorithms**

#### *Appendix A.1. Harris Hawks Optimization Algorithm*

Heidari and et al. [36] proposed a new nature-inspired optimization algorithm called Harris Hawks Optimizer. They were inspired by the cooperative behavior and chasing style of Harris hawks. The modeling of this technique is based firstly on an exploration phase; afterwards, the transition from exploration to exploitation, then the

exploitation phase and, finally, the soft besiege. The modeling is taken on for all strategies for exploring a prey, surprise pounce and different attacking methods of Harris hawks. The pseudo-code of the HHO algorithm is proposed below.

**Algorithm A1**: Pseudo code of HHO

Initialize the population size and max iteration (*Kmax*) Initialize a set random rabbit location, within the limits *Xi min* ≤ *Xi rabbit* ≤ *Xi max*. Evaluate the objective function for all rabbits While (k < *Kmax*) Calculate the fitness of hawks Set *xrabbit* in the best location for each hawk do Update the initial energy *E*0, energy E and jump strength J; *E*0 = 2rand () − 1, *E* = 2*E*0 1 − *t T* , J = 2(1 − rand ()) if (|E| ≥ 1) then Exploration phase if (|E| < 1) then Exploitation phase if (r ≥ 0.5 and |E| ≥ 0.5) then Soft besiege else if (r ≥ 0.5 and |E| < 0.5) then Hard besiege else if (r < 0.5 and |E| ≥ 0.5) then Soft besiege with progressive rapid dives else if (r < 0.5 and |E| < 0.5) then Hard besiege with progressive rapid dives Return *xrabbit*

*Appendix A.2. Artificial Electric Field Algorithm*

Anita and Yadav [37] were inspired by Coulomb's law of electrostatic force to create a novel artificial electric field algorithm. The concepts of electric field and charged particles provide us a strong theory for the working force of attraction or repulsion between two charged particles. The pseudo code of the AEFA algorithm is proposed in Algorithm A2.

**Algorithm A2**: Pseudo code of AEFA

Initialize a set of random population *Xi B* = *X*1 *B*, *X*2 *B*, ... , *XN B* of N size, within the limits *Xi min* ≤ *Xi B* ≤ *Xi max*. Initialize the velocity to a random value Evaluate the fitness of whole population Set the iteration to zero Reproduction and Updating While criteria not satisfied do Calculate K (t), best (t) and worst (t) for i = 1: N do Evaluate the fitness values Calculate the total force in each direction Calculate the acceleration V*i* (t + 1) = rand () × V*i* (t) + a*i* (t) X*i* (t + 1) = X*i* (t) + V*i* (t + 1) end for end while

Appendix A.2.1. Invasive Weed Optimization Algorithm

Invasive weed optimization is a numerical stochastic optimization algorithm inspired by colonizing weeds, which was introduced in 2006 by Mehrabian and Lucas [38]. In IWO, a certain number of weeds make up the

whole population, and each weed comprises a set of decision variables. Weeds are a serious threat to desirable plants because they are plants that are invasive and hardy.

Weeds are plants which are vigorous and invasive; they pose a serious threat to desirable, cultivated plants in agriculture. Weeds have shown to be very robust and adaptive to change in the environment. The IWO optimization algorithm has been modeled based on four steps: initialization, reproduction, spatial dispersal and competitive exclusion.

#### • *Initialization and Production*

Firstly, the population is spread over the research space randomly; afterwards, each plant is allowed to produce seeds depending on its own fitness; the production of seeds is not only allowed for the better plants' fitness as in the other evolutionary algorithms, but the reproduction step of IWO is also proposed to give a chance to infeasible individuals to survive and reproduce similar to the mechanism which occurs in nature. The weeds producing seeds can be formulated as follows:

$$\text{Weed}\_n = \frac{f - f\_{\text{min}}}{f\_{\text{max}} - f\_{\text{min}}} (s\_{\text{max}} - s\_{\text{min}}) + s\_{\text{min}} \tag{A1}$$

where in each iteration, *f* is the current weed's fitness. *fmax* and *fmin* represent the max and min fitness values, respectively. *smax* and *smin* represent the max and min values of the weeds, respectively.

• *Spatial Dispersal*

The generated seeds are being randomly distributed over the search space such that they abode near the parent plant. However, the standard deviation (σ) of the random function will be reduced in every iteration, the nonlinear alteration equation of the standard deviation at each iteration is presented as follows:

$$
\sigma\_{intcr} = \frac{(iter\_{max} - iter)^n}{\left(iter\_{max}\right)^n} (\sigma\_{initial} - \sigma\_{final}) + \sigma\_{final} \tag{A2}
$$

where *itermax* is the maximum iteration, *n* is the nonlinear modulation index, σ*initial* and σ*final* are the initial and final values of standard deviation, respectively.

• *Competitive Exclusion*

In a colony, the maximum number allowed of plants is limited; for that, competitive exclusion is applied. The plant that leave no offspring would go extinct; otherwise, they can survive. After some iterations, the number of plants in a colony will reach its maximum through the reproduction step, the seeds and their parents are ranked together, and all plants in the research space are considered as weeds; afterwards, weeds with lower fitness are eliminated.

The overall steps of the IWO algorithm are presented in Algorithm A3.

**Algorithm A3:** Pseudo code of IWO

Initialize a set of random weeds, *weedi B* = *weed*1 *B*, *weed*2 *B*, ... , *weedN B* within the limits *weedi min* ≤ *weedi B* ≤ *weedi max*. Set the IWO's parameters Evaluate the objective function for all weeds While (*iter* < *itermax*) Calculate the best and worst fitness in the colony Calculate the σ for each weed in the colony Calculate the number of seeds following the fitness of each weed Add the seeds to their parents in the colony if *Sizemax* ≤ *Nbpopulation* Sort the new population according to their fitness Eliminate the worst fitness in order to achieve the *Sizemax* allowed end if end for Update iteration *iter* = *iter* + 1 end while Return the final best solution

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Economic Optimal Scheduling of Wind–Photovoltaic-Storage with Electric Vehicle Microgrid Based on Quantum Mayfly Algorithm**

**Ximu Liu, Mi Zhao \*, Zihan Wei and Min Lu**

College of Mechanical and Electrical Engineering, Shihezi University, Shihezi 832003, China **\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** The effectiveness of energy management systems is a great concern for wind–photovoltaicstorage electric vehicle systems, which coordinate operation optimization and flexible scheduling with the power grid. In order to save system operation cost and reduce the energy waste caused by wind and light abandonment, a time-sharing scheduling strategy based on the state of charge (SOC) and flexible equipment is proposed, and a quantum mayfly algorithm (QMA) is innovatively designed to implement the strategy. Firstly, a scheduling strategy is produced according to the SOC of the battery and electric vehicle (EV), as well as the output power of wind–photovoltaic generation. In addition, the minimum objective function of the comprehensive operation cost is established by considering the cost of each unit's operation and electricity market sale price. Secondly, QMA is creatively developed, including its optimization rule, whose performance evaluation is further carried out by comparisons with other typical bionics algorithms. The advantages of QMA in solving the low-power multivariable functions established in this paper are verified in the optimization results. Finally, using the empirical value of the power generation and loads collected in enterprise as the initial data, the mayfly algorithm (MA) and QMA are executed in MATLAB to solve the objective function. The scheduling results show that the time-sharing scheduling strategy can reduce the system's cost by 60%, and the method decreases energy waste compared with ordinary scheduling methods, especially when using QMA to solve the function.

**Keywords:** microgrid; economic scheduling; clean energy; quantum mayfly algorithm (QMA)

#### **1. Introduction**

With the rapid development of the national economy, demand for fossil fuels is increasing. However, the output of traditional energy is limited, and the utilization of renewable energy has become the general trend. As new energy sources, wind and light are widely distributed and can be permanently used. They are characterized by low cost and strong environmental protection. Nevertheless, there are considerable random fluctuations of output power due to natural conditions, which affect the security and stability of the power system after grid connection [1–3]. Therefore, designing a stable and low-cost scheduling strategy based on the structure of the wind–photovoltaic-storage electric vehicle complementary power generation system has become a key issue. By coordinating wind–photovoltaic and flexible devices, we can ensure the reliability of the system, improve the economy, and increase environmental friendliness [4,5].

In the process of modeling, parameter uncertainty becomes a significant problem that threatens system security with the increasing scale of the microgrid. In order to solve this problem, robustness becomes one of the parameters worthy of study [6,7]. The improvement of robustness can make the system maintain its balance when it is disturbed, thereby reducing the loss caused by disturbances [8]. Meanwhile, some other power consumption units are introduced, such as electric vehicles (EVs) and fuel cells, on the basis of wind–photovoltaic-storage microgrid architecture, which can ensure the system maintains

**Citation:** Liu, X.; Zhao, M.; Wei, Z.; Lu, M. Economic Optimal Scheduling of Wind–Photovoltaic-Storage for Electric Vehicle Microgrid Based on Quantum Mayfly Algorithm. *Appl. Sci.* **2022**, *12*, 8778. https://doi.org/ 10.3390/app12178778

Academic Editors: S. M. Muyeen and Mohamed Benbouzid

Received: 7 August 2022 Accepted: 28 August 2022 Published: 31 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

a stable operation by interacting power with other components when it is impacted [9,10]. In terms of the objective function, some microgrid systems, including fossil fuels, need to be considered regarding carbon emissions to reduce the environmental pollution caused by them [11]. Therefore, multi-objective functions combined with economic objectives and carbon emissions are established in recent studies. Furthermore, in order to ensure the normal operation of the microgrid system, it is necessary to ensure power conservation is met, and power flow constraints, including power flow calculation, should be added to the scheduling study [12,13]. Additionally, equipment parameters are important constraints to enable the system to operate under the superior performance of the equipment, thereby ensuring equipment safety and power quality.

In order to realize the scheduling strategy and accurately find the optimal value of the objective function, a variety of high-performance intelligent algorithm schemes based on biological habits are proposed, and verified in specific cases [14–16]. Bionics algorithms, such as the flying sparrow search algorithm, flower pollination algorithm, and other heuristic algorithms, are proposed and applied to practical problems [17–19]. Among them, as a new bionics algorithm, the mayfly algorithm (MA) has better optimization performance than the traditional intelligent algorithms. The integration of the advantages of various algorithms has significantly improved the optimization process in terms of speed and accuracy [20]. However, when it comes to high-dimensional complex problems, it is still difficult to jump out of local optimal regions simply by relying on its own mechanism, and the convergence accuracy of the algorithm is not high enough because the step size of the mayfly is too short when it moves. Based on the above shortcomings, logistic mapping is used to improve the optimization ability of mayfly individuals [21,22], or the learning factor is changed to improve the degree of self-adaptation [23].

In order to describe the research status of microgrid scheduling and its characteristics more clearly, the studies conducted for our paper, and other related studies, are summarized in Table 1, which is represented as follows:

**Table 1.** Comparison of relevant studies.

Through the above table, it can be seen that the robustness and stability of the system are studied by some scholars. In these studies, robust model predictive control (RMPC) and other methods are used to deal with the uncertainties of renewable energy. In the microgrid system containing fossil energy, reducing carbon emissions is also an important issue that needs to be solved to improve the environmental protections of the system. More studies focus on the microgrid system containing EVs, namely the EVs charging and discharging states, as well as the stability of the EVs themselves.However, it is also worth studying how to interact with the microgrid while ensuring its stable operation. Meanwhile, some studies are aimed at improved algorithms, which have been greatly modified in optimization rate, convergence, and escape from local deadlock; however, there is still room for further improvement. Therefore, the problems that are concentrated on in this paper are represented as follow:

(1) The EV is taken as one of the dispatching objects, and its interaction state with the microgrid is judged according to its SOC. The microgrid system is able to make full use of power in this measure.

(2) The QMA algorithm raised in this paper is compared with the moth–flame optimization (MFO), grey wolf optimize (GWO), particle swarm optimization (PSO), whale optimization algorithm (WOA), sine cosine algorithm (SCA), and MA to verify its advantages, and is applied to the solution of the objective function.

The rest of the paper is organized as follows: The structure model of the microgrid is established in Section 2. Meanwhile, the scheduling strategy of each unit, as well as the objective functions with constraint conditions of economic dispatching regarding the microgrid, are established. In Section 3, the operation rules of the QMA are analyzed, and the comparison between various intelligent algorithms is realized. The algorithms are combined with the objective function in Section 4, and the economic scheduling of the microgrid in a specific area is carried out according to the actual situation. Section 5 summarizes the content and contributions of this paper, and puts forward prospects for future studies.

#### **2. Microgrid Structure and Dispatching Strategy**

#### *2.1. Microgrid Structure*

The model of the wind–photovoltaic-storage electric vehicle microgrid system is shown in Figure 1. Power output can be achieved by wind turbines and photovoltaic generations, while the load can only consume power. As flexible equipment, the consumption and production of electric energy can be realized by EVs and batteries. Meanwhile, energy interaction can be implemented by the accession of the microgrid to the main grid .

**Figure 1.** The model of wind–solar-storage electric vehicle microgrid system.

#### *2.2. Scheduling Strategy Research*

In this paper, in order to minimize the overall cost by dispatching each generation unit, the following requirements should be met in the scheduling strategy to ensure power supply reliability:

In order to achieve the above requirements, the dispatching strategy aiming at the system economy is designed as follow: The output of wind power and photovoltaic generations are coordinated, and participate in the dispatching according to the state of the battery and EV on charge. Finally the power balance is realized through the interaction process with the main grid.

a. Wind and photovoltaic generation

The empirical output value of wind and photovoltaic generation is taken as the actual output power, and the scheduling is carried out after increasing the cost of wind and solar abandoning in order to reduce energy waste. The output power of wind and photovoltaic generation after dispatching is compared with the empirical value. If the output power is lower than the empirical value, the output power after dispatching can be used; otherwise, the empirical value has to be maintained.

b. Battery generation

Since photovoltaic and wind generation is affected by daily environmental changes, the charging and discharging process of batteries can meet the demand. According to power load trends within 24 h, power supply states are divided into peak, normal, and valley periods. Charging and discharging strategies are adopted based on the current period and the battery state of charge (SOC). The specific process is shown in Figure 2.

**Figure 2.** The scheduling strategy of battery generation.

According to Figure 2, charging and discharging can be carried out only if the battery SOC is within a certain range. When the charging and discharging standard is reached by the battery, the interactive state is judged according to its SOC value and the time period. Thereby, the charging and discharging power is calculated. The specific calculation method is provided as follow:

$$P\_{\text{loss}}(t) = SO\mathbb{C}\_x(t)\mathbb{C}\_x \tag{1}$$

where the charging and discharging power of battery at time *T* is described as *Pbess*(*t*). The SOC of the battery at time *T* is shown as *SOCx*(*t*), which is the scheduling target of the battery. The rated capacity of the battery collected from empirical data is defined as *Cx*.

c. Electric vehicle

The output of an EV can be judged according to its SOC. If the SOC can maintain the normal operation of the vehicle, the power is outputted to the system; otherwise, the system's power is consumed. The flow is expressed in Figure 3.

**Figure 3.** The scheduling strategy of EV.

Similar to battery scheduling, the output state of EVs can be judged according to their SOC. The output of the SOC after scheduling, according to the idea shown in Figure 3, can be used for further calculations, and the specific expression is represented as follows:

$$P\_{cv}(t) = SOC\_c(t)C\_{c^\*} \tag{2}$$

where interactive power between EVs and the microgrid at time *T* is given as *Pev*(*t*). The SOC of the EV at time *T* is shown as *SOCe*(*t*), and the fixed capacity of EV is expressed as *Ce*.

#### *2.3. Expression Construction*

#### 2.3.1. Objective Function

In this paper, objective functions are established for the economy of the wind–photovoltaicstorage electric vehicle microgrid system under a grid-connected operation. Consider the operation, maintenance, and investment cost of each unit in a system, and the cost of grid sales and purchases. Meanwhile, taking into account the operating costs of EVs and the cost policies issued by the government for new energy, the optimization goal can be achieved through relevant scheduling. The expression is as follows:

$$Minimize\ F(\mathbf{x}) = F\_{wv} + F\_\mathbf{x} - F\_\mathbf{G} - F\_{EV} + F\_{BT} \tag{3}$$

where the operation and maintenance costs of wind and photovoltaic generation are represented as *Fwv*. The power generation and investment costs of the battery are determined as *Fx*. The power grid sales and purchase costs are simplified as *FG*. The operating costs of EVs and policy expenses related to new energy power generation are given as *FEV* and *FBT*. The decision variable of the objective function is the output power of each unit, which is the real number. Moreover, the specific coupling relationships between decision variables and objective functions (3) are shown in subsections a–e.

a. Wind and photovoltaic generation

The equipment loss, manual repair, and operating reserve storage costs are included in *Fwv*, which can be expressed as follows:

$$F\_{uvv} = \sum\_{t=1}^{T} \mathbb{C}\_{\text{v}} P\_{wf}(t) + \sum\_{t=1}^{T} \mathbb{C}\_{\text{v}} P\_{vf}(t) + F\_{RE} \tag{4}$$

where the operation and maintenance cost coefficients of wind and photovoltaic generation are defined as *Cw* and *Cv*, respectively. These coefficients are based on the ratio of the total operation and investment costs of the equipment to the total annual generated power within a year obtained in enterprise. The output power of wind and photovoltaic power generation is computed as *Pw f* and *Pv f* , respectively. The operating reserve storage cost is simplified as *FRE*, which can be expressed as follows:

$$F\_{RE} = \mathbb{C}\_{RE}[P\_{vf}(t) + P\_{wf}(t)]\tag{5}$$

where the operating reserve storage cost coefficients of wind and photovoltaic generation are given as *CRE*.

b. Battery generation

In Equation (6), *Fx* consists two parts, which are the acquisition and maintenance costs of batteries. The acquisition cost is converted into each work process, and the specific expression is calculated as follows:

$$F\_{\rm x} = \frac{F\_{\rm cap}C}{P\_{\rm bess}T\_{\rm n}} |\sum\_{t=1}^{T} P\_{\rm bess}(t)| + \frac{K\_{o}P\_{\rm bess,max}}{365} + K\_{M}\frac{1}{\Delta t}P\_{\rm bess}(t)\Delta t\tag{6}$$

where the acquisition cost and depreciation coefficient are defined as *Fcap* and *C*, respectively. The annual operating cost coefficient and hours are represented as *Ko* and *Tn*, respectively. *Pbess* and *Pbess*,*max* are battery-rated power and maximum charge and discharge power, respectively. Meanwhile, the maintenance cost coefficient is described as *KM*, and the duration from the beginning to the end of charging and discharging the battery is shown as *T*.

c. Interaction with main grid

The microgrid system can interact with the main grid when the power is vacant or redundant. When the interactive power is positive, the power is purchased from the main grid; otherwise, the interaction cost *FG* is the sum of the two. The equation is simplified as follows:

$$F\_G = \sum\_{t=1}^{24} P\_{grid}(t)C\_{price,grid} \tag{7}$$

where *Pgrid*(*t*) refers to interactive power with the main grid, and the time-sharing sale and purchase price of the main grid is expressed as *Cprice*,*grid*.

d. Policy of government

Subsidies are paid for wind and photovoltaic generation, which encourage companies to adopt new sources of electricity. Meanwhile, fees are charged for energy waste to reduce the energy loss caused by wind and solar discarding. The formula is computed as follows:

$$F\_{BT} = \sum\_{t=1}^{T} \mathbb{C}\_{q} \left[ P\_{\text{wy}}(t) + P\_{\text{vy}}(t) - P\_{\text{wy}}(t) - P\_{\text{vy}}(t) \right] - S\_{bt} \tag{8}$$

where the wind and solar discarding coefficient is given as *Cq* and *Sbt*, which represents the government subsidy. The actual utilization power of wind and photovoltaic generation is determined as *Pwy*(*t*) and *Pvy*(*t*), respectively.

e. Electric vehicles

As a client device, the cost of selling and purchasing power to the microgrid is mainly included in the cost of EVs. The specific expression is given as follows:

$$F\_{EV} = \sum\_{t=1}^{24} P\_{EV}(t)\mathbb{C}\_{price, EV} \tag{9}$$

where the electricity purchase price of an EV is defined as *Cprice*,*EV*.

#### 2.3.2. Constraint Condition

In order to ensure the safe and stable operation of the grid system, the power balance should be regarded as the basic constraint condition, and the output power and SOC of flexible devices should be limited to ensure the normal operation of equipment.

a. Power balance

In any period of time, all output, load, flexible device charging and discharging, and interaction powers between the microgrid system and the main grid should be kept in balance, and the constraint equation is given as follows:

$$P\_{wy}(t) + P\_{vy}(t) + P\_{grid}(t) + P\_x(t) + P\_{EV}(t) = P\_{load}(t) \tag{10}$$

where *Px* is the charging and discharging power of the battery, and the total load is expressed as *Pload*.

b. Wind and photovoltaic generation

Due to the limitations of the equipment parameters, the output power of wind and photovoltaic generation is within a certain range. Meanwhile, in order to ensure the normal operation of the system when the power required surges, it is necessary to add operating reserves. In detail, 30% of power generated by wind and photovoltaic power generation is stored, and the rest is the maximum power supplied to the system. Furthermore, the change of power in each time period should also be constrained in a certain range. Therefore, the relevant constraints are set as follows:

$$\begin{cases} \begin{aligned} 0 \le P\_{wy}(t) \le 0.7P\_{wf}(t) \\ 0 \le P\_{vy}(t) \le 0.7P\_{vf}(t) \\ R\_{w - down}\delta t \le \delta P\_w \le R\_{w - up}\delta t \\ R\_{v - down}\delta t \le \delta P\_v \le R\_{v - up}\delta t \end{aligned} \end{cases} \tag{11}$$

where the adjustable power amplitude of wind and photovoltaic generation in a certain period is represented as *δPw* and *δPv*, respectively. The downward climbing rate of wind and photovoltaic generation during adjustment is expressed as *Rw*−*down* and *Rv*−*down*, respectively. Meanwhile, *Rw*−*up* and *Rv*−*up* give the upward climbing rate for wind and photovoltaic generation during adjustment, respectively. The adjustment time is represented by *δt*.

c. Flexible generation

Batteries are not only constrained by power, their safety performance and service life are affected by their SOC , and the operation of an EV is also affected by its SOC. Specific conditions are simplified as follows:

$$\begin{cases} SOC\_{\textit{cmin}} \le SOC\_{\epsilon}(t) \le SOC\_{\textit{cmax}}\\ SOC\_{\textit{xmin}} \le SOC\_{\pi}(t) \le SOC\_{\textit{xmax}}\\ \quad 0 \le P\_{\textit{bess},t}(t) \le P\_{\textit{bess},\textit{max}} \end{cases} \tag{12}$$

*SOCemin*, *SOCxmin*, *SOCemax*, and *SOCxmax* are the upper and lower limits of their SOC.

In addition, to ensure the safe and stable operation of batteries, the charging and discharging rates should be controlled within a certain range. The specific constraints are shown as follows:

$$
\delta SOC\_x \le 0.2 \tag{13}
$$

where the adjustable range of a battery's SOC in each period is expressed as *δSOCx*.

#### **3. Improved Mayfly Algorithm**

#### *3.1. Traditional Mayfly Algorithm*

The MA is a bionics algorithm derived from the social behavior of the mayfly. Inspired by the movement mode and reproduction process of female and male populations, the optimal and suboptimal individuals in each population are selected. Meanwhile, the optimal offspring generation is obtained through mating between the optimal male and female individuals, and the suboptimal offspring generation is obtained in the same way. The direction of movement of each mayfly is influenced by the dynamics of individual and collective optimal positions, and female mayflies target male mayflies towards their positions.

a. The movement of male mayfly

The flight mode of male mayflies is similar to the movement mode of birds in a particle swarm algorithm, and the direction and distance of male mayflies are adjusted according to their own flight experiences and that of individuals around them. The specific method is shown as follows:

$$\mathbf{x}\_{i}^{u+1} = \mathbf{x}\_{i}^{u} + \boldsymbol{\upsilon}\_{i}^{u+1} \tag{14}$$

where *xn i* and *vn i* are the current position and speed, respectively, of the male mayfly *i* on the *n*th search. The values are calculated as follows:

$$\boldsymbol{\upsilon}\_{ij}^{n+1} = \boldsymbol{\upsilon}\_{ij}^{n} + \boldsymbol{a}\_{1}\boldsymbol{e}^{-\beta\hat{\boldsymbol{\alpha}}\_{p}^{2}}(\boldsymbol{p}\_{best\_{i}^{j}} - \boldsymbol{X}\_{ij}^{n}) + \boldsymbol{a}\_{2}\boldsymbol{e}^{-\beta\hat{\boldsymbol{\alpha}}\_{p}^{2}}(\boldsymbol{g}\_{best\_{i}^{j}} - \boldsymbol{X}\_{ij}^{n}) \tag{15}$$

Because male mayflies perform a dance on the surface of water to attract females, the position of the male mayflies is constantly changing, meaning they do not build a high speed. *vn ij* is the speed of the *n*th search of the mayfly *i* at *j* dimension, and *xn ij* is the position at that time. *α*1 and *α*2 are estimated from the positive attraction coefficients of social interaction, and *β* is the visibility coefficient of the mayfly. Meanwhile, the optimal locations of the individual and collective mayflies are expressed as *pbestj i* and *gbestj i* , respectively. Additionally, the distances from current position to *pbestj i* and *gbestj i* are defined as *lp* and *lg*, respectively, and are calculated as follows:

$$\|\mathbf{x}\_i - \mathbf{X}\_i\| = \sqrt[2]{\sum\_{j=1}^n (\mathbf{x}\_{ij} - \mathbf{X}\_{ij})^2} \tag{16}$$

For the best mayfly in the population, a fixed dance pattern is needed to be maintained. Meanwhile, a random element is introduced to ensure that the speed is constantly changing. The calculation is described in this case as follows:

$$
\sigma\_{ij}^{n+1} = \sigma\_{ij}^{n} + d \times r \tag{17}
$$

where the dance coefficient is simplified as *d*, and *r* is the random natural number within [−1, 1].

b. The movement of female mayfly

The movement of the female mayfly depends on the attraction of the male mayfly, and their position renewal depends on the increase of speed, which can be expressed as follows:

$$y\_i^{n+1} = y\_i^n + \upsilon\_i^{n+1} \tag{18}$$

The speed update is a certain process, which means that in order to ensure the quality of offspring, the best female needs to be attracted by the best male, the second-best female by the second-best male, and so on. The speed update expression is given as follows:

$$w\_{ij}^{n+1} = \begin{cases} \left. \upsilon\_{ij}^{\text{n}} + \alpha\_2 e^{-\beta l\_f^2} (x\_{ij}^{\text{n}} - y\_{ij}^{\text{n}}) \right. \\ \left. \upsilon\_{ij}^{\text{n}} + \text{g} \times r \end{cases} \tag{19}$$

where the position of the female mayfly is expressed as *yn ij*, the random walk coefficient of the female mayfly is represented as *g*, and *lf* is determined by the distance between the male and female mayflies.

#### c. The mating of male and female mayflies

In the process of parent mating, the optimal and suboptimal individuals in the male and female populations should be selected for mating and reproduction based on their fitness functions. The results of interbreeding, which produces the optimal and suboptimal offspring, are calculated as follows:

$$\begin{cases} \
of fsring1 = L \times m + (1 - L) \times f\_m \\
of fsring2 = L \times f\_m + (1 - L) \times m \end{cases} \tag{20}$$

where the male and female in the parent generation are represented as *m* and *fm*, respectively, and *L* is a random natural number within a specific range.

#### *3.2. Ma with Quantum Idea*

The traditional MA can find the optimal value in a single-peak function accurately by taking advantage of the characteristics in mayfly reproduction. However, with regard to the large population and complicated process, the search speed is slow and the convergence is not fine. Meanwhile, it is easy to fall into local deadlock when dealing with multi-peak functions. Therefore, the quantum idea is introduced on the basis of a traditional MA, thereby forming the QMA. Because the position and velocity of the mayfly cannot be determined simultaneously in quantum space, the wave function is used to represent the position of the mayfly, and the Monte Carlo method is used to solve the problem. The particle update expression is shown as follows:

$$\begin{cases} m\_{\text{best}^n} = \frac{1}{N} \sum\_{i=1}^q P\_{\text{best}^n\_i}(i = 1 \dots n) \\\ P\_i^t = \gamma \times P\_{\text{best}^n\_i} + (1 - \gamma) \times g\_{\text{best}}^u \\\ X\_i^{n+1} = P\_i^n \pm \varepsilon |m\_{\text{best}^n} - x\_i^n| \log \frac{1}{a} \end{cases} \tag{21}$$

where the numbers of individuals and iterations are represented as *N* and *n*, respectively. *mbestn* is the average historical optimal position of the male mayfly, and *Pn i* is obtained from the updated position of the *i*th male mayfly in the *n*th iteration. *r* and *a* are uniformly distributed values within (0,1), and *c* is the final random motion parameter.

The execution steps of the QMA are given as follows:

Step 1: Initialize the positions of the female and male mayflies in the space.

Step 2: Calculate the average optimal location *mbest* of the male mayflies according to the first equation of (21).

Step 3: Calculate the fitness value and sorting according to Formula (3), and compare with the previous iterative value. If the current function value is less than the previous iteration, the current mayfly position is updated for the individual optimal position, otherwise it keeps the previous iteration. Hence, the optimal male individual *pbest* and collective locations *gbest* are obtained.

Step 4: Calculate the new positions of the female and male mayflies according to Formula (18) and the second equation of (21), respectively, and mate in sequence.

Step 5: Calculate the fitness function and update *pbest* and *gbest*.

Step 6: Repeat Step 2 to 5 until the stop condition is met.

The QMA flow chart is shown in Figure 4.

**Figure 4.** The flow chart of QMA.

#### *3.3. Performance Analysis of Qma*

In order to verify the superiority of the QMA, the single-peak and multi-peak functions are tested by MFO, GWO, PSO, WOA, SCA, MA, and QMA in this section. We combine the biological habits and existing literature on each algorithm, and the initialization parameters and test functions are shown in Tables 2 and 3, respectively. Meanwhile, to ensure the effectiveness of the comparison results, the number of biological populations in all bionics algorithms is set as 10. The function expression is represented by *Fun*, and the expression dimension is expressed as *Dim*. *Lb* and *Ub* are the upper and lower boundaries of the variables, respectively.

**Table 2.** Algorithm parameters.

Sphere function is represented as F1, which is a typical single-peak function. The optimization time and local search ability of the algorithm can be effectively verified in this type of function. The Rastrigrin function is simplified as F2, and the validity of breaking away from local deadlock in the algorithm can be proved in this multi-peak function. The performances of different optimization algorithms in solving these two functions are expressed in Table 4. The iteration curves of these functions run once are shown in Figures 5 and 6.

**Table 4.** Comparison results in different algorithms.

**Figure 5.** The optimization of algorithms using Sphere function.

**Figure 6.** The optimization of algorithms using Rastrigrin function.

As can be seen from the Figures and Table, the average value, variance, and standard deviation of the different algorithms after 30 iterations in solving functions are different. The optimal value can be found quickly in the WOA, MFO, and GWO algorithms once they get rid of local deadlock when solving multi-variable functions. However, it is easy for them to fall into a local optimum, which is greatly affected by the number of iterations and randomness. Hence, it is not suitable for these algorithms to solve such functions. Compared with the SCA, PSO, and MA algorithms, although the convergence speed is not fast as the above algorithms, the rate is stable, and it is not easy to get into local deadlock. They have good performance in solving multi-variable lower-power functions, and MA has the most stable optimization rate among them. On the basis of MA, the stability of optimization is improved further in the QMA, and the optimal solution of the objective function can be found better in it. In the function where the optimal value is zero, the optimal value found by the QMA in the optimization process is lower than other algorithms, and the performance is more stable with little variance. In the scheduling problem of the microgrid, many variables are obtained in an economic objective function, and the stability of optimization is needed to ensure this. Hence, it is suitable for the QMA to find the optimal value of the economic scheduling function based on the microgrid established in this paper.

#### **4. Simulation Analysis of Microgrid**

shown in Figure 8.

Entering into cooperation with an enterprise, a typical winter day in the region where the enterprise is located is taken as an example to make economic scheduling. The actual data, such as wind–photovoltaic power and residential electricity consumption obtained from enterprise, are input as empirical values. Based on the economic scheduling objective function of the wind–photovoltaic-storage electric vehicle microgrid established in Section 2, the QMA and MA algorithms are used to solve it. The algorithms are invoked, combining with the microgrid economic objective function, and the number of the population in the algorithms is set as 100 after several simulations. In this condition, the objective function the algoriths can commendably solve the problems regarding optimization speed and convergence. Hence, the output power of each unit within 24 h is obtained. The output empirical values of wind and photovoltaic generation collected from enterprise are expressed in Figure 7, and the empirical values of residential electricity consumption are

**Figure 7.** Empirical values of wind and photovoltaic generation.

**Figure 8.** Empirical values of fixed load.

Under the load conditions, the optimization algorithms are used in this paper to find the minimum output cost under the condition of maintaining the security and stability of the system. Combined with empirical data obtained from the enterprise, the status of various components in the microgrid and the equivalence coefficient are shown in Tables 5 and 6. The related parameters of the time-of-use price are shown in Table 7.

**Table 5.** Status of various components in the microgrid.

**Table 6.** Equivalence coefficient.

**Table 7.** Market parameters.

When MA and QMA are used to solve the objective function, the SOC of the battery and EV are shown in Figures 9 and 10, respectively . The recharging of the battery is represented as an increase in SOC; otherwise, the battery delivers power to the grid. When the SOC of the battery is between 0.3 and 0.7, it can maintain its operation in the best condition. Taking Figure 10 as an example, when QMA is used to solve the system objective function, the battery is charged in three time periods: 1:00–2:00, 15:00–16:00, and 23:00–24:00. Meanwhile, discharge is performed in four time periods: 10:00–12:00, 14:00–15:00, 18:00–19:00, and 20:00–21:00. The rest of time, the battery remains neither charged nor discharged. Furthermore, if there is more than one EV in the area, their SOC is not continuous with great fluctuations, and they can operate normally when the SOC is within the range of 0.3 to 0.7.

**Figure 9.** The battery and EV SOC under the objective function solving by MA.

**Figure 10.** The battery and EV SOC under the objective function solving by QMA.

The output of each unit when MA and QMA solve the objective function is represented in Figures 11 and 12, respectively. The load power curve is the absolute value of consumed power.

Abundant power is generated by wind turbines within 0:00 to 4:00, and the microgrid consumes redundant power by interacting with the main grid. From 4:00 to 11:00, the output power of wind–photovoltaic generation and EVs supplies power to load; meanwhile, the redundant power is consumed by the main grid. From 11:00 to 16:00 is the peak period of electricity consumption. During this time, almost all power generation is used to maintain the power required by load. The load demand gradually decreases from 16:00 to 18:00 and reaches the peak again from 18:00 to 21:00. In this period, wind generation and EVs are in a power output state, while photovoltaic generation occasionally outputs power, and the interaction with the main grid absorbs or transfers power to the main grid as the

load requires. Finally, load demand drops from 21:00 to 24:00, and wind generation and EVs mainly send out power and feedback to the main grid.

**Figure 11.** The output of each unit under the objective function solving by MA.

**Figure 12.** The output of each unit under the objective function solving by QMA.

When the output power is at 10% of its maximum capacity it is stored by the battery, and the wind–photovoltaic generation produces power according to the empirical value. The operation cost of the microgrid is shown in Figure 13; furthermore, the operation cost of the microgrid obtained by solving the objective function using the scheduling strategy in Section 2 is shown in Figure 14.

The results of the value and standard deviation of MA and QMA after 100 and 500 iterations are shown in Table 8. In terms of these results, the operating costs under the empirical value are higher than the costs after participating in scheduling according to the strategy. Moreover, the optimization rate of QMA is higher than that of MA, and the optimal solution found by the improved algorithm after 500 iterations is obviously better than the solution obtained before the improvement.

**Figure 13.** Microgrid operation cost utilizing regular scheduling mode.

**Figure 14.** Microgrid operation cost after time-sharing scheduling.

**Table 8.** The cost of microgrid.

In summary, the scheduling strategy provided in this paper can effectively improve the economy of the microgrid system, and the improved QMA has a high level of superiority when solving the objective function.

#### **5. Conclusions**

According to the operation state of a wind–photovoltaic-storage electric vehicle microgrid connected to the grid in a certain area, a time-sharing scheduling strategy is designed in this paper; furthermore we established the minimum cost objective function according to economic requirements and constraints, such as wind and light abandonment. Meanwhile, a QMA algorithm is innovatively proposed and applied to the solution of the objective function. In summary, the main conclusions of this paper are expressed as follows:

In future research, new energy devices need to be taken into account, and the robustness of the system should also be paid attention to because of the increasing complexity of the system. Meanwhile, in order to make managers more intuitive, namely to understand the operation state of the system, the connection between information management and the physical layer should be strengthened. It is necessary to study the use of hierarchical scheduling, source and charge interactions, and other methods to meet the needs of different managers.

**Author Contributions:** X.L. wrote the manuscript, designed the scheduling strategy, drew the flow chart and table, and performed the simulation experiment and analysis. M.Z. supervised the study, helped obtain data, and checked the contents of the manuscript. Z.W. verified the simulation data and contributed to the writing and typesetting of the manuscript. M.L. provided practical suggestions and evaluated the feasibility of the application. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was Supported by the Natural Science Foundation of China under Grant No. 61563045.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Thanks to the electric power enterprises in Liaoning Province for providing data support for this research.

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

#### **Acronyms**

The following acronyms are used in this paper:

#### **References**

### *Article* **Smart Cooperative Energy Supply Strategy to Increase Reliability in Residential Stand-Alone Photovoltaic Systems**

**Carlos Roldán-Blay 1,\*, Carlos Roldán-Porta 1, Eduardo Quiles 2 and Guillermo Escrivá-Escrivá 1**

**Featured Application: The proposed study can be applied to optimize energy consumption, generation and storage in neighbour communities with distributed energy resources.**

**Abstract:** In reliability studies of isolated energy supply systems for residential buildings, supply failures due to insufficient generation are generally analysed. Recent studies conclude that this kind of analysis makes it possible to optimally design the sizes of the elements of the generation system. However, in isolated communities or rural areas, it is common to find groups of dwellings in which micro-renewable sources, such as photovoltaic (PV) systems, can be installed. In this situation, the generation and storage of several houses can be considered as an interconnected system forming a cooperative microgrid (CoMG). This work analyses the benefits that sharing two autonomous installations can bring to each one, from the point of view of reliability. The method consists of the application of a random sequential Monte Carlo (SMC) simulation to the CoMG to evaluate the impact of a simple cooperative strategy on the reliability of the set. The study considers random failures in the generation systems. The results show that the reliability of the system increases when cooperation is allowed. Additionally, at the design stage, this allows more cost-effective solutions than single sizing with a similar level of reliability.

**Keywords:** microgrid; cooperation; energy; reliability; Monte Carlo; smart energy control

#### **1. Introduction**

The increase in energy demand and dependence on fossil fuels [1] has boosted the integration of renewable energies [2]. This has favoured the development of distributed energy resources (DERs) [3], especially in rural or isolated areas [4].

Small size DERs in Spain are usually implemented through the so-called stand-alone photovoltaic (SAPV) energy systems [5–7], which generally include batteries to store energy and supply night demand. These types of systems are currently very attractive in rural areas thanks to the new regulations that have regulated self-consumption since 2019 [8]. In a general case, the energy storage system (ESS) could also be based on renewable energies, as in the case of micro hydroelectric plants or green hydrogen generators [9,10]. The great disadvantage of this type of system is the reliability of supply, since a profitable design usually implies a large number of hours per year of energy not supplied (ENS), increasing indices such as the loss of load probability (LOLP) to inappropriate levels [11]. Due to the high price of ESSs, especially those based on renewable energies, an optimized design of the generation system is of vital importance. In this sense, there are studies that focus on improving the performance of these systems [12–14]. On the other hand, other works are focused on improving the management of generation, demand and storage resources to obtain maximum reliability with the lowest possible cost [15,16].

**Citation:** Roldán-Blay, C.; Roldán-Porta, C.; Quiles, E.; Escrivá-Escrivá, G. Smart Cooperative Energy Supply Strategy to Increase Reliability in Residential Stand-Alone Photovoltaic Systems. *Appl. Sci.* **2021**, *11*, 11723. https:// doi.org/10.3390/app112411723

Academic Editor: Mohamed Benbouzid

Received: 10 November 2021 Accepted: 7 December 2021 Published: 10 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).