## 3 Methodology of Solution
### Slave Stage  
The slave optimization stage in master-slave optimization algorithms corresponds to the core (i.e., heart) of the metaheuristic since this explores and exploits the solution space [22]. In the case of the optimization applied to the electrical networks, the slave stage is typically associated with the recursive evaluation of the power flow problem for each individual \(x_{j}^{t}\). In this research, we employ the power flow approach known as the successive approximation method initially reported in [23] for AC distribution networks, with the main advantage that its convergence can be guarantee through the application of the Banach fixed-point theorem [24]. The recursive power flow formula to solve power flow problems in AC distribution networks takes the structure presented in Equation (26) [25]:  
\[\mathbb{V}_{d,h}^{m+1}=\mathbb{V}_{dd}^{-1}\Big{[}\mathbf{diag}^{-1}\left( \mathbb{V}_{d,h}^{\star,m}\right)\left(\mathbb{S}_{pp,h}^{\star}-\mathbb{S}_{d,h }^{\star}\right)-\mathbb{Y}_{ds}\mathbb{V}_{s,h}\Big{]}, \tag{26}\]  
where \(m\) represents the counter of iterations, \(\mathbb{V}_{d,h}\) is the vector that has all the voltage variables in the demand nodes in complex form (note that for \(t=0\), we consider the plane voltages, i.e., \(\mathbb{V}_{d,h}^{t+1}=1\angle 0^{\circ}\)). \(\mathbb{Y}_{dd}\) and \(\mathbb{Y}_{ds}\) are sub-components of the general nodal admittance matrix that connect the demand and slack nodes among them. \(\mathbb{S}_{pp,h}^{\star}\) is a vector that contains all the complex power injections in the PV sources (where it is assumed that they are designed with unity power factor); \(\mathbb{V}_{s,h}\) represents the vector with the complex voltage outputs in the slack sources (i.e., known voltage values). It is worth mentioning that the operator \(\mathbf{diag}(\theta)\) transforms the vector \(\theta\) into a diagonal matrix; and \(\theta^{\star}\) is the complex conjugate value of the vector \(\theta\). The recursive evaluation of the power flow formula (26) ends when the convergence error is met between two consecutive voltage iterations. This criterion is defined in (27):  
\[\max\Big{\{}\Big{|}\Big{|}\mathbb{V}_{d,h}^{t+1}\Big{|}-\Big{|}\mathbb{V}_{d,h }^{t}\Big{|}\Big{|}\Big{\}}\leq\zeta, \tag{27}\]  
where \(\zeta\) represents the error of convergence which is assigned as \(1\times 10^{-10}\) based on the recommendation of [24]. Note that the NMA sends the combination of the output powers for the PV sources with their locations to the slave stage, i.e., \(\mathbb{S}_{pp}^{\star}\), which are considered inputs for the power flow formula (26). With this information, the fitness function that guides the exploration and exploitation process in the master optimization stage is calculated. Note that, to calculate the annualized energy purchasing costs at the substation buses, all the voltages provided by the solution of the power flow problem were evaluated in (28), i.e., \(\mathbb{V}_{d,h}\), which allowing to compute the total slack complex power generation as defined in Equation (28). \[\mathbb{S}_{cg,h}^{\star}=\mathbf{diag}\Big{(}\mathbb{V}_{s,h}^{\star}\Big{)} \left[\mathbb{Y}_{sd}\mathbb{V}_{d,h}+\mathbb{Y}_{ss}\mathbb{V}_{s,h}\right]. \tag{28}\]  
The evolution through the solution space for the NMA is guided by the value of the fitness objective function value, which corresponds to the summation of the objective function of the studied problem with some penalties associated with the constraints to ensure the solutions' feasibility.