A Wolf Pack Algorithm for Active and Reactive Power Coordinated Optimization in Active Distribution Network

This paper presents an active and reactive power dynamic optimization model for active distribution network (ADN), whose control variables include the output of distributed generations (DGs), charge or discharge power of energy storage system (ESS) and reactive power from capacitor banks. To solve the high-dimension nonlinear optimization model, a new heuristic swarm intelligent method, namely wolf pack algorithm (WPA) with better global convergence and computational robustness, is adapted so that the network loss minimization can be achieved. In this paper, the IEEE33-bus system is used to show the effectiveness of WPA technique compared with other techniques. Numerical tests on the modified IEEE 33-bus system show that WPA for active and reactive multi-period optimization of ADN is exact and effective.


Introduction
With increasing penetration of decentralized generation, distributed energy storage devices (ESD), control equipments and advanced communication networks, the distribution system is gradually becoming to the active one from traditionally passive network [1]. But high penetration of DGs access to the power grid may lead to voltage fluctuation or overvoltage, which seriously restrict the consumptive capacity of ADN for renewable energy generations (RESs). Therefore, it is necessary that the active and reactive power of dispatchable DGs need to be optimally dispatched by distribution network operators (DNOs). So the active and reactive power coordinated optimization (ARPCO) has become a research topic.
In the past years, many researchers have also engaged in optimal programming and operation. Numbers of programming methods or intelligent techniques were applied to optimize reactive power dispatch (OPRD). For example, particle swarm optimization (PSO) [2,3], artificial bee colony (ABC) [4] and the hybrid differential evolution with ant system applied by C.-M. Huang [5] to solve OPRD. On the compelling problem of ADNs control, Reference [6] has proposed a new bidirectional converter as the interface of combined generating and electric storage systems with the grid. In reference [7], a broadcast-based unified control algorithm was brought forward to provide reactive power support for the grid by a seamless control of heterogeneous energy resources such as distributed storage systems and demand-responsive loads. In [8], a real-time distributed generation control scheme and a real-time decentralized load control scheme have been designed for power networks with tap-changing transformers.
Distribution networks are characterized by reduced line lengths with a non-negligible resistance over reactance (R/X) ratio, so the ARPCO of ADNs is more significant. In [9], a two-stage optimal schedule model for ADN has been proposed to dispatch active and reactive power of DGs. Marco Bronzini has presented a centralized control model for distributed resources in ADN and solved the problem using quasi-Newton method [10]. Levron et al. [11] presented a dynamic search algorithm based on a load flow for the OPF problem in microgrids with ESDs and RESs, but the method worked well only for systems with a small scale system. A solution for ARPCO in ADNs was explicated by Reference [12], but the charging and discharging periods of ESS were fixed so that the solutions were low-quality. In [13], a mixed-integer second-order cone programming (MISOCP) model has been proposed for the optimal operation problem in ADN. This method can convergence to the optimal value by using existing optimization software, but the solution is not accurate enough. Graditi G. [14] models the 24h behavior of appliances by means of a real-valued function and then extended the GSO algorithm to solve multi-objective energy management optimization problems for smart grids with direct control over shiftable loads. This method was verified comparable or even better than NSGA-II [15]. In [16], a comprehensive multi-objective optimal technique was proposed for optimal operation of ADN. The method was included the internal procedure and the external one. The former applied the Nondominated Sorting Genetic Algorithm II for the optimal management, and the latter was responsible for choosing the design features such as ratings and/or types. This paper proposed the ARP coordinated dynamic multi-period optimization model, and applies a new swarm intelligent optimization algorithm-WPA to solve ARPCO problem. WPA has been proposed by [17] and is especially suitable for solving high-dimension optimization problems. This paper is organized as follows: Section 2 discusses the ARPCO model, and section 3 briefly describes WPA technique. Then, WPA implementation in solving ARPCO model is explicated in section 4. Section 5 presents the simulation results and discussion. At last, section 6 states the conclusion of this paper. The main contribution of this paper includes the realistic ARPCO model considering ESDs and application of the WPA to solve it, whose main advantages are: 1) a realistic and precise model; 2) computational robustness; and 3) good convergence performance.

