## 3 Implementation of the Harmony Search Based Optimization Algorithm  
We use a recently developed meta-heuristic method based on the improvisation process of the jazz musicians: harmony search algorithm. We integrated the developed harmony search based optimization module into the experiment model.  
The implementation of the harmony search optimization algorithm is given as follows.  
1. The first step includes selection of the parameters of the harmony search algorithm: harmony memory size (HMS), harmony memory consideration rate (HMCR), pitch adjusting rate (PAR). Note that, HMS represents the solution candidates in each iteration. HMCR and PAR are equivalent to crossover rate and mutation rate in genetic algorithms.
2. Randomly create initial solution candidates in the allowed ranges of the variables (active and reactive power limits). A sample HM structure is shown as given in Equation (2). Note that the HMS is \(m\), and the number of variables to be optimized is \(n\) for this specific example. \[HM=\begin{bmatrix}P_{\text{C. Device}_{1}}^{1}&Q_{\text{C. Device}_{1}}^{1}&P_{ \text{C. Device}_{2}}^{1}&Q_{\text{C. Device}_{2}}^{1}&\cdots&P_{\text{C. Device}_{n}}^{1}&P_{\text{C. Device}_{n}}^{1}\\ P_{\text{C. Device}_{1}}^{2}&Q_{\text{C. Device}_{1}}^{2}&P_{\text{C. Device}_{2}}^{2}&Q_{\text{C. Device}_{2}}^{2}&\cdots&P_{\text{C. Device}_{n}}^{2}&P_{\text{C. Device}_{n}}^{2}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ P_{\text{C. Device}_{1}}^{m}&Q_{\text{C. Device}_{1}}^{m}&P_{\text{C. Device}_{2}}^{m}&Q_{\text{C. Device}_{2}}^{m}&\cdots&P_{\text{C. Device}_{n}}^{m}&P_{\text{C. Device}_{n}}^{m}\end{bmatrix}\] (2)
3. Input each row of this matrix to ladder iterative technique based power flow module. Use the voltage magnitude outputs to calculate the objective function values one by one. Save those values.
4. Perform improvisation and create a new solution candidate by applying the following procedure for all variables. * Compare HMCR with a randomly generated number in [0,1] interval. * If this number is smaller than HMCR, construct a new member of the solution candidate vector by randomly picking an element from the candidate solution matrix. For example, for the first variable, one can randomly select the 5\({}^{th}\) solution candidate however for the next variable another solution candidate, say the 7\({}^{th}\) solution candidate may be selected. Moreover, a mutation like operation may be performed. For this aim, another random number in between 0 and 1 is created and compared to the value of PAR. If it is bigger than PAR nothing needs to be done, for the opposite case then the value of the variable is perturbed either in positive or negative direction randomly. * If this number is not smaller than HMCR, then a random number in the range of the corresponding variable is created. Note that all the variables considered for this specific problem are either the active or reactive power outputs of DERs. These limits are defined according to the condition of the test system.
5. Find the numerical value of the objective function of the new solution candidate vector and compare it to the worst one in HM. If it is bigger nothing needs to be done otherwise replace the row related to the worst one with the new solution candidate vector.  
Figure 2: Sample 4 node system.  
If a predefined maximum number of iterations are run, or other predefined stopping criteria are reached stop, otherwise continue from step 3.