Optimal distributed generation and battery energy storage units integration in distribution systems considering power generation uncertainty

This paper proposes an application of the recent metaheuristic rider optimization algorithm (ROA) for determining the optimal size and location of renewable energy sources (RES) including wind turbine (WT), photovoltaic (PV), and biomass-based Distributed Generation (DG) units in distribution systems (DS). The main objective function is to minimize the total power and energy losses. Power loss-sensitivity factor (PLSF) is used with the ROA to determine the suitable candidate buses and accelerate the solution process. The Weibull and Beta probability distribution functions (PDF) are employed to characterize the variability of wind speed and solar radiation, respectively. The high penetration of intermittent renewable resource together with demand variations has introduced many challenges to distribution systems such as power ﬂuctuations, voltage rise, high losses, and low voltage stability, therefore battery energy storage (BES) and dispatchable Biomass are considered to smooth out the ﬂuctuations and improve supply continuity. The standard 33 and 69-bus test systems are used to verify the effectiveness of the proposed technique compared with other well-known optimization techniques. The results show that the developed approach accelerates to the near-optimal solution seamlessly, and in steady convergence characteristics compared with other techniques.

Standard (RPS) [2]. By accepting an RPS, it is an obligation to produce a certain percentage of total electricity generation from RES during a specific date. However, the energy production of wind and photovoltaic generation is determined by local weather conditions, meaning that its production is variable and difficult to predict. In case of high levels of RES, large and sudden changes in the power generated by RES can threaten the provision of a continuous power supply. Further, the amount of WT and PV energy that can be absorbed by the power system may be significantly limited due to the inability of conventional generating units to counteract the fluctuations in WT and PV. To encourage the use of these energy sources and support the stability of the system, BES and dispatchable biomassbased DG can be used to avoid violation of system constraints and decrease power fluctuations in DS [3]. Several studies have shown that the optimal allocation of DGs such as PV and WT can improve the reliability, security and efficiency of the system, by mitigating their variable nature problems. Indeed, an unplanned distribution of RES can have various negative effects on system operation and stability [4].
To date, many previous studies have been focused on determining the optimal allocation and sizing of RES in distribution networks, in particular considering hybrid planning. Recent trends are to use metaheuristic optimization methods such as the Genetic Algorithm (GA) [5], the Bat Algorithm (BA) and Particle Swarm Optimization (PSO) [6], An efficient analytical [7], Stochastic Fractal Search Algorithm (SFSA) [8]. In [9], Electrostatic Discharge Algorithm (ESDA) has been used to solve the DG unit's allocation problem with the aim of minimizing the power losses and improving the voltage stability. Mansur et al. have introduced an Artificial Ecosystem-based Optimization method for optimal allocation of DGs in the DS in terms of minimizing annual power losses in [10].
An analytical method was presented in [11] for optimal sizing autonomous hybrid PV+WT units, taking into account the time factor and system costs for the given load. In [12,13], an effective optimization sizing model has been developed for hybrid solar/wind power generation system. The mixed-integer linear programming (MILP) approach has been used in [14] to determine the optimal placement and sizing of DGs in radial DSs. Koutroulis et al. [15] presented the genetic algorithm (GA) for finding the optimal size of PV+WT units with BES in distribution system. In [16] An improved Harmony Search algorithm (HSA) approach have been presented, and the proposed work is aimed at evaluating the optimization solution to find the optimal size of the PV/battery hybrid unit to deliver and achieve the required power demand of the case study. In [17] Geographic Information System module-based structure has been proposed to optimally size the stand-alone photovoltaic (PV) -diesel units in rural areas. This structure determines the best location based on technical, economic, reliability, social and environmental criteria. In [18] hybrid structure has been proposed to minimize the value of the life cycle and employed to assess the viability of the grid independent PV unit with hydrogen storage considering for the reliability of the hybrid system (HS). In [19] hybrid fuzzy logic controller (FLC) and the harmony search algorithm (HSA) has been proposed. This work focuses on sizing a hybrid energy system including the PV units and WT with/without considering the backup system (ESS) to reduce economic costs and increase reliability. In [20] probabilistic-based planning approach has been presented for identifying the optimal fuel mix of different types of renewable DG units in order to minimize the annual energy losses in the DS without violating the system constraints.
In [21] weighted aggregation particle swarm optimization (WAPSO) method has been proposed for optimal location and sizing of PV and WT DGs in DS by considering power loss minimization, voltage stability and network security improvement. The stochastic nature of solar irradiance and wind speed are accounted using suitable probabilistic models. In [22][23]an analytical method has been proposed to determine the optimal size and power factor of individual DG unit, simultaneously, at each location for minimizing power losses by considering the time-varying load demand combined with the possible characteristics of DG units. In [24] presents a concept for connecting a hybrid PV and BES system in a commercial DS with the aim to minimize energy loss and improvement voltage stability. In this concept, the BES unit is considered as a dispatchable source, whereas the PV unit is considered as a non-dispatchable source. In [25], a hybrid energy system consisting of two subsystems with two scenarios has been developed. In the first one, these two subsystems have been connected to allow load flow between them while there is no connection between the subsystems in the second scenario. Determining the optimal size of PV and WT with BES, as well as the best combination in HS has been investigated in [26]. Power flow and related power losses due to the influence of HSs with the grid have been also examined. Mansur et al. have been introduced the tree growth algorithm for optimal allocation of DGs in the DS for minimizing annual power losses [27]. A charge and discharge strategy of BES units have been presented to mitigate the unexpected changes in PV outputs [28]. Efficient voltage regulation in DSs by managing the BES units' output on the consumer side with high PV penetration has been introduced in [29]. The optimal size planning of BES units and PV-based DG units for minimizing energy loss has been studied in [30]. In [31] proposed the Crisscross Optimization Algorithm (COA) and Monte Carlo Simulation (MCS) for allocation of DGs in the DS for minimizing the power losses and the cost. In [32], a probabilistic planning technique has been suggested based on mixed integer nonlinear programming (MINLP) and has been implemented for energy loss minimization with optimal integration of RESs in a rural DSs. A simplified analytical method has been proposed in [33] for optimal DG integration to reduce power loss in DS, and the results of power loss reduction were analysed using the other four methods for one, two and three DGs. In [34], the exact loss equation-based analytical method has been used to find the optimal size and the location of a single DG in different DSs. In [35] the mixed-integer nonlinear programming (MINLP) method have been used to determine optimal location and size of DGs with the aim of power loss minimization. In [36], the classical Kalman filter algorithm is used to compute the size and location of the DG. In [37] proposed Mixed Integer Non-Linear Programming (MINLP) based DG planning model for cost minimization. A Linear Programming (LP) based model for maximizing DG injection in deregulated environment is presented in [38]. In [39], a multi-objective MINLP model for a trade-off between operating cost and DG energy is proposed for rural microgrid. In [40] formulated optimal sizing and siting problem as Mixed Integer Linear Programming (MILP) problem. Since the optimal placement DG unit is MINLP formulation, linear model might affect the power flows in some cases. The minimization of energy losses for stochastic renewable generation is formulated in [41].
Overall, the above review shows that significant work has been done on the optimal size and placement of DG units of various types with BESs in grid-connected and gridindependent HSs. While most research has focused on off-grid PV&WT generation systems, few studies have analysed gridconnected hybrid power systems. Only a few studies have been reported on the allocation of DG units taking into account timevarying load demand to minimize energy losses [20][21][22][23]. In addition, most of these approaches assumed that the DG unit operates at fixed power factors. In these studies, only location and size were considered, while the optimal power factor of a single DG unit, which would be an important part of minimizing energy losses, was ignored. Several papers have been discussed on the optimal power factor to reduce power losses [24][25]; however, the optimal power factor of an individual DG unit for each location was not considered. In addition, a "controlled" power factor was proposed for each period in a given range to minimize energy losses [10] and to minimize reactive support from the transmission network [12]. However, most of the work presented assumes that the DG output is dispatchable and controllable. Most of the available methods cannot simulate the stochastic nature of the power output of DG units. This article addresses some of the shortcomings of existing methods.
This paper first proposes an application of the recent heuristic rider optimization algorithm (ROA) for optimal allocation of PV and WT-based DGs in DSs. The optimal scheduling BES and Biomass units is discussed with the aim of converting the system to be dispatchable and minimizing the total energy losses. Power loss-sensitivity factors (PLSF) are used to determine suitable candidate buses and accelerate the solution process. The voltage stability index (VSI) is calculated after the optimization process to measure the performance of the proposed method on the stability of the DS. The Weibull and Beta probability distribution functions (PDF) are employed to characterize the variability of wind speed and solar radiation. The performance of the developed technique is tested on the 33 bus and 69 bus DSs and for a fair compression compared with Harris Hawks Optimizer (HHO) [42], Henry gas solubility optimization (HGSO) [43] and other techniques in the literature that solved the same problem. In addition, comparison of statistical indicators of ROA vs HHO, and HGSO has been done to check random nature of proposed algorithm.
This paper is structured as follows: PLSF is described in Section 2. Probabilistic power generation of PV and PV based DGs, Biomass, BES, and load modelling are presented in Section 3. Section 4 proposes the mathematical formulation and optimal PV, WT, Biomass, and BES sizing. The solution process of the proposed technique is presented in Section 5. Numerical results are described in Section 6 and conclusions of the paper are mentioned in Section 7.