ARPCO Model
The equipments such as dispatchable distributed energy resources (DER) and switchable shunt capacitor banks (SCB) are usually controlled for optimal operation of ADN. So the active power of DER and reactive power of SCB are chosen as control variables, ie, u= [P DG In general, the minimum energy loss in the scheduling period is selected as the objective function for the optimal operation of ADN, f expressed as follow: Where T is the period number of a complete dispatching cycle, N l is the number of transmission lines, V t and θ t , as state variables, are respectively the bus voltage amplitude and phase angel in the periodt. G i(m,n) is the conductance of the branch-i connecting node-m and node-n.
The equality constraint equations include load flow suggested in [8] and the operation equality constraints of distributed generator (DG) and ESS, as follows: 2) Where P G may be the active power from DG, ESS or substation. For the ESS, its charge power is positive and discharge power is negative.
Distributed generator (DG): Where P DG i , Q DG i , φ are the active power , reactive power and power factor angle of DG-i respectively. ESS: t is the remaining energy of the ESS-i in the period-t, N s is the number of ESS. In addition, the inequality constraints include operation constraints of ADN, as follow: DG constraints: P DGi min ≤ P Dgi t ≤P DGi max (5) Where P DGi min and P DGi max are the upper and lower limits of active power generation respectively. ESS constraints: Charge or discharge power and energy are restricted by their upper and lower limits, as follow: For longer life of the ESS, the energy stored in ESS is generally from 20% to 90% of E j max in fact [13]. Constraints of power from substation: In order to restrain the influence of the AND power fluctuation on the transmission network, the constraint of power from substation is considered as follow: 2) Shunt capacitors restricted by their limits: t is the node-i voltage amplitude in the period-t; I ij t is the current of the buanch-ij in the period-t; V i min and V i max are the voltage limits, I ij max is the max branch current.

WPA
Based on the cooperative hunting characteristics of wolves, Wu Husheng, etc put forward a new intelligent algorithm-Wolf pack algorithm (WPA) [17]. It has been proven to has better global convergence and computational robustness, and especially suitable for solving high-dimension and multimodal function optimization problems with other classical intelligent algorithm such as PSO, fish swarm algorithm, genetic algorithm (GA) and so on. In nature, the wolf belongs to candidate family and lives in a pack consisting of 5-12 wolves on average. In general, the common wolf pack can be divided into three parts: a lead wolf, some scout wolves and ferocious wolves. The lead wolf, as a leader, is always the smartest and most ferocious one. Scout wolves hunt around for prey. Ferocious ones are responsible for rounding up the prey.
WPA is different from previous swarm intelligence optimization algorithm. It has three different search abilities: hunting behavior, calling behavior, siege behavior, and two intelligent rules: a winnertakes-all productive rule for lead wolf and a randomly renewable mechanism named survival of the stronger for the wolf pack.
The parameters such as X i , N, k max , T max , S , L near , α, β and the objective function value Y have been defined explicitly and all the components have been discussed in detail in [8], so the main computation steps are only described below.
Step1: Initialization. For a D dimension of the search space, the follow parameters are initialized: Step2: The wolf with the maxium Y value (ie, Y lead ) is selected as lead wolf, and the rest N sw (an inter value between N/(α+1) and N/α) wolves with better Y value, as scount wolves, begin to scout in D dimension solution space according to equation (10)  Step3: The rest N fw (=N-N sw -1) wolves take calling behavior according to equation (11). If Y i ≥Y lead , go to step2; otherwise the wolf-i continue running until L (i, l) ≤ L near , then go to step 4; Where x id k+1 is the position of ferocious-wolf-i at the (k+1)th iteration; step b is the summoning step length, g d k is the position of lead wolf in the dth variable space at the kth iteration. L(i, l) is the distance between the ith wolf and the lead wolf and expressed as Manhattan distance.
Step4: The position of siege wolves is updated according to equation (12); Where λ is a random number uniformly distributed at the interval [-1, 1], step c is the siege step length. The relationship between the step length step a , step b and step c in dth variable space is as follows: (13) Step5: Update the lead wolf to the one with maximum objective function value, and then delete R wolves with worst objective function value from the wolf pack and meanwhile produce R wolves according to equation (14) to replace them; Where R is an integer and randomly selected at the interval [n/(2β),n/β] .
Step6: If the program reaches the precision requirement or the maximum number of iterations, the position and Y lead of lead wolf, the problem optimal solution, will be outputted, otherwise go to step2 to continue iteration until termination condition is met.
In the above steps, step 2 is fine search, reflecting the local optimal solution search precision; step 3 is a rough search, reflecting the local optimal solution search efficiency; step 4 is gradually refined search, which reflects the accuracy that the local optimal solution is also global one; step 5 produces a new generation of wolves, which not only retains the excellence of the local optimal solution founded by original wolves, but also increases the probability of reaching the optimal solution to guarantee the global optimality of the algorithm.

