Multi-Period Optimal Reactive Power Dispatch Using a Mean-Variance Mapping Optimization Algorithm

: Optimal reactive power dispatch plays a key role in the safe operation of electric power systems. It consists of the optimal management of the reactive power sources within the system, usually with the aim of reducing system power losses. This paper presents both a new model and a solution approach for the multi-period version of the optimal reactive power dispatch. The main feature of a multi-period approach lies on the incorporation of inter-temporal constraints that allow the number of switching operations in transformer taps and capacitor banks to be limited in order to preserve their lifetime and avoid maintenance cost overruns. The main contribution of the paper is the constraint handling approach which consists of a multiplication of sub-functions which act as penalization and allow simultaneous consideration of both the feasibility and optimality of a given candidate solution. The multi-period optimal reactive power dispatch is an inherently nonconvex and nonlinear problem; therefore, it was solved using the metaheuristic mean-variance mapping optimization algorithm. The IEEE 30-bus and IEEE 57-bus test systems were used to validate the model and solution approach. The results allow concluding that the proposed model guarantees an adequate reactive power management that meets the objective of minimizing power losses and keeping the transformer taps and capacitor bank movements within limits that allow guaranteeing their useful life.


Introduction
The optimal reactive power dispatch (ORPD) consists of the management of reactive generation and voltage control resources of a power system usually with the aim of minimizing power loses. The ORPD is carried out after the scheduling of the active power dispatch and must be performed without altering it; with the exception of the reference bus that must compensate power losses [1]. Although the main objective of the ORPD is power loss reduction, other objectives-such as improving voltage profile, minimizing voltage deviation from a given value, and improving voltage stability-may be considered [2].
From the standpoint of mathematical optimization, the ORPD is a mixed integer nonlinear programming problem which combines discrete control variables such as tap setting transformers and reactive compensators, and continuous variables such as generator bus voltage set points [3]. The first attempts to solve the ORPD resorted to classical optimization techniques, such as linear programming [4], nonlinear programming [5], quadratic programming [6,7], interior point method [8], Newton-based method [9,10], gradient-based algorithm [11], and decomposition approach [12]. Despite their optimal results in other fields, when solving the ORPD, these techniques are easily trapped in local optimal solutions due to the nonlinear and nonconvex nature of the ORPD problem. Furthermore, classic optimization techniques are unable to tackle problems with discontinuous and The multi-period version of the ORPD (MP-ORPD) encompasses several time intervals and is usually planed for 24 h. The MP-ORPD is a more challenging problem than its singleperiod counterpart since it implies constraints regarding the total number of maneuvers in transformer taps, capacitor, and reactor banks throughout the time interval considered. The MP-OPRP is carried out after the execution of a unit commitment (UC) process [39,40]. The UC can be described as a multi-period optimal power dispatch where active power is adjudicated to generators as function of their price offers (bids) and limits. The UC, executed prior to the MP-ORPD, may take into account the presence of renewable generation [41,42]. In this case, the expected active power injections are inputs to the MP-ORPD, which is envisaged not to alter the previous active dispatch found by the UC-with the exception the slack bus that must account for power losses.
Being a more challenging optimization problem, there are few studies that approach the MP-ORPD. Some of the methodologies used to solve this problem include interior point methods [43][44][45], genetic algorithms [46], artificial electric field algorithm [47], dynamic programming [48,49], mathematical programming approaches [50][51][52][53], and decomposition techniques [54]. The main modeling aspects considered in the aforementioned studies are summarized in Table 1. Note that the proposed approach complements other works presented in the specialized literature based on the fact that that the stopping criterion is not only a given number of evaluations but also a goal accomplished. Furthermore, daily and inter-hour constraints are considered for both transformer and reactive compensation banks. As main contributions, this paper proposes a new model and solution approach for the MP-ORPD. One of the special features of the proposed approach is the constraint handling approach that complements the one proposed in [55] initially proposed for the single-period ORPD. In this case, the objective function along with the constraints consist on a multiplications of sub-functions. A goal on power loss reduction is set at the beginning of the algorithm. If a given candidate solution complies with the goal of power losses and with all constraints, then the result of the multiplication is one (meaning that the problem is both feasible and optimal); otherwise the multiplication is less than one. This feature facilitates the optimization process and allows us to easily verify the quality of a given candidate solution. Furthermore, the number of switching operations in transformer taps and capacitor banks is limited according to data recommended by the manufacturer in order to preserve the lifetime of these devices and avoid maintenance costs overruns. A mean-variance mapping optimization (MVMO) is used to solve the proposed MP-ORPD model. It is worth mentioning that the effectiveness and applicability of the MVMO algorithm to solve the single-period ORPD was proven by the authors in [56] where it was compared with other metaheuristic approaches that include artificial bee colony (ABC), firefly algorithm (FA), differential evolution (DE), comprehensive learning particle swarm optimization (CLPSO), biogeography-based optimization (BBO), and improved gravitational search algorithm (IGSA). The rest of the document is organized as follows: Section 2 presents the mathematical modeling of the MP-ORPD; Section 3 describes the solution approach based on a MVMO technique; Section 4 presents the tests and results in two benchmark power systems; finally, Section 5 presents the conclusions of the work.

