#### Abstract

In this paper, an improved slime mould algorithm (ISMA) is proposed for finding the optimal location and sizing of photovoltaic distributed generation units (PVDGUs) in unbalanced distribution systems (UDSs). The proposed method is developed by changing the control variable update mechanism of the original slime mould algorithm (SMA). Total power losses on distribution lines and voltage deviation index of all buses are minimized under the consideration of various constraints. The location and sizing of PVDGUs found by ISMA are added into the cosimulation between MATLAB and OpenDSS for reaching other remaining parameters of UDSs. In addition, sunflower optimization (SFO), social-ski drive (SSD), cuckoo search algorithm (CSA), salp swarm algorithm (SSA), bonobo optimizer (BO), and SMA are also run for finding PVDGUs placement solution on the unbalanced three-phase IEEE 123-bus test feeder. As a result, the proposed ISMA can reduce the power loss up to 78.88% and cut the voltage deviation up to 1.4779 pu while that of others is from 69.10% to 78.87% and from 1.5759 to 1.4996 pu. Thus, ISMA should be used to place PVDGUs in UDSs meanwhile the cosimulation between MATLAB and OpenDSS should be applied as a power flow calculation tool for UDSs.

#### 1. Introduction

##### 1.1. Importance of DGUs in the Distribution System

In recent years, conventional power sources, including nuclear power plants, hydropower plants, and thermal power plants, brought terrible consequences to people and the Earth such as radiation exposure, floods, and global warming. So, the reduction of the use of these sources was implemented while renewable power sources with smaller capacity, which were known as distributed generation units (DGUs), were encouraged to be installed, especially in distribution systems.

The integration of DGUs in general and PVDGUs into distribution power networks is growing strongly with many benefits gained. Many studies have shown various benefits thanks to the integration of DGUs in the distribution system. Because of that, many countries have invested significantly in research to integrate DGUs into the electricity distribution systems to gain benefits like flexibility, reliability, environmental friendliness, economic benefits, etc [1]. Numerous studies have shown that selecting and connecting the suitable DGUs to the distribution system can achieve positive results such as reduced power losses, enhanced voltage profile, improved power quality, reduced fuel cost, and increased reliability and security of the system [2]. However, to maximize benefits, the determination of position and capacity for DGUs must be taken seriously [3]. Once the location and sizing of DGUs is predetermined incorrectly, it will cause many unwanted problems on the electrical system such as power losses, voltage sag, voltage flicker, and fault current [4]. Therefore, the determination of the optimal location and capacity of DGUs is necessary for maximizing both economic and technical benefits.

##### 1.2. Related Work

Most studies have applied popular optimization algorithms such as genetic algorithm (GA) and particle swarm optimization algorithm (PSO) to find the location and capacity of DGUs before connecting them to the distribution system. The main purpose of this is to reduce power loss on branches and improve the voltage profile [5–9]. The authors in [5, 6] suggested GA and [7–9] applied PSO for determining the integration of DGUs with objective functions of improving voltage stability, reducing losses with the lowest investment cost. Through the application of technology-economic factors, the suitable solution of DGUs is found with the smallest generated cost. GA and PSO are fairly simple methods in operation and are popular in solving optimal problems. However, it is easily trapped in the local search space as well as the low convergence rate. In addition, another optimization tool is called the artificial bee colony (ABC) algorithm is also proposed for finding the suitable DGUs with a major role in enhancing voltage stability and decreasing losses [10, 11]. This method has a quite stable convergence and it is built based on the principle of operation of the bee colony in finding food sources. Due to going through all three phases such as employed bee, onlooker bee, and scout bee phases, the convergence speed of this method is relatively slow. Similarly, another effective method, biogeography-based optimization (BBO), is also applied in searching for the optimal installation of DGUs. BBO is created based on the study of the distribution of biological species in habitats. Reference [12] used BBO to solve the optimal problem for the connection of DGUs with considering the multiobjective function. The results showed that BBO has many advantages over compared methods. Recently, a new nature-inspired optimization method is introduced. This method is developed based on sunflower’s motion and it is called the sunflower optimization algorithm (SFO) [13]. SFO is built on the natural phenomenon of the movement of sunflowers towards sunlight under consideration of pollination between adjacent sunflowers and this algorithm is stimulated by inverse square law radiation. The intensity of the radiation will depend on the distance between the sunflowers and the sun [14]. The effectiveness of the SFO method is clearly shown in [15]. The authors applied the SFO method for determining the suitable place and sizing of DGUs with the goal of decreasing power losses and increasing the voltage profile in the distribution systems. In addition, the salp swarm algorithm (SSA) is also receiving much attention recently. SSA is a metaheuristic optimization algorithm which is developed by the swarming behavior of salps in the oceans [16]. The authors in [17–19] used the SSA method for optimizing the installation of DGUs in the integrated systems. Thanks to the suitable integration of DGUs, total operating costs are minimized while power quality is significantly improved in the grid. Besides, a novel algorithm which is called social ski-driver (SSD) has been introduced. The characteristic of SSD is stochastic exploration like the roads where ski-drivers slip downhill. In this algorithm, there are four main parameters including the position of the agents, previous best position, mean global solution, and velocity of the agents. The number of parameters should be determined optimally in order to be most effective in finding the optimal solutions [20]. This method promises to be a powerful tool in solving the optimization problems. In addition, another optimization algorithm, which is named bonobo optimizer (BO), is recently published [21]. This algorithm is developed based on the social behavior and reproductive strategies of bonobos. The features of bonobos have been artificially modeled in mathematical formulas to solve optimization problems with high efficiency. Similarly, a new method for stochastic optimization is slime mould algorithm (SMA), is also introduced in [22]. SMA is inspired from the oscillation mode of slime mould in the nature. The properties of the slime mold have been constructed as a mathematical model with using adaptive weights to simulate the process of generating positive and negative feedback on a slime mold oscillation. This oscillation has been modeled to form an active method of searching and discovering food sources optimally [22, 23]. Several authors in [23] have also demonstrated the superior effectiveness of SMA through finding the optimal parameters for the solar photovoltaic system. The obtained results showed SMA is an efficient method with low data processing time in the solving optimization problems.

##### 1.3. Proposals, Novelties and Contributions