WPA for solving ARPCO problem
In order to apply the WPA to find the optimal control variables of the ARPCO problem , the vector of polulation can be expressed as follows: The objective function Y is evaluated according to equation (16) as follow: Y = -f = -P LOS (16) The initial position of each artificial wolf can be formed as follow: Where ω is distance determinant factor. For the ARPCO problem, the iteration termination condistion is expresse as follow: Where is the convergence control constant number. In order to calculate the objective function, position variables from each wolf are firstly mapped into the power flow equations and the flow program is executed to obtain the network loss by MATPOWER. Y is evaluated according to equation (16). Then the simulation is processed according to the WPA solving procedure from step1 to setp6 described above. During the process, these variables out of bound are tagged at their boundaries. The simulation is repeated until meeting the termination condition. The implementation of WPA for ARPCO problem is depicted in Figure 1. Or T > T max ? Figure 1. Flowchart of propsed WPA for solving ARPCO

Simulations
In order to apply the simulations for solving ARPCO problem by WPA are implemented using Matlab R2012b on a Windows 7 Professional Intel i5-3210M CPU 2.5GHz 8GB RAM. In this paper, the IEEE 33 node radial distribution system which PV, SVC, ESS, CB are integrated into is selected as the example system, as shown in figure 2. The network with 37 branches is radially operated. The total load power: P load = 3635kW, Q load = 2265kvar; for each PV: P PV max = 300kW, cosφ = 0.95; for each shunt capacitor (SC): Q SC min = 25kvar, Q SC max = 100kvar; for each ESS: P ES max = 240 kW, P ES min = -200 kW, E max = 1200kWh.

Figure 2. Modified IEEE 33-bus test system
The load curve and the sunshine intensity curve of the system at the moment of 24 are shown as figure 3. For simplifying the analysis, all of the PV use the same sunshine curve and load node use the same load curve.

Figure 3. Load and sunshine profile
It is seen from figure 3 that the peak value of the load and the peak point of the light does not appear at the same time. If there are no control equipments such as energy storage, PV will abandon light due to the light load so that the load absorb power only from the substation during the heavy load period , and then the heavy load node voltage will drop and the voltage difference further rise, which causes a large number of network loss.
To test the robustness of this method under the large volume data condition, MATLAB cubic spline interpolation method is used to interpolate the above 24 hour data as 96 data points so that the interval acquiring the data for the day-ahead scheduling is 15 minutes.
WPA parameters were specified as according to reference [17]: N = 50, K max = 1000, S = 200, T max = 15, α = 4, β = 5, ε = 10 -3 . After 68 iteration steps (shown as fig. 4), the global optimal solution were found and the mean square error is about 0.0005, the total time is 20. 517s. The ARPCO results are shown as figure 5(a) ~ 5(e). PV is almost in max power state because of ESS. Figure 5. (a) shows that ESS start to charge at load trough period and to discharge at load peak period so as to reduce the total loss in the scheduling period. Again, energy is stored when output of PV increases and injected into the system when load power increases, shown as figure 5. (b). Network loss is smaller because of smaller voltage difference at load trough period and is larger because of larger voltage difference at load peak period. As a result, the total network loss reduces compared with no optimization, shown as figure 5. (e) with less total power loss than the loss in [13]. Figure 5. (c) shows that SCB operation occurs when the load demand and PV generation increase, and their operation numbers are strictly restricted in the allowable range. Figure. 5(d) shows that the distribution system absorbs less active and reactive power from the substation after ARPCO process in order to improve the penetration of PV. From table 1, it can be seen that the energy loss obtained by WPA is the lowest and the WPA can converge to the optimal solution of the problem in a reasonable amount of time. Therefore the WPA for solving ARPCO problem is a competitive algorithm.

Conclusions
This paper proposed the ARPCO model considering the output of DG, charge or discharge power of ESS and reactive power from SBs as control variables. And minimizing the energy loss of the distribution network was taken as the optimal objective. The use of WPA technique guarantees convergence to optimality in a reasonable amount of time.
The modified IEEE 33-bus test system was used to demonstrate the accuracy of the ARPCO mathematical models and the efficiency of the WPA for solving it. When compared with PSO and SOCP technique used in [13], WPA shows more effective because of minimum total energy loss and reasonable convergence speed.
In future, the stochastic of renewable resources can be considered in solving ARPCO problem to demonstrate the robustness of WPA. In addition, an advanced WPA with better performance will be proposed in future.