Mathematical Modeling
The mathematical formulation of the MP-ORPD presented in this paper is an extension of the formulation presented in [55]. In this case, new constraints are added to account for maximum number of daily maneuvers in transformers and shunt elements; furthermore, the objective function is also modified to consider the enforcement of inter-temporal constraints.

Objective Function
The objective function given by Equation (1) is a multiplication of several sub-functions over a time horizon of 24 h. One of the terms is a daily loss reduction target indicated as f Lossd . The other terms allow an alternative handling of constraints related to voltage limits at system nodes, as well as apparent power flows and number of maneuvers on transformers taps. In this case, f V N (i, t) and f CL (j, t) are sub-functions that indicate limits on voltage magnitudes at node i, and apparent power flows in line j for time period t, respectively. f CTh (i, t) is the sub-function that indicates limits on number of maneuvers in transformer k at time period t.

Enforcement of Voltage Limits
The mathematical expression of the sub-function that controls voltage limits is given by Equation (2). In this case V min i and V max i are the minimum and maximum voltage limits at bus i, respectively; λv is is a penalty constant for voltage deviations, and V i,t is the voltage of bus i at time interval t.

Enforcement of Power Flow Limits
The sub-function that enforces apparent power flow limits is given by Equation (3). In this case, λb is a penalty factor for power flows out of limits, Load Rj max is the maximum apparent power flow allowed by transmission line j, and Load Rj,t is the apparent power flow of line j at time interval t.

Active Power Losses Target
The goal of active power losses, defined by the system operator, is given by Equation (4). In this case, λl is a penalty factor for power losses greater than the target defined by the system operator, Loss R represents the daily total power losses of the system, and P lossd is the target of daily power losses.
f lossd = e λl(Loss R −P lossd ), P lossd > Loss R 1, P lossd ≤ Loss R (4) Figure 1 illustrates the sub-functions used to evaluate limits on voltages, apparent power flows, and the goal on daily power losses. It is important to note that the subfunctions described below for transformer taps switching constraints and for daily capacitor switching limits feature a similar shape as those illustrated in Figure 1.

Limits on Transformer Taps Maneuvers
Equations (5)-(8) enforce the limits on number of transformer maneuvers. In this case, λ CT , C Max h,T,k and C OT (k, t) are the exponential penalty factor for deviations from the maximum allowed transformer maneuvers, the maximum cost of hourly maneuvers allowed in the k-th transformer and the hourly operating cost for maneuvers performed in the k-th transformer in period t, respectively. The last two expressions are given by Equations (6) and (7), respectively; where CT k is the cost of operating the k-th transformer, T max h,k is the maximum number of operations allowed for the k-th transformer, NT k,t is the number of operations performed by the k-th transformer in the period t, and P base is the active base power. Equation (8) indicates the number tap positions moved within two consecutive time periods of a given transformer.

Daily Limit of Transformer Tap Operations
A maximum number of daily operations is defined for each transformer in order to preserve their lifetime and avoid maintenance cost overruns. Equations (9) and (10) enforce the maximum daily switching operations of transformers. In this case, λ CT is an exponential penalty factor for number of tap operations out of limits, C Max d,T,k is the maximum daily switching cost allowed for the k-th transformer, and C OT (k, t) is the switching cost of transformer k in time interval t. In Equation (10), CT k is the cost of operating the k-th transformer, and Tmax d,k is the maximum number of daily operations allowed for the k-th transformer.

