Dynamic Optimal Power Flow in FSWGs Integrated Power System

INTRODUCTION

Wind energy is the world’s fastest growing renewable energy source. With the increasing levels of wind generator penetration in modern power systems, one of major challenges in the present and coming years is the optimization control, such as optimal power flow including wind farms (Li and Chen, 2008).

The well-known advantages of FSWGs (fixed speed wind generators) are it is robust, easy and relatively cheap for mass production, there are still widely use in power system with FSWGs. Since FSWGs always draw reactive power from the grid, in most cases, capacitors are connected in parallel to the generator to compensate for the reactive power consumption, the P-bus power flow model of FSWGs bus is proposed in this study.

In this study, the problems of DOPF (Dynamic Optimal Power Flow) including FSWGs are researched. The expectation model of wind generators' active power outputs is adopted. DOPF is a method to schedule the online generator outputs with the predicted load demands over a certain period of time so as to operate an electric power system most economically. Normally, it is solved by dividing the entire dispatch period into a number of small time intervals, then a static economic dispatch has been employed to solve the problem in each interval (Pothiya et al., 2008; Uturbey and Costa, 2007). In this study, DOPF model, which takes all conventional units cost minimum as the objective function and takes the whole time and the inherent relations of different intervals into account in FSWGs integrated system, is established.

DOPF is a typical complex multi-constrained non- convex non-linear programming problem when considering the valve-point effect of conventional generators. Moreover, it is mixed integer when considering the discreteness of FSWGs reactive compensation devices. However, both lambda-iterative and gradient technique methods in conventional approaches to the problems are calculus-based techniques and require a smooth and convex cost function and strict continuity of the search space.

On intelligent algorithms, easy to handle discrete variables, no special requirements for the objective function, no problems of dimension disaster, strong global search capability and the greater probability in finding out the global optimum solution are the prominent advantages of these algorithms. However, the computational accuracy and convergence rate are not high for most intelligent algorithm.

In current study, we see a myriad of improved algorithms. Many of these algorithms incorporate some improvement strategies in the algorithms itself while some improved algorithms are the result of hybridizing two algorithms (He et al., 2008; Yuan et al., 2009). The Taguchi method is an established approach for robust design, applying the idea from statistical experiment design for evaluating and improvements in products, processes and equipment (Taguchi et al., 2000, 2005; Ross, 1989). The key to this concept is to improve the quality of a product by minimizing the effect of the causes of variation without the elimination of the relevant causes. The two major tools used in the Taguchi method are: (1) Signal-to-Noise Ratio (SNR) which measures quality and (2) orthogonal arrays which are used to study many design parameters simultaneously (Taguchi et al., 2000, 2005; Ross, 1989). This study adopts the concept of Taguchi method into PSO algorithm to improve its convergence which is so called Taguchi-PSO (TPSO). By applying the algorithm to dealing with IEEE 39-bus system, compared with PSO algorithm, the experimental results show that the algorithm is an effective way to solving DOPF problems in FSWGs integrated power system because of higher computational accuracy, better robustness and convergence performance.

IEEE 39-bus system is used to illustrate the effectiveness of the proposed method compared with those obtained from PSO algorithm. The test results show that the proposed method is effective and has a certain practicality.

DOPF MODEL IN FSWGS INTEGRATED POWER SYSTEM

Due to the random variation of the wind velocities and load demands, it is difficult to research the DOPF in the power system including wind farms. For simplifying this problem, the dividing-stage strategy is adopted in this study. According to the wind velocity forecasting curves and the load forecasting curves in the planning horizon, the expectations of wind generators' outputs and the load demands at dispatch interval can be calculated.

Constraints: Constraints include equality and inequality constraints. The equation constraint is the power flow formulation constraint while inequality constraints including generator power output, ramp rate and bus voltage are as in Eq. 1-4. For the FSWGs bus, absorbed reactive power from the power system is bounded as in Eq. 4. The constraints of real power generation limit and the ramp rate are taken into account as in Eq. 1.

(1)

(2)

(3)

(4)