POWER LOSS SENSITIVITY FACTOR
In this paper, PLSF is used to determine the candidate buses to reduce the search space. In order to explain the factor, Figure 1 shows the equivalent grid-connected two buses radial DS [1,9]. where, P i and P k are active power flow out of bus i and k, Q i and Q k are reactive power flow out of bus i and k, V i and V k are the voltage magnitudes. R ik and X ik are distribution line resistance and reactance between bus i and k. P PV +BES ,k and Q PV +BES ,k are the active and reactive power output of PV+BES at bus k, respectively; P d,i and P d,k are the active power load at bus i and k, Q d,i and Q d,k are reactive power load at bus i and k, respectively. The active power loss between buses i and k is calculated as [1,44,45]: The total energy losses E loss in 24 h are expressed as follows: where, Δt is the time duration, in this study (1 h). The PLSF is calculated by (3) [33,34]. Figure 2 and Figure 3 show the PLSFs for the 33-bus and the 69-bus test systems, respectively. After calculating the PLSFs of all buses, they are sorted in descending order. The buses with higher PLSF are considered more suitable for installing PV, WT, Biomass, and BESs units (up to 50% of system buses) [44][45][46].

Power generation by PV and WT model
Power generation from PV and WT is highly dependent on weather conditions, such as solar radiation, ambient temperature, and wind speed. Hence, wind speed and solar radiation at a specified location must be appropriate modelled in order to obtain good solutions. The standard deviation (SD) and mean of hourly wind speed and solar radiation per day are calculated using collected historical data. To incorporate the output power of the solar DG and wind-based DG units as multistate variables in the planning formulation, the continuous probability distribution function (PDF) of each has been divided into states (periods), in each of which the solar irradiance and wind speed are within specific limits. In other words, for each time segment, there are number of states for the solar irradiance and wind speed [20]. In this paper, the stage for solar radiation and wind speed are 0.05 kW/m 2 and 1 m/s, respectively. The average value of each