In terms of choosing the test system, most previous studies have ignored the consideration in unbalanced distribution systems. Determining the location and capacity of DGUs in an unbalanced distribution system is more complex than a balanced distribution system due to their effects on power flow analysis [24, 25]. To further clarify this issue, the authors in [26] conducted experimental simulations on both balanced and unbalanced distribution systems with and without DGUs. The obtained results showed that the current magnitude and voltage magnitude varied according to load levels and DGU’s penetration levels. Besides, the balanced and unbalanced distribution system’s power flows may differ significantly. Therefore, the choice of the system for simulation and analysis is very important in connecting the device/machine into the power grid. Moreover, the authors [27] also analyzed the effects of connecting DGUs in a three-phase unbalanced secondary distribution system by using the novel sensitivity analysis (SA) method. The results have shown that a dramatic change in the integrated distribution system comes with tremendous economic and technical benefits. For the same problem with [27], the authors in [28] proposed the PSO as an optimization method and the cosimulation framework as a calculation and data processing tool. In that paper, the authors have also mentioned the superiority of using the OpenDSS analysis support tool [29] in improving calculation speed and enhancing accuracy. Overall, the previous studies [30–32] mostly focused on the consideration of the multiobjective function. However, the effectiveness of the multiobjective function has not been clearly and accurately evaluated. The problem of trade-off in the multiobjective functions always happens and this can cause controversy in determining the best solution that meets all the stated criteria. Therefore, considering a single objective function will have more advantages in evaluation. In this paper, the single objective function is proposed to fairly evaluate the effectiveness of the proposed method and compared methods in determining the optimal location and sizing of DGUs, specifically photovoltaic distributed generation units (PVDGUs). As shown in [33], power losses depend greatly on factors such as the type of used conductors, transformers, and generators. in the distribution system. Hence, losses of the system can only be minimized by various techniques that cannot be completely removed and it is a leading factor in considering the economic and technical efficiency of the grid. Thus, total power loss minimization becomes one of the objective functions in this study. Besides, the enhancement of voltage stability is also paid much attention to a distribution system [34]. For solving this issue, the voltage improvement is considered as a single objective function in the network. In short, in this paper, the power loss reduction and the voltage improvement are proposed to be the two single objective functions under consideration of constraints such as voltage limits, current limits, and generated power limits of each PVDGU.

Generally, the optimization issues, especially the PVDGUs optimization, have been solved by the optimization methods. However, these methods are not effective and have low stability. Therefore, in this paper, an improved stochastic optimizer, called an improved slime mould algorithm (ISMA), is proposed to solve the optimization problem.

In summary, the main contributions in this paper can be listed as follows:(1)Improved slime mould algorithm (ISMA) with the application of two proposed modifications is developed. The first modification is implemented to improve the effectiveness of the second formula of the new solution updating process; meanwhile the second modification abandons the old selection condition of SMA and applies a new condition based on a fitness function. Thanks to the simultaneous combination of the two modifications, ISMA is superior to SMA in terms of search speed and solution quality.(2)Test system is recommended: Balanced distribution system and an unbalanced distribution system may have differences in voltage and current magnitudes. Therefore, to ensure the accuracy of power flow analysis, a complex unbalanced distribution system is proposed to be a test system, IEEE 123-bus test feeder.(3)Cosimulation tool is suggested: conventionally, the three-phase power flow can be solved using the full rectangular coordinate and the Newton–Raphson method. However, the computational speed is relatively slow. In addition to the Newton–Raphson method, forward backward sweep method is also a popular power flow method for distribution systems. The same disadvantage of the two methods when solving unbalanced systems is to run power flow three times for three single phases. This advantage leads to the time-consuming manner for the two methods. On the contrary, the computing speed of OpenDSS is very fast. In fact, OpenDSS does not rebuild the *Y* matrix during iterations and it is run one time only for three unbalanced phases simultaneously. In addition, convergence and accuracy in the OpenDSS are also great. Thus, MATLAB and OpenDSS are combined as a cosimulation tool to solve the optimal problem for improving accuracy and data processing speed. This cosimulation promises to be the best cosimulation tool for solving power flow problem.(4)Benefits of integrating the suitable PVDGUs into the distribution system are demonstrated: when PVDGUs are accordingly integrated into the system; it can help to reduce losses, improve voltage profile, enhance voltage stability, etc. Thus, the optimal location and capacity of the PVDGUs is proposed to maximize economic and technical benefits.

The structure in this paper is divided into 5 sections as follows: the objective functions and constraints are presented in problem formulation, Section 2. Next, two modifications in the second method and the condition for selecting the method in the equations of generating new solutions with simulation tools are introduced and analyzed in detail in the proposed method and simulation tools, Section 3. Then, the application of ISMA in finding suitable PVDGUs in an unbalanced distribution system is presented in the proposal of ISMA in optimizing the position and capacity of PVDGUs, Section 4. Besides, all obtained results and analysis as well as the expanded research fields are presented in simulation results, Section 5. Lastly, the summary and conclusions for the whole paper as well as exciting future works which continue this research are shown in the conclusions, Section 6.

#### 2. Problem Formulation

In this study, the location and capacity of PVDGUs is considered to connect to an unbalanced distribution system optimally. Losses and voltage profile are the two main factors, which commonly used to evaluate the power quality as well as the economic benefits for a distribution system. Thus, minimizing losses and improving voltage with strict constraints have become the two main goals of this paper. These constraints are in place to ensure that the proposed solution meets the technical criteria in a distribution system.

##### 2.1. Unbalanced Distribution Systems

For clarifying the characteristic of unbalanced distribution systems, a typical 123-bus system in Figure 1 is used for the analysis of loads, impedance, and voltage. The system is comprised of 123 buses, 124 lines, and 91 loads. Among the 91 loads, there are only 3 three-phase loads and 3 two-phase loads while the number of single-phase loads (*N*_{spl}) accounts for the highest rate with 85 loads. There are 2 balanced and 1 unbalanced three-phase loads among 3 three-phase loads and there are 2 unbalanced two-phase loads and 1 balanced two-phase loads among the 3 two-phase loads. The view on the number of loads indicates that there are three-phase, two-phase and one-phase distribution lines in the unbalanced system and the load of each phase at the same bus can be equal or totally different. The detail of the distribution lines and the loads at different buses can be clarified in the section.

Normally, a distribution line connecting Bus *a* and Bus *b* as shown in Figure 2 has three conductors with the same impedance and loads of phase *A*, phase *B*, and phase *C* at the receiving Bus *b*. The three series impedances are represented as in which each impedance has two parts, resistance and reactance shown in equation (1). Similarly, the load of phases is represented as at Bus *b*, and mathematically formulated in equation (2).

