PSO-Based Multi-Criteria Economic Dispatch Considering Wind Power

Document Sample
PSO-Based Multi-Criteria Economic Dispatch Considering Wind Power Powered By Docstoc
					      PSO-Based Multi-Criteria Economic Dispatch
     Considering Wind Power Penetration Subject to
                 Dispatcher’s Attitude
                                               Lingfeng Wang and Chanan Singh
                                       Department of Electrical and Computer Engineering
                                                     Texas A&M University
                                                   College Station, TX 77843

   Abstract— Significant attention has been paid to the renewable     desired. Particle swarm optimization (PSO) is a salient meta-
energy resources such as wind power in recent years. It has          heuristics which has drawn much attention recently due to
potential benefits in curbing emissions and reducing the con-         its ability to solve highly nonlinear and complex engineering
sumption of irreplaceable fuel reserves. However, the penetration
of wind power into traditional fuel-based generation systems         optimization problems as well as its outstanding convergence
will also have some implications such as security concerns due       performance. In this study, a multi-objective particle swarm
to its unpredictable nature. Thus, in economic power dispatch        optimization (MOPSO) algorithm is developed to derive the
with wind power penetration, a reasonable tradeoff between           tradeoff solutions for economic dispatch accounting for wind
system risk and operational cost is desired. In this paper, a bi-    penetration. Furthermore, considering the different attitudes of
objective economic dispatch problem considering wind penetra-
tion is formulated, which treats economic and security impacts       dispatchers towards the wind power integration, we used the
as conflicting objectives. Different fuzzy membership functions       fuzzy constraints to indicate the security level in terms of wind
are used to reflect the dispatcher’s attitude toward the wind         penetration and wind power cost. Different fuzzy constraints
power penetration. A modified multi-objective particle swarm          including linear and quadratic functions can be used to reflect
optimization (MOPSO) algorithm is adopted to develop a power         the dispatcher’s optimistic, neutral, or pessimistic attitude
dispatch scheme which is able to achieve compromise between
economic and security requirements. Numerical simulations in-        toward wind penetration.
cluding sensitivity analysis are reported based on a typical IEEE       The remainder of the paper is organized as follows. Section
test power system to show the validity and applicability of the      2 presents the wind penetration model. Section 3 formulates
proposed approach.                                                   the target economic dispatch problem including its multiple
                                                                     objectives as well as design constraints. Section 4 introduces
                                                                     the inner working of particle swarm optimization algorithms.
                      I. I NTRODUCTION                               The proposed MOPSO algorithm is discussed in Section 4.
   Economic dispatch (ED) is a crucial function in the mod-          Simulation results and analysis are presented in Section 5. Fi-
ern Energy Management System (EMS). Its objective is to              nally, conclusions are drawn and future research is suggested.
schedule the power generation properly in order to minimize
the total operational cost [7-10, 14, 15]. More recently, wind                 II. W IND P OWER P ENETRATION M ODEL
power has drawn much attention as a promising renewable                 Because of its unpredictable and variable characteristics,
energy resource, which has shown some prospects in cutting           the integration of wind power into the traditional thermal
fuel consumption and reducing the pollutants emission [1-6,          generation systems will inevitably bring about the operator’s
11, 17, 18, 20]. However, accurate forecasting of the expected       concern on system security. Fuzzy definition regarding wind
generation from a wind park is nearly impossible primarily           penetration is a viable way to represent the penetration level
due to the stochastic nature of the wind as well as the highly       of the wind power, since it is usually difficult to determine the
nonlinear transform from wind speed to electrical energy. This       optimal amount of wind power that should be integrated into
inevitably brings about the security concerns about incorpo-         the conventional power grids [16].
rating wind power in the traditional power system, and risk             As shown in Figure 1, a membership function µ regarding
mitigation issues should be properly addressed. For instance,        the wind penetration is defined to indicate the system security
the dynamic stability may be violated due to wind fluctuations.       level. It can be mathematically expressed in the following [16]:
Therefore, in order to achieve the compromise between system
                                                                               1,                        W ≤ W (PD )min
risk and total running cost, it is highly desirable to investigate
                                                                                      W (PD )max −W
how to dispatch the power in an appropriate manner for the              µ=                             , Wmin ≤ W ≤ Wmax
                                                                               W (PD )max −W (PD )min