where, P_i,min and P_i,max are the maximum and minimum limits of the power generation of unit i, P_i^t is the real power output of unit i at the tth interval, P_i^t-1 is the real power output of unit i at the t-1th interval; U_Ri is the up-ramp limit of the ith generator (in units of MW/time-period) and D_Ri is the down-ramp limit of the ith generator (in units of MW/time-period); ΔT is time interval, N_gen is the number of conventional generating units and N is the number of system buses (excluding slack bus); V_i^t is the voltage magnitude output of bus i at the tth interval; Q^t_Gi is the reactive power output of conventional generating unit i at the tth interval; Q_ab^t is the absorbed reactive power of wind farm bus at the tth interval; max is the maximum value of the variable, min is the minimum value of the variable.

After calculating the power flow, the state variables, power loss and real power output of the slack bus generator corresponding to the current control variables are available. The real power output of the slack bus generator will be set to the limit if it violates the limit. After handling overlimit of the real power output of the slack bus generator, the system power balance constraints as in Eq. 5 must meet, otherwise adding Eq. 5 as penalty terms to the objective function to form a generalized objective function. Details of the generalized objective function used in this study are given in Eq. 8.

(5)

where, ΔP^t is the unbalance of the real power at the tth interval, N_gen-1 represents the number of conventional generating units excluding the slack bus, P_sl^t is the real power output of the slack bus generator after handling its overlimit at the tth interval, P_ls^t is the total power loss at the tth interval, P_ld^t is the total load expectation at the tth interval, P^t_w,av is the expectation of wind generators' real power outputs at the tth interval.

Objective function: Due to the fact that wind generation does not consume the fuel, the utility must purchase all the energy produced by wind generating units. Consequently, the objective is to minimize the following total incremental fuel cost function F associated to N_gen dispatchable units for T intervals in the given time horizon, subject to the above-mentioned equality and inequality constraints.

(6)

The inclusion of valve-point loading effects makes the modeling of the fuel cost function of the unit more practical. This increases the non-linearity and local optima in the solution space. Also, the solution procedure can easily trap in the local optima in the vicinity of optimal value. The fuel cost function of the ith unit F (P_i^t) with valve-point loadings are represented as follows:

(7)

where a_i, b_i and c_i are cost coefficients and e_i, f_i are constants from the valve-point effect of the ith generating unit.

Evaluation function: We must define the evaluation function for evaluating the fitness of each individual in the population. In the most of the nonlinear optimization problems, the constraints are considered by generalizing the objective function using penalty terms.

FSWGs always draw reactive power from the grid. According to the principle of the nearest reactive power compensation, the required reactive power of wind farm doesn’t absorb from the system as much as possible, but mainly from the wind farm reactive compensation devices. Because reactive power compensation devices are group switched capacitors in wind farms, the maximum of absorbed reactive power of wind farm from the grid is set to one set of group switched capacitors which is considered by generalizing the objective function using penalty terms. Considering efficiency maximization of wind farm reactive power compensation devices, the injected reactive power of wind farm bus is not constrained.

To sum up, the above problems are generalized as follows:

(8)

where, K_V, K_Q, K_C and K_D are variable overlimit penalty coefficients, V_i^t is the voltage magnitude of bus i at the tth interval (excluding the slack bus and PV bus ); Q^t_Gi is the reactive power output of generator i at the tth interval; Q_ab^t is the absorbed reactive power of wind farm bus at the tth interval; V_i^lim and Q_Gi^lim denote the violated upper or lower limits; Q_ab^lim denotes the violated lower limit that is defined as:

(9)

where, Q_ab^max is the upper limit value of the absorbed reactive power of wind farm bus which is set to 1Mvar in the study.

In this study, K_V, K_Q and K_C are set to 1,1,1 respectively. Because the unbalance of the real power ΔP is hard to meet, an adaptive penalty function to handle penalty coefficient K_D is adopted, where, k is the algorithm’s current iteration number; β is a relative violated value of the constraints, γ is a multi-stage assignment value, a is the power of the penalty value.