Solar radiation modelling
It is considered that the probabilistic nature of solar radiation follows the Beta PDF [46,47]. The Beta PDF of solar radiations (kW/m 2 ) in the time interval 't'is defined as Equation (4): where, f b (s t ) is the Beta PDF of s t , α t and β t are the shape rates of Beta PDF and Г depict Gamma function.
The shape rates of Beta PDF can be found using mean (μ) and SD (σ) of radiation for a suitable time interval as follows [46,47]:

PV array power generation
The PV array hourly average power output corresponds to an exact time interval 't' (P t PV ) expressed as Equation (6) [47]. A typical day for three years is generated in p.u., as shown in Figure 4.
where 'g' denotes a stage factor and n s is the solar radiation discrete stage number. S t g is the gth stage of solar radiation at tth time interval. Solar radiation and ambient temperature are the basic dominant factors that affect the PV array power output. The PV power generation with average solar radiation (s ag ) for the gth stage is estimated as [27,47]: where N PVmod is PV modules total number. The volt-ampere characteristic of the PV module can be definite for a given stage using the following equations [22,27]:

Wind speed modelling
Weibull PDF was chosen to evaluate the stochastic behaviour of wind speed at a predetermined duration of time [21,48]. Weibull PDF for wind speed v t (m/s) at the tth time interval can be calculated as: The shaping rate (k t ) and scale rate (c t ) at t th time interval are expressed as [21,36]: where, t v and t are mean and SD of wind speed at time interval 't'.

WT power generation
The hourly WT average output power corresponds to a specific time interval 't' (P t WT ) can be expressed as Equation (11). A typical day for three years is generated in p.u., as shown in Figure 4.
where 'g' denotes a stage factor and n s is the number of wind speed discrete stage. v t g is the gth stage of wind speed at tth time interval.
The WT power generation [21] with an average wind speed (v ag ) for stage "g" is expressed as: where P r is the nominal power rate that WT can be generated; v cout is cut-out; cut-in (v cin ) and nominal (v r ) wind speed, respectively, constants A and B are achieved as [21]:

Biomass modelling
Biomass has become one of the world's largest energy sources (REF). Despite the availability of various biomass resources, it was not optimally used to contribute to improving the energy sector's production. The share of the agricultural sector in waste production is higher. Recently, these types of wastes have been used as a source of biogas to serve rural areas. In the development of conversion technology, as in many countries, it can be widely used to generate electricity [49]. Currently, the most important biomass processing technologies are pyrolysis, gasification, incineration, and digestion. The anaerobic digestion process is used for raw materials such as bulgur, while other processes include the use of dry feedstocks. The advantage of using this type of renewable energy, biomass-based DG units produce little emissions [50].

Dispatchable source
The biomass-based DG units used as a synchronous generator. Besides, the biomass-based DG output can be distributed accordingly to the load curve. The Biomass-based DG unit's capacity factor (CF) is determined as the ratio of the area under the power output curve in p.u. to the total period [22].

Nondispatchable source
The output of biomass-based DG can be also considered nonintermittent (constant) at its nominal capacity for 24 h a day. Its CF is equal to one.

BES modelling
Over the last few years, several types of energy storage (ES) have been intensively studied. These include super-capacitors, electrochemical battery, compressed air ES, flywheel ES, and superconducting magnetic ES. For this paper lithium-ion (Li-ion) BES is selected with a roundtrip efficiency of 77%. This is the most popular type of BES for today's portable electronics, with one of the best energies and weight ratio, no memory effect, and slow charge loss on non-use. Incorrect handling may cause a Li-ion battery explosion. Lithium-ion BES is becoming increasingly popular with defence, aerospace, and automotive applications due to its high energy intensity [51]. The Charging and discharging of the BES are given in Equation (15) [51].
where, P disch BE S k is discharged power by the BES for a t time duration. P ch BE S k is the power charged by the PV to the BES, that is, the BES is being charged. E BE S k (t) is the energy stored in the BES at the time (t). Δt is the time duration of each segment. d and c are respectively, the discharging and charging efficiency.
BES must be satisfying the restrictions from Equations (16) to (18).
Power limits [46,51]: Stored energy limits [39]: Starting and ending energy limits [51]: where P disch,max BE S k is the maximum discharge power rate;P ch,max BE S k is the maximum charging power rate; E BE S k min and E BE S k max are the minimum and maximum energy limits of the BES; E BE S k (0) is the initial energy of the BES; E BE S k s is the BES unit initial limit of stored energy. For the energy balance, the stored energy E BE S k (T) is set the same as the initial stored energy. In this paper, it is assumed that the BES unit's minimum and maximum capacity limits are 20% and 90%, respectively [52].

Load modelling
A 24-hour daily load curve with a peak of 1 p.u. is used in the proposed model, as shown in Figure 4 [53]. The load factor (LF) can be determined as Equations (19) follows [17]: A time-varying voltage-dependent load demand model depends on the time and voltage. Consequently, the voltagedependent load demand model in [54], which included variable loads at time t, can be calculated as: (20) where P k and Q k are actual and reactive power which is at bus k; P ok and Q ok are load actual and reactive power at bus k; V k represents the voltage at buss k, and n p = 1.51 and n q = 3.4 are the actual and reactive load voltage indexes [54].

Objective function
The considered objective function in this study is the minimization of the total power and energy losses in grid connected HS. In this study the PV, WT, WT+Biomass, and PV+BES units are optimally integrated in the HSs to achieve the intended goal, taking into account the operational constraints. The schematic diagram of the grid connected HS is shown in Figure 5. The objective function can be described as follows [20][21][22][23]: where, E loss is the total energy loss. The VSI is used for indicating the voltage stability of DS. High VSI of any bus identifies the less sensitive bus to the voltage collapse. The VSI of bus k is calculated as follows [9]: The objective function is subjected to following constraints such as DG size, bus voltage, power factor, and transmission line total power limitation.

Equality constraints
The generated power must meet the following constraints [9]: where, P slack and Q slack represent the active and reactive power of slack bus. N DG is the installed DG units' number, L is the total branch number.

Voltage constraint
Bus voltages must meet the requirements between its lower V min and upper V max limits:

Integrated DG units power constraint
The installed DGs size must be limited as [9]: where, P max DG and P min DG represent the maximum and minimum DG unit's active power, respectively. Q max DG and Q min DG represent the maximum and minimum DG units' reactive power, respectively.

DG's power factor limits
The DG unit's PF must be within PF DG ,min and PF DG ,max as [9]:

Branch capacity limitation
The maximum capacity of branch must meet the following limits:

