By the end of 2013, the installed wind capacity in China has been the largest in the world, and the wind has become the thirdlargest electrical power source in China [1]. On the other hand, limited predictability and controllability of wind generation due to its intermittent characteristics still affect the quality of its power supply [2]. To enhance the reliability of power system operation under fluctuant load and wind power, spinning reserve (SR) provided by thermal generators needs to be scheduled. As the penetration of wind generation increases rapidly, it is of great importance to explore more reliable and economic strategies to address wind power uncertainty in reserve scheduling problem [3].
The reserve scheduling with wind generation integration has been well studied in many literatures, in which those methodologies can be classified into three categories: the deterministic reserve margin methods (DRM), the stochastic optimization (SO) and the robust optimization (RO). As of DRM, the reserve capability is usually evaluated based on the possible largest forecast error of realtime wind power output [4], [5]. However, it is very difficult to determine the representative operation condition for DRMbased decisionmaking. As an improvement, SO can consider more complicated operating conditions than DRM with certain assumption of wind power probability distribution given. For instance, in [6] a Beta function is used to fit the distribution of realtime output of wind farms; in [7] a chanceconstrained optimization model is proposed to compute the reserve requirement; Weibull distribution of wind speed is used to obtain the probability distribution of wind generation output in [8]. However, there exist two SO's disadvantages: 1) SO problem above usually is a nonconvex optimization problem which poses optimal solution difficulty [7]; 2) accurate evaluation of probability distribution is hard in some real situations [9].
From another aspect, in recent years the RObased approaches have attracted considerable attention for their probability distributionfree advantages [10][25]. As a branch of RO methods, the robust linear optimization (RLO) is a convex optimization and can well control the degree of solution's conservatism [10][15]. The method of RLO has experienced several evolutions since Soyster first developed the RO method in 1973 [12]. To address the RO's overconservatism problem, BenTal and Nemirovski introduced a robust linear optimization, where the rows of the constraint matrix belong to some ellipsoid uncertainty sets and the robust counterpart is a secondorder cone programming [13]. To overcome the numerical calculation difficulty, an improved RLO method was proposed in [14], in which a linear robust counterpart is developed. However, all above RLO methods considering distributional information assume that the uncertainty of random variables is symmetrical with respect to the mean value. To handle the uncertainty asymmetry, Kang further proposed an RLO method called Kang's robust linear optimization (KRLO) [15].
To our best knowledge, the existing RLObased reserve scheduling methods have not included the research on asymmetrical wind power with respect to its mean value, e.g., [26] uses wind power interval without mean value; [27] considers symmetrical wind power distributional information, which indicates the upspinning reserves and the downspinning reserves may have the same distribution interval. However, wind power may be asymmetrical distributed in some cases [6][9]. In this paper, a robust reserve scheduling method is introduced to pursue the minimum spinning reserve cost, in which spinning reserve allocation as an additional AGC (automatic generation control) functionality is achieved by a linearization decision rule. This KRLObased method applies to asymmetrical wind power distribution, and can control the degree of solution's conservatism by an adjustable robustness budget. Meanwhile, the wind power uncertainty is represented by an asymmetrical uncertainty set of KRLO with given mean value and maximum/minimum limits. The robust counterpart handling of KRLO, which is based on dual theory, can convert the proposed model into deterministic formulations.
The rest of the paper is organized as follows. In Section Ⅱ, the mathematical formulation of reserve scheduling is introduced. The KRLO method and the solution framework are presented in Section Ⅲ. Numerical examples are provided in Section Ⅳ. And Section Ⅴ concludes the paper.
Ⅱ. MATHEMATICAL FORMULATION A. Problem DescriptionThe following assumptions are made for problem formulation: 1) only wind power uncertainty is considered and bus loads are regarded as wellforecasted; 2) each wind generator is connected to the transmission network from a single bus; 3) wind} power can be asymmetrically distributed; 4) a standard DC power flow is used [28]; 5) unit commitment solution has been solved in advance. The first and second assumptions are common in the reserve scheduling problem. The third assumption is more practical compared with those symmetricaldistribution cases. The fourth is based on a fact that in a transmission network the active power is mainly phase angle related and reactive power is voltage magnitude related. Under normal condition, voltage magnitude variation is small and a DC power flow can be considered in transmission network [28]. The last assumption ensures that the state of units has been determined before reserve scheduling.
Spinning reserve serving to AGC is allocated to each participating generator according to some distribution vector,
Download:


Fig. 1 Schematic diagram illustrating the AGC functionality required for the proposed reserve scheduling method. 
The reserve scheduling aims to reduce the procurement cost of generation and reserves. The following linear cost objective function is adopted [26],
$ \begin{equation} \min\quad {\pmb C}_{\rm g}^T {\pmb P}_{\rm G} + {\pmb C}_{\rm up}^T {\pmb R}_{\rm up}+ {\pmb C}_{\rm down}^T {\pmb R}_{\rm down} \label{eq1} \end{equation} $  (1) 
where
Constraints are as follows:
1). Nodal power balance:$ \begin{equation} {\pmb B}{\pmb \theta} + {\pmb g}_{\rm s} + ({\pmb g}_{\rm f} + \Delta {\pmb g}_{\rm f} )+ {\pmb p}_{\rm w} = {\pmb p}_{\rm d} \label{eq2} \end{equation} $  (2) 
where
$ \begin{equation} {\pmb P}_{\rm w}^{\rm l} \le {\pmb P}_{\rm w}\le {\pmb P}_{\rm w}^{\rm u} \label{eq3} \end{equation} $  (3) 
where
Responding to the time scale and practical demand in engineering, the random variable ${\pmb P}_{\rm w}$ can be divided into two parts: the forecasted mean value and the uncertain deviation:
$ \begin{equation} {\pmb P}_{\rm w}=\pmb{\mu}_{\rm w} + \Delta{\pmb p}_{\rm w} \label{eq4} \end{equation} $  (4) 
where
$ \begin{equation} \left\{ \begin{aligned} \pmb{\omega}^{\rm F} = {\pmb P}_{\rm w}^{\rm u}  \pmb{\mu}_{\rm w}\\ \pmb{\omega}^{\rm B} = \pmb{\mu}_{\rm w}{\pmb P}_{\rm w}^{\rm l}\\ \end{aligned} \right. \label{eq5} \end{equation} $  (5) 
where
$ {\pmb g}_{\rm f}^{\min}\le {\pmb g}_{\rm f} + \Delta {\pmb g}_{\rm f}\le {\pmb g}_{\rm f}^{\max} $  (6a) 
$ {\pmb g}_{\rm s}^{\min}\le {\pmb g}_{\rm s} \le {\pmb g}_{\rm s}^{\max} $  (6b) 
where
$ \Delta{\pmb g}_{\rm f} = {\pmb d}({\pmb l}^T \Delta {\pmb p}_{\rm w}) $  (7a) 
$ {\pmb R}_{\rm down}\le \Delta {\pmb g}_{\rm f}\le {\pmb R}_{\rm up} $  (7b) 
$ {\pmb R}_{\rm down}, {\pmb R}_{\rm up} \le \eta({\pmb g}_{\rm f}^{\max}  {\pmb g}_{\rm f}^{\min}) $  (7c) 
$ {\pmb 1}^T {\pmb d} = 1 $  (7d) 
where
$ {\pmb F} = {\pmb T}{\pmb \theta} $  (8) 
$ {\pmb f}_{\max} \le {\pmb F} \le {\pmb f}_{\max} $  (9) 
where
In the model, decision variables are
Considering SR allocation and transmission capability constraints, the reserve scheduling with uncertain wind power, as an additional AGC functionality [30], is provided. The operating controller may monitor the deviation of the wind power and use the distribution vector
The reserve scheduling problem in (1)(9) is a special linear programming problem with wind power uncertainty embedded, which will be handled by the tractable reformulation below.
Ⅲ. TRACTABLE REFORMULATION A. KRLO With Asymmetric Data UncertaintyFor a general linear optimization model considering uncertain data:
$ \left\{\begin{aligned} &\min \quad{\pmb cx}\\ &{\rm s.t.} \quad {A}{\pmb x} \le {\pmb b}\\ &\qquad ~~{\pmb l} \le {\pmb x} \le {\pmb u} \end{aligned} \right. $  (10) 
where
RLO was proposed to find a solution of (10) immune to the uncertain data in [14]. RLO methods are based on the following notation. The elements of
Studies in [14] and [26] are both based on the assumption that
To handle the asymmetrical uncertainty where
$ \begin{equation} \Re(\Gamma_i) = \left\{ {\pmb a_i} \left \begin{aligned} &a_{ik} \in [\bar{a}_{ik}\beta_{ik}t_{ik}^{\rm B}, \bar{a}_{ik} + \beta_{ik}t_{ik}^{\rm F}]\\ &0\le \beta_{ik}\le 1; \sum\limits_{k\in J_i}\beta_k\le \Gamma_i; \forall k\in J_i\\ \end{aligned} \right. \right\} \label{eq11} \end{equation} $  (11) 
where
With dual theory, the robust counterpart of (10) and (11) can be given in (12):
$ \begin{equation} \left\{ \begin{aligned} &\min \quad{\pmb {cx}}\\ &{\rm s.t.} \sum\limits_{j=1}^n \bar{a}_{ij}x_j + \Gamma_i z_i + \sum\limits_{k\in J_i} p_{ik} \le b_i, ~ i=q, \ldots, m\\ &z_i + p_{ik}\ge \max (t_{ik}^{\rm F} x_k, t_{ik}^{\rm B} x_k);~ i=1, \ldots, m;~ \forall k\in J_i\\ &z_i \ge 0;~ p_{ik} \ge 0;~ i=1, \ldots, m;~ \forall k\in J_i\\ &{\pmb l} \le {\pmb x}\le {\pmb u}\\ \end{aligned} \right. \label{eq12} \end{equation} $  (12) 
where
By introducing the robust counterpart, the nominal linear programming is converted into a deterministic formulation. The robust counterpart (12) has the same optimal solution with the nominal linear programming (10) when the solution is immune to the uncertainty of (11), which is proved in [15].
B. Robust Reserve SchedulingThis section serves mainly to achieve the decoupling between the decision variables and uncertain variables in the equality constraints (2), and convert the optimization model into the general form (10).
1). Variables SimplificationNodal phase angle vector
$ \begin{equation} { S\pmb{F}} + {\pmb g}_{\rm s} + ({\pmb g}_{\rm f} + \Delta {\pmb g}_{\rm f}) + {\pmb p}_{\rm w} = {\pmb p}_{\rm d} \label{eq13} \end{equation} $  (13) 
where
The power balance constraint covered in equality (13) can be formulated as
$ \begin{equation} \sum\limits_{i\in {IB}} {\pmb g}_{\rm si} + \sum\limits_{u\in {UB}}({\pmb g}_{\rm fu} + \Delta {\pmb g}_{\rm fu}) + \sum\limits_{j\in {WB}}({\pmb \mu}_{\rm wj} + \Delta {\pmb p}_{\rm wj}) = \sum\limits_{n\in { LB}} {\pmb p}_{\rm dn} \label{eq14} \end{equation} $  (14) 
where IB, UB, WB and LB respectively denote the sets of outputfixed unit buses, spinning reserve unit buses, wind generator buses and load buses, respectively.
Further, by substituting equality constraints (13) to the transmission limit constraints (9), inequality constraints with uncertain wind power are obtained:
$ \begin{equation} {\pmb f}_{\max}\le {\pmb S}^{1} [{\pmb p}_{\rm d}{\pmb g}_{\rm s}({\pmb g}_{\rm f} + \Delta{\pmb g}_{\rm f}){\pmb p}_{\rm w}]\le {\pmb f}_{\max}. \label{eq15} \end{equation} $  (15) 
So far, the equality constraint with uncertain variables (2) is converted into the equality (14) and inequality (15) without nodal phase angle vector
Notice that there are no equalities with uncertain data in (12) and there are uncertain variables in equality (14). Therefore, to eliminate the uncertain wind power in (14), we substitute constraints (4), (7a) and (7d) to (14). And then, the equality with respect to the system power balance under wind power forecast is obtained:
$ \begin{equation} \sum\limits_{i\in {IB}}\! {\pmb g}_{\rm si} \!+\! \sum\limits_{u\in {UB}}\!{\pmb g}_{\rm fu} \!+\! \sum\limits_{j\in {WB}}\!{\pmb \mu}_{\rm wj} \!= \!\sum\limits_{n\in {LB}}\!{\pmb p}_{\rm dn}. \label{eq16} \end{equation} $  (16) 
By integrating constraints (6a)(6b), (7a)(7d) and (15), we get the following inequality constraints with the only uncertain vector
$ \left\{ \begin{aligned} {\pmb S}^{1}[{\pmb p}_{\rm d}{\pmb g}_{\rm s}{\pmb \mu}_{\rm w}{\pmb d{\rm l}^T}\Delta{\pmb p}_{\rm w}]&\le {\pmb f}_{\max}\\ {\pmb S}^{1}[{\pmb p}_{\rm d}{\pmb g}_{\rm s}{\pmb \mu}_{\rm w} {\pmb d{\rm l}^T}\Delta{\pmb p}_{\rm w}]&\le {\pmb f}_{\max}\\ \end{aligned} \right. $  (17) 
$ \left\{ \begin{aligned} {\pmb g}_{\rm f}  {\pmb d}({\rm l}^T\Delta{\pmb p}_{\rm w})&\le {\pmb g}_{\rm f}^{\max}\\ {\pmb g}_{\rm f}  {\pmb d}({\rm l}^T\Delta{\pmb p}_{\rm w})&\le {\pmb g}_{\rm f}^{\min}\\ \end{aligned} \right. $  (18) 
$ \left\{ \begin{aligned} {\pmb d}({\rm l}^T\Delta{\pmb p})&\le {\pmb R}_{\rm up}\\ {\pmb d}({\rm l}^T\Delta{\pmb p})&\le {\pmb R}_{\rm down}.\\ \end{aligned} \right. $  (19) 
The uncertainty set of deviation
$ \begin{equation} \Re (\Gamma) = \left\{ \Delta {\pmb p}_{\rm w}\left \begin{aligned} &\Delta p_{wi}\in [\beta_i w_i^{\rm B}, \beta_i w_i^{\rm F}]\\ &0\le \beta_i \le 1; \sum\limits_{i=1}^{WB}\beta_i \le \Gamma\\ \end{aligned} \right. \right\} \label{eq20} \end{equation} $  (20) 
where the robustness budget
Thus, the reserve scheduling model (1)(9) can be further converted into a deterministic programming problem that consists of linear equality/inequality constraints.
C. Algorithm Complexity AnalysisThe optimization problem can be transformed into the robust counterpart (12) automatically via the MATLAB interface YALMIP [30]. The algorithm for the problem could be any common linear programming algorithm, which is embedded in many solvers, such as the CPLEX [31]. Therefore, the algorithm in this paper is of polynomial time complexity and fast enough when wind generation changes.
Ⅳ. CASE STUDOESIn this section, a simulation analysis of varying robustness budget and ramp rate will be conducted with asymmetrical wind power distribution. In addition, this section also presents a comparison between KRLO method and an adjustable robust method with symmetrical wind power distribution.
A. Revised Garver's 6bus Test SystemThe revised Garver's 6bus system [10] is used to verify the effectiveness of the introduced approach. As shown in Fig. 2, four wind farms are connected to the power grid from different buses. Parameters of the load (L_{1}L_{5}), generators, wind farms, branches, and the cost coefficients are shown in Tables ⅠⅢ, respectively. All of the thermal generators are engaged in the reserve allocation;
Download:


Fig. 2 Configuration of modified Garver's network. 
The optimization is solved by the solver CPLEX in the interface YALMIP of MATLAB.
1). Effect of Robustness BudgetThis case investigates the effect of the robustness budget. Fixed ramp rate parameter
Download:


Fig. 3 Downspinning reserves provided by each unit ( 
Download:


Fig. 4 Upspinning reserves provided by each unit ( 
As is seen in Table Ⅳ, the total cost of the generation and reserves with a fixed
This case considers the effect of the ramp rate. Fig. 5 and Fig. 6 show the up/down spinning reserves provided by each participating unit with
Download:


Fig. 5 Downspinning reserves provided by each unit ( 
Download:


Fig. 6 Upspinning reserves provided by each unit ( 
It can be seen from the comparison between $\eta = 1/6$ and $\eta = 1/4$ that a bigger value of $\eta$, which corresponds to a larger ramp rate provided by participating units, will reduce the total cost in Table Ⅳ, and increase the reserves provided by unit 1 in Fig. 5 and Fig. 6. The reason is that a larger ramp rate ensures the power mismatch can be compensated more by the optimal reserves purchased from unit 1, which is an affordable approach in this test system.
In addition, the total up/down spinning reserves with
This case investigates the differences between the KRLO method and the adjustable robust optimization (ARO) method of [26] with symmetric wind power. Based on the method in [26], each deviation of the wind generation belongs to an uncertain interval without robustness budget, i.e.,
By setting the
The IEEE 39bus system is used to test the scalability of the method. Without loss of generality, it is assumed that three wind farms are connected with the test system at bus 16, 23 and 26, respectively. The data of wind farms are listed in Table Ⅵ. The generators at bus 30, 31, 35 and 38 are assumed to provide SR and the power outputs of the other units are fixed. All other system data can be referred in [32].
The total generation and reserve cost, total up spinning reserves and down spinning reserves are shown in Table Ⅶ. It can be seen that the cost and reserves increase with parameter
The allocations of up/down spinning reserves are shown in Fig. 7 and Fig. 8. In Fig. 7, unit 1 provides the maximum down spinning reserves, and the reserves increase with robustness budget being larger. The same situation occurs in Fig. 8.
Download:


Fig. 7 Downspinning reserves provided by each unit ( 
Download:


Fig. 8 Upspinning reserves provided by each unit ( 
By comparison between the revised Garver's 6bus system and IEEE 39bus system, the scalability of the method is illustrated.
Ⅴ. CONCLUSIONSA robust method is proposed in this paper, as an analysis tool for reserve scheduling considering asymmetrical wind power distribution. This method can serve as an additional AGC function to allocate the SR with varying solution's conservatism, and ramp rate. Simulation results verify that the approach is effective.
Besides, storage system will be helpful for the reserve scheduling problem, which will be our future work. On the other hand, this work is based on linear dual theory that has strong duality property, therefore, the robust counterpart is equal to the original problem. If the problem has nonaffine constraints, the robust counterpart is hard to obtain, which will serve as another future work of us.
[1]  China's wind power installed capacity in 2013, China National Renewable Energy Centre, Beijing, 2013. 
[2]  F. Bouffard and F D Galiana, "Stochastic security for operations planning with significant wind power generation, " IEEE Trans. Power Syst. , vol. 23, no. 2, pp. 306316, May 2008. http://ieeexplore.ieee.org/document/4470561/ 
[3]  N. Menemenlis, M. Huneault, and A. Robitaille, "Computation of dynamic operating balancing reserve for wind power integration for the timehorizon 148 hours, " IEEE Trans. Sust. Energy, vol. 3, no. 4, pp. 692702, Oct. 2012. http://ieeexplore.ieee.org/document/6202398/ 
[4]  J. M. Morales, A. J. Conejo, and J. PerezRuiz, "Economic valuation of reserves in power systems with high penetration of wind power, " IEEE Trans. Power Syst. , vol. 24, no. 2, pp. 900910, May 2009. http://ieeexplore.ieee.org/document/4814478/ 
[5]  M. A. OrtegaVazquez and D. S. Kirschen, "Estimating the spinning reserve requirements in systems with significant wind power generation penetration, " IEEE Trans. Power Syst. , vol. 24, no. 1, pp. 114124, Feb. 2009. 
[6]  S. Tewari, C. J. Geyer, and N. Mohan, "A statistical model for wind power forecast error and its application to the estimation of penalties in liberalized markets, " IEEE Trans. Power Syst. , vol. 26, no. 4, pp. 20312039, Nov. 2011. http://ieeexplore.ieee.org/document/5765544/ 
[7]  Y. C. Guan, Y. S. Wang, and Z. Y. Tan, "Environmental protection and security considered dynamic economic dispatch for wind farm integrated systems, " in Proc. 2012 AsiaPacific Power and Energy Engineering Conf. , Shanghai, China Proc., pp. 14. 
[8]  J. Hetzer, D. C. Yu, and K. Bhattarai, "An economic dispatch model incorporating wind power, " IEEE Trans. Energy Convers. , vol. 23, no. 2, pp. 603611, Jun. 2008. http://ieeexplore.ieee.org/document/4505391/ 
[9]  W. C. Meng, Q. M. Yang, and Y. X. Sun, "Guaranteed performance control of DFIG variablespeed wind turbines, " IEEE Trans. Control Syst. Technol. , vol. 24, no. 6, pp. 22152223, Nov. 2016. http://ieeexplore.ieee.org/document/7419243/ 
[10]  Y. Chen, J. Y. Wen, and S. J. Cheng, "Minimum loadcurtailment in transmission network planning considering integrated wind farms, " Proc. CSEE, vol. 31, no. 34, pp. 2027, Dec. 2011. 
[11]  X. N. Han, J. M. Li, J. Y. Wen, X. M. Ai, J. H. Li, and W. H. Luo, "Optimization for robust energy storage allocation in power system with multiple wind farms integrated, " Proc. CSEE, vol. 35, no. 9, pp. 21202127, May 2015. 
[12]  A. L. Soyster, "Convex programming with setinclusive constraints and applications to inexact linear programming, " Operat. Res. , vol. 21, no. 5, pp. 11541157, Oct. 1973. 
[13]  A. BenTal and A. Nemirovski, "Robust solutions of linear programming problems contaminated with uncertain data, " Mathem. Progr. , vol. 88, no. 3, pp. 411424, Sep. 2000. 
[14]  D. Bertsimasl and M. Sim, "The price of robustness, " Operat. Res. , vol. 52, no. 1, pp. 3553, Feb. 2004. 
[15]  S. C. Kang, "Robust linear optimization using distributional information, " Ph. D. dissertation, Boston Univ., Boston, USA, 2008. 
[16]  A. BenTal, A. Goryashko, E. Guslitzer, and A. Nemirovski, "Adjustable robust solutions of uncertain linear programs, " Mathem. Progr. , vol. 99, no. 2, pp. 351376, Mar. 2004. http://www.sciencedirect.com/science/article/pii/S0167637799000164 
[17]  R. W. Jiang, J. H. Wang, and Y. P. Guan, "Robust unit commitment with wind power and pumped storage hydro, " IEEE Trans. Power Syst. , vol. 27, no. 2, pp. 800810, May 2012. http://ieeexplore.ieee.org/document/6069602/ 
[18]  H. X. Ye and Z. Y. Li, "Robust securityconstrained unit commitment and dispatch with recourse cost requirement, " IEEE Trans. Power Syst. , vol. 31, no. 5, pp. 35273536, Sep. 2016. http://ieeexplore.ieee.org/document/7331341/ 
[19]  Y. An and B. Zeng, "Exploring the modeling capacity of twostage robust optimization: variants of robust unit commitment model, " IEEE Trans. Power Syst. , vol. 30, no. 1, pp. 109122, Jan. 2015. http://ieeexplore.ieee.org/document/6812211/ 
[20]  Q. M. Yang, S. Jagannathan, and Y. X. Sun, "Robust integral of neural network and error sign control of MIMO nonlinear systems, " IEEE Trans. Neural Netw. Learn. Syst. , vol. 26, no. 12, pp. 32783286, Dec. 2015. http://ieeexplore.ieee.org/document/7234927/ 
[21]  E. H. Delage, "Distributionally robust optimization in context of datadriven problems, " Ph. D. dissertation, Stanford Univ., Stanford, California, USA, 2009. 
[22]  R. Jiang, M. Zhang, G. Li, and Y. Guan, "Benders decomposition for the twostage security constrained robust unit commitment problem, " Ph. D. dissertation, Univ. of Florida, USA, 2011. 
[23]  W. C. Meng, Q. M. Yang, Y. Ying, Y. Sun, Z. Y. Yang, and Y. X. Sun, "Adaptive power capture control of variablespeed wind energy conversion systems with guaranteed transient and steadystate performance, " IEEE Trans. Energy Convers. , vol. 28, no. 3, pp. 716725, Sep. 2013. http://ieeexplore.ieee.org/document/6570509/ 
[24]  B. P. G. V. Parys, D. Kuhn, P. J. Goulart, and M. Morari, "Distributionally robust control of constrained stochastic systems". IEEE Trans. Autom. Control , vol.61, no.2, 2015. 
[25]  S. Zymler, D. Kuhn, and B. Rustem, "Worstcase value at risk of nonlinear portfolios, " Management Science, vol. 59, no. 1, pp. 172188, Jan. 2013. 
[26]  Z. G. Li, W. C. Wu, B. M. Zhang, and B. Wang, "Adjustable robust realtime power dispatch with largescale wind power integration, " IEEE Trans. Sust. Energy, vol. 6, no. 2, pp. 357368, Apr. 2015. http://ieeexplore.ieee.org/document/7001653/ 
[27]  W. Wei, F. Liu, S. W. Mei, and Y. H. Hou, "Robust energy and reserve dispatch under variable renewable generation, " IEEE Trans. Smart Grid, vol. 6, no. 1, pp. 369380, Jan. 2015. http://ieeexplore.ieee.org/document/6862932/ 
[28]  G. Andersson, "Modeling and Analysis of Electric Power Systems", Zurich, Switzerland: ETH Zurich, 2008, Lecture Notes. 
[29]  M. Vrakopoulou, K. Margellos, J. Lygeros, and G. Andersson, "A probabilistic framework for reserve scheduling and N1 security assessment of systems with high wind power penetration, " IEEE Trans. Power Syst. , vol. 28, no. 4, pp. 38853896, Nov. 2013. http://ieeexplore.ieee.org/document/6570751/ 
[30]  J. Lofberg, "YALMIP: A toolbox for modeling and optimization in MATLAB, " in Proc. 2004 IEEE Int. Symp. Computer Aided Control Systems Design, Taipei, China, pp. 284289. 
[31]  CPLEX11. 0 Users Manual, ILOG. SA., Gentilly, France, 2008. 
[32]  M. A. Pai, Energy Function Analysis for Power System Stability. US: Springer, 1989. 