For unbalanced three-phase loads, there are some cases of the unbalanced power such as three unbalanced phases, and two balanced phases with one unbalanced phase. For the sake of easy understanding, the unbalanced three-phase loads at Bus *b* can be mathematically formulated as the following equation:

In practical, the demands of three phases are not requested for some cases while one-phase loads and two-phase loads need to be supplied. Consequently, the unbalanced distribution can be modeled in Figure 3.

**(a)**

**(b)**

As shown in Figure 3, Figure 3(a) shows a one-phase distribution line and Figure 3(b) shows a two-phase distribution line in which the one-phase distribution line supplies electricity to a single-phase load at phase A and the two-phase distribution line supplied electricity to two loads at phase A and phase B. The two figures can be mathematically formulated as follows:

Equation (4) shows the case with load at phase *A*. So, the loads at phase *B* and phase *C* are zero. On the contrary, Equation (5) can present the characteristic of both the balanced two-phase loads and the unbalanced two-phase loads. Loads at phase *A* and phase *B* are the same for the case with balanced two-phase loads whilst the loads are different for the case with unbalanced two-phase loads. About the impedance, the distribution lines can be represented as the following formulas:

In summary, the unbalanced distribution systems in practical are very complicated as shown in the figures and equations.

##### 2.2. Objective Function

As mentioned, this paper focuses on two objective functions including total power loss (TPL) and voltage deviation index (VDI). The applied mathematical equations for describing those two objective functions are shown as follows.

###### 2.2.1. Total Power Losses (TPL)

Power loss reduction in a distribution system brings tremendous benefits economically and technically. Therefore, connecting appropriate PVDGUs will be a great idea for minimizing losses. The mathematical equation of total power losses (TPL) is shown as follows:

###### 2.2.2. Voltage Deviation Index (VDI)

Connecting proper PVDGUs can improve voltage quality and enhance power grid reliability. So, the voltage deviation index should be considered like a main goal and its mathematical equation is shown as follows:

In equation (9), in pu is the average voltage of phases at the *b*^{th} bus. In the unbalanced distribution systems as explained in Section 2.1, some distribution lines have three phases but other ones may have one or two phases. So, the voltage of the end bus in the line should be the average value of the voltage of phases. For three cases with three phases, one phase, and two phases, the voltage average is calculated by the following equation:where and are voltage of the first phase and the second phase at Bus *b* in which the first phase and the second phase are two different phases among phases *A*, *B,* and *C* of distribution systems.

##### 2.3. Constraints

###### 2.3.1. The Power Flow Balance Constraints

In this study, losses at a higher frequency are ignored due to a very tiny current flow. Therefore, losses mainly occur at the fundamental frequency. All components are calculated and considered with the same the fundamental frequency and the power flow balance constraint has the following form [12]:

###### 2.3.2. The Voltage Limits

According to IEC Std. 50160, the lower and upper bounds should be kept from 0.90 pu to 1.10 pu, respectively. However, narrowing the voltage limits is better in enhancing voltage quality for low and medium voltage distribution systems, and the best range for applying the voltage limits is suggested to be from 0.95 pu to 1.05 pu [12]. Consequently, the lower bound and upper bound of voltage, which are represented by *V*_{min} and *V*_{max}, are selected to be 0.95 and 1.05 pu in the study. As suggested in [35, 36], the voltage of buses should be constrained by the following equation:

However, it should be noted that in the unbalanced distribution system, an example bus can have three values for voltage because of the presence of unbalanced loads. Hence, the constraint is modified as follows:

Here and are determined by the following equation:

###### 2.3.3. Generation Limits of Photovoltaic Distributed Generation Units

The total capacity of proposed PVDGUs should not exceed the total load demand and the capacity limit of each PVDGU is constrained as follows [28]:

###### 2.3.4. The Branch Current Limits

The thermal capacity must not exceed the allowable limit of each conductor. Thus, the current limit is defined in the following model [9]:

Solving the power flow problem is one of the important steps in the process of optimizing location and size of PVDGUs in unbalance distribution networks. To evaluate the fitness function of each solution, the power flow solutions need to be evaluated to obtain parameters such as the branches’ power loss, the buses’ voltage, and the branches’ current magnitude.

#### 3. The Proposed Method

In 2020, Li et al. developed Slime mold algorithm (SMA) and showed its outstanding performance over many previous metaheuristic algorithms [22]. SMA reached results as expected for a set of complicated benchmark functions such as global optimums, stability, and convergence speed. However, these functions are not difficult enough because they were not comprised of complicated constraints and different search spaces for different control variables. Thus, SMA still has big limitations for such complicated problems and the placement of DGUs in unbalance distribution systems is one example. In the section, the existing shortcomings of SMA are analyzed and then modified to form an improved version, called the improved slime mold algorithm (ISMA). To implement ISMA for the problem, it should be simulated by using OpenDSS and MATLAB software. So, the connection between OpenDSS and MATLAB is also expressed in the section.

##### 3.1. Conventional Slime Mould Algorithm

Slime mould algorithm (SMA) has been developed recently and proved to be more effective than a huge number of popular and famous metaheuristic algorithms [22]. Each particle in the population of SMA is represented as a solution Sol_{m} where *m* is from 1 to *N*_{ps}. Sol_{m} in the population is within the lower bound solution (LB) and upper bound solution (UB), and randomized as follows:

SMA has one solution update technique and it updates the solution once for each iteration based on the following equation:

In the formula, Sol_{rd1} and Sol_{rd2} are the two randomly selected solutions in the current population. *V*_{A}, *V*_{B}, and Weight are vectors with the same dimension as each solution *m*. *V*_{A}, *V*_{B}, and Weight are formulated as follows [22]:

In the computational formula of Weight, the Condition indicates that represents the ranking of the first half of the population. *rd*_{1} and *rd*_{2} are vectors with the same dimension as solution *m* in which *rd*_{1} and *rd*_{2} denote the random numbers in the interval of [0, 1] and [−1, 0], respectively.

In addition, CD in equation (19) is obtained by the following equation:

After producing new solution by using (19), the new solution is verified and fixed by using the following equation:

Fitness function of the *m* new solution is calculated and represented by . Among values (where *m* = 1, …, *N*_{ps}), the lowest fitness value is set to *Ft*_{best} and the highest fitness value is set to *Ft*_{worst}. Then, the selection technique is implemented to retain better solutions and abandon worse solutions. The process of selection is done by the following model:

In the last step, the so-fart best fitness is determined by using the lowest fitness value of (*m* = 1, …, *N*_{ps}) and the so-fart best solution Sol_{sfbest} is obtained accordingly.

In summary, the application of the applied SMA can be described in the following steps. Step 1: Select values for *N*_{ps} and *M*_{Iter} Step 2: Employ (18) to get the initial population Step 3: Calculate fitness function *Ft*_{m} Step 4: Determine *Ft*_{best} and *Ft*_{worst}, Set *C*_{Iter} = 1 Step 5: Employ (23) and (24) to obtained *A* and *B* Step 6: Randomly produce *V*_{A} and *V*_{B} satisfying (20) and (21) Step 7: Employ (22) to obtain Weight Step 8: Employ (25) to calculate CD Step 9: Employ (19) to generate new solution Step 10: Use (26) to verify and fix Step 11: Calculate and determine *Ft*_{best} and *Ft*_{worst} from Step 12: Use (27) and (28) to get and Step 13: Select the lowest fitness value from *Ft*_{m} and determine Sol_{sfbest} Step 14: if *C*_{Iter} and *M*_{Iter} are equal, stop implementing SMA and reporting results. Otherwise, set *C*_{Iter} to *C*_{Iter} + 1 and back to Step 5.

###### 3.1.1. Improved Slime Mould Algorithm

As shown in equation (19), SMA uses two different ways for producing new solutions. The formula of the equation is called the first method and the following one is named the second method. The advantages and disadvantages of the two methods are analyzed as follows:(1)The first method searches around the so-far best solution with a distance by using . The first method is the effective creation of SMA through applying the fitness function (i.e., Weight factor in equation (19)) and the creation supported SMA to have a better search method than other metaheuristic algorithms.(2)The second method modifies the old solutions by using vector *V*_{B} and the second method using () for updating the new solution is not a potential search strategy. In fact, *V*_{B} is a randomly produced vector within −*B* and *B* as shown in equation (21) while *B* has a range from 0 (for the case ) to 1 for the case as shown in equation (24). From equations (21) and (24), *V*_{B} is a vector with terms from −1 to 1. As a result, () can be from (−Sol_{m}) to (Sol_{m}) and the new solution is equal to from (−Sol_{m}) to (Sol_{m}). Clearly, if is around Sol_{m}, it may be a potential new solution. But for the case that is close to 0 or close to (−Sol_{m}), the new solution is certainly low quality solution. Furthermore, the new solution is always smaller than the old one Sol_{m} and there are not any demonstrations that promising solutions have to be smaller than the old one. The noted issue in the second method is that B is approaching 0 when the current iteration is approaching the maximum iteration (i.e., => ) and *V*_{B} is also approaching 0 as a result. Thus, the new solution will become a Zero solution with 0 value for all control variables. In [22], benchmark functions were solved for reaching the minimum value and the global solutions for approximately all functions are Zero solution. The global solutions and the produced solution by the second method are unintentionally the same and SMA becomes a very strong method.

In this paper, control variables have a large range and the global optimum is certainly not the zero solution. This fact indicates the second method using () to update new solutions is ineffective and it should be replaced with another method. In this section, a proposed method is developed to replace the second method in equation (19) and it is presented as follows:

Equation (29) is the first modification of the proposed ISMA. Currently, ISMA has two different methods for producing new solutions in which the first method is taken from the conventional SMA and the second method is taken from the first modification. However, the condition for using either the first method or the second method should not be taken from SMA. As shown in equation (19), the condition is where CD is calculated by using equation (25), which is . For the sake of simplicity, the result of () is represented by *x* and the result of () is represented by *y*. As showing the values of *x*, CD, and *y*, *x* is set the range of [0, 5] for a better view of the three factors. It is noted that can be much higher than 5; however, we limit the maximum value at 5 only because CD at *x* = 5 is approximately equal to 1 and CD is still one if *x* > 5 is set. On the contrary to *x* and CD, the values of *y* is either 0 (for ) or 1 (for ). After producing the matrix with random values rdn*, y* is obtained by comparing rdn and CD. As a result, Figure 4 is plotted to show the values of *x*, CD, and *y*.

Observing Figure 4 seen that most values of *y* are 1 and few values are 0. For a better view of the influence of on *y*, is selected to be from 5 to 1000 and Figure 5 is plotted. It is seen that *y* is equal to 1 for all values of *x*. Hence, it concludes that *y* is equal to 1 for approximately all values of , and the first method in equation (19) is used for approximately all solutions over the search process with *M*_{Iter} iterations whereas the second method is hardly ever used. In the second modification, we propose to change the condition of using the comparison between rdn and CD by using a measured condition. The first method is to search around the best solution, so it should be applied for the old solutions with low quality. The second method is to exploit around the current old solutions, so the around searched solutions should be promising solutions and the search space around the solutions are possible to be the presence of higher quality. So the condition for selecting the first or the second method is as follows:where *Ft*_{mean} is the mean of fitness values of the whole population.

As a result, the technique for producing the new solutions of ISMA is summarized as follows:

Generally, the first method in equation (31) just makes a jump for finding new solutions around the best solution (). However, when the quality classification condition of solutions is fulfilled (), then the second method is applied. This second method contains () to produce one more jumping step. So, two generated jumps in the second method are to search for new solutions far from current good quality solutions (). This facilitates the minimization of falling into the local optimal region and opens opportunities to discover new solutions in the larger possible search area. In summary, in order to strike a balance between exploiting new solutions around the best solution and finding new solutions in the feasible search space, Equation (31) should be applied under the selection condition for using the first method or the second method. This combination with reasonably selected condition promises to create the great opportunities in searching the better potential solutions than the original method.

##### 3.2. OpenDSS Software and Cosimulation between OpenDSS and MATLAB

###### 3.2.1. OpenDSS Software

OpenDSS is the distribution system simulator (DSS), which is an open source version used to support the calculation and simulation for distribution systems in the frequency domain. This is developed by the Electric Power Research Institute (EPRI) to address the gaps of other distribution system analysis tools. Besides, it is built to create a flexible and reliable research platform that aims to serve distribution system analysis applications, especially the distribution system with the integration of distributed generation sources as a synchronous generator, photovoltaic/wind turbine generator, capacitor, etc. The general structure of this software can be briefly presented as shown in Figure 6 [37].