Capacitor Banks Daily Switching Limit
The number of maneuvers in capacitor banks must also be limited to avoid affecting its useful life and reduce maintenance costs. The enforcement of switching limits on capacitor banks are given by Equations (11)- (14). In this case, λ CC is an exponential penalty factor for number of capacitor maneuvers out of limits, C Max d,C,l is the maximum number of daily maneuvers in a capacitor bank, C OC (l, t) is the cost of hourly maneuvers performed by the l-th capacitor bank. In Equation (12), CC l is is the cost of operating the l-th capacitor bank and Cmax d,l is the maximum daily number of switching operations allowed for the l-th capacitor bank. In Equation (13), NC l,t is the number of hourly switching operations performed by the l-th capacitor bank in the period t, and in Equation (14), C l,t , C l,t−1 , and ∆C are reactive power capacity injected by the l-th capacitor bank in time t and time t − 1, and the reactive power injection step of the l-th capacitor bank, respectively.

Solution Approach
The MP-ORPD is a nonlinear and nonconvex optimization problem. These type of problems are better handled by metaheuristic techniques than by classical mathematical optimization. In this case, the MVMO was selected since it has proven to be effective for solving the single-period ORPD problem as reported in [57,58]. Furthermore, the MVMO algorithm features faster convergence speed than others metaheuristic techniques [57].
The MVMO approach shares several similarities with other metaheuristic algorithms; however, its main feature is a special mapping function applied to transform the variable domain for mutating the offspring, which is based on the mean and variance of the n-best population. The swarm version of this algorithm includes a population of candidate solutions [58,59].
The basic MVMO operates over a set of candidate solutions, its objective is to execute a correct and accurate optimization with a minimum quantity of objective function evaluations. The inner search space of all variables within MVMO is restricted to [0, 1]. Thus, the minimum and maximum limits of the variables must be normalized in this range. During each iteration, it is not possible that any vector solution component violates the corresponding limits. To achieve this objective, a special mapping function (h-function) was developed. The inputs of this function are the means and variance of the best solutions that the MVMO has found so far. Note that the output of this mapping function will always be within the range [0, 1], avoiding violation of variable limits during the search process. The shape and mapping curves are adjusted according to the progress in the search process and the MVMO updates the candidate solution around the best solution in each step of the iteration [60]. MVMO is able to search around the best local solution with a small possibility of being trapped in a local optimal solution. The basic considerations of the MVMO metaheuristic are detailed below.

Fitness Evaluation and Constraint Handling
The chi-squared test is applied for each candidate solution. The feasibility of the solution is checked and a fitness valued is assigned. The static penalization approach is used to handle constraints. Since the control variables are self-limited, all depending variables are restricted by applying the fitness function, as indicated in Equation (15), where f is the original objective function, n is the number of constrains, β is the order of the penalization term, v i is the penalization coefficient of the ith constraint, and g i represents the inequality constraints.

Enhanced Mapping
The mapping function transforms a variable x * i varied randomly with unity distribution to another variable x i , which is concentrated around the mean value. The new value of the ith variable is determined by Equation (16). Where h x , h 1 , and h 0 are the outputs of the mapping function based on the different inputs given by Equation (17). The mapping function (h-function) is parameterized as indicated in Equation (18), where s 1 and s 2 are shape factors that allow asymmetric variations of the mapping function. Finally, the shape factor is calculated as indicated in Equation (19). In this case, f s is a scale factor that allows controlling the search process during each iteration, whilex i and v i are the mean and variance of the solution archive, respectively. A more detailed description on how to adjust shape and scale is presented in [61].

Solution Archive
This archive constitutes the base knowledge of the algorithm to guide the search direction. Therefore, the n best individuals that the MVMO has found so far are saved in this archive. The fitness value for each individual is also saved. The following rules are adjusted to compare the generated individuals in each iteration and the existing archived solutions aiming to avoid losing high-quality solutions: (i) any feasible solution is preferred over any unfeasible solution; (ii) between two feasible solutions, the one with the best fitness value is preferred; (iii) between two unfeasible solutions, the one with the lowest fitness value is preferred. The update is only done if the new individual is better than those that are currently in the archive. Feasible solutions are located in the upper part of the archive. These solutions are organized according to their fitness value. Unfeasible solutions are organized according to their fitness and located in the lower part of the archive. Once the archive is complete with n feasible solutions, any unfeasible solution candidate would not have the chance to be saved in the archive.

Offspring Generation and Stopping Criterion
The first positioned solution in the archive (the best so far), named as xbest, is assigned as the parent. For offspring generation a variable selection is performed. MVMO searches the mean value saved in the solution archive for the best solution only in m selected directions. This means that only these dimensions of the offspring will be updated while the D − m remaining dimensions take the corresponding values of xbest, being D the dimension of the problem (number of control variables). Then, mutation is carried out for each m-selected dimension. Finally, the search process of the MVMO stops after a predetermined number of fitness evaluations is carried out.