4.2.1
PV sizing Figure 6 [22,23] shows the proposed grid-connected PV+BES HS conceptual model. This HS is designed for installation on the roofs of commercial buildings. The concept is to convert every non-dispatchable PV unit output to a dispatchable PV unit with a combination of PV and BES unit, to minimize the total power system loss for a given load level. This strategy can generate a daily amount of dispatchable energy, E (PV +BES ) [22,23]. The E PV is PV unit generating energy, during the 24 h cycle of the day. E Grid PV is the transferred energy to the grid. The PV unit's excess energy is used to charge the BES unit, E ch BES instead of cutting it when the PV output power is bigger than load during the day. When the output power of the PV is low or zero overnight, this stored energy is then discharged to the grid, E disch BES . PV and BES units are located on the same bus to avoid energy losses in the charge state of BES. The daily energy amount of PV+BES unit for all the duration (T = 24 h) on bus k is expressed as [22,23,46]: where, P (PV +BES ) (t ) is the combination of PV and BES unit's active power output on bus k during a given day time interval t. The daily charge and discharge energy at bus k, shown in Figure 6, is calculated by Equation (29).
The PV+BES and the PV based DG unit total output energy on the given k bus are expressed as [23]: where, E Grid PV k is the PV unit energy amount transferring to the grid on the bus k. The BES unit charge and discharge energy on the bus k with a roundtrip efficiency ( BES = c × d ) is given as follows [22,23]: The PV unit total output energy on the bus k is given as [23]: The maximum PV unit power generation during a specific time duration over a 24-hour cycle of the day is used to indicate the PV unit nominal power or optimal size on the bus k [22].
where, CF unit PV = P unit PV E unit PV , P unit PV is PV unit maximum power output, and E unit PV is the PV unit generated energy during 24 h.
is then taken from Figure 6. When BES is less than a unit, P PV k increases. Anyway, as shown by the simulation results, E Grid PV k increases slightly compared to E Grid ′ PV k . Therefore, the optimal size of PV is obtained from Equations (32) and (33) as [22,23]:

BES sizing
The optimal BES unit size at bus k is calculated based on the nominal power (P BE S k ) and the energy capacity (E BE S k ) as well as it can hold all the excess a PV unit energy, which reduces the power losses for each time interval at the minimum. Maximum charge and discharge over a specific duration during the 24-hour cycle of the day is used to determine BES unit nominal power [22,23,46]. The maximum energy of charge during this period is used to identify BES unit energy capacity.

Sizing a combination of WT and biomass (WT+Biomass)
Here, it is assumed that the WT based DG is non-dispatchable, it belongs to the developers of the DG and it is controlled by a utility company. The biomass-based DG is dispatchable, owned, and controlled by a utility company. Figure 7. [24] shows an example of the proposed WT and Biomass hybrid system (WT+Biomass) power output curve. The concept is a combination of each non-dispatchable WT unit output and dispatchable Biomass units, to minimize the total power system loss for a given load level [24]. The daily energy amount of WT and Biomass unit for all the duration at bus k can be expressed as [24]: where, P (WT +Biomass) (t ) is the output power of the combination of WT and Biomass units at bus k during a given day time interval t.
The proposed procedure of the combination WT+Bomass is described as: Stage 1: Run the power flow, determine the optimal size, location, and PF of the WT+Biomass considering the total minimum power loss over all the time interval. Stage 2: Find each non-dispatchable WT unit optimal size or maximum power output (P max WT ) over all the time interval on the condition that, WT unit power output is no bigger than the output of WT+Biomass unit as determined in Stage 1, at each time interval. Stage 3: Find the optimal output power of each WT unit for time interval t as expressed in Equation (36) [24].
where, p.u. DG WT output(t) is the WT output in p.u. at time interval t.
Stage 4: Calculate the output power for each dispatchable Biomass-based DG unit. This is equal to the power output of the combination of WT and Biomass unit in Stage 1 minus the power output of the WT unit in Stage 3, for each time interval t, then find the dispatchable Biomassbased DG unit optimal size or maximum power output overall time interval.
In Figure 7, the WT output power follows the expected WT power generation curve in Figure 4. In this case, the wind power penetration is at its maximum. The biomass-based DG units are used as an additional dispatch source to fill a portion of the energy that the WT based DG units cannot use.