Meanwhile, several experiments have been done in order to obtain the penalty parameters. In this study, if β≤1 then α = 1, otherwise α = 2. Furthermore, if β≤0.001, then γ = 1, else, if β≤0.01 then γ = 10, else, if β≤0.1 then γ = 30, else, if β≤1 then γ = 100, otherwise γ = 300.

THE CONCEPT OF TAGUCHI METHOD

Taguchi methods were developed by Genichi Taguchi as an aid to improving the implementation of total quality control in Japanese industry. He described his methods as 'quality engineering and, in fact, the methods take an engineering approach to the understanding of process information. The methods are based on the design of experiments to provide near optimal quality characteristics for a specific objective. Taguchi defines quality in a negative manner as the loss imparted to society from the time the product is shipped (Taguchi et al., 2000).

The method is based on several statistical concepts which have proven to be valuable tool in the subject of quality improvement. Many Japanese manufacturers have applied this approach for improving product and process quality with unprecedented success. Taguchi essentially utilizes the conventional statistical tools, which has been simplified by identifying a set of stringent guidelines for experiment layout and the analysis of results. Recently Western industries have begun to recognize Taguchi’s method as simple but highly effective approach in improving product and process quality. Taguchi method has a wide range of applications in other areas (Cornejo-Mazon et al., 2008; Nava-Arenas et al., 2008). Two major tools used in the Taguchi method are the orthogonal array and the SNR. This study incorporates the two-level orthogonal arrays (Table 1), whereby the two levels are denoted by number 1 and 2 in the Table 1. The basic concept of the structure and usage of two-level orthogonal arrays are briefly described here. This is adequate in the context of describing the development of the hybrid TPSO algorithm.

Table 1:	L8 (27) orthogonal array

The general form of two-level standard orthogonal arrays can be represented by the following mathematical term:

(10)

where, n = 2^k, n is the number of experimental runs, k is a positive integer which is greater than 1, 2 is the number of levels for each factor and n-1is the number of columns in the orthogonal array. The size of the orthogonal array to be used depends on the problem instance. At a time only two particles (two factors) are executed at two different levels for n number of experiments. The two-level standard orthogonal arrays most often used in practice are L₄(2³), L₈(2⁷), L₁₆(2¹⁵), L₃₂(2³¹), L₆₄(2⁶³), L₁₂₈(2¹²⁷) and so on.

Orthogonal arrays are readily composed and are available from many texts (Ross, 1989). The way that they are constructed is to have (1) each level of every factor appear the same number of times in every column of the array and (2) each combination of factors between any two columns appears the same number of times. For any given values of (number of levels) and (number of factors), (number of tests) is determined as the smallest number that satisfies (1) and (2) above.

In this study, a two-level orthogonal array is used. There are m factors, where m is the number of design factors (variables) and each factor has two levels. If, m < n-1, only the first columns are used, while the other n-1- m columns are ignored. For example, there are 6 factors with two levels for each factor. We only need 6 columns to allocate these factors and L₈(2⁷) is enough for this purpose because it has 7 columns. Table 1 shows an orthogonal array L₈(2⁷). The number on the left of each row is called the run number or experiment number and runs from 1 to 8. In fact, each row represents a particle with only first 6 columns. In our methodology, we randomly choose two particles in every iteration. For this example, perform 8 experiments based on Table 1 whereby the number (either 1 or 2) of each column is the dimension from first and second particle, respectively.

Many designed experiments use matrices called orthogonal arrays for each experimental run and for analyzing the data. The array is called orthogonal because all columns can be evaluated independently of one another, as shown in Table 1. As one column represents a particle, therefore each column has its own fitness value (f₁, f₂, … f₈). The primary goal in conducting this matrix experiment is to determine the best or the optimal level for each factor. In taguchi based PSO, two tools of taguchi method i.e. orthogonal array and SNR are used. The basic idea is to identify the one best particle in an orthogonal array experiment, which consists of a set of experiments conducted at two different levels.

Besides two-level orthogonal arrays, there exists three-level and four-level orthogonal arrays which can provide better precision to the specified process. The detailed description of how these orthogonal arrays are formed can be found in the books written by Ross (1989), Roy (1990) and Taguchi et al. (2000).

OVERVIEW OF PARTICLE SWARM OPTIMIZATION

Kennedy and Eberhart (1995) first introduced the Particle Swarm Optimization (PSO) method. It is one of the optimization method categorized in the family of evolutionary computation. The method has been found to be robust in solving real-world problems featuring non-differentiability, high dimension, multiple optima and non-linearity which is widely used in various fields (Heng et al., 2006; Sutha and Kamaraj, 2008; Lai-Jun et al., 2009; Yang et al., 2009; Zakermoshfegh et al., 2008). PSO algorithm is a model that mimics the movement of individuals (fishes, birds, or insects) within a group (school, flock and swarm). Similar to GA, a PSO consists of a population refining its knowledge of the given search space. PSO is inspired by particles moving around in the search space. The individuals in a PSO thus have their own positions and velocities. Each particle moves in the search space with velocity which is dynamically adjusted and balanced based on its own best movement (pBest) and the best movement of the group (gBest).

Instead of using evolutionary operators such as selection, mutation and crossover, each particle in the population moves in the search space with velocity which is dynamically adjusted and all particles are assumed to be of no volume. Each particle keeps track of its coordinates in the search space, which are associated with the best solution it has achieved so far. This value is known as pBest. Another best value that is tracked by the global version of the particle swarm optimizer is the overall best value or the best solution in the population is called gBest.

The PSO concept consists of, at each time step, changing the velocity of each particle toward its pBest and gBest solutions. The movement is weighted by a random term, with separate random numbers being generated toward pBest and gBest values. For example, the ith particle consisting d dimensions is represented as X_i = (X_i,1, X_i,2, X_i,3, …, X_i,d). The same notation applied to the velocity, V_i = (V_i,1, V_i,2, V_i,3, …, V_i,d). The best previous position of the ith particle is recorded and represented as pBest_i = (pBest_i,1, pBest_i,2, pBest_i,3,… pBest_i,d). In the case of minimization that we consider in this study, the value of pBest_i with lowest fitness is known as gBest. The modification of velocity and position can be calculated using the current velocity and the distance from pBest_i to gBest as shown in the following formulas

(11)

(12)

where, i ∈ 1...N, j ∈ 1...d, t ∈ 1...T with N is the number of population size, d is the number of dimension and T is the number of maximum generation.

The position, X of each particle is updated for every dimension for all particles in each iteration. This is done by adding the velocity vector to the position vector, as described in Eq. 12 above. In Eq. 11, w is known as the inertia weight. This parameter was introduced by Shi and Eberhart (1998) to accelerate the convergence of PSO. Suitable selection of w provides a balance between global and local explorations, thus requiring less iteration on average to find sufficiently optimal solution. Low values result in particles moving in the region far from the optimum value before being tugged back. On the other hand, high values result in abrupt movement toward target regions.

The parameters ρ₁ and ρ₂ are set to constant values, which are normally given as 2.0 whereas r₁ and r₂ are two random values, uniformly distributed in [0,1]. The constants, ρ₁ and ρ₂ represent the weighting of the stochastic acceleration terms that pull each particle toward pBest and gBest positions.

TAGUCHI-PARTICLE SWARM OPTIMIZATION

Here, we illustrate with a simple example the concept of Taguchi-PSO. The Taguchi method is performed after the entire population is updated via PSO equations. Following this Taguchi method, further elaborations are given in a list of steps below. All the pseudocodes are written in programming style for L₈(2⁷),other two-level standard orthogonal arrays such as L₄(2³), L₁₆(2¹⁵), L₃₂(2³¹), L₆₄(2⁶³) and L₁₂₈(2¹²⁷) etc., can be analogized.

For instance, for the case of factor = B and level = 1, we have:

The same calculation is applied to the similar factor with level = 2:

Step 6 above is important as a better dimension value is assigned to optimal particle P^o from two randomly chosen particles. This technique provides an interesting philosophy of vertical optimization which is very useful for many real world optimization problems. For better performance in regards to convergence, the two particles in step 1 of Taguchi method are randomly chosen from pBest vector. As of such, the Taguchi method will therefore improve the dimensions value among pBest, an effective way leading to faster convergence and better accuracy. We employ the similar methodology in this work for better performance.

IMPLEMENTS

For DOPF problem including FSWGs, there are T dispatches by N_gen-1 conventional generating units and many wind farms (for simplicity, one wind farm is adopted, many wind farms may be analogized). A particle array of control variable vectors is:

(13)

(14)

where, P is particle vector, g is the number of population particles, P_n^t is the real power output of nth generating unit at the tth interval, Q_c^t is the number of group switched capacitors in wind farm at the tth interval:

NUMERICAL RESULTS

To verify the effectiveness and efficiency of the adopted TPSO for DOPF problems including FSWGs, IEEE 39-bus power system is used as the test systems. The procedure has been implemented in MATLAB 7.0 programming language and numerical tests are carried on a Pentium 4 2.4 G computer. The wind farm including 60 wind generators with the same type, the rating power of which reaches 36 MW, is connected to the system at the bus 14. For the single FSWG, its rating voltage is 690 V, the excitation reactance x_m is 2.2059O, the sum of the stator reactance x₁ and rotor reactance x₂ are 0.1998O. For simplifying the analysis, the load size is considered invariable in the planning horizon. The planning horizon is divided into 12 intervals and every interval is 0.5 h. The wind generators' outputs are shown in Table 2. IEEE 39-bus system data are given by Zimmierman et al. (2009) and Victoire and Jeyakumar (2005). The parameters of the conventional generating units are shown in Table 3 and 4.

The maximum capacity of the reactive power compensation is 50 Mvar which constitutes of 50 set of group switched capacitors with per unit 1 Mvar. So, the upper limit value of the absorbed reactive power of wind farm bus is set to1Mvar in the study.

To demonstrate the superiority of the proposed approach for DOPF problems, simulation results have been compared with PSO method. Owing to the randomness in intelligent algorithms, two algorithms are executed 20 times when applied to the test system.

Table 2:	The wind farm data in different periods

Table 3:	The parameters of conventional generating units

Table 4:	The parameters and limits of conventional generating units

Table 5:	Best solution obtained using TPSO method
Total production cost : 457039.8613 $/h

Table 6:	Best solution obtained using PSO method
Total production cost: 458852.7267 $/h

For DOPF problem including FSWGs, Table 5 and 6 list the best control variables found by TPSO and PSO algorithm for IEEE 39-bus system, respectively. It is clearly shown that, by using TPSO, the total production cost savings of 1812.8654$/h is obtained compared with PSO algorithm. Hence, it is justified that TPSO approach gives the exact minimum dispatch solution. From Table 7, the best, worst and average cost values are 457039.8613 $/h, 458177.5291 $/h, 457518.5863 $/h and 458852.7267 $/h, 460137.6551 $/h, 459394.2806 $/h, respectively with TPSO and PSO after 20 independent trials. From the results, the superiority of TPSO strategies over PSO can be noticed. The difference between the best and worst solutions are 1137.6678 $/h with TPSO. At the same time, the difference between the best and worst solutions is 1284.9284$/h with PSO. Moreover, the best and worst solutions obtained by TPSO are very close to the average value, which proves that TPSO is more robust and consistent. In conclusion, it is clearly shown that TPSO is the most accurate and gives the exact minimum dispatch solution.

Table 7:	Comparison of best, worst and average cost values in the IEEE 39-bus system

Table 8:	Average execution time comparison in the IEEE 39-bus system

The average execution time taken to complete the fixed number of iterations (T_fix) and the average execution time taken to converge into the lower solution range (T_low) for 20 trials are shown in Table 8 for IEEE 39-bus system.

For IEEE 39-bus system, PSO takes an average execution time of 4235.83 sec to complete 150 iterations. PSO converges faster than TPSO by reason of the small sub-memeplex generation number of TPSO. In comparison to PSO, TPSO has additional components, i.e., the Taguchi method. This extra burdens increase the execution time of TPSO. TPSO takes 4317.79 sec more than PSO to complete 150 iterations. Nevertheless, TPSO takes only 2140.45 sec to converge into the lower solution range (457039-457107 $/h), PSO are not able to converge into the lower solution range.