For this OpenDSS’s structure, it allows for importing user-written code into the software. The environment in OpenDSS contains script files for circuit definition with the transformer specification, line data, load data, etc. These script files are appropriately defined by the user. Besides, OpenDSS is developed to allow interaction with other popular software such as MATLAB, PYTHON, VBA, C+, and C#. This interaction is considered as coordinated simulation and this is implemented through the windows component objective model (COM). In other words, OpenDSS allows script development or modification from other software to control the OpenDSS model via this COM server.

###### 3.2.2. Cosimulation between OpenDSS and MATLAB

OpenDSSS can be driven by other software via the COM server as mentioned above. In this study, OpenDSS has interacted with MATLAB for coordinate simulation to perform the analysis of PVDGU impacts in the distribution system. The coordination process between OpenDSS and MATLAB is briefly summarized as Figure 7 [38].

In this study, a hybrid interface is built for necessary data transmission between two software by using the COM server. OpenDSS is limited in solving the problem that involves changing the control variable. Therefore, in this work, OpenDSS needs to coordinate with MATLAB to solve the considered problem in the distribution system. OpenDSS is a powerful tool for calculating the power flow in large distributed systems with fairly fast simulation times and using MATLAB to control the operations of OpenDSS is a great hybrid interface. For this particular case, the grid data is described in the OpenDSS. The power flow is calculated and the obtained results are sent out by OpenDSS. The required values are transferred to the MATLAB through the COM server. Besides, the written optimization algorithm in MATLAB is responsible for computing and proposing the suitable solution in distributed resource integration into the considered grid. This proposed solution is transferred to OpenDSS via the COM server for recalculating the power flow. This process is implemented at each cycle until the loop requirement is satisfied.

#### 4. The Implementation of the Applied Method for the Problem

##### 4.1. Improved Slime Mould Algorithm

In this paper, the implemented methods are in charge of finding the most suitable location and the best generation for installed photovoltaic systems. Hence, each solution *m* among the set of solutions contains the two main factors, which are represented by and (where *k* = 1, …, *N*_{DG} and *m* = 1, …, *N*_{ps}). The solution *m* Sol_{m} is formulated mathematically and randomly produced as follows:where LB and UB are the minimum solution and maximum solution containing the minimum values and maximum values of the location and power generation of distribution systems. The two solutions should be defined before implementing the search process of the applied methods as the two following formulas:

As shown in the boundary solution, and are the minimum location and maximum location in the unbalanced distribution system corresponding to Bus 2 and Bus *N*^{b}. and are the minimum rated power and maximum rated power of each installed PVDGU. Besides, the dependent variables of the problem are comprised of , , and , and , , and which are the results from running OpenDSS software.

##### 4.2. OpenDSS for Solving the Power Flow of Unbalanced Distribution System

This simulation tool is designed to simulate most analysis types of the distribution planning. It is capable of quickly and accurately solving problems related to power system analysis such as power flow, dynamic, harmonics, and short circuit calculations. The application of DSS to solve the problem of power flow analysis can be summarized briefly based on the process of building and calculating the matrix for the elements in the electrical system as Figure 8 [37].