RIDER OPTIMIZATION ALGORITHM (ROA)
ROA is a recent metaheuristic optimization technique. It is inspired by a group of riders [55]. The group number is divided into four groups. The first group name is bypass rider (BR), the second group is follower rider (FR), the third group is overtaker rider (OR) and the fourth group is attacker rider (AR). The group numbers have equal riders BR = F = O = A = R∕4. The groups have different strategies to reach the final solution.
The basic stages of the proposed algorithm can be summarized as follows: Stage 1: Randomly initialize the position of the group, and parameters of rider, such as steering angle, brake, acceleration, and gear [55].
where, R is the number of total riders, it is also equal to the number of groups (G). D is a dimension of coordinates.X t ( j , k) is the position of the jth rider at time t.
T t (i, j ) is steering angle of the ith rider vehicle at time t and can be found as follows [55]: where, i is position of the ith rider vehicle, is coordinate angle and they can be expressed as: The gear, accelerator and brake are can be initialized as follows [55]: where, E i is the gear of the ith rider's vehicle and at time t it takes from the set {0, 1, 2, 3, 4}. The initial value of the gear is zero. e i is the accelerator of the ith rider's vehicle and ranges from 0 to 1. The initial value of the gear is zero. K i is the brake of the ith rider's vehicle and its taking value is within 0-1. The initial value of the brake is one. The riders can adjust the speed of the vehicle. Therefore, the maximum speed cane be expressed as [55]: where, X i U and X i L are maximum and minimum value in i th rider's position, respectively. T OFF is off time of race. V i max is maximum speed of the ith rider.
The speed limit of the gear for ith rider is can be expressed as: where |E| is number of the gears. Stage 2: The rider's success rate (SR) is calculated as follows [55].
where X j is jth rider position. L T is the position of the target. Stage 3: The leading rider can be found as follows. The nearest rider from the goal position is considered the leading rider. Stage 4: The riders' positions can be updated as follows: 1. The BR can achieve the goal by bypassing the leading path.
This refers to the BR does not follow the LR.
BR's position can be updated as follows: where, and are random [0,1] numbers, and are random [1, R] numbers.
2. The FR is one who follows or depends on the LR on most of the axis.
FR's position can be updated as follows: where, i is coordinate selector, X L is the LR position, L is the LR index, T t j ,i is angle of steering in the ith coordinate, d t j is the distance to be travelled of the ith rider.
3. The OR follows his position to achieve the goal, in conformity with the nearest position of LR.
OR's position can be updated as follows: where X t ( j , i ) is the position of the rider in the ith coordinate, Q I t is indicator of rider direction at time t.
4. The AR is an aggressive player. It takes the rider's position to achieve the goal using full speed.
AR's position can be updated as follows: (50) where, X L (L, k) is the position of the rider in kth coordinate. Stage 5: New LR can be found after the position update process.
Stage 6: The riders' parameters updating. The parameters of rider use the activity counter for updating their value, which is updated in conformity with the SR.
1. The activity counter update can be calculated as follows: 3. Gear update can be calculated as follows [43]: (53) 4. Accelerator update can be calculated as follows [55]: 5. Brake update can be calculated as follows [55]: Stage 7: Riding off time. At the end of the race, the LR becomes the winner. Table 1 presents the parameters of the proposed algorithm and operational limitation of the study [1,9,55].
The overall procedure of proposed algorithm to solve optimization problem in this study is shown in Figure 8.

NUMERICAL RESULTS
In this section, the proposed approach is examined on the standard 33-bus and 69-bus test systems. It is simulated using MAT-LAB R2018b.
To prove the effectiveness of proposed approach, the following scenarios are considered: Scenario #1: In this scenario, the optimal allocation of multiple PV, WT units at unity and optimal power factors for minimizing power losses without voltage-dependent time-varying load demand is considered. The single line diagram of these systems is shown in Figures 9  and 10 [1,9]. Base kV and MVA of the system are 12.66 kW and 100 MVA for 33-69-bus DSs. The active and reactive power loads of the 33-bus test system are 3,715 kW and 2,3 kVAr, respectively. Initial active power loss and the minimum voltage Other information can be found in [56]. The active and reactive power loads of the 69-bus test system are 3801.5 kW and 2694.6 kVAr, respectively. Initial active power loss and the minimum voltage of the test system are 224.98 kW and 0.90919 p.u, respectively. Other information about the 69-bus test system can be found in [57]. The mean and St for each hour per day are calculated using historical data of solar radiation and wind speed per hour and given in [21].

