A Dynamic Economic Dispatch Model Incorporating Wind Power Based on Chance Constrained Programming

In order to maintain the stability and security of the power system, the uncertainty and intermittency of wind power must be taken into account in economic dispatch (ED) problems. In this paper, a dynamic economic dispatch (DED) model based on chance constrained programming is presented and an improved particle swarm optimization (PSO) approach is proposed to solve the problem. Wind power is regarded as a random variable and is included in the chance constraint. New formulation of up and down spinning reserve constraints are presented under expectation meaning. The improved PSO algorithm combines a feasible region adjustment strategy with a hill climbing search operation based on the basic PSO. Simulations are performed under three distinct test systems with different generators. Results show that both the proposed DED model and the improved PSO approach are effective.


Introduction
Dynamic economic dispatch (DED), which determines the optimal generation scheme to meet the predicted load demand over a time horizon satisfying the constraint such as ramp-rate limits of generators between time intervals, is crucial for power system operation [1][2][3].Prior to the widespread OPEN ACCESS use of alternate sources (solar, wind) of energy, the DED problem involved only conventional generators.In recent years, wind power has experienced an explosive growth, and has shown great potential in fuel savings and environmental protection.However, its uncertainty and intermittency also makes it challenging to find a proper dispatch scheme for a wind penetrated power system.
The key issue associated with the incorporation of wind power is how to deal with its fluctuation and intermittence considering the required reliability and security of power systems [2].Up to now, different authors have proposed models to solve economic dispatch (or unit commitment) problems for wind-penetrated power systems [4][5][6][7][8][9][10][11][12][13][14][15][16].Chen [4] proposed a method to incorporate wind generators into the generation scheduling problem.Special reserve constraints are established to operate the power system within the required stability margin, which is also adopted by Zhou et al. [5] and Jiang et al. in [6].These models are deterministic and cannot accurately describe the effect of wind fluctuations on system operation.
To effectively and safely use wind power, researchers usually forecast the wind speed or wind generation over a time horizon in advance, and then obtain the statistic distribution of wind speed or wind generation.Based on the distribution function and estimated loads, researchers can determine the dispatch scheme.Hetzer et al. presented an ED model incorporating wind power with the stochastic wind speed characterization based on the Weibull probability density function [7].Penalty costs for overestimation and underestimation of available wind energy are also considered.In [8], a modified ED optimization model with wind power penetration is presented, and risk-based up and down spinning reserve constraints are presented considering the uncertainties of available wind power, load forecast error and also generator outage rates.It has been discussed in [9][10][11] that forecasted wind generation usually follows a beta distribution function, and following this, the authors of [12] and [13] proposed ED models incorporating the impact of wind variability.Miranda and Hang developed in [14] an ED model including wind generators using concepts from the fuzzy set theory, and they added a penalty cost factor for not using the available wind power capacity.It is noted clearly in references [7,8,[12][13][14] that the constraints of ramp rates are not taken into account.In [15] and [16], a scenario-based approach is utilized to model the uncertainty of wind power and stochastic models are proposed to solve the optimal scheduling of the generators in wind penetrated power system.This scenario approach requires a large number of scenarios to ensure the quality of solution, and usually suffers huge computation cost.
Chance constrained programming (CCP) is a kind of stochastic optimization approach [17].It is suitable for solving optimization problems with random variables either included in constraints or in the objective function.CCP has been studied to solve the transmission planning problem in [17,18] and the stochastic unit commitment problem with uncertain wind power output in [19].In [17] and [18], the same formulation of chance constraint is applied to transmission planning, and it is in the form that the not-overload probability for transmission line is required to be more than a specified confidence level.In [19], the chance constraint is applied to describe policies to ensure the utilization of wind power, and a two stage stochastic optimization is presented to solve the unit commitment problem.In this paper, the chance constraint is used to describe that the probability that actual wind generation is greater than or equal to the scheduled wind power is more than a given confidence level.
In a manner somewhat similar to the wind model in [12] and [13], this paper uses a beta distribution to characterize the actual wind generation for each individual period.Up spinning reserve (USR) and down spinning reserve (DSR) are required to deal with the forecasting errors in load and the sudden fluctuations in wind generation [15].In general, the reserve requirement with respect to forecast errors in load is defined to be a fixed percentage [5] (e.g., 5%) of the load demand.As to the reserve requirements caused by variation of wind generation, up reserve requirement (URR) induced by the sudden fall in wind generation and down reserve requirement (DRR) induced by the sudden raise in wind generation are both taken into account in our work.The conditional expectation when actual wind generation is smaller than the scheduled wind is calculated, and URR is defined as the difference between the scheduled wind power and the conditional expectation of wind power.Similarly, DRR is defined as the difference between the conditional expectation when actual wind generation is greater than or equal to the scheduled wind power and the scheduled wind power.
This paper proposes a stochastic DED model based on chance constraints in wind power penetrated systems.New formulations of the spinning reserve constraints are considered under the expected meaning.An improved particle swarm optimization (PSO) approach is proposed to optimize the model.In order to demonstrate the efficacy of the proposed model and the PSO approach, various comparisons are performed under two different test systems.Results show that the above-mentioned model is effective and the proposed PSO approach is able to solve such model.
The rest of this paper is organized as follows: in Section 2, the stochastic DED model is formulated.Section 3 introduces the beta distribution as the basis of the equivalent transformation of the proposed DED model, and a deterministic DED model is obtained.In Section 4, an improved PSO approach including feasible region adjustment (FRA) strategy and hill climbing search operation (HCSO) is applied to solve the deterministic model.Section 5 presents a discussion of the numerical results.Finally, conclusions are drawn in Section 6.

DED Problem Formulation
In this section, the objective function of DED problem is described and the DED model based on CCP is formulated.

Objective Function
In wind penetrated power systems, wind production is usually regarded as zero cost, the DED model with the valve point effect usually takes the following form [20][21][22][23][24]: where: cost f is the total generation cost over the whole time horizon; T is the number of periods; I is the number of thermal units; , i t p is the power output (MW) of the i th unit corresponding to time period t ; , ( ) C p is the generation cost of the i th unit corresponding to time period t ; , ( ) E p is the valve point loading effect of the i th unit corresponding to time period t ; For the thermal units, the generation cost can be approximated by a quadratic function of the power output, which is practical for most of the cases, and is given by: 2 , , , ( ) where i a , i b and i c are cost coefficients for the i th unit.
, ( ) i it E p is expressed as follows: where i e and i f are coefficients related to valve point effect of the ith unit.min i p is the minimum generation limit of unit i.

System and Unit Constraints
This DED problem is subjected to a variety of system and unit constraints, which include power balance constraints, generation limits of units, ramp rate limits and spinning reserve constraints.These constraints are discussed below.

Power Balance Constraints
Total power generation must equal the load demand pd,t in all time period: , , , 1 where pw,t is the scheduled wind power of wind farm at time t.

Generation Limits of Thermal Units and Wind Farm
The output of each thermal unit and wind farm must lie in between a lower and an upper bound.These constraints are represented as follows: , m a x 0 where max i p is the maximum generation limit of thermal unit i, and wmax is the installed capacity of wind farm.

Chance Constraint on Wind Power
Due to the stochastic nature of wind power, the schedule for wind may not be realized on a scheduled day.So we introduce the following chance constraint: where wt is a random variable representing the wind power generation at time t.ρ is the confidence level.Equation (7) defines the probability that the scheduled wind power can be realized is greater than or equal to ρ.Furthermore, Equation ( 7) sets a reasonable upper bound for wind power generation, and the probability that this upper bound can be realized is no less than ρ.The larger the confidence level is, the less stochastic the wind power is, and hence the more reliable the power system is.Especially, when ρ = 1, there is no wind power in system, and the DED problem is changed into a deterministic one.

Ramp Rate Limits of Thermal Units
The ramp rate limits restrict the operating range of all the units for adjusting the generation between two periods.The generation may adjust the dispatch level as: where Δ are the upper and lower ramp rate limits, respectively.T60 is the operating period, i.e., 1 h.

Spinning Reserve Constraints
Due to the penetration of wind power and its intermittent nature, additional reserve needs to be provided to support the wind generation variation [15].The USR supports the forecast errors in load and a sudden decrease in wind power, and the DSR contributes to the sudden raise in wind power.Both the USR and DSR are supplied by the ramping capacity of thermal units, and are formulated by Equations ( 9) and (10) The conditional expectation is obtained by the following equation: where f w is the probability density function of wind power at time t.
Similarly, , r is calculated using the two equations as mentioned in Equations ( 13) and ( 14): x E w w p = < Taking a closer look at the optimization model ( 1)-( 14), it contains chance constraint and thus is a stochastic DED model.In the next section we will discuss the equivalent transformation of the stochastic model to a deterministic one.

Beta Distribution of Wind Power
In the real word, the actual wind generation is a function of wind speed that randomly changes all the time [13].In our study, normalized wind power in each period is regarded as a random variable which follows beta distribution [9][10][11][12][13].We denote the probability density function (PDF) by ( ) t X t f x as follows: where max / t t x w w = , and [0,1] and ( , ) B α β is the well-known beta function defined as: Parameters α and β can be determined by following equations: where t μ and t σ are the mean value and standard deviation of the predicted wind generation at time t on the scheduling day, respectively.The cumulative density function (CDF) of X can be easily obtained as shown below:

Equivalent Transformation
Putting the aforementioned discussion of the beta distribution, we can get the following CDF ( ( ) Then, Equations ( 12) and ( 14) are transformed into the following equations as expressed in Equations ( 20) and ( 21 Both Equations ( 20) and ( 21) will be solved using numerical integration methods because of the complexity.In addition, we rewrite the Equation (7) Given all that, the equivalent transformation of the stochastic DED model is mathematically represented as:

Improved PSO Approach
The problem in Equation ( 24) is a high-dimension and non-convex optimization problem, so it is very difficult to find analytical solutions.In recent decades, many salient approaches have been developed to solve such problems, such as genetic algorithm [20,21], differential evolution [22], evolutionary programming [25,26], and PSO [6,[27][28][29][30][31][32].PSO, first introduced by Kennedy and Eberhart, is a population-based optimization technique, and conducts its search using a population of particles [27,28].Each particle is a candidate solution to the problem and is moved toward the optimal point by adding a velocity with its position.The position and the velocity of the jth particle in the D dimensional search space can be expressed as [ , , , ] [ , , , ] respectively.Each particle has its own best position ( , 1,2,  ) corresponding to the personal best objective value obtained at generation k.The global best particle is denoted by k gbest , which represents the best particle found so far at generation k among the whole population.The new velocity and position of each particle at generation k + 1 are calculated as shown below [6,29]: where: J is the population size; ω( ) k is the dynamic inertia weight factor, and can be dynamically set with the following equation [6]: where max ω and min ω are initial and final inertia weight factors and set to 0.9 and 0.4 respectively.K is the maximum iteration number. 1 ( ) k ϕ and 2 ( ) k ϕ are time-varying acceleration coefficients corresponding to cognitive and social behavior [6], and are set with the following equations shown in Equations ( 28) and ( 29): where 1i ϕ , 2i ϕ are the initial values of 1 ( ) k ϕ and 2 ( ) k ϕ , and are set to 2.5 and 0.5 respectively; , and are set to 0.5 and 2.5 respectively.
It is not always effective to solve the problem with equality and inequality constraints using basic PSO, and the particles (solutions) which satisfy the inequality constraints usually violate the equality constraints.In this section, we propose a FRA strategy over these particles which violate equality constraints.In addition, in order to improve the quality of the best solution, HCSO is applied to update the global best particle along with the iteration.In subsections 4.1 and 4.2, we will separately discuss the FRA strategy and the HCSO in detail.

Feasible Region Adjustment Strategy
Let us rewrite the position of the jth particle as the following matrix: where M is the number of generators including thermal units and one wind farm, i.e., M = I + 1, and than the demand, the FRA strategy will be applied to the jth particle.The FRA strategy is described as follows: ( In Equation ( 31    From above discussion, we can learn that the FRA strategy ensures that the particles after adjusted still satisfy the generation limit constraints.

Hill Climbing Search Operation
Hill climbing is an optimization technique which belongs to the local search family.It is an iterative algorithm that starts with an arbitrary solution, and then attempts to find a better solution by incrementally changing a single element of the solution.In our work, in order to obtain a better solution, we adopt HCSO to update the global best particle along with the iteration.When the fitness of the global best particle does not change continuously for a fixed number of times, HCSO can be used as explained in Equations ( 35 Δ and n Δ are the upper ramp rates of unit m and n.In this paper, we adopt the upper ramp rates in simulation.However, the lower ramp rates can be used either.min( , ) m n Δ Δ ⋅ε is the linear decrease step size of hill climbing operation and ε is listed in Table 2 in Section 5.Both Equations ( 35) and (36) have the same step size and the opposite direction, which ensures that the global best particle will not violate equality constraints after hill climbing operation.

DED Constraints Handling Using PSO-HCSO
The fitness function of particle is the generation cost bounded with the penalty functions as shown in Equation (37): where λ is the penalty factor corresponding to the constraints.Here, we use uniform penalty factor for , g t PF , i.e., 1 8 e λ = + ; , g t PF is the penalty function.It is noted that the solution should not contain any penalty for the constraint violation.

The Procedure of Improved PSO Approach
The procedures for implementing the PSO approach are given as the following steps: Step 1: Initialize the parameters, such as population size J, the maximum iteration number K, and the maximum number of hill climbing search operation H. Set the sequence number of iteration K = 1 and the number of times that the global best particle does not change continuously 1 hcount = .
Step 2: Create a swarm of particles as the initial population, including random position and velocity.For any particle which violates the equality constraints, FAR strategy is utilized.Evaluate the fitness of particles and obtain the initial 0 gbest and 0 j pbest , 1, 2, , j J =  .
Step 3: Calculate ( ) , and then update the position and velocity of each particle among the population according to Equations ( 25) and (26).FAR strategy is applied to any particle which violates the equality constraints. Step

Simulation Results and Discussions
In order to verify the effectiveness of the proposed DED model with wind power, three distinct test systems (i.e., system 1, system 2 and system 3) are employed in this paper.System 1 contains six thermal units and one wind farm (WF) which is also used in system 2 and 3. System 2 contains 15 thermal units.The characteristics of the 6 and 15 thermal units can be obtained from [27].System 3 has 26 thermal units derived for a IEEE 24-bus system.The wind farm is located in the Inner Mongolia Autonomous Region in China, and contains 132 wind turbines.The total installed capacity is 198 MW.Table 1 gives the hourly expected value and standard deviation of wind power forecast on the scheduled day, and also the corresponding beta parameters which are calculated by Equation (17).
The parameters of improved PSO are shown in Table 2.In the table, h is the sequence number of Hill climbing search operation.The improved PSO approach has been implemented on a personal computer with four processors running at 3.2 GHz and equipped with 4 GB of RAM memory using Matlab 7.9.0.

Comparisons Among the Three Cases of the DED Model
In each test system, we perform 50 trials using the improved PSO approach considering three cases as follows: Case (1): the DED model without considering wind power; Case (2): the DED model without considering wind effect in the reserve constraint; Case (3): the proposed DED model in this paper.
Both casea (1) and ( 2) can be derived from case (3).The first case can be obtained by setting the confidence level (ρ) of case (3) to 1, and the latter is obtained by removing the URFW and DRRW in the reserve constraint.Comparison results among three cases mentioned above are shown in Table 3.It is evident that the average generation cost in case ( 1) is the highest among the three systems; this is because that no wind power is utilized in case (1).Wind power in cases ( 2) and (3) affords a certain proportion of the load and thus reduces the fuel consumption.Average generation cost decreases gradually as the confidence level decreases for cases (2) and (3) among the three systems.The smaller the confidence is, the lower the average generation cost is.Taking system 1 for example, the average cost in cases ( 2) and (3) decreases by 5.02% and 4.93%, respectively, compared to case (1) when ρ = 0.9.When ρ decreases to 0.1, the average generation cost decreases by 9.00% and 7.78%, respectively, compared to case (1).
Figures 2-4 show the average cost change trend for different confidence level scenarios for cases (2) and (3) among the three systems.As seen in the figures, the average cost in the three systems increases gradually while the confidence level increases.There is not much difference in the average cost between case (3) and case (2), however, case (3) considers the URFW and DRRW in the reserve constraints, and this makes the systems more reliable.We choose the solutions with the minimum generation cost from 50 trials as the optimal solutions by varying ρ from 0.1 to 0.9.The optimal solutions (power output during each period for each unit and wind farm) for system1 and 2 under ρ = 0.9 are shown in Figures 5 and 6.Power outputs for period 1, 5, 9, 13, 17, 21 in system 3 are listed in      ) under the optimal solutions in each system are shown in Figures 7-9.As shown, wind penetration in each system under ρ = 0.9 is the lowest, because under that situation each system requires the highest reliability and security.As the confidence level decreases, wind penetration almost increases and does not exceed the wind penetration limits in the three systems.

Discussion about the Optimal Reserve Allocation
Due to the stochastic and intermittent nature of wind power, the USR and DSR are used to ensure the reliability and security of systems with wind farms.According to Inequations ( 9) and ( 10), the USR, DSR, URFW and DRRW for each period are calculated based on the optimal solutions when ρ = 0.9 serves as an example.Figures 10 and 11 show the optimal USR and DSR allocation for 24 h in system 1 respectively.Figures 12 and 13 give the optimal USR and DSR allocation for 24 h in system 2. It can be learned that the USR and DSR provided by thermal units in both systems can effectively cover the sudden fall and increase in wind power.The same conclusion can also be drawn for system 3, so the corresponding figures for USR and DSR are not displayed.
Figure 10.The optimal USR allocation and URR plus reserve for load forecast error in system 1 with ρ = 0.9.

Figure 11.
The optimal DSR allocation and DRR for 24 h in system 1 with ρ = 0.9.The optimal USR allocation and URR plus reserve for load forecast error for 24 h in system 2 with ρ = 0.9.

Figure 13.
The optimal DSR allocation and DRR for 24 h in system 2 with ρ = 0.9.

Comparisons between the PSO-HCSO and PSO without HCSO
In order to investigate the effect of the parameter H on the improved PSO algorithm, we set H as 100, 200, 300, 400 and perform the algorithm on system 1 with ρ = 0.9 for 50 trials.The average generation cost and average time are listed in Table 5.The average generation costs decrease gradually and average time increases as H increases.

Conclusions
Wind power provides energy savings and environmental protection benefits.However, the intermittency and uncertainty of wind power generation require that conventional units provide additional reserve to ensure the stability and reliability of a wind power-penetrated power system.This paper formulates a DED model incorporating wind power based on CCP.The uncertain nature of wind generation is represented by a beta distribution function.In order to ensure the reliability of the power system, a chance constraint is included, and conditional expectation is presented to calculate the up and down spinning reserves.The proposed DED model is then numerically solved using an improved PSO approach in three different test systems.Results show that the proposed model can effectively respond to sudden wind power falls or raises.The improved PSO approach is fit to solve the DED model.
The results also show that the average generation cost and wind penetration are dependent on the confidence level.If the confidence level is increased, the wind penetration will be reduced, which results in higher reliability of the power system, and the average generation cost will be increased.Conversely, if the confidence level is decreased, less reliability allows more wind power to be incorporated in the power system.The average generation cost will be reduced.

Author Contributions
All the authors contributed to this paper.Haifeng Zhang formulated the DED model, performed experiments and edited the manuscript.Wushan Cheng reviewed and revised the manuscript.

Figure 1 .
Figure 1.Concept illustration of , u w t r and , d w t r .

Figure 2 .
Figure 2. Average generation cost under different confidence level in system 1. 0

Figure 3 .
Figure 3. Average generation cost under different confidence level in system 2.

Figure 4 .
Figure 4. Average generation cost under different confidence level in system 3.

Figure 5 .
Figure 5. Optimal generation output of units and wind farm in system 1 with ρ = 0.9.

Figure 6 .
Figure 6.Generation output of units and wind farm in system 2 with ρ = 0.9.

Figure 7 .
Figure 7. Wind penetration and its limit under different confidence levels in System 1.

Figure 8 .
Figure 8. Wind penetration and its limit under different confidence levels in System 2.

Figure 9 .
Figure9.Wind penetration and its limit under different confidence levels in system 3.

Figure 12 .
Figure12.The optimal USR allocation and URR plus reserve for load forecast error for 24 h in system 2 with ρ = 0.9. ) p is the reserve level to support forecast error in demand.T10 is 10 min.maxrare the URFW and the DRRW to follow the sudden decrease and increase in wind power at time t.Figure1shows the conceptual illustration of , r is calculated by the difference between the scheduled wind power and conditional expectation of wind power when the actual wind generation is less than scheduled, and is shown in Equation (11):, , , as follows: The element of matrix ymt is the output of the mth generator at time t.The row vectors represent the output of the individual generator each hour one day.The sum of the output of individual generator (i.e., power supply) each hour must be equal to the demand, i.e., We use the following equation to adjust the output of the thermal units and the wind farm as shown in Equation (31): )

Table 1 .
Hourly expected value and standard deviation of wind power forecast and beta parameters.

Table 2 .
Parameters for improved PSO algorithm.

Table 3 .
Comparisons among the three cases in the three systems.

Table 4 .
Optimal output of units and wind farm (WF) in system 3 with ρ = 0.9.

Table 5 .
Effect of H on the improved PSO algorithm with ρ = 0.9.
Power output of thermal unit i at time t Scheduled wind power of wind farm at time t Coefficients related to valve point effect of thermal unit i Minimum and Maximum generation limits of thermal unit i Acceleration coefficients corresponding to cognitive and social behavior Global best position at generation k in the whole population Maximum number of the HCSO hcountThe number of times that gbest does not change continuously t w Actual wind generation, a random variable max w k gbest