power system with wind penetration. In this paper, the problem                    0,                      W ≥ W (PD )max
is formulated as a bi-objective optimization problem through                                                                    (II.1)
simultaneous minimization of both risk level and operational         where W is the wind power incorporated in economic dis-
cost. For this purpose, an effective optimization algorithm is       patch; W (PD )min is the lower bound of wind power penetra-
         µ                                                                            µ

         1                                                                             1


         0                                                                             0

                         W ( PD ) min    W ( PD ) max          W                                  W ( PD ) min W1    W2    W3 W ( PD ) max        W
                        WC ( PD ) min   WC ( PD ) max         WC

Fig. 1. Fuzzy linear representation of the security level in terms of wind   Fig. 2. Fuzzy quadratic representation of the security level in terms of wind
penetration and wind power cost.                                             power penetration.

tion, below which the system is deemed secure; W (PD )max is                 shape is shown in Figure 3.
the upper bound of wind power penetration, above which the                        
system is considered as insecure due to the wind perturbations.                    1,                            W C ≤ W C(PD )min
                                                                             µ=      ac W C 2 + bc W C + cc , W Cmin ≤ W C ≤ W Cmax
Both W (PD )min and W (PD )max are dependent on the total
                                                                                     0,                           W C ≥ W C(PD )max
load demand in the power dispatch.
   The above defined membership function can also be repre-                   where ac , bc , and cc determine the curve shape of the quadratic
sented in terms of the operational cost for incorporating wind               function defined in terms of the running cost of wind power.
         1,                     W C ≤ W C(PD )min
              W Cmax −W C                                                              1
   µ=        W Cmax −W Cmin ,    W Cmin ≤ W C ≤ W Cmax
           0,                    W C ≥ W C(PD )max
where W C is the running cost of wind power in the power
dispatch; W C(PD )min is the lower bound cost for producing
                                                                                                 WC ( PD ) min WC1   WC2   WC3    WC ( PD ) max   WC
wind power, below which the system is seen as secure;
W (PD )max is the upper bound cost for including wind power,
above which the system is considered as insecure due to the
wind intermittency. In a similar fashion, both W C(PD )min                   Fig. 3. Fuzzy quadratic representation of the security level in terms of wind
                                                                             power cost.
and W C(PD )max are dependent on the total load demand in
the power dispatch. In this study, sensitivity studies are also
carried out to illustrate the impact of different allowable ranges
of wind power penetration as well as different running costs                                    III. P ROBLEM F ORMULATION
of wind power on the final solutions obtained.                                   The problem of economic power dispatch with wind pen-
   To reflect dispatcher’s different attitudes toward wind power              etration consideration can be formulated as a bi-criteria op-
penetration, a quadratic membership function can be defined                   timization model. The two conflicting objectives, i.e., total
as follows [16]:                                                             operational cost and system risk level, should be minimized
                                                                             simultaneously while fulfilling certain system constraints. This
        1,                          W ≤ W (PD )min                          bi-objective optimization problem is formulated mathemati-
    µ=   aw W 2 + bw W + cw ,        Wmin ≤ W ≤ Wmax                         cally in this section.
         0,                          W ≥ W (PD )max
                                                                             A. Problem objectives
where aw , bw , and cw are the coefficients of the quadratic
function, which determine its curve shape reflecting the dis-                   There are two objectives that should be minimized simulta-
patcher’s attitude toward wind power. As shown in Figure 2,                  neously, that is, system security level and the total operational
by selecting different coefficients aw , bw , and cw , different              cost.
curve shapes of the quadratic function can be defined. For                      • Objective 1: Minimization of system risk level
the identical risk level µ0 , the wind power costs for different                   From the security level function defined in (II.1) and (II.
defined functions w1 < w2 < w3 . The curves corresponding                           2), we know that the larger the value of membership
to these three values reflect the pessimistic, neutral, and                         function µ is, the more secure the system will become.
optimistic attitudes of the dispatcher toward the wind power                       If the wind penetration is restricted under a certain level,
integration, respectively.                                                         the system can be considered as secure. On the contrary,
   In a similar fashion, the security level can also be defined                     if excessive wind penetration is introduced into the power