33-bus system
The obtained results of proposed method for optimal size, location and power factor of single and multiple DG units in the 33-bus system are listed in Table 2 [58], PSO [7], Exhaustive Load Flow (ELF) [36], Improved Analytical (IA) [36], Mixed Integer Non-Linear Programming (MINLP) [36], Novel heuristic approach (NHA) [59], HHO and HGSO. Besides, the convergence curves of the three methods (i.e. ROA, HHO, and HGSO) are shown in Figure 11. The results show that the ROA accelerates to the near-optimal solution seamlessly, and in steady convergence characteristics compared with HHO and HGSO algorithms. In addition, the integration of WT-based DG provides better results in terms of power loss reduction, voltage profile and VSI compared with PV-based DG due to its capacity to generate reactive power.

69-bus system
The results of proposed method for optimal size, location and power factor of single and multiple DG units in the 69-bus system are listed in Table 5     minimum power loss in one and several PV-based and WTbased DG units compared with those obtained by Bat Algorithm [6] Particle Swarm Optimization (PSO) [6], Hybrid Salp Swarm Algorithm (HSSA) [58], EA [7], Exhaustive Load Flow (ELF) [36], Improved Analytical (IA) [36], Mixed Integer Non-Linear Programming (MINLP) [36], Novel heuristic approach (NHA) [59], combined Non-Linear programming with Power Loss Sensitivity (NLP&PLS) [60], HHO and HGSO. Besides, the convergence curves of the three methods (i.e. ROA, HHO, and HGSO) are shown in Figure 12. The results show that the ROA accelerates to the near-optimal solution seamlessly, and in steady convergence characteristics compared with HHO and HGSO algorithms. In addition, the integration of WT-based DG provides better results in terms of power loss reduction, voltage profile and VSI compared with PV-based DG due to its capacity to generate reactive power.

FIGURE 12
Convergence curve of the 69-bus system

Power loss (KW)
On the whole, the heuristic techniques are distinguished by its randomness. Therefore, many tests have been performed to prove the robustness of the ROA with 15 independent runs. The optimization objective function's convergence curves of the three methods (i.e. ROA, HHO, and HGSO) of the P loss are shown in Figures 11 and 12, respectively. In comparison with HHO, and HGSO algorithms, the results show that the ROA accelerates to the near optimal solution seamlessly, and in steady convergence characteristics. Also, the efficacy and robustness of the proposed ROA over HHO, and HGSO algorithms are verified based on statistical factors, where many trials have been made. The minimum, maximum, mean, the standard deviation (SD) of the objective function after 15 runs for standard 69node DS with integrating of three PVs with PLSF as presented in Table 8 and without PLSF as presented in Table 9. These compressions used to check if the near-optimal solution located inside this limited search space or not. From these compressions, we can see the PLSF is suitable for the optimization problem. In addition, the proposed ROA algorithm has more robust statistical indicators than the other algorithms, such as

Scenario 2 considers the minimization of energy loss as the main objective function
In this scenario, uncertainty model of power generation is applied using the optimal sizes and locations obtained in Scenario 1, the time-varying daily load demand is used, as shown in Figure 4. The effect of load variation on the voltage profile over 24 h for 33-bus and 69-bus are shown in Figure 21.
Three PV and WT units are integrated at buses 13, 24, and 30 for the 33-bus system and 18, 11, and 61 for the 69-bus system with 24-hour power generated, as shown in Figures 13 and 14.
A significant improvement occurs in the voltage profiles and loss reduction when the power of PV and WT is available from 6:00 to 19:00 h and from 1:00 to 19:00 h, as shown in Figure 21(b,c) and Figure 22(b,c). respectively. However, the lack of integration of PV and WT is investigated after from 7:00 to 9:00 h and from 13:00 to 14:00 h, when the voltage is still under the lower limit of 0.95 p.u., due to the very small probabilistic depends of PV power and an increase in load demand. When the   16 Daily PV+BES outputs for the 69-bus system voltage is over the upper limit of 1.05 p.u., due to the high probabilistic depends of WT power and decrease in load demand.