MVMO Swarm Variant
In its swarm variant, MVMO is initiated with n p particles. In this case, each particle or candidate solution has its own solution archive and mapping function. In the process, each candidate solution executes m steps to identify an optimal set of independent solutions. Later, the particles exchange information. In some cases, some particles are very close to each other, which means that there is information redundancy. This redundancy is eliminated by discarding redundant particles. As with particle swarm optimization, a local and global best (gbest) solution are defined. Furthermore, the normalized distance between each particle to the local best and global best solutions are calculated. The i-th particle is discarded from the process if its normalized distance is less than a certain predefined threshold. If such threshold is zero, all particles are taken into account in the whole process; otherwise, if the threshold is one all particles are discarded except from the global best. In this case, after (m * n p + n p ) fitness evaluations only one particle,the gbest, remains.
Intermediate threshold values entail better adaptation to any optimization problem [58].
After independent evaluation, and if the particle is further considered, its searching will be directed toward the global solution by assigning global best instead of local best solution, as parent for the particle's offspring. The remaining steps are identical with those of the classical MVMO. Figure 2 illustrates the flowchart of the implemented methodology. A detailed description of swarm MVMO variant is presented in [58,62].

Tests and Results
To show the effectiveness and applicability of the proposed approach several tests were carried out with two benchmark power systems: namely, the IEEE 30 and 57 bus test systems. This section details the application of the proposed model and solution approach.

Description of the Test Systems
Several simulations were performed on the IEEE 30-bus [63] and 57-bus test systems [64]. The IEEE 30-bus test system consists of 41 branches, 24 PQ buses, and 19 control variables; these include 6 generation units, 4 transformers, and 9 capacitors. The IEEE 57-bus test system consists of 80 branches, 50 PQ buses, and 25 control variables; these are 7 generation units, 15 transformers, and 3 capacitors. For both systems, a time horizon of 24 h is defined, for which active and reactive power demands are established at each node. The per unit demand curve shown in Figure 3 illustrate the load curves of both systems. It should be noted that three types of users, namely residential, commercial, and industrial, were considered.  Table 2 presents the power losses for the base case and the goal aimed at for each test system. It is worth mentioning that this goal must be defined by the power system operator based on its knowledge of the system. In this case, the goal of power losses for each system was obtained through a sensitivity analysis. One of the features that distinguish the proposed approach from other methods reported in the specialized literature is the handling of constraints. The sub-functions used to enforce constraints must be parameterized with several lambdas. In this case, the penalty factors (λ v , λ b , λ ct , λ cc , λ cr ) for each of the sub-functions that make up the total fitness function are set to 0.05, except for λ l associated with the sub-function that enforces the maximum power flow through transmission lines which is set to 0.005, making this consideration less restrictive. These values of lambdas were established through sensitivity tests. Likewise, the cost per switching of transformer tap changers ( CT k ) is defined as USD 6 per maneuver, and the cost per switching for capacitor banks (CC l ) is defined as USD 4 per maneuver. The maximum number of evaluations was set to 40,000, and 40 particles were used as initial candidate solutions.
The maximum number of daily operation allowed for the k-th transformer and Tmax d,k is established according to technical data. The daily limit is set to 114 maneuvers, taking into account 500,000 operations in 12 years [65]. As regards capacitor banks the maximum number of maneuvers recommended by the manufacturer is 100,000 in 20 years (to avoid premature maintenance interventions), which yields a maximum of 13 daily maneuvers [66].

Results with the IEEE 30-Bus Test System
The IEEE 30-bus test system used in this paper is illustrated in Figure 4, to which 9 capacitor banks (each of 4 MVAr) proposed in [67] were incorporated. This results in 19 control variables which correspond to 6 controllable voltages at generation buses 1, 2, 22, 27, 23, and 13; 4 transformers with tap changers at branches 6-9, 6-10, 4-12, and 4-12; and 9 capacitor banks located at nodes 10, 12, 15, 17, 20, 21, 23, 24, and 29. In this case, the maximum demand is 181.73 MW in hour 14. Furthermore, voltage limits were set within 0.95-1.1 pu. After running the MVMO algorithm, several feasible solutions were found for the MP-ORPD, which comply with the objective function and constraints. All tests were carried out on a laptop AMD Ryzen 7 with 20 G of RAM memory. Table 3 indicates the statistical results for 25 independent runs of the algorithm. The second column indicates the power loses achieved by the algorithm. It should be noted that despite the variety of solutions, most of them were significantly close to the goal of 59.9 MW of power losses set for the IEEE 30-bus test system. The third column indicates the processing time of the algorithm.  Figure 5 illustrates the convergence of the proposed algorithm for three independent runs. As the problem is set as a minimization the optimal solution is −1. Note that although 40,000 evaluations are set as the limit, the convergence is achieved at approximately 22,000 evaluations.  Figure 6 illustrates the voltage profile at generation bus 2 for three randomly selected runs of the algorithm. In this case, most voltage magnitudes are improved with respect to the base case. It is worth mentioning that the algorithm does not look for voltage improvement as such, but instead for a multi-period operation that guarantees power system limits and minimizes active power losses. In this case, all voltage magnitudes are within the permissible operating range.