in terms of the operational cost of wind power. Its function                       dispatch, the system may become insecure. Thus, here we
      define an objective function which should be minimized                       For normal system operations, real power output of each
      in order to ensure system security:                                         generator is restricted by lower and upper bounds as
                                          1                                                           min            max
                                R(µ) =                            (III.1)                          PGi ≤ PGi ≤ PGi                 (III.8)
                                                                                           min         max
                                                                                  where PGi and PGi are the minimum and maximum
  •   Objective 2: Minimization of operational cost
                                                                                  power from generator i, respectively.
      The cost curves of different generators are represented
                                                                              •   Constraint 2: Power balance constraint
      by quadratic functions with sine components. The super-
                                                                                  The total power generation and the wind power must
      imposed sine components represent the rippling effects
                                                                                  cover the total demand PD and the real power loss
      produced by the steam admission valve openings. The
                                                                                  in transmission lines PL . For the linear membership
      total $/h fuel cost F C(PG ) can be represented as follows:
                                                                                  function, this relation can be represented by
               F C(PG ) =              ai + bi PGi +       2
                                                       ci PGi                                PGi + Wmax − µ ∗ ∆W = PD + PL                 (III.9)
                                           min                                    For the quadratic membership function, the relation can
                            + |di sin[ei (PGi − PGi )]|           (III.2)
                                                                                  be expressed by
      where M is the number of generators committed to the
                                                                                   M                                   b2
      operating system, ai , bi , ci , di , ei are the cost coefficients                       bw         µ − (cw −    4aw )

      of the i-th generator, and PGi is the real power output of                       PGi −     ±                            = PD +PL (III.10)
                                                                                             2aw               aw
      the ith generator. PG is the vector of real power outputs
      of generators and defined as                                                 The sign of the last term in (III.10) is determined by
                                                                                  the curve shape of the defined quadratic function. The
                      PG = [PG1 , PG2 , . . . , PGM ]             (III.3)         transmission losses can be calculated based on the Kron’s
      The running cost of wind power can be represented in                        loss formula as follows:
      terms of the value of membership function µ which indi-                                M   M                   M
      cates the system security level. For the linear membership                   PL =              PGi Bij PGj +         B0i PGi + B00 (III.11)
      function case,                                                                       i=1 j=1                   i=1

                                                             M                    where Bij , B0i , B00 are the transmission network power
       W C(PG , µ)      = Cw (Wav − (PD + PL −                   PGi ))           loss B-coefficients. It should be noted that the transfer
                                                             i                    loss of the wind power is not considered in this study.
                        − µ ∗ ∆W C + W Cmax                       (III.4)     •   Constraint 3: Available wind power constraint
                                                                                  The wind power used for dispatch should not exceed the
      where Wav is the available wind power from the wind
                                                                                  available wind power from the wind park:
      farm, Cw is the coefficient of penalty cost for not using
      all the available wind power, PD is the load demand, PL                                                 M

      is the transmission loss, and                                                          0 ≤ PD + PL −         PGi ≤ Wav              (III.12)
                     ∆W C = W Cmax − W Cmin .                     (III.5)     •   Constraint 4: Security level constraint
      For the quadratic membership function case,                                 From the definition of membership function shown from
                                                                                  (II.1) to (II.4), the values of µ should be within the
                                                                                  interval of [0, 1]:
       W C(PG , µ)      = Cw (Wav − (PD + PL −                   PGi ))
                                                                                                          0 ≤ µ ≤ 1.             (III.13)
                              bc         µ − (cc −   4ac )
                                                                            C. Problem statement
                        −        ±                                (III.6)
                             2ac               ac
                                                                               In summary, the objective of economic power dispatch op-
      The sign of the last term in (III.6) is determined by the             timization considering wind penetration is to minimize R(µ)
      curve shape of the defined quadratic function. Thus, the               and T OC(PG , µ) simultaneously subject to the constraints
      total operational cost T OC can be calculated as                      (III.7)–(III.11).
           T OC(PG , µ) = F C(PG ) + W C(PG , µ)                  (III.7)
                                                                              IV. M ECHANISM OF PARTICLE S WARM O PTIMIZATION
                                                                               Particle swarm optimization (PSO) is a population-based
B. Problem constraints
                                                                            stochastic optimization technique, which was inspired by the
   Due to the physical or operational limits in practical sys-              movement pattern in a bird flock or fish school [12, 13]. In
tems, there is a set of constraints that should be satisfied                 PSO, individuals (i.e., particles) move around in a multidi-
throughout the system operations for a feasible solution.                   mensional search space to approach the optima, where each
   • Constraint 1: Generation capacity constraint                           point represents a solution to the target problem. Initially a