6.3
Scenario 3 optimal planning of the PV+BES power generation for minimizing the total energy losses and improving the voltage profile As mentioned above, BESs are installed at the PV locations and the proposed technique is used at every hour to optimally determine the power of the PV+BES with the aim of minimizing the energy losses. Using the power of PV+BES, the optimal size for the PV and BES can be achieved using the proposed model in Section 4.2. Figures 15 and 16, shows the PV+BES units output power curves at buses 13,24 30 for 33 bus systems and 11, 28,61 for 69 bus system. The power output curves of combination model units follow the load demand curve in Figure 4 because PV+BES units are considered dispatchable sources. Figures 17  FIGURE 17 PV and PV+ BES power outputs curves at bus 13 of 33-bus system FIGURE 18 PV and PV+ BES power outputs curves at bus 11 of 69-bus system and 18 show the hourly output power of PV and PV + BES at buses 13 and 11 for 33 bus and 69 bus systems, respectively. This PV output curve exactly follows the expected PV power output curve shown in Figure 4. The PV unit's maximum power output is determined at 12 h. and it specifies the optimal size of PV-based DG. From Figures 17 and 18, it can be observed that the PV unit generates the energy against PV+BES unit energy accommodated by the system to keep the minimum power loss. The hourly differences between the two schemes determine the charge and discharge energy of the BES unit. Using the power output of PV + BES, the optimal sizes for PV and BES can be achieved using the PV + BES model.
The maximum power output difference is found at 13 h., which provides the maximum charging power or the nominal power rating of the BES unit (P BES ). Similarly, for the remaining units at buses 24, 30, and 18, 61, the results are obtained.
Tables 10 and 11 present the results for the PV unit's sizes as well as the nominal power and power output of the BES units for each location for 33-bus and 69-bus test systems, respectively. Scenario 3 optimal planning of the WT+Biomass power generation for minimizing the total energy losses and improving the voltage profile Figure 12 shows the optimal hourly output curves for WT and Biomass units at bus 18. The total generation of DG units of WT + biomass in each time interval in Figure 12 is also follows to the load demand curve in Figure 3. At the same time,   the wind power output corresponds to the wind power output curve in Figure 3, provided that the wind power penetration is maximum. As noted above, biomass-based DG units are used as additional supplies to fill some of the energy that WT-based DG units cannot afford. The WT based DG units in Figure 12 cannot afford the load demand of Figure 3. The maximum power output difference found at hour 13, gives the Biomass-based DGs maximum power output or nominal power rating (P Biomass ) of Biomass-based DG unit. Similar results are obtained for the WT+Biomass units at buses 18 and 61. Also,  Figure 13(d) shows the effects of combined DGs installation on voltage profiles. From these figures it is observed that the voltage profile is significantly improved in the case of locating the DGs in their optimal allocation and optimal power factor using the proposed method. Besides, this improvement is increased in the case of using WT+Biomass-based DG compared with PV+BES-based DG units due to its ability to generate reactive power and at optimal PF. Table 12 and 13 summarize the simulation results obtained in Scenario 4 for 33 bus and 69 bus test systems, respectively. These tables include DG unit type, size, location, and PF for each type.

Energy losses analysis: 33-bus and 69-bus test systems
A noticeable decrease in energy loss in scenarios 2-4 is observed compared with the base case. The decline in the highest energy loss is determined in scenario 3-4, while the smallest energy loss is realized in Scenario 1. Total energy losses and its loss reduction for a given day are presented in Tables 14 and 15 for 33 bus and 69 bus systems, respectively. The maximum reduction in energy loss is 94.04 % and 98 % for 33 bus and 69 bus systems, respectively.

CONCLUSION
In this paper, a recent metaheuristic rider optimization algorithm (ROA) technique with power loss-sensitivity factor  (PLSF) has been proposed to solve the optimal allocation of single and multiple photovoltaic (PV), wind turbine (WT), Biomass, PV+Battery energy storage (BES), and WT+Biomass units in distribution systems, taking into account the timevarying load demand, optimal power factor, and probabilistic power generation. The voltage stability index (VSI) is calculated after the optimization process to measure the performance of the proposed method on the stability of the distribution systems. The Beta and Weibull probability distribution functions (PDF) models have been employed to describe the stochastic nature of solar radiation and wind speed. The objective function has been formulated as the minimization of active power loss of distribution system. The 33-bus and 69bus test systems have been used to demonstrate the effectiveness of the proposed approach. Three PV, WT PV+BES, and WT+Biomass units have been considered to check the capability of the proposed approach for minimizing energy loss and enhancing voltage profile. The performance of the proposed approach has been compared with those based on Harris Hawks Optimizer (HHO), Henry gas solubility optimization (HGSO) and other optimization techniques in literature. The total energy loss has been significantly decreased in the case of integrating WT+Biomass units, due to generating reactive power. Besides, it observed that the proposed approach has better convergence characteristics than other techniques. In addition, according to the results integration of the PV, WT, PV+BES and WT+Biomass units into 33-bus system can reduce energy losses 49 %, 63%,64% and 94% compared to base case, respectively. For 69-bus system reduced to 52 %, 66%, 68% and 98% compared to base case, respectively.