With the initial bus data and line data from the power distribution system, a base matrix *Y* (or 1/*Z*) is created for each element. After each element’s matrix is built, the overall system admittance matrix (*Y*) is generated by aggregating all the matrices of the element and solving them through a sparse matrix solver. The relationship between the phase angle and voltage magnitude is appropriately established after all power delivery elements are maintained interconnection to the considered power system. Power conversion elements inject current and these currents are added appropriately to their place in the vector matrix (*I*_{inj}). The bus voltage and the feeder current vectors for each phase in the distribution system can be presented as follows:

Subsequently, the sparse matrix will be solved to find out the phase angle and amplitude values of all voltages at each phase in the system. The voltage value at each bus is converged to the most accurate value through each iteration and stops this if the tolerance is satisfied. This convergence is processed quite quickly because the *Y*-matrix system is built only once in the early stage and it is not rebuilt during processing. This has contributed to reduce the processing time in the power flow analysis [38]. Besides, many researchers have also analyzed and mentioned the great accuracy of OpenDSS in solving problems in distributed systems [39]. Thereby, OpenDSS is really the best tool in solving power flow problems.

##### 4.3. Fitness Function

The fitness function of each solution *m* is represented by *FF*_{m} and it is used to measure the effectiveness of the determined location and rated power of PVDGUs installed in the unbalanced systems. Before calculating the fitness value, objective function and penalty terms must be obtained. The study considers two single objectives, including total power loss and voltage deviation index, which are obtained by employing equations (7) and (9), respectively. On the other hand, penalty terms can be calculated by comparing dependent variables and their limits as follows:where and are the penalty terms for the violation of current in the *f*^{th} feeder and the violation of voltage at Bus *b*.

Finally, fitness functions for reduction of total power loss and reduction of voltage deviation index are, respectively, determined by the following equation:where and are the total power loss and voltage deviation index of the *m*^{th} solution, respectively.

##### 4.4. Correction of Control Variables

The proposed method generates new solutions where each new solution is comprised of the location and capacity of each PVDGU. To ensure that the new solution is always within the search limits, the position and capacity of each PVDGU is checked and modified appropriately by applying two equations as follows:

##### 4.5. Termination Condition for Iterative Algorithm

To ensure objective comparability between the proposed method and other methods, 50 trial runs are implemented for all methods. The best solution in each trial run is saved for comparing and evaluating the effectiveness of each method. In each trial run, the maximum number of iterations is considered and appropriately selected to ensure the execution methods have converged completely.

##### 4.6. The Cosimulation between OpenDSS and MATLAB

The implemented cosimulation between OpenDSS and MATLAB for determining the optimal position and capacity of PVDGUs in an unbalanced distribution system is briefly presented by the flowchart as Figure 9.

This implementation is carried out in two different environments. In the OpenDSS, a program is written by the user for describing the entire distribution system. It includes main circuit script and element scripts such as line modeling, load modeling, and transformer modeling. OpenDSS will be run from the main circuit script for data processing and computation to solve the power flow problem. The results obtained from this calculation (phase angle, voltage magnitude, power losses on each branch, etc) will be transferred to MATLAB via COM server for processing as MATLAB input data. MATLAB will conduct analysis and generate the possible solutions of location and sizing of PVDGUs in the search area. Proposed solutions will be sent back to OpenDSS for integrating and processing additional data as OpenDSS input data at the first iteration. The received solution will be properly added to the system’s script and the power flow is resolved by using a power flow solver in OpenDSS with the integration of PVDGUs in the distribution system. Obtained results as OpenDSS output data will be transmitted through the COM server to MATLAB to process, calculate and evaluate solution quality. The optimization algorithm that is written in MATLAB will rely on the quality evaluation of existing solutions to suggest better solutions. These new solutions will also be passed over to OpenDSS to handle data at the next iteration. This process is repeated until the optimization algorithm’s criterion is satisfied, i.e., the current iteration is equal to the maximum iteration. After each loop is finished, the next solution with better quality than the previous one will be updated as the best current solution and convergence will occur. As a result, the most suitable solution for the position and capacity of PVDGUs in the distribution system will be found.

#### 5. Simulation Results

In this paper, the proposed method (ISMA) and six other methods (SFO, SSD, CSA, SSA, BO, and SMA) are implemented for finding the optimal location and capacity of three PVDGUs in IEEE 123-bus test feeder as Figure 1. The used data of the system is taken from IEEE PES [40]. For an accurate and objective evaluation, all methods are simulated in 50 trial runs with randomly generated initial solutions. MATLAB and OpenDSS are two used software and this simulation is done on a personal computer with a 4.0 GB RAM, 1.80 GHz CPU, and Intel Core i5. For all implemented methods, the main parameters have been surveyed and selected appropriately to have a fair assessment from the obtained simulation results. The total created generation numbers at each loop are surveyed from 20 to 40 with a step size of 5 and the suitable result is 30. Besides, the total iteration numbers are also chosen enough to guarantee the complete convergence of algorithms with 200 loops. Other remaining parameters of SFO and SSA, SSD, CSA, and BO are taken from [20, 21, 41, 42], respectively. Moreover, for building the unbalanced distribution systems, the model of circuit elements in OpenDSS has been described in detail in [37]. As mentioned, OpenDSS is applied as software to solve the power flow problem in the unbalanced system. This software is interacted with MATLAB software through the Component Object Model server Dynamic-link library for cosimulation [43, 44]. In addition, two single objective functions are considered independently, including of minimizing the total power losses on transmission lines and improving the voltage deviation index in the distribution system while the allowable limits of the load voltage and branch current are satisfied.

##### 5.1. Case 1: Minimizing the Total Power Losses

Power loss is one of the important factors that strongly influence the evaluation of the system quality as well as the economic benefits. Therefore, the reduction of losses on the transmission lines is a necessary matter that should be studied carefully. In this case, the total power losses are considered for reducing to a minimum under consideration of load voltage and branch current limits. Besides, in order to ensure the objectivity of the proposed method and the compared methods, the results of the 50 trial runs are collected for analysis and evaluation.

The values of the best fitness, the worst fitness, and the average fitness of the proposed method are compared in detail with others such as SFO, SSD, CSA, SSA, BO, and SMA as shown in Table 1. For the assessment between the improved method and the original method, the best fitness, the worst fitness, and the average fitness values for ISMA are 20.2221, 22.4978, and 20.7153 while for SMA are 20.2352, 23.4059, and 20.9881, respectively. These values indicate that the performance of ISMA is better than SMA. In other words, the modifications on the new control variable updating process of ISMA have reached a positive effect in enhancing performance. To clearly justify the effectiveness, the improvement level (IL) was calculated and the last three columns show the IL of the proposed ISMA over others for the mean fitness, worst fitness, and the best fitness. The three values are, respectively, 1.2998%, 3.8798%, and 0.0647% as comparing results with SMA. Similarly, ISMA is also compared with other popular methods. As obtained results, the best fitness and the worst fitness values of ISMA are better than the five remaining methods from 20.2340 to 21.7190 and from 22.6942 to 25.5408 with the levels of improvement from 0.0588% to 6.8921% and 0.8654% to 11.9143%. This shows that the solutions with the best quality and the worst quality proposed within the 50 random runs from ISMA are always more optimal and efficient than other solutions. Moreover, in order to demonstrate the high stability of the proposed method, the calculation of the average value of the fitness function is performed. The average fitness value of 50 trial runs from ISMA is the smallest as compared to other popular methods from 20.9122 to 23.2920 with the corresponding improvement level from 0.9416% to 11.0626%. This shows that the proposed method is more stable than the other methods in finding the optimal solution in the same search area.

Figure 10 shows the fitness values of all implemented methods during the 50 trial runs with the initial solutions are randomly generated. For convenience in reviewing and evaluating the stability of the implemented methods, all fitness values are rearranged from low to high. As a result, the curve of ISMA has proved the effectiveness and stability of this method compared to other methods with the same iteration number and created generation number in each iteration. In addition, the convergence characteristic of ISMA is also better than other methods like Figure 11 presented. During the 200 iterations of the best solution, ISMA exhibited convergence fairly quickly; it found the better quality solution than others and the best solution at the 138^{th} iteration. This shows the outstanding properties of the proposed method in the process of finding the feasible solution for solving the optimization problem.

Like presented in Figure 12, by using the proposed solution from ISMA, the power losses on transmission lines are reduced as compared with the original system. Specifically, the total losses on the system are significantly reduced from 95.77 kW to 20.22 kW, corresponding to 78.88% in the loss reduction after installed suitable PVDGUs and the obtained loss of ISMA is smallest as compared to other methods. This shows that the proper installation of PVDGUs has had a strong impact on reducing losses on the transmission lines. In other words, the benefit from loss reduction depends on the quality of the optimal solution that has been proposed. All of the above-given arguments prove that ISMA is really an efficient and suitable method for solving the problem of the optimal place and sizing of PVDGUs in the unbalanced multiphase distribution system.

The best solutions of PVDGUs placement found by previous and implemented methods are presented in Table 2. As demonstrated from simulation results, the power loss reduction entirely depends on the optimal solutions that the implemented methods have proposed. Three methods (ISMA, SMA, and CSA) proposed Buses 47, 65, and 72, and they could achieve the relatively same good results for loss reduction, which is 78.88% for ISMA, 78.87% for SMA, and 78.80% for CSA. Besides, another option consisting of Buses 47, 65, and 76 found by BO could also lead to good results like the three above-given locations. BO achieved the loss reduction of 78.87%, which is approximately equal to that of the best method, 78.88%. Thereby, it shows that the two key locations to maximize the benefits in reducing power loss are Buses 47 and 65, and one remaining position is Bus 72 or Bus 76. In addition, PSO suggested the installation locations at Buses 47, 67, and 72, but it got the worst loss reduction with only 69.10%. The cause is that PSO did not properly determine good capacity for each PVDGU. In fact, BO, SMA, and ISMA proposed approximately the same total capacity, which are respectively 2809.0 kW, 2807.7 kW, and 2808.9 kW, and the capacity of each PVDGU from the three algorithms are nearly equal. On the contrary, the total capacity found by PSO was 2940.0 kW and the capacity for each PVDGU is also different from that of BO, SMA, and ISMA. Clearly, the capacity of each PVDGU obtained by PSO was not effective as others. Therefore, determining the suitable capacity for each position is extremely important and it greatly contributes to high loss reduction to distribution systems. In summary, the three best locations for this case should be Buses 47, 65, and 72, and the best capacities for them are, respectively, 907.2 kW, 341.1 kW, and 1560.6 kW.

##### 5.2. Case 2: Minimizing the Voltage Deviation Index

In this case, the enhancement of the voltage deviation index is considered as an objective function with under the constraints of the load voltage profile and the current on the branches in the test system. In evaluating a distribution system, the voltage deviation index becomes an important criterion and it is the basis for analysis and evaluation. The consideration for improving this index will drastically enhance the quality of the voltage for the entire system. As mentioned, the proposed method is compared to six other effective methods with 50 trial runs.

The best, the worst, and the average values of the fitness function from the proposed method and other popular methods are compared in detail in Table 3. The best fitness, the worst fitness, and the average fitness values of ISMA are 1.4779, 1.6322, and 1.5523, respectively, meanwhile those values of SMA are higher and equal to 1.5036, 1.6749, and 1.5726. From the three values, it is obtained that the levels of improvement are 1.7092% for the best fitness, 2.5494% for the worst fitness, and 1.2909% for the mean fitness. The best fitness and the worst fitness of the other five remaining are, respectively, from 1.4974 to 1.5759, and 1.6327 to 1.7031. The proposed ISMA can reach the levels of improvement from 1.3023% to 6.2187% for the best fitness and from 0.0306% to 4.1630% for the worst fitness. In addition, the average fitness over 50 trial runs indicates that the proposed ISMA can reach better stability from 0.4745% to 5.0058%. In summary, the proposed ISMA is very successful in finding solutions for the problem, and it has high improvement levels over six other implemented algorithms.

Figure 13 shows the fitness function values of the proposed and compared methods in 50 trial runs, and all fitness values are rearranged from low to high for the convenience of evaluation. Most of the points of fitness values on the ISMA curve are lower than those of the original SMA and other compared methods. This represents the stability level of the methods at random trial runs. The obtained results from Figure 13 also showed that ISMA have better solution quality and stability as compared with the implemented methods. Besides, the convergence of seven methods in the best solution is also clearly presented in Figure 14. ISMA found better quality solutions than others at the 100^{nd} iteration and the best solution is also determined at 123^{th} iteration. Thereby, ISMA is a method with quite good convergence and it is an effective method in finding the optimal solution.

From obtained results, the voltage values of all phases before and after installation of PVDGUs for Case 1 and Case 2 are plotted according to distance from the substation as shown in Figures 15–17, respectively. In other words, the *y*-axis is the values of voltage amplitude (pu) in phases such as phase *A*, phase *B* and phase *C*. These voltage values are distributed in the *x*-axis by distance (km) calculated from the substation. Since the considered system is a multiphase system, the voltage values of the phases for each bus are different and they are clearly denoted in Figures 15–17. Due to the characteristics of the transmission line connection and the distribution of the load concentrated in phase *A* of the system, the system has a higher number of bus voltages in phase *A* than in phase *B* and phase *C*. As shown, thanks to the application of the proposed solution from ISMA, the minimum value of rms voltage is significantly improved from 0.979 pu to 0.989 pu for Case 1 and from 0.979 pu to 0.990 pu for Case 2. Not only that, the maximum value of rms voltage before and after the installing PVDGUs also changed positively from 1.037 pu to 1.030 pu approximately for Case 2 and that value has no improvement after connecting PVDGUs in Case 1.

Besides, the voltage deviation of the system has been greatly enhanced as Figure 18 after optimally integrating the position and capacity of the PVDGUs into the system. The smaller this deviation, the voltage quality of the system will be more stable. Specifically, the sum of voltage deviation values for all buses without PVDGUs connection is 2.2035 pu. However, after the integration of three suitable PVDGUs into this system with the two single goals of reducing the total power loss and minimizing the voltage deviation index, the total voltage deviation values at the buses drop to 1.8523 pu and 1.4779 pu, respectively. As a result, Figure 18 shows the concrete evidence of the positive change of voltage deviation at each bus in the distribution system. Generally, the voltage deviation of the buses with suitable PVDGU’s connection in the second objective function is better than in the other two cases of no PVDGU connection and PVDGU’s connection in the first objective function. This shows that, when properly connected PVDGUs into the system, the voltage deviation values are significantly reduced and this contributes to enhance voltage quality for an unbalanced multiphase distribution system. In other words, the well-being of the voltage depends strongly on the optimization of the solution as well as on the efficiency of the proposed method. Like all things proved, ISMA provides the optimal solution for the place and sizing of PVDGUs in the unbalanced multiphase distribution system with higher quality and more stability than other methods under the same conditions.

The most appropriate position and capacity of PVDGUs of all implemented methods are presented in Table 4. As shown in the table, ISMA is the best method in improving the voltage profile with the lowest VDI of 1.4779 by proposing Buses 46, 67, and 106, whereas the second-best method, BO reaches the VDI of 1.4974 by proposing Buses 22, 46 and 67. Although the two methods did not use the same locations, the total capacity of the three PVDGUs is approximately the same, 2767.1 kW for BO and 2767.0 for ISMA. Focusing on the locations and total capacity of the third-best method (SMA) and the worst method (SSA) indicates that different locations and different capacity were used. SMA used Buses 38, 45, and 67 and the total capacity of 2777.4 kW. SSA used Buses 48, 69, and 97, and capacity of 2735.0 kW. Among the three proposed locations of other remaining methods, only one or two buses are identical to those of ISMA. On the other hand, the total capacity from these methods is also much different from that of ISMA. In summary, the best position and the best capacity of three PVDGUs for reaching the best voltage profile are Buses 46, 67, and 106, and 813.8 kW, 1900.0 kW, and 53.2 kW.

##### 5.3. Discussion on the Objective Functions and Expanded Research Fields

In this study, we suggested the use of two single objective functions in the two cases for determining the optimal location and sizing of photovoltaic distributed generation units. The reduction of the total power loss on the branches is an important criterion regarding the economic problem of the distribution system. The reduced loss will contribute greatly in saving the operating cost of the system [45]. Thus, Case 1 focuses on minimizing power losses (TPL). In addition, the phase voltages at each bus of the system in this case have also been constrained within the allowable limits to avoid overvoltage by high penetration of PVDGUs. In addition, to improve the voltage quality, another important factor, which is the voltage deviation index (VDI), is also considered [46]. Hence, this index minimization is assigned as the goal of Case 2.

Figures 15–17 indicate the voltage profiles of the based system without PVDGUs and hybrid systems with PVDGUs for Case 1 and Case 2, respectively. Although all cases have voltages within the permissible limits of 0.95 pu and 1.05 pu, the voltage range between the cases differs significantly. In detail, the phase voltages at buses are in the range of [0.979, 1.037] pu for base system without PVDGUs, while this range is [0.989 1.037] pu for Case 1 and [0.990, 1.030] pu for Case 2. Obviously, Case 2 has a better voltage profile than others because all values are close to 1.0 pu. On the other hand, Figure 18 sees the voltage deviation of Case 2 is the smallest. On the other hand, the sum of voltage deviation values for all buses for the base case is 2.2035 pu, but the result is much smaller in Case 1 and Case 2, which are, respectively, 1.8523 pu and 1.4779 pu. In summary, different objective functions can result in different achievements and we can choose the most appropriate objective function for our purpose. If we are interested in the economic aspect and ignore voltage quality improvement, Case 1 will be a more suitable choice than the remaining cases. On the contrary, if we need a very good voltage profile, Case 2 should be applied. This indicates the advantage of the single objective function in facilitating the appropriate selection by purpose.

On the other hand, the two single objectives can be reached simultaneously by using a multiobjective function [47], but it cannot reach a solution with the same total power loss and the same voltage profile as Case 1 and Case 2. Normally, there are two cases for the solution with the two objectives. In the first case, total power loss is better, but voltage is worse. In the second case, the results are the opposite. To solve the multiobjective problem, two weight factors are associated with each objective and higher values of the weight factor are set for a more important objective [2]. The two weight factors are constrained to be within 0 and 1 and their sum should be always 1.0. For the problem, the application of ISMA or other algorithms is like the study with two single objectives.

#### 6. Conclusions

This paper proposes ISMA for determining the optimal location of PVDGUs in the unbalanced distribution system, IEEE 123-bus test feeder for minimizing the power loss and improving the voltage quality. Besides, the study also proposed a cosimulation combining MATLAB and OpenDSS for solving the considered problem. This combination offers many advantages in improving the speed of data processing and enhancing the accuracy of power flow computation in the unbalanced distribution system. Confidently, it is the best cosimulation tool for solving the power flow problems. In addition, by applying the solutions for optimal installation of PVDGUs, the power loss and voltage deviation index in the system have dropped significantly to 20.222 kW (corresponding to 78.88%) and 1.4779 pu for ISMA, while these values are 20.235 kW (corresponding to 78.87%) and 1.5036 pu for SMA, respectively. Thereby, the improvement in ISMA is more efficient than the original SMA. The proposed method is also better than other implemented methods such as SFO, SSD, CSA, SSA, and BO, and other published methods such as PSO and JAYA. Like the obtained results, the solution of ISMA is better than others from 69.10% to 78.87% in the power loss reduction and from 1.5759 pu to 1.4996 pu in the minimization of voltage deviation index. In short, ISMA is really effective in maximizing economic and technical welfare.

As shown previously, the study has significant contributions to unbalanced distribution systems for loss reduction, voltage improvement, and voltage deviation reduction between phases; however, there are still serious problems existing in the unbalanced distribution systems such as surplus power of PVDGUs based on renewable energies for the cases of low load demand and a high deviation between real and estimated wind speed and solar radiation. So, the placement of the battery energy store system (BESS) and smart inverters in the system can be the upcoming directions of the study. The energy storage function of BESS and the smart function of the inverter can lead to high benefits without causing adverse effects. The electric components can be located at the most suitable nodes by using the proposed ISMA or a newly developed algorithm with higher performance.

#### Nomenclature

TPL: | Total power losses |

N_{f} : | The number of feeders in the system |

: | The active power loss of the f^{th} feeder in the system |

VDI: | Voltage deviation index |

N^{b}: | The number of buses in the system |

: | The i^{th} power load in the system |

: | The capacity of the k^{th} PVDGU |

: | The power is supplied by grid |

: | The maximum number of iteration and the current iteration |

N_{ps}: | The population size or number of solutions in the sole solution set |

V^{b}: | The b^{th} bus voltage in the system |

: | The rated current magnitude of the f^{th} feeder |

: | The current magnitude of the f^{th} feeder |

: | The randomly generated integer numbers in the range from 1 to number of search agents |

: | The randomly generated integer numbers in the range from 1 to number of search agents |

V_{max}, V_{min}: | The minimum and maximum voltage of the system |

: | The number of control variables |

N_{L}: | The number of loads |

: | The current magnitude of phase A, B and C at the f^{th} feeder |

: | The resistance of phase A, B and C at the f^{th} feeder |

N_{spl}: | The number of single phase loads |

N_{DG}: | The number of PVDGUs |

: | The minimum and maximum voltage of the b^{th} bus |

: | The minimum and maximum of PVDGU’s capacity |

: | The resistance of phase A, B and C |

: | The reactance of phase A, B and C |

: | The active power of phase A, B and C at the b^{th} bus |

: | The reactive power of phase A, B and C at the b^{th} bus |

: | The location and capacity of the k^{th} PVDGU corresponding to the m^{th} solution |

: | The real part of voltage of phase A, B and C at the b^{th} bus |

: | The imaginary part of voltage of phase A, B and C at the b^{th} bus |

FF_{A}, FF_{B}: | The fitness value of objective functions |

Sol_{sfbest}: | The best solution |

Ft_{best}, Ft_{worst}, Ft_{m}: | Fitness function of the best, worst and the m^{th} solutions |

rdn: | A random number within 0 and 1. |

#### Data Availability

The data used to support the study are included in the paper.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by Foundation for Science and Technology Development of Ton Duc Thang University (FOSTECT), website: https://fostect.tdtu.edu.vn, under Grant no. FOSTECT.2022.10.