bunch of particles are randomly created and set into motion                                     V. T HE P ROPOSED A PPROACH
through this space. In their movement, each particle adjusts                     The standard PSO algorithm is not suited to resolve multi-
its position based on its own experience as well as the expe-                 objective optimization problems in that no absolute global op-
rience of a neighboring particle by utilizing the best position               timum exists there, but rather a set of non-dominated solutions.
encountered by itself and its neighbors. At each generation,                  Thus, to render the PSO algorithm capable of dealing with MO
they observe the “fitness” of themselves and their neighbors                   problems, some modifications become necessary. When using
and move toward those with a better position. In this way, PSO                stochastic search based algorithms to optimize multi-objective
combines both local and global search methods together in                     problems, two key issues usually arise in the algorithm design.
order to improve its search effectiveness and efficiency. Unlike               First, the fitness evaluation should be suitably designed to
other evolutionary computation algorithms including Genetic                   guide the search toward the set of Pareto-optimal solutions.
Algorithms (GA), PSO has no evolution operators such as                       Second, the diversity of the population should be maintained
crossover and mutation. The optima is obtained via following                  by refraining the search from premature convergence. In this
the current optimum particles by the potential particles. This                study, the classic PSO algorithm is revised accordingly to
simple algorithm turns out to be highly effective in a diverse                facilitate a multi-objective optimization approach, i.e., multi-
set of optimization problems.                                                 objective particle swarm optimization (MOPSO). The Pareto-
   Let x and v denote a particle coordinates (posi-                           dominance concept is used to evaluate the fitness of each
tion) and its corresponding flight speed (velocity) in the                     particle and thus determine which particles should be selected
search space. Therefore, the i-th particle is represented as                  to store in the archive of non-dominated solutions. Somehow
xi = [xi1 , xi1 , . . . , xid , . . . , xiM , xi,M +1 ] in the (M + 1)-       similar to the elitism used in evolutionary algorithms and the
dimensional space. Each particle keeps track of its coor-                     tabu list used in tabu searches, the best historical solutions
dinates in the solution space which are associated with                       found by the population are recorded continuously in the
the best solution it has achieved so far. This fitness                         archive in order to serve as the non-dominated solutions
value is called pbest. The best previous position of the                      generated in the past. Furthermore, due to the global attraction
i-th particle is recorded and represented as pbesti =                         mechanism in PSO, the historical archive of previously found
[pbesti1 , pbesti2 , . . . , pbestid , . . . , pbestiM , pbesti,M +1 ]. An-   non-dominated solutions would make the search converge
other “best” value that is tracked by the particle swarm                      toward globally non-dominated solutions highly possible.
optimizer is the best value obtained so far by any particle
in the neighbors of the particle. When a particle takes all the
                                                                              A. Constrained PSO
population as its topological neighbors, the best value is a
global best and is called gbest. The index of the best particle                  Because the standard PSO does not take into account how
among all the particles in the group is represented by the                    to deal with the constraints, the constraints handling mech-
gbestd . The rate of the velocity for particle i is represented               anism should be added to ensure the solution feasibility in
as vi = (vi1 , vi2 , . . . , vid , . . . , viM , vi,M +1 ). The modified       constrained optimization problems such as power dispatch. In
velocity and position of each particle can be calculated using                the proposed MOPSO, a simple constraint checking procedure
the current velocity and the distance from pbestid to gbestd                  called rejecting strategy is incorporated. When an individual is
as shown in the following formulas:                                           evaluated, the constraints are first checked to determine if it is
                                                                              a feasible candidate solution. If it satisfies all of the constraints,
                                                                              it is then compared with the non-dominated solutions in the
       (t+1)               (t)                                   (t)
      vid        = w ∗ vid + c1 ∗ rand() ∗ (pbestid − xid )                   archive. Here the concept of Pareto dominance is applied
                 + c2 ∗ Rand() ∗ (gbestd − xid ),
                                                                  (IV.14)     to determine if it is eligible to be chosen to store in the
                                                                              archive of non-dominated solutions. The constraint satisfaction
                                                                              checking scheme used in the proposed algorithm proves to be
                                                                              quite effective in ensuring the feasibility of the non-dominated
 (t+1)         (t)     (t+1)
xid      = xid +χ∗vid     , i = 1, 2, . . . , N, d = 1, 2, . . . , M +1.      solutions.
where N is the number of particles in a group, M + 1                          B. Individual (particle) representation
is the number of members in a particle, t is the pointer
                                                                                 It is crucial to appropriately encode the individuals of
of generations, χ ∈ [0, 1] is the constriction factor which
                                                                              the population in PSO for handling the economic dispatch
controls the velocity magnitude, w is the inertia weight factor,
                                                                              application. The power output of each generating unit and
c1 and c2 are acceleration constants, rand() and Rand() are
                                                   (t)                        the value of membership function are chosen to represent
uniform random values in a range [0, 1], vi is the velocity
                                         (t)                                  the particle position in each dimension, and positions in
of particle i at generation t, and xi is the current position
                                                                              different dimensions constitute an individual (particle), which