Results with the IEEE 57-Bus Test System
The IEEE 57-bus test system, illustrated in Figure 10, features 80 branches and 25 control variables: 7 generators with voltage control located at buses 1, 2, 3, 6, 8, 9, and 12; 15 on load tap changing transformers at branches 4-18 (2 transformers), 21-20, 24-26, 7-29, 34-32, 11-41, 15-45, 14-46, 10-51, 13-49, 11-43, 40-56, 39-57, and 9-55; and 3 capacitor banks located at nodes 18, 25, and 53. Note that this system has more transformers than capacitor banks. For this reason, a single-step reactive power injection of 10 MVAr is considered for each capacitor bank. After running the MVMO algorithm, several feasible solutions were found for the MP-ORPD, which comply with the objective function and constraints. All tests were carried out on a laptop AMD Ryzen 7 with 20 G of RAM memory. Table 4 indicates the statistical results for 25 independent runs of the algorithm. The second column presents the power losses achieved by the algorithm. Note that despite the variety of solutions, most of them were significantly close to the goal of 219.5 MW of power losses set for the IEEE 57-bus test system. The third column indicates the computational time of the algorithm.  Figure 11 illustrates the convergence of the proposed algorithm for three independent runs. The maximum number of evaluations was set to 40,000; nonetheless, it can be observed that the algorithm achieves convergence at approximately 25,000 evaluations. Figure 11. Convergence of the MVMO algorithm for the IEEE 57-bus test system. Figure 12 illustrates the voltage setpoints of generator located at bus 3 for three independent runs of the algorithm. It should be noted that all voltage setpoints are higher than the ones given in the base case. Figure 13 illustrates the voltage at load bus 9 for three runs of the algorithm. Note that the MP-ORPD allows having higher voltages in load buses; nonetheless none of them is out of the operative limits as indicted by the constraints. Figure 14 illustrates the taps movement for the transformer located at branch 7-29. It should be noted that tap variations correspond to the limits previously established.    Figure 15 illustrates the reactive compensation provided by the capacitor bank located at bus 18. Note that this compensation is used almost all the periods in run 2; nonetheless, it is not used in certain periods of time for runs 1 and 3. It is worth mentioning that the compensation does not need to be fully active all the time, since it is used when needed, to comply with the constraints and objective function.

Conclusions
This paper presented a new mathematical model for the multi-period version of the ORPD problem. Besides traditional limits on voltage magnitudes and power flows, the proposed MP-ORPD model incorporates inter-temporal constraints related to the number of maneuvers of transformer taps and capacitor banks. Such constraints allow to maintain the useful life of these elements. A key aspect of the proposed modeling is the design of the objective function. This consists of a multiplication of sub-functions in which penalization due to constraint violations is included. The main advantage of this objective function is the fact that both feasibility and optimality are easily verified.
Several tests carried out on two benchmark power systems evidenced the effectiveness and applicability of the proposed approach. For both power systems is was found that the algorithm was capable of finding high-quality solutions within reasonable computational time.
The MVMO algorithm proved to be an effective tool for solving the proposed MP-ORPD model. This technique presented adequate convergence as evidenced in the results. It is worth mentioning that for all the solutions found by the MVMO, the operative constraints and goal on power losses were achieved. Furthermore, there are notable improvements in the voltage profiles and in the reduction in transformer taps and capacitor bank maneuvers.
Although the MVMO proved to be an effective tool for solving the MP-ORPD, its main limitation lies in the high computational time required to find a solution. Nonetheless, such drawback is shared by all metaheuristic techniques. It is worth mentioning that the computational time is not a major inconvenience, since the MP-ORPD as proposed in this paper is not envisaged to be carried out online. A solution to the this issue and a future work is to explore parallel computation alternatives to reduce the computational time of the MP-ORPD.