of particle i at generation t. The particle swarm optimization
                                                                              is a candidate solution for the target problem. The position in
concept consists of, at each time step, changing the velocity
                                                                              each dimension is real-coded. The i-th individual PGi can be
of (accelerating) each particle toward its pbest and gbest
                                                                              represented as follows:
locations. Acceleration is weighted by a random term, with
separate random numbers being generated for acceleration                      PGi = [PGi1 , PGi2 , . . . , PGid , . . . , PGiM , µi ], i = 1, 2, . . . , N
toward pbest and gbest locations.                                                                                                                 (V.16)
where M is the number of generators and N is the population                          For the value of membership function µ,
size; PGid is the power generated by the d-th unit in i-th                                             (t+1)      (t+1)         (t+1)
individual; and µi is the value of the membership function in                                        µi        = µi       + RT µi       (V.22)
i-th individual. Thus, the dimension of a population is N ×                          where RT is the turbulence factor. The turbulence term
(M + 1).                                                                             is used here to enhance the diversity of solutions.
                                                                                 •   Step 11: Update the archive which stores non-dominated
C. Algorithm steps                                                                   solutions according to Pareto-optimality based selection
   The proposed MOPSO is applied to the constrained power                            criteria [19].
dispatch problem in order to derive out the optimal or near-                     •   Step 12: If the current individual is dominated by the
optimal solutions. Its computational steps include:                                  pbest in the memory, then keep the pbest in the memory;
                                                                                     Otherwise, replace the pbest in the memory with the
   • Step 1: Specify the lower and upper bounds of generation
                                                                                     current individual.
     power of each unit as well as the range of security level.
                                                                                 •   Step 13: If the maximum iterations are reached, then go
   • Step 2: Randomly initialize the individuals of the popu-
                                                                                     to Step 14. Otherwise, go to Step 7.
     lation. Note that the speed and position of each particle
                                                                                 •   Step 14: Output a set of Pareto-optimal solutions from
     should be initialized such that each candidate solution
                                                                                     the archive as the final solutions.
     (particle) locates within the feasible search space.
   • Step 3: For each individual PGi of the population, the
     transmission loss PLi is calculated based on B-coefficient                  VI. S IMULATION AND E VALUATION OF THE P ROPOSED
     loss formula.                                                                                 A PPROACH
   • Step 4: Evaluate each individual PGi in the population                      In this study, a typical IEEE 30-bus test system with 6-
     based on the concept of Pareto-dominance.                                generators [19] is used to investigate the effectiveness of
   • Step 5: Store the non-dominated members found thus far                   the proposed MOPSO approach. The system configuration is
     in the archive.                                                          shown in Figure 4. The system parameters including fuel cost
   • Step 6: Initialize the memory of each particle in which a                coefficients and generator capacities are listed in Table I. The
     single local best pbest is stored. The memory is contained               sinusoidal term in (III. 2) is not considered in this study due
     in another archive.                                                      to its relatively minor impact on the total fuel costs. The B-
   • Step 7: Increment iteration counter.                                     coefficients are shown in (VI.23). The load demand PD used
   • Step 8: Choose the personal best position pbest for                      in the simulations is 2.834 p.u., the available wind power Wav
     each particle based on the memory record; Choose the                     is 0.5668 p.u, and the coefficient of penalty cost Cw is set 20
     global best gbest from the fuzzified region using binary                  $/p.u.
     tournament selection [19]. The niching and fitness sharing
     mechanism is also applied throughout this process for
     enhancing the diversity of solutions.
   • Step 9: Update the member velocity v of each individual
     PGi based on (IV.14). For the output of each generator,
        (t+1)             (t)                                     (t)
       vid      = w ∗ vi + c1 ∗ rand() ∗ (pbestid − PGid )
                + c2 ∗ Rand() ∗ (gbestd − PGid ),
                  i = 1, . . . , N ; d = 1, . . . , M             (V.17)
      where N is the population size, and M is the number of
      generating units. For the value of membership function
       (t+1)              (t)                                           (t)
      vi,M +1   = w ∗ vi + c1 ∗ rand() ∗ (pbesti,M +1 − µi )
                + c2 ∗ Rand() ∗ (gbestM +1 − µi ),                 (V.18)
  •   Step 10: Modify the member position of each individual
      PGi based on (IV.15). For the output of each generator,                 Fig. 4.   IEEE 30-bus test power system
                      (t+1)          (t)          (t+1)
                     PGid       =   PGid    +   χvid              (V.19)
      For the value of membership function µ,
                                                                              A. Comparison
                       (t+1)         (t)         (t+1)
                      µi        =   µi     +   χvi,M +1           (V.20)         Since PSO algorithms are sometimes quite sensitive to
      Following this, add the turbulence factor into the current              certain parameters, the simulation parameters should be appro-
      position. For the output of each generator,                             priately chosen. In the simulations, both the population size
                                                                              and archive size are set to 100, and the number of generations
                     (t+1)          (t+1)           (t+1)
                   PGid      = PGid         + RT PGid             (V.21)      is set to 500. The acceleration constants c1 and c2 are chosen
                                                                                                                                  
                                         0.1382   −0.0299   0.0044      −0.0022                        −0.0010      −0.0008
                                       −0.0299   0.0487    −0.0025     0.0004                          0.0016       0.0041        
                                        0.0044   −0.0025   0.0182      −0.0070                        −0.0066      −0.0066        
                             Bij =                                                                                                                                (VI.23)
                                       −0.0022   0.0004    −0.0070     0.0137                          0.0050       0.0033        
                                       −0.0010   0.0016    −0.0066     0.0050                          0.0109       0.0005        
                                        −0.0008   0.0041    −0.0066     0.0033                          0.0005       0.0244

                               TABLE I
       F UEL COST COEFFICIENTS AND GENERATOR CAPACITIES                                         1
                                                                                                                                 Quadratic function (Pessimistic)
                                                                                               0.9                               Linear function (Neutral)
          Generator i   ai    bi    ci      min
                                           PGi     max
                                                  PGi                                                                            Quadratic function (Optimistic)
          G1            10    200   100    0.05   0.50                                         0.8

          G2            10    150   120    0.05   0.60                                         0.7
          G3            20    180   40     0.05   1.00
          G4            10    100   60     0.05   1.20

                                                                              Security level
          G5            20    180   40     0.05   1.00
          G6            10    150   100    0.05   0.60


as 1. Both turbulence factor and niche radius are set to 0.02.                                 0.2
The inertia weight factor w decreases when the number of
generations increases:
                         wmax − wmin                                                            0.25   0.3   0.35   0.4     0.45     0.5     0.55      0.6      0.65
          w = wmax −                     × iter          (VI.24)                                                      Wind power (p.u.)
where itermax is the number of generations and iter is the
                                                                    Fig. 5.                    Different curve shapes of membership functions.
current number of iterations. From the above equation, we
can appreciate that the value of w will decrease as the iteration                                        TABLE II
number increases. In the search process, the most efficient way                      E XAMPLE SOLUTIONS FOR DIFFERENT DESIGN SCENARIOS .
to locate the optimal or near-optimal solutions in a complex
large search space is first to move to the smaller solution                          Generators/objectives            Pessimistic        Linear         Optimistic
                                                                                    PG1                                0.0678          0.0500            0.0605
space as promptly as possible, and then seek out the desired                        PG2                                0.2463          0.2430            0.2425
solution in this space via thorough search. The parameter w is                      PG3                                0.3863          0.4029            0.4221
defined to regulate the size of search step of each particle.                        PG4                                0.9164          0.9300            0.9354
                                                                                    PG5                                0.4034          0.3990            0.3490
At first, the value of w is set relatively large in order to                         PG6                                0.3166          0.2929            0.2897
drive the particle to the solution area quickly. Then, when                         W                                  0.5043          0.5232            0.5408
the particle approaches the desired solution, the size of each                      Cost ($/hour)                     518.893          515.436          512.147
                                                                                    Risk level                         6.5864          6.49894          6.31094
search step becomes smaller in order to prevent the particle
from flying past the target position during the flight. In this
way, the desired solutions can be sought out through gradual
refinement. The simulation program is coded using C++ and            fronts evolved using the proposed MOPSO are shown in
executed in a 2.20 GHz Pentium-4 processor with 512 MB of           Figure 6.
RAM. In simulations, the minimum and maximum allowable                 As shown in the figure, the Pareto-optimal solutions are
wind power penetrations are set as 10% and 20% of the total         widely distributed on the tradeoff surface due to the diver-
load demand, respectively. The running cost of wind power           sity preserving mechanisms used in the proposed MOPSO
is calculated based on its linear relationship with the amount      algorithm. Unlike the single-objective optimization, in multi-
of wind power integrated, i.e., W C = σW . The coefficient           objective optimization the decision maker can choose a suit-
σ indicating the running cost of wind power is set 50 $/p.u.        able solution based on his/her preference from a pool of non-
in the simulation. The parameters used in the simulations are       dominated solutions. We can also appreciate that for the same
listed below and different function curves are shown in Figure      risk level calculated from different membership functions,
5.                                                                  the optimistic design has the lowest operational cost since it
                                                                    includes the largest amount of wind power among all of the
   • Quadratic representation (optimistic design): aw          =
                                                                    three designs.
      −9.9607, bw = 4.94, cw = 0.4;
   • Linear representation (neutral design): Wmin =0.2834,
      Wmax =0.5668;                                                 B. Sensitivity analysis
   • Quadratic representation (pessimistic design): aw =               A study on sensitivity analysis is carried out in order to
      4.9803, bw = −7.7629, cw = 2.8.                               illustrate the impacts of different allowable ranges of wind
   The illustrative non-dominated solutions derived in different    power penetration as well as different running costs of wind
design scenarios are listed in Table II and the Pareto-optimal      power on the final non-dominated solutions derived. Here the
               20                                                                                         20
                                                      Quadratic function (Pessimistic)                                                                                       [6%,16%]
               18                                     Linear function (Neutral)                           18                                                                 [8%,18%]
                                                      Quadratic function (Optimistic)                                                                                        [10%,20%]
               16                                                                                         16

               14                                                                                         14

               12                                                                                         12
  Risk level

                                                                                             Risk level
               10                                                                                         10

                8                                                                                          8

                6                                                                                          6

                4                                                                                          4

                2                                                                                          2

                0                                                                                          0
                500        510        520        530          540        550         560                   500     510     520     530   540      550      560         570     580       590
                                       Operational cost ($/hour)                                                                    Operational cost ($/hour)

Fig. 6.             Pareto fronts obtained based on different membership functions.        Fig. 7.             Pareto fronts obtained for different wind penetration ranges.

linear membership function is used. As indicated in (II.1), for                                           20
                                                                                                                                                             Running cost = 40 $/p.u.
the same range of wind penetration, we can see that the more                                              18                                                 Running cost = 50 $/p.u.
wind power is integrated, the more insecure the system will                                                                                                  Running cost = 60 $/p.u.
become. We herein quantify the impact of wind penetration
through numerical simulations by changing the permissible                                                 14

ranges of wind power penetration [W (PD )min , W (PD )max ].                                              12
                                                                                             Risk level

In the simulations for determining the impact of different
allowable wind penetration ranges, the running cost of wind
power keeps unchanged, i.e., σ = 50$/p.u. In a similar                                                     8
fashion, in the simulations for examining the impact of running                                            6
costs of wind power, the penetration range of wind power is
fixed, i.e., [W (PD )min , W (PD )max ] = [10%∗PD , 20%∗PD ].                                               4

The derived Pareto fronts are shown in Figure 7 and Figure 8,                                              2
respectively. From the figures, we can appreciate that the non-
dominated solutions vary with the different ranges of allowable                                            500       510         520       530        540        550         560         570
                                                                                                                                       Operational cost ($/hour)
wind penetration as well as the running costs of wind power.
In Figure 7, at the same risk level, the design scenario with the
largest value of maximum allowable wind penetration Wmax                                   Fig. 8.             Pareto fronts obtained for different running costs of wind power.
has the lowest cost since the most portion of wind power is
integrated. In Figure 8, at the identical risk level, the scenario
with the lowest running cost of wind power results in the                                  wind power penetration level and operational costs are used
lowest overall cost since the same amount of wind power is                                 in the construction of economic dispatch models. A multi-
integrated with the lowest cost.                                                           objective particle swarm optimization (MOPSO) algorithm
                                                                                           is developed to derive the optimal tradeoff solutions with
                             VII. C ONCLUDING R EMARKS                                     respect to the two specified design objectives. Different design
   The paper is primarily intended to investigate the integration                          scenarios can be formulated according to dispatcher’s attitudes
of wind power into conventional power networks and its                                     toward wind power integration in terms of risk and cost. A nu-
impact on the generation resource management due to its                                    merical application example is used to demonstrate the validity
nondispatchable characteristic. Wind power is environmentally                              and applicability of the proposed optimization procedure. In
attractive since it can help to spur the reductions in fossil fuel                         the further investigations, probabilistic methods may also be
and natural gas consumption. Wind power needs less opera-                                  adopted to handle various operational and planning problems
tional cost since it does not consume fossil fuels and natural                             in a more effective fashion due to the intermittency of wind
gases. However, due to the intermittent and variable nature                                power.
of the wind power, it is usually quite difficult to determine
how much wind power should be used in order to guarantee                                                                         ACKNOWLEDGMENT
both power system security and operational cost reduction. In                                The authors would like to thank the financial support from
this paper, fuzzy representations of system security in terms of                           NSF for this research (Grant No. ECS0406794: Exploring the
Future of Distributed Generation and Micro-Grid Networks).

                              R EFERENCES
[1] Ahlstrom, M., Jones, L., Zavadil, R., and Grant, W. (2005). The future of
    wind forecasting and utility operations, IEEE Power & Energy Magazine,
    November/December, pp. 57–64.
[2] Bagen and Billinton, R. (2004). Evaluation of different operating strate-
    gies in small stand-alone power systems. IEEE Transactions on Power
    Systems, Vol. 20, No. 3, pp. 654–660.
[3] Billinton, R. and Bai, G. (2004). Generating capacity adequacy associated
    with wind energy. IEEE Transactions on Power Systems, Vol. 19, No. 3,
    pp. 641–646.
[4] DeMeo, E. A., Grant, W., Milligan, M. R., and Schuerger, M. J. (2005).
    Wind plant integration: costs, status, and issues, IEEE Power & Energy
    Magazine, November/December, pp. 38–46.
[5] Douglas, J. (2006). Putting wind on the grid, EPRI Journal, Spring 2006,
    pp. 6–15.
[6] Eriksen, P. B., Ackermann, T., Abildgaard, H., Smith, P., Winter, W., and
    Garcia, R. (2005). System operation with high wind penetration, IEEE
    Power & Energy Magazine, November/December, pp. 65–74.
[7] Farag, A., Al-Baiyat, S., and Cheng, T. C. (1995). Economic load dis-
    patch multiobjective optimization procedures using linear programming
    techniques, IEEE Transactions on Power Systems, Vol. 10, pp. 731–738.
[8] Gaing, Z.-L. (2003). Particle swarm optimization to solving the economic
    dispatch considering the generator constraints, IEEE Transactions on
    Power Systems, Vol. 18, No. 3, pp. 1187–1195.
[9] Hota, P. K. and Dash, S. K. (2004). Multiobjective generation dispatch
    through a neuro-fuzzy technique, Electric Power Components and Sys-
    tems, 32: 1191–1206.
[10] Jayabarathi, T., Jayabarathi, K., Jeyakumar, D. N., et al., (2005).
    Evolutionary programming techniques for different kinds of economic
    dispatch problems, Electric Power System Research 73, pp. 169–176.
[11] Jenkins, N., et al. (2000). Embedded Generation, IEE Power and Energy
    Series 31, The Institution of Electrical Engineers.
[12] Kennedy, J. and Eberhart, R. (1995). Particle swarm optimization, IEEE
    Proceedings of the International Conference on Neural Networks, Perth,
    Australia, pp. 1942–1948.
[13] Kennedy, J. and Eberhart, R. (2001). Swarm Intelligence, Morgan
    Kaufmann Publishers, San Francisco.
[14] Lee, K.Y., Yome, A.S., and Park, J.H. (1998). Adaptive Hopfield neural
    networks for economic load dispatch, IEEE Transactions on Power
    Systems, Vol. 13, No. 2, pp. 519-526.
[15] Lin, C.E. and Viviani, G.L. (1984). Hierarchical economic dispatch
    for piecewise quadratic cost functions, IEEE Transactions on Power
    Apparatus and Systems, Vol. 103, No. 6, pp. 1170-1175.
[16] Miranda, V. and Hang, P. S. (2005). Economic dispatch model with fuzzy
    constraints and attitudes of dispatchers, IEEE Transactions on Power
    Systems, Vol. 20, No. 4, Nov., pp. 2143–2145.
[17] Piwko, R., Osborn, D., Gramlich, R., Jordan, G., Hawkins, D., and
    Porter, K. (2005). Wind energy delivery issues, IEEE Power & Energy
    Magazine, November/December, pp. 67–56.
[18] Smith, J. C., (2005). Wind of change: issues in utility wind integration,
    IEEE Power & Energy Magazine, November/December, pp. 20–25.
[19] Wang, L. F. and Singh, C. Multi-objective stochastic power dispatch
    through a modified particle swarm optimization algorithm, Special Ses-
    sion on Applications of Swarm Intelligence to Power Systems, Proceed-
    ings of IEEE Swarm Intelligence Symposium, Indianapolis, May, 2006,
    pp. 127–135.
[20] Zavadil, R., Miller, N., Ellis, A., and Muljadi, E. (2005). Making
    connections: wind generation challenges and progress, IEEE Power &
    Energy Magazine, November/December, pp. 26–37.