# Practical Robust Optimization Method for Unit Commitment of a System with Integrated Wind Resource

Authors: Yuanchao Yang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Unit commitment, one of the significant tasks in power system operations, faces new challenges as the system uncertainty increases dramatically due to the integration of time-varying resources, such as wind. To address these challenges, we propose the formulation and solution of a generalized unit commitment problem for a system with integrated wind resources. Given the prespecified interval information acquired from real central wind forecasting system for uncertainty representation of nodal wind injections with their correlation information, the proposed unit commitment problem solution is computationally tractable and robust against all uncertain wind power injection realizations. We provide a solution approach to tackle this problem with complex mathematical basics and illustrate the capabilities of the proposed mixed integer solution approach on the large-scale power system of the Northwest China Grid. The numerical results demonstrate that the approach is realistic and not overly conservative in terms of the resulting dispatch cost outcomes.


---

## Body

## 1. Introduction

Unit commitment (UC) is a critically important function in operation scheduling, as it provides the linkage from the maintenance and production scheduling to economic dispatch [1]. In today’s competitive electricity market environment, UC is a basic tool used by an independent system operator (ISO) or regional transmission organization (RTO) to clear the day-ahead markets and by a resource manager to optimize its offering strategy [2]. The UC solution provides the ISO with the optimal hour-by-hour schedule for start-up and shut-down of generation units taking into account various constraints, including system security under critical component outages and load uncertainty and important physical limitations and obligations to meet the base-load and cycling demands.The conventional grid is subject to a wide range of uncertainties in the loads, component outages, and market behaviors. The system operator uses the reserve requirements to specify the capacity to manage the impacts of all the sources of uncertainties [3]. Since the unexpected events may be predicted with high confidence, regional reliability councils or system operators can make certain rules allocating reserve to various units, typically a given percentage of forecasted peak demand, to trade-off between reliability and system economy.The significant penetration of wind generation in the grid raises a complex challenge to system operators caused by the highly uncertain and variable wind power generation pattern which is difficult to predict. The conventional “reserve requirements” can no longer cope with system uncertainties under the new environment. Furthermore, the spatial correlations between wind farms governed by atmospheric conditions complicate the difficulty on system operation and security; thus it is important to quantify the uncertainties and correlations of the wind power outputs at the geographically distributed wind sites.Many research efforts have been made to deal with wind power uncertainty in UC [4–7]. Much of the work has focused on the solution of the optimization problem with the explicit representation of uncertainty. There are two principal approaches in the literature: stochastic optimization and robust mathematical programming. The first approach explicitly makes use of probability distributions to represent various sources of the uncertainties and solves the problem in terms of expectation with discrete scenarios of future uncertain events [4–6]. This approach is simple and effective but has to address the challenges including determining the probabilities of the scenarios that can adequately capture a broad range of uncertainties and developing tractable numerical methods for dealing with huge number of scenarios with application to large-scale systems.The robust optimization techniques recently gained attentions for optimization problems under parameter uncertainties [8, 9]. There are also some technical publications introducing the robust optimization methods into the realm for power system operation [7, 9–11]. Le Cadre et al. [10] presents a robust approach for the single-bus contingency constrained UC problem. A two-stage adaptive robust optimization model is proposed for unit commitment problem with load and wind forecasting uncertainties [7, 9]. Since the wind power uncertainty is rather difficult to characterize in terms of analytic probability distributions, a common way to describe the forecasted wind power output is to give its intervals an associated confidence. Furthermore, Benders’ decomposition method is applied to solve the min-max optimization problem to achieve feasibility in the “worst case” with a nonlinear cost function. Wang et al. [11] propose a Linear Decision Rules approach to solve unit commitment problem with parameter uncertainty, but no known numerical results on the performance of UC based on this approach are provided. In general, the robust optimization approach can provide a more reasonable way to quantify the uncertain information in power systems, especially the renewable resources, such as the wind; it also requires less information than the scenario tree approach and thus can give a schedule results for system operators in a more robust and reliable manner. One of the main challenges of this approach currently is that it needs major computational effort for the decomposition and coordination process to converge the above two-stage robust model. Furthermore, in the above publications, the spatial correlation between nodal injections caused by distributed wind sites, which would affect nodal wind power injections and thus real power flow through transmission lines, is not explicitly considered in the UC problem formulation.To effectively address the above challenges, we propose a new formulation and method for solving the UC problem for a system integrated with wind resource based on robust optimization schemes. The wind uncertainty is modeled as an interval with an associated confidence depending on the individual risk tolerance [12–14]. However, it does not mean the wind farm power output can be independently selected within the above intervals; in other words, for example, if wind farms 1 and 2 have a positive correlation between their power output, then wind farm 2 will have a very low probability to create a small power output under the fact that wind farm 1 has a large power output. To consider the grid line power flow constraints, we explicitly represent the correlations between the uncertain wind power outputs based on the principal component analysis (PCA) techniques. The PCA techniques can capture effectively the correlation relationship between the interrelated random variables and its based transformation allows us to convert a large number of interrelated variables into uncorrelated principal components (PCs) [15, 16]. In this way, the multivariate statistical wind power interdependency can be represented by a series of uncorrelated variables, PCs; this will make it easier to devise a solution approach for the extended UC problem. If PCA-based transformation is not used, the data we collected from the wind power forecasting cannot be directly used in the day-ahead scheduling UC problem since only the data of intervals does not completely represent the information for the wind power output. Only after transforming the correlated wind power outputs by their corresponding uncorrelated principal components, we can use the proposed solution approach to obtain the final scheduling results. To conclude, the PCA technique is the key step for preprocessing the data for the correlated wind power outputs. Finally, it is worth noting that the PCA technique can only perform an orthogonal transformation so as to obtain an uncorrelated data. If the distribution considered is multivariate Gaussian, resulting PCs will be independent; otherwise, PCs will be uncorrelated but still dependent. In fact, wind power output is a non-Gaussian stochastic process and thus the diagonal covariance matrix of PCs only implies that they are uncorrelated. In this paper, we adopt the preprocessing and transformation techniques which have been proposed in the previous publications [17, 18] to obtain approximately Gaussian datasets of wind power output. After that, we can thus treat the PCs as approximately mutually independent.To solve the above two-stage mixed integer nonlinear optimization problem with uncertain parameters, we develop a practical solution approach computationally tractable for the large power systems. We first treat the uncertainty intervals as the worst-case realization and cast the problem into a deterministic formulation. Then, we convert the deterministic counterpart into a single-stage mixed integer linear programming (MILP) form, to avoid the two-stage iteration and convergence issue that is dealt with within decomposition approaches [7, 9]. To reach that goal, an approximate analytical solution of the inner-stage problem is obtained. The analytical solution obtained is aggregated as linear constraints into the deterministic reformulation. The restated single-stage problem has a MILP form and thus can be easily implemented under the software environment of advanced commercial optimization packages which have been widely used in the system operation. Numerical testing results demonstrate that the new approach is computationally tractable and thus could be useful for the operation of real and large-scale power systems.The rest of the paper is organized as follows: we present the formulation of the extended UC problem with wind generation sources at geographically distributed locations in Section2; we discuss the proposed solution approach in Section 3; in Section 4, we illustrate the application of the methodology to a large-scale power system; we conclude the paper in Section 5.


## 2. The UC Problem Formulation with Wind Generation

We devote this section to the development of the explicit representation of the uncertainty in wind resource outputs with the correlations of wind farm outputs at different locations. We consider a central UC problem with 24-hour commitment horizon. We account for the impacts of contingency cases, such as generation and transmission line outage contingencies, in the specification of the reserves requirements.As wind generation penetration deepens, better utilization of probabilistic data is significant and many system operators are involved in collecting probabilistic wind power information [19, 20]. It first collects the measurements of the random wind power output which are used to construct an approximation of the mass density function of the wind output random variables. It identifies the interval containing a subset of the samples of the random variables and calls this the interval of uncertainty. The sum of the probabilities of the samples in the subset determines the associated probability. In this paper, we describe the day-ahead wind power output randomness in terms of an uncertainty interval with an associated confidence probability. Thus, wts, the realization of power output at the wind site s in period t, is within the interval [w_ts,w-ts]; wt, the realization of system wind power output in period t, is within the interval [w_t,w¯t].Furthermore, the correlations between the uncertain wind power outputs among geographically distributed wind farms are described based on PCA techniques [15]. For a given set of observations of power outputs from S wind farms through eigenvalue analysis of the wind power covariance matrix, the PCA-based transformation allows us to recast the interrelated variables into the uncorrelated PCs. Thus, the correlated realization of wind power output wts within forecast interval [w_ts,w-ts] at t can be expressed by the realization of the uncorrelated PCs vtμ:(1)wt1,…,wtST=Φtvt1,…,vtST,where Φt is derived by PCA.Thus, the realization of total system wind powerwt=∑s=1Swts can be expressed by a linear combination [15] of the realization of PCs vtμ as ∑μ=1Sξtμvtμ. Here, ξtμ is simply the sum of the values on the column μ of Φt. As such, the impact of the S wind output realizations on the real power flow on line l is ∑s=1Sψln~swts, which we restate as the linear combination of the realization of PCs vtμ:∑μ=1SHl,tμvtμ. Here, Hl,tμ is the inner product of [ψln~1,ψln~2,…,ψln~S]T and column μ of Φt. We use the uncorrelated representation mentioned above in the following paper.The UC formulation also needs to explicitly account for the change on the objective function after the problem considers the uncertain wind issues. The objective function of conventional central UC problem formulation is to minimize the sum of start-up cost of the controllable generation resources in a power system and the fuel cost or economic dispatch (ED) cost. As the penetration of wind resources deepens, the uncertainty in the wind power outputs requires the ISOs to provide UC schedules that are robust for all possible realization of such random behavior. The decision making of a UC problem for system operators is thus divided into a two-stage stochastic program, where the first stage of the problem represents a day-ahead UC, and the second stage represents an hour-ahead ED of the entire system, given the fixed day-ahead schedule of controllable units and one of the realization of uncertain information. Thus, an effective quantification of the system ED costs in the new UC problem formulation is a “worst” value of all ED costs under a given UC decision and all realizations of uncertainties within the interval of uncertainty. The objective function of the new UC problem, stated in (2), is to minimize the sum of two terms: day-ahead start-up costs and a hour-ahead ED cost which is then maximized over the uncertainty interval with an associated confidence probability.(2)minpit,uit,rit,uiUPt,uiDNt⁡∑i=1Iσiuit+∑t=1Tχ-t.Here, we define∑i=1Iσi[ui(t)] as the start-up costs for the entire commitment horizon, where σi[·] is a function of ui(t) [21]. The “worst” ED cost at period t is defined as(3)χ-t=maxdt∈d_t,d¯t⁡χt,where χ(t), defined in (4), is the optimal solution of ED problem under one particular realization of system net loads:(4)χt=minpit∑i=1ICipits.t.Eq.5,6withgivendt∈d_t,d¯tandEq.7and8.Here,Ci[·] is the fuel cost of unit i and is usually an increasing convex quadratic function of the power output pi(t). The constraints in the optimization problem stated in (4) are the basic constraints of a standard ED problem; considering they are also the part of the key operating constraints in UC formulation, we will give their statement in the following part for describing all the related system and single unit operating constraints, including the constraints in (4): (5)∑i=1Ipit=dt,∀t=1,2,…,T,dt∈d_t,d¯t,where d_t=∑n=1Nd^tn-w¯t, d¯t=∑n=1Nd^tn-w_t.Equation (5) is the energy balance equation that matches the system supply by controllable units and net load at each time period. The system net load is defined by load demand minus system total wind power output. Since the day-ahead wind power output randomness is quantified by uncertainty interval, the net load dt will be within an interval [d_t,d¯t], as cast in (5). The energy balance constraints in (4) are represented by (5) with given dt.(6)flm≤∑i=1Iψln^ipit-∑n=1Nψlnd^tn+∑μ=1SHl,tμvtμ≤flM∑μ=1Sξtμvtμ=∑n=1Nd^tn-dt,dt∈d_t,d¯t,∀t=1,2,…,T,∀l=1,2,…,L.Equation (6) is the DC approximation of the nonlinear AC flow equations and has the form that the real power line flow is within the maximum and the minimum real power flow allowed [3]. The use of the nonlinear AC power flow modeling approach in the UC problem may not be practical and is computationally demanding. In the literature, it has been standard to simplify the modeling under some reasonable assumptions and making use of the so-called injection shift factors (ISFs) to replace the nonlinear AC real power flow equations by a set of linear equation, which is stated as (6) [22, 23]. The ISF ψln of line l is the (approximate) sensitivity of the change in the line l real power flow with respect to a change in the injection at node n and the withdrawal of an amount equal to the injection change at the slack bus. In (6), as mentioned above, we use the uncorrelated formulation ∑μ=1SHl,tμvtμ to represent ∑s=1Sψln~swts; the realization of total system wind power wt=∑s=1Swts, which is equal to total system load ∑n=1Nd^tn minus net system load dt, is represented by the uncorrelated formulation ∑μ=1Sξtμvtμ. The line power flow constraints in (4) is represented by (6) with given dt.(7)uit·pim≤pit≤uit·piM-rit,∀t=1,2,…,T,∀i=1,2,…,I.Equation (7) is the constraint that the generation level of unit i should be larger than minimum generation limits for unit i; the sum of generation level and spinning reserve should be smaller than maximum limits for unit i.(8)pit-1-Δi·uit≤pit≤pit-1+Δi·uit,∀t=1,2,…,T,∀i=1,2,…,I.Equation (8) is the ramp rate constraint; that is, the speed at which a unit can increase or decrease its generation level is bounded in a range.(9)∑i=1Irit≥r^with0≤rit≤uit·riM,∀t=1,2,…,T,∀i=1,2,…,I,where riM is the maximum reserves contribution provided by unit i within the specified response time and ramping capability constraint.The feasible region of the discrete decision variables of controllable units is formulated by the following constraints:(10)uit-uit-1=uiUPt-uiDNtuiUPt+uiDNt≤1,∀t=1,2,…,T,∀i=1,2,…,I,(11)uiUPt+∑t′=t+1minT,t+κ¯i-1uiDNt′≤1,uiDNt+∑t′=t+1minT,t+κ_i-1uiUPt′≤1,∀t=1,2,…,T,∀i=1,2,…,I.The specific form of the constraints is similar to [24]. Equations (10) are logic constraints between on and off status and the start-up and shut-down actions. Equations (11) include the minimum time a unit must operate in once it has been started up and the minimum shut-down time once a unit is decommitted.We may view the UC formulation above stated as (2)–(11) have a two-stage mathematical structure. The outer stage in (2) and (5)–(11) is a statement of the complex optimization problem for the UC problem with mixed integer/continuous decision variables and with the explicit representation of the confidence intervals. The inner stage in (3) and (4) is to get the “worst” ED charges for a fixed commitment and one particular realization of system net loads, which is an optimization problem with a max-min form. In the optimization problem stated by (4), the related constraints are (5) with given dt, (6) with given dt, and (7)-(8). For a realistic system, this is a large-scale nonlinear optimization problem.


## 3. Thrust of the Proposed Solution Approach

Solving the above two-stage optimization problem is challenging since there is the representation of the uncertainty by interval sets in (5) and (6) at the outer-stage problem and thus it needs to equivalently recast the original formulation into a deterministic counterpart, especially requiring handling the associated correlation of the wind power injection, before we address any approach to solve this problem; the max-min structure in (3) and (4) at the inner-stage problem makes the overall problem in general NP-hard and very difficult to solve.The solution approach we propose in this paper is as follows: at first, we equivalently recast the formulation of (2)–(11) into a deterministic form. To find the equivalently deterministic reformulation, we restate (5) and (6) in the outer-stage problem which have the uncertain sets [d_t,d¯t].After the reformulation, we present a new approach by converting the deterministic counterpart into a single-stage MILP form, to avoid the iteration and convergence issue that is dealt with within the decomposition approaches [7, 9] as mentioned in Section 1. To reach that goal, an approximate analytical solution of the inner-stage problem is obtained. The analytical solution obtained is aggregated as linear constraints into the deterministic reformulation without applying decomposition solution methods. The restated single-stage problem has a MILP form and thus can be easily implemented under the software environment of advanced commercial optimization packages.


### 3.1. Deterministic Reformulation of Problem by (2)–(11)

Reformulation of (5). We notice that the representation of (5) indicates that, for every possible realization of dt∈[d_t,d¯t], there is at least one pi(t), i=1,2,…,I, with the ability to maintain the energy balance for t=1,2,…,T. To make it clear, the notations are introduced as follows:p¯i(t): actual upper level of power output of unit i in period tp_i(t): actual lower level of power output of unit i in period tIn other words,p¯i(t) and p_i(t) are represented by the following equations:(12)p¯it=min⁡piM·uit-rit,pit-1+Δi·uit,p_it=max⁡pim·uit,pit-1-Δi·uit,∀t=1,2,…,T,∀i=1,2,…,I.Thus, (7)-(8) can be equivalently restated by (13): (13)p_it≤pit≤p¯it,∀t=1,2,…,T,∀i=1,2,…,I.Based on the above aided notations, the following theorem will give a proof that (5) can be equivalently represented by (14):(14)d_t≥∑i=1Ip_it,d¯t≤∑i=1Ip¯it,∀t=1,2,…,T.Theorem 1.
Equations (5), (7), and (8) in Section 2 can be equally represented by (12)–(14).Proof of Theorem1.
(  1) Sufficiency. If (5) and (7)-(8) hold, from (7)-(8), we can obtain (13) based on the definition of p¯i(t) and p_i(t) shown in (12); then from (13), we can have ∑i=1Ip_i(t)≤∑i=1Ipi(t)≤∑i=1Ip¯i(t) for t=1,2,…,T. Thus, from (5), we obtain (14).
(  2) Necessity. If (12)–(14) hold, from (12)-(13), we obtain (7)-(8) based on the definition of p¯i(t) and p_i(t); from (12)-(13), we have ∑i=1Ip_i(t)≤∑i=1Ipi(t)≤∑i=1Ip¯i(t) for t=1,2,…,T. Equation (14) indicates the interval [∑i=1Ip_i(t),∑i=1Ip¯i(t)] can cover the interval [d_t,d¯t] for t=1,2,…,T; thus, for every possible realization of dt∈[d_t,d¯t], there is at least one ∑i=1Ipi(t)∈[∑i=1Ip_i(t),∑i=1Ip¯i(t)], with the ability to maintain ∑i=1Ipi(t)=dt for t=1,2,…,T. Thus, (5) holds.Reformulation of (6). We notice that the representation of (6) indicates that, for every possible realization of dt∈[d_t,d¯t], there is at least one pi(t), i=1,2,…,I, which is restricted by (5) with given dt and (7)-(8), with the ability to maintain the real power flow limits of line l within the limits.We notice that, for a givendt∈[d_t,d¯t], this corresponds to multiple values of ∑μ=1SHl,tμvtμ under different realization of power output at the S wind sites in period t. Thus, to make it clear, the notations are introduced as follows:γ¯l,t: upper bound of ∑μ=1SHl,tμvtμ for a given dtγ_l,t: lower bound of ∑μ=1SHl,tμvtμ for a given dtFor a givendt, γ¯l,t and γ_l,t can be stated as the optimal solution of the following optimization problems:(15)γ¯l,t=maxvtμ,μ=1,2,…,S∑μ=1SHl,tμvtμs.t.∑μ=1Sξtμvtμ=∑n=1Nd^tn-dtw_ts≤vt1,vt2,…,vtSΦts≤w¯ts∀s=1,2,…,S,(16)γ_l,t=minvtμ,μ=1,2,…,S∑μ=1SHl,tμvtμs.t.∑μ=1Sξtμvtμ=∑n=1Nd^tn-dtw_ts≤vt1,vt2,…,vtSΦts≤w¯ts∀s=1,2,…,S.Equations (15) and (16) are intended separately to find the upper and lower bound of the impact of the S wind output realizations, for a givendt, on the real power flow of line l, where we state it as ∑μ=1SHl,tμvtμ. The constraints of the optimization problems in (15) and (16) are as follows: the realization of total system wind power ∑μ=1Sξtμvtμ is equal to total system load ∑n=1Nd^tn minus net system load dt; the realization of the system wind power output in period t, which is represented by [vt1,vt2,…,vtS]Φts, is within the interval [w_ts,w-ts].From (15)-(16), we notice that, for any given dt∈[d_t,d¯t], to maintain the power flow limits on line l, it has the following constraints: (17)flm≤∑i=1Iψln^ipit+γ_l,t-∑n=1Nψlnd^tn,∑i=1Iψln^ipit+γ¯l,t-∑n=1Nψlnd^tn≤flM,∀t=1,2,…,T,∀l=1,2,…,L.In other words, for any givendt∈[d_t,d¯t], (6) can be equivalently recast by (17).From (15)-(16), we notice that γ¯l,t and γ_l,t are a function of dt. Furthermore, the value of pi(t) in (17) is restricted by (5) with a given dt∈[d_t,d¯t] and (7)-(8); thus, ∑i=1Iψln^ipi(t) is also a function of dt.We introduce other notations as follows:k~1l(t): minimum value of ∑i=1Iψln^ipi(t) subject to (7)-(8) and (5) with a given dt∈[d_t,d¯t]k~2l(t): maximum value of ∑i=1Iψln^ipi(t) subject to (7)-(8) and (5) with a given dt∈[d_t,d¯t]k1l(t): upper value of k~1l(t)+γ¯l,t within [d_t,d¯t]k2l(t): lower value of k~2l(t)+γ_l,t within [d_t,d¯t]Based on the above definition, we havek1l(t), k2l(t), k~1l(t), and k~2l(t):(18)k1lt=maxdt∈d_t,d¯tk~1lt+γ¯l,t,(19)k2lt=mindt∈d_t,d¯tk~2lt+γ_l,t,where k~1l(t) and k~2l(t) are the optimal solution of optimization problems stated by (20) and (21) with given dt∈[d_t,d¯t]. Here (7)-(8) have been represented by (13) based on Theorem 1:(20)minpit∑i=1Iψln^ipit,s.t.∑i=1Ipit=dt,p_it≤pit≤p¯it,∀i=1,2,…,I,(21)maxpit∑i=1Iψln^ipit,s. t.∑i=1Ipit=dt,p_it≤pit≤p¯it,∀i=1,2,…,I.From the above notations, the following theorem will give a proof that (6) can be equivalently recast by (22): (22)flm≤k2lt-∑n=1Nψlnd^tn,k1lt-∑n=1Nψlnd^tn≤flM,∀t=1,2,…,T,∀l=1,2,…,L.Theorem 2.
Equation (17) can be equivalently recast by (22).Proof of Theorem2.
(  1) Sufficiency. If (17) holds, then, for any given dt∈[d_t,d¯t], we have flm≤∑i=1Iψln^ipi(t)+γ_l,t-∑n=1Nψlnd^tn. Then we have flm≤k~2l(t)+γ_l,t-∑n=1Nψlnd^tn for any given dt∈[d_t,d¯t] based on the definition of k~2l(t). For that matter, we can have flm≤k2l(t)-∑n=1Nψlnd^tn based on the definition of k2l(t). Similarly, if (17) holds, then, for any given dt∈[d_t,d¯t], we have ∑i=1Iψln^ipi(t)+γ¯l,t-∑n=1Nψlnd^tn≤flM. Then we have k~1l(t)+γ¯l,t-∑n=1Nψlnd^tn≤flM for any given dt∈[d_t,d¯t] based on the definition of k~1l(t). For that matter, we can have k1l(t)-∑n=1Nψlnd^tn≤flM based on the definition of k1l(t). To conclude, (22) holds.
(  2) Necessity. If (22) holds, then, based on the above notations, we have flm≤k2l(t)-∑n=1Nψlnd^tn and k1l(t)-∑n=1Nψlnd^tn≤flM. Thus, we can have flm≤k~2l(t)+γ_l,t-∑n=1Nψlnd^tn and k~1l(t)+γ¯l,t-∑n=1Nψlnd^tn≤flM for any given dt∈[d_t,d¯t] based on the definition of k~1l(t) and k~2l(t). For that matter, for any given dt∈[d_t,d¯t], there at least exists one pi1(t), i=1,2,…,I, which is the optimal solution of the optimization problem in (21) with the ability to make flm≤∑i=1Iψln^ipi1(t)+γ_l,t-∑n=1Nψlnd^tn and at least exists one pi2(t), i=1,2,…,I, which is the optimal solution of the optimization problem in (20) with the ability to make ∑i=1Iψln^ipi2(t)+γ¯l,t-∑n=1Nψlnd^tn≤flM. Under the large penetration level of wind power for a power system, it usually has ∑i=1Iψln^ipi1(t)+γ_l,t≤∑i=1Iψln^ipi2(t)+γ¯l,t; thus, for any line l, we can always find at least one pi(t), i=1,2,…,I, to make (17) feasible. Furthermore, from the real engineering practice and the definition of the injection shift factor of line l at node n [23], for any line l, we usually have the same result to make an order for (ψl1,ψl2,…,ψlN) from the smallest value to the largest value. In other words, for any line l, the optimal solution for (21) or (20) is the same and thus there usually exists a common one, pi(t), which is identical for different transmission lines. To conclude, (17) holds.
To sum up, through the deterministic reformulation, the deterministic counterpart is an optimization problem stated by (2)–(4), (9)–(14), and (22). The inner-stage problem is represented by (3)-(4), where the constraints in the inner-stage problem have been recast as (12)–(14) and (22); the outer-stage problem is stated by (2), (9)–(14), and (22).


### 3.2. Recast the Deterministic Counterpart into a Single-Level MILP Form

The deterministic counterpart has a two-stage structure: the outer-stage problem is represented by (2), (9)–(14), and (22) where k1l(t) and k2l(t) in (22) are optimization problems given by (18)–(21); the inner-stage problem χ-(t) represented by (3) and (4) is an optimization problem with a max-min form in which χ(t) is usually a quadratic programming problem and is to determine the dispatch cost for a given wind power realization, which is then maximized over [d_t,d¯t]. To avoid the iteration and convergence issue that is dealt with within decomposition approaches, an exact analytical solution of k1l(t) and k2l(t) in the outer-stage problem and an approximate analytical solution of χ-(t) in the inner-stage problem are obtained. The analytical solutions obtained are aggregated as linear constraints into the deterministic reformulation.An Approximate Analytical Solution of χ - ( t ). To make the problem represented by (3)-(4) easy to solve, we make the following assumptions: use linear function representing Ci[·], not considering (22). The second assumption means the inner-stage optimization problem does not consider the line power flow limits. However, the outer-stage optimization problem still keeps line power flow limits as a critical constraint which is represented by (22).The reason that the above assumptions are reasonable and involving no loss of generality are as follows:(1) the coefficient of the quadratic term of Ci[·] is usually small; thus the cost error by the linear approximation is small [25]; (2) χ-(t) is a recourse cost which represents the “worst” possible ED cost but is not to be the real ED cost; the real ED cost will be determined with the final realization of all the wind power outputs. In the UC problem formulation, χ-(t) is to provide a reference to system operator that χ-(t) is the largest possible fuel cost under a given commitment result. Thus, doing an approximation of χ-(t) only means that it does an approximation of the “worst” possible ED cost; it does not mean the real ED cost will be largely changed after the approximation.Thus, we assumeCi[pi(t)] in (4) as bi·pi(t)+ci, where bi and ci are coefficients. Equation (4) is stated as follows for given dt:(23)χt=minpit∑i=1Ibi·pit+ci·uit,s.t.∑i=1Ipit=dt,p_it≤pit≤p¯it,∀i=1,2,…,I.Equation (23) is a linear programming (LP) problem, where the related constraints are as follows: (5) with given dt; (7)-(8) which are recast by (13).We find that the LP problem in (23) has the same structure as that of (19) in [26] and thus can obtain an analytical condition of χ(t) directly without solving (23). Let all controllable units permute on 0≤b1≤b2≤⋯≤bI; based on Theorem 4 in [26], we obtain the optimal solution to the LP problem in (23) as the following equations:(24)pm∗t=p¯mtifm≤k′-1,pk′∗t=dt-∑m=1k′-1p¯mt-∑m=k′+1Ip_mt,pm∗t=p_mtifm>k′.The value of the optimal objective function,χ(t), is given as(25)χt=∑i=1Ici·uit+∑m=1τ-1bm-bτp¯mt+∑m=τ+1Ibm-bτp_mt+bτdt∀t=1,2,…,T,where bτ are permuted on nondecreasing order. The integer number τ, τ=1,2,…,I, satisfies(26)ητt=∑m=1τp¯mt-p_mt+∑m=1Ip_mt,ητ-1t≤dt≤ητt.We notice that (25) is a piecewise linear nondecreasing function of dt∈[d_t,d¯t], since bτ≥0, τ=1,2,…,I, and is on a nondecreasing order. Thus, the optimal solution to the optimization problem in (3) is given by dt=d¯t. The value of the optimal objective function, χ-(t), is given as(27)χ¯t=∑i=1Ici·uit+∑m=1τ-1bm-bτp¯mt+∑m=τ+1Ibm-bτp_mt+bτdt,where τ satisfies ητ-1(t)≤d¯t≤ητ(t).Thus, (27) is the approximate analytical solution of the inner-stage optimization problem; however, it is represented as a nonlinear form, since it is nontrivial to get what τ in 1,2,…,I satisfies ητ-1(t)≤d¯t≤ητ(t).The following theorem gives a proof that (27) can be equivalently recast by series of linear constraints and thus make the inner-stage problem be aggregated as a MILP form. To aid analysis, we introduce the following variables ςτ(t), ∀τ=1,2,…,I, where bτ are permuted on a nondecreasing order:(28)ςτt=∑i=1Ici·uit+∑m=1τ-1bm-bτp¯mt+∑m=τ+1Ibm-bτp_mt+bτdt∀t=1,2,…,T.Theorem 3.
Equation (27) can be equivalently represented by a series of linear constraints.Proof of Theorem3.
From (28), ∀τ=1,2,…,I-1, ςτ+1(t)-ςτ(t) can be restated by(29)ςτ+1t-ςτt=bτ-bτ+1ητt-d¯t∀t=1,2,…,T.
From (29), we notice that when τ′ is an integer number (1≤τ′≤I) such that ητ′-1t≤d¯t≤ητ′(t), with τ increasing from 1 to τ′-1, the values of ςτ+1(t)-ςτ(t) are nonnegative; with τ increasing from τ′ to I-1, the values of ςτ+1(t)-ςτ(t) are nonpositive. Thus, the maximum value of ςτ(t) for τ=1,2,…,I is ςτ′(t) and we have χ¯(t)=ςτ′(t) from (27). In (2), the objective function is to minimize the total cost; thus (27) can be equivalently recast in terms of the following linear constraints:(30)χ¯t≥ςτt,∀t=1,2,…,T,∀τ=1,2,…,I.
Thus, the inner-stage problem in (3)-(4) is recast by the linear constraints (28) and (30).Exact Analytical Solution of k 1 l ( t ). From (18), (20), and (15), we notice that k1l(t) has a form where, in (20), k~1l(t) is to determine the minimum line power flow injected from all controllable units and, in (15), γ¯l,t is to determine the maximum line power flow injected from all wind sites; in (18), k1l(t) is to determine the maximum line power flow of the sum of k~1l(t) and γ¯l,t over [d_t,d¯t].Similar to the way we get the analytical solution of (23), based on Theorem 4 in [26], we obtain the value of the optimal objective function of the LP problem in (20), k~1l(t), as follows:(31)k~1lt=∑m=1j-1ψln^m-ψln^jp¯mt+∑m=j+1Iψln^m-ψln^jp_mt+ψln^jdt∀t=1,2,…,T,where the units i are permuted on ψln^1≤ψln^2≤⋯≤ψln^I. The integer number j, j=1,2,…,I, satisfies(32)ρjt=∑m=1jp¯mt-p_mt+∑m=1Ip_mt,ρj-1t≤dt≤ρjt.To obtain the analytical solution of (15), we firstly state (15) by its dual form with a similar structure to (23). With the Strong Duality theorem, γ¯l,t is also the optimal solution of the dual problem. We obtain the optimal solution of the dual problem of (15) (described in Theorem 4) based on [26, Theorem 4], as follows:(33)γ¯l,t=∑m=1δ-1ψln~m-ψln~δw¯tm-w_tm+∑n=1Nd^tn-dtψln~δ+∑s=1Sw_tsψln~s-ψln~δ∀t=1,2,…,T,where the wind sites are permuted on ψln~1≥ψln~2≥⋯≥ψln~S.The integer numberδ, δ=1,2,…,S, satisfies(34)ωδ,t=∑n=1Nd^tn-∑m=1δw¯tm-w_tm-∑s=1Sw_ts,ωδ,t≤dt≤ωδ-1,t.We notice from (31) and (33) that k~1l(t) and γ¯l,t are both piecewise linear functions of dt. Thus, all the values of k~1l(t) and γ¯l,t over [d_t,d¯t] can be represented by a finitely affine function of their extreme points, which separately are the values of k~1l(t) and γ¯l,t when dt=ρj(t), j=1,2,…,I, and dt=ωδ,t, δ=1,2,…,S. Similar to Theorem 3, we give a proof in [27] that the maximum values of k~1l(t) and γ¯l,t over [d_t,d¯t] are reached at dt=d_t or d¯t. Based on the definition of k1l(t), we obtain that the value of k1l(t) is reached at dt=d_t or d¯t. Furthermore, k1l(t) can be also recast as a linear form [27]. The linear representation of k2l(t) can also be stated in a similar way.Theorem 4.
λ l , t +, αl,ts, and βl,tss=1,2,…,S are the dual variables for the constraints in (15). We permute wind sites on ψln~1≥ψln~2≥⋯≥ψln~S; the dual problem is stated as follows:(35)minλl,t+,αl,ts,βl,ts∑s=1Sw¯tsαl,ts-∑s=1Sw_tsβl,ts+∑n=1Nd^tn-dtλl,t+(36)s.t.ΦtT⋮-ΦtT⋮∑s=1SΦtsαl,t1⋮αl,tSβl,t1⋮βl,tSλl,t+=Hl,t1⋮Hl,tS(37)αl,ts≥0,βl,ts≥0,∀s=1,2,…,S.
From PCA, we knowΦtT is invertible, and (36) is recast as(38)αl,t1-βl,t1⋮αl,tS-βl,tS=ψln~1⋮ψln~S-λl,t+⋮λl,t+.
The following is a proof that the optimal solution of the dual problem stated by (35)–(37) is(39)αl,t1=ψln~1-ψln~δ,⋮αl,tδ-1=ψln~δ-1-ψln~δ,αl,tδ=0,⋮αl,tS=0,λl,t+=ψln~δ,βl,t1=0,⋮βl,tδ=0,βl,tδ+1=ψln~δ-ψln~δ+1,⋮βl,tS=ψln~δ-ψln~S,where ωδ,t≤dt≤ωδ-1,t.Proof of Theorem4.
We restate (35) as (40) based on (38):(40)∑s=1Sw¯ts-w_tsαl,ts+∑s=1Sw_tsψln~s+∑n=1Nd^tn-dt-∑s=1Sw_tsλl,t+.
Ifλl,t+ is +∞ or -∞, the value in (40) is +∞; if λl,t+ decreases, the third term of (40) decreases while the first term increases. If λl,t+>ψln~δ, the increasing rate of the first term is less than the decreasing rate of the third term; thus, the whole value in (40) decreases. Otherwise, the whole value increases. Thus, the optimal solution of the dual problem is the value represented by (39).Based on the above proof,γ¯l,t can be recast by (33).


## 3.1. Deterministic Reformulation of Problem by (2)–(11)

Reformulation of (5). We notice that the representation of (5) indicates that, for every possible realization of dt∈[d_t,d¯t], there is at least one pi(t), i=1,2,…,I, with the ability to maintain the energy balance for t=1,2,…,T. To make it clear, the notations are introduced as follows:p¯i(t): actual upper level of power output of unit i in period tp_i(t): actual lower level of power output of unit i in period tIn other words,p¯i(t) and p_i(t) are represented by the following equations:(12)p¯it=min⁡piM·uit-rit,pit-1+Δi·uit,p_it=max⁡pim·uit,pit-1-Δi·uit,∀t=1,2,…,T,∀i=1,2,…,I.Thus, (7)-(8) can be equivalently restated by (13): (13)p_it≤pit≤p¯it,∀t=1,2,…,T,∀i=1,2,…,I.Based on the above aided notations, the following theorem will give a proof that (5) can be equivalently represented by (14):(14)d_t≥∑i=1Ip_it,d¯t≤∑i=1Ip¯it,∀t=1,2,…,T.Theorem 1.
Equations (5), (7), and (8) in Section 2 can be equally represented by (12)–(14).Proof of Theorem1.
(  1) Sufficiency. If (5) and (7)-(8) hold, from (7)-(8), we can obtain (13) based on the definition of p¯i(t) and p_i(t) shown in (12); then from (13), we can have ∑i=1Ip_i(t)≤∑i=1Ipi(t)≤∑i=1Ip¯i(t) for t=1,2,…,T. Thus, from (5), we obtain (14).
(  2) Necessity. If (12)–(14) hold, from (12)-(13), we obtain (7)-(8) based on the definition of p¯i(t) and p_i(t); from (12)-(13), we have ∑i=1Ip_i(t)≤∑i=1Ipi(t)≤∑i=1Ip¯i(t) for t=1,2,…,T. Equation (14) indicates the interval [∑i=1Ip_i(t),∑i=1Ip¯i(t)] can cover the interval [d_t,d¯t] for t=1,2,…,T; thus, for every possible realization of dt∈[d_t,d¯t], there is at least one ∑i=1Ipi(t)∈[∑i=1Ip_i(t),∑i=1Ip¯i(t)], with the ability to maintain ∑i=1Ipi(t)=dt for t=1,2,…,T. Thus, (5) holds.Reformulation of (6). We notice that the representation of (6) indicates that, for every possible realization of dt∈[d_t,d¯t], there is at least one pi(t), i=1,2,…,I, which is restricted by (5) with given dt and (7)-(8), with the ability to maintain the real power flow limits of line l within the limits.We notice that, for a givendt∈[d_t,d¯t], this corresponds to multiple values of ∑μ=1SHl,tμvtμ under different realization of power output at the S wind sites in period t. Thus, to make it clear, the notations are introduced as follows:γ¯l,t: upper bound of ∑μ=1SHl,tμvtμ for a given dtγ_l,t: lower bound of ∑μ=1SHl,tμvtμ for a given dtFor a givendt, γ¯l,t and γ_l,t can be stated as the optimal solution of the following optimization problems:(15)γ¯l,t=maxvtμ,μ=1,2,…,S∑μ=1SHl,tμvtμs.t.∑μ=1Sξtμvtμ=∑n=1Nd^tn-dtw_ts≤vt1,vt2,…,vtSΦts≤w¯ts∀s=1,2,…,S,(16)γ_l,t=minvtμ,μ=1,2,…,S∑μ=1SHl,tμvtμs.t.∑μ=1Sξtμvtμ=∑n=1Nd^tn-dtw_ts≤vt1,vt2,…,vtSΦts≤w¯ts∀s=1,2,…,S.Equations (15) and (16) are intended separately to find the upper and lower bound of the impact of the S wind output realizations, for a givendt, on the real power flow of line l, where we state it as ∑μ=1SHl,tμvtμ. The constraints of the optimization problems in (15) and (16) are as follows: the realization of total system wind power ∑μ=1Sξtμvtμ is equal to total system load ∑n=1Nd^tn minus net system load dt; the realization of the system wind power output in period t, which is represented by [vt1,vt2,…,vtS]Φts, is within the interval [w_ts,w-ts].From (15)-(16), we notice that, for any given dt∈[d_t,d¯t], to maintain the power flow limits on line l, it has the following constraints: (17)flm≤∑i=1Iψln^ipit+γ_l,t-∑n=1Nψlnd^tn,∑i=1Iψln^ipit+γ¯l,t-∑n=1Nψlnd^tn≤flM,∀t=1,2,…,T,∀l=1,2,…,L.In other words, for any givendt∈[d_t,d¯t], (6) can be equivalently recast by (17).From (15)-(16), we notice that γ¯l,t and γ_l,t are a function of dt. Furthermore, the value of pi(t) in (17) is restricted by (5) with a given dt∈[d_t,d¯t] and (7)-(8); thus, ∑i=1Iψln^ipi(t) is also a function of dt.We introduce other notations as follows:k~1l(t): minimum value of ∑i=1Iψln^ipi(t) subject to (7)-(8) and (5) with a given dt∈[d_t,d¯t]k~2l(t): maximum value of ∑i=1Iψln^ipi(t) subject to (7)-(8) and (5) with a given dt∈[d_t,d¯t]k1l(t): upper value of k~1l(t)+γ¯l,t within [d_t,d¯t]k2l(t): lower value of k~2l(t)+γ_l,t within [d_t,d¯t]Based on the above definition, we havek1l(t), k2l(t), k~1l(t), and k~2l(t):(18)k1lt=maxdt∈d_t,d¯tk~1lt+γ¯l,t,(19)k2lt=mindt∈d_t,d¯tk~2lt+γ_l,t,where k~1l(t) and k~2l(t) are the optimal solution of optimization problems stated by (20) and (21) with given dt∈[d_t,d¯t]. Here (7)-(8) have been represented by (13) based on Theorem 1:(20)minpit∑i=1Iψln^ipit,s.t.∑i=1Ipit=dt,p_it≤pit≤p¯it,∀i=1,2,…,I,(21)maxpit∑i=1Iψln^ipit,s. t.∑i=1Ipit=dt,p_it≤pit≤p¯it,∀i=1,2,…,I.From the above notations, the following theorem will give a proof that (6) can be equivalently recast by (22): (22)flm≤k2lt-∑n=1Nψlnd^tn,k1lt-∑n=1Nψlnd^tn≤flM,∀t=1,2,…,T,∀l=1,2,…,L.Theorem 2.
Equation (17) can be equivalently recast by (22).Proof of Theorem2.
(  1) Sufficiency. If (17) holds, then, for any given dt∈[d_t,d¯t], we have flm≤∑i=1Iψln^ipi(t)+γ_l,t-∑n=1Nψlnd^tn. Then we have flm≤k~2l(t)+γ_l,t-∑n=1Nψlnd^tn for any given dt∈[d_t,d¯t] based on the definition of k~2l(t). For that matter, we can have flm≤k2l(t)-∑n=1Nψlnd^tn based on the definition of k2l(t). Similarly, if (17) holds, then, for any given dt∈[d_t,d¯t], we have ∑i=1Iψln^ipi(t)+γ¯l,t-∑n=1Nψlnd^tn≤flM. Then we have k~1l(t)+γ¯l,t-∑n=1Nψlnd^tn≤flM for any given dt∈[d_t,d¯t] based on the definition of k~1l(t). For that matter, we can have k1l(t)-∑n=1Nψlnd^tn≤flM based on the definition of k1l(t). To conclude, (22) holds.
(  2) Necessity. If (22) holds, then, based on the above notations, we have flm≤k2l(t)-∑n=1Nψlnd^tn and k1l(t)-∑n=1Nψlnd^tn≤flM. Thus, we can have flm≤k~2l(t)+γ_l,t-∑n=1Nψlnd^tn and k~1l(t)+γ¯l,t-∑n=1Nψlnd^tn≤flM for any given dt∈[d_t,d¯t] based on the definition of k~1l(t) and k~2l(t). For that matter, for any given dt∈[d_t,d¯t], there at least exists one pi1(t), i=1,2,…,I, which is the optimal solution of the optimization problem in (21) with the ability to make flm≤∑i=1Iψln^ipi1(t)+γ_l,t-∑n=1Nψlnd^tn and at least exists one pi2(t), i=1,2,…,I, which is the optimal solution of the optimization problem in (20) with the ability to make ∑i=1Iψln^ipi2(t)+γ¯l,t-∑n=1Nψlnd^tn≤flM. Under the large penetration level of wind power for a power system, it usually has ∑i=1Iψln^ipi1(t)+γ_l,t≤∑i=1Iψln^ipi2(t)+γ¯l,t; thus, for any line l, we can always find at least one pi(t), i=1,2,…,I, to make (17) feasible. Furthermore, from the real engineering practice and the definition of the injection shift factor of line l at node n [23], for any line l, we usually have the same result to make an order for (ψl1,ψl2,…,ψlN) from the smallest value to the largest value. In other words, for any line l, the optimal solution for (21) or (20) is the same and thus there usually exists a common one, pi(t), which is identical for different transmission lines. To conclude, (17) holds.
To sum up, through the deterministic reformulation, the deterministic counterpart is an optimization problem stated by (2)–(4), (9)–(14), and (22). The inner-stage problem is represented by (3)-(4), where the constraints in the inner-stage problem have been recast as (12)–(14) and (22); the outer-stage problem is stated by (2), (9)–(14), and (22).


## 3.2. Recast the Deterministic Counterpart into a Single-Level MILP Form

The deterministic counterpart has a two-stage structure: the outer-stage problem is represented by (2), (9)–(14), and (22) where k1l(t) and k2l(t) in (22) are optimization problems given by (18)–(21); the inner-stage problem χ-(t) represented by (3) and (4) is an optimization problem with a max-min form in which χ(t) is usually a quadratic programming problem and is to determine the dispatch cost for a given wind power realization, which is then maximized over [d_t,d¯t]. To avoid the iteration and convergence issue that is dealt with within decomposition approaches, an exact analytical solution of k1l(t) and k2l(t) in the outer-stage problem and an approximate analytical solution of χ-(t) in the inner-stage problem are obtained. The analytical solutions obtained are aggregated as linear constraints into the deterministic reformulation.An Approximate Analytical Solution of χ - ( t ). To make the problem represented by (3)-(4) easy to solve, we make the following assumptions: use linear function representing Ci[·], not considering (22). The second assumption means the inner-stage optimization problem does not consider the line power flow limits. However, the outer-stage optimization problem still keeps line power flow limits as a critical constraint which is represented by (22).The reason that the above assumptions are reasonable and involving no loss of generality are as follows:(1) the coefficient of the quadratic term of Ci[·] is usually small; thus the cost error by the linear approximation is small [25]; (2) χ-(t) is a recourse cost which represents the “worst” possible ED cost but is not to be the real ED cost; the real ED cost will be determined with the final realization of all the wind power outputs. In the UC problem formulation, χ-(t) is to provide a reference to system operator that χ-(t) is the largest possible fuel cost under a given commitment result. Thus, doing an approximation of χ-(t) only means that it does an approximation of the “worst” possible ED cost; it does not mean the real ED cost will be largely changed after the approximation.Thus, we assumeCi[pi(t)] in (4) as bi·pi(t)+ci, where bi and ci are coefficients. Equation (4) is stated as follows for given dt:(23)χt=minpit∑i=1Ibi·pit+ci·uit,s.t.∑i=1Ipit=dt,p_it≤pit≤p¯it,∀i=1,2,…,I.Equation (23) is a linear programming (LP) problem, where the related constraints are as follows: (5) with given dt; (7)-(8) which are recast by (13).We find that the LP problem in (23) has the same structure as that of (19) in [26] and thus can obtain an analytical condition of χ(t) directly without solving (23). Let all controllable units permute on 0≤b1≤b2≤⋯≤bI; based on Theorem 4 in [26], we obtain the optimal solution to the LP problem in (23) as the following equations:(24)pm∗t=p¯mtifm≤k′-1,pk′∗t=dt-∑m=1k′-1p¯mt-∑m=k′+1Ip_mt,pm∗t=p_mtifm>k′.The value of the optimal objective function,χ(t), is given as(25)χt=∑i=1Ici·uit+∑m=1τ-1bm-bτp¯mt+∑m=τ+1Ibm-bτp_mt+bτdt∀t=1,2,…,T,where bτ are permuted on nondecreasing order. The integer number τ, τ=1,2,…,I, satisfies(26)ητt=∑m=1τp¯mt-p_mt+∑m=1Ip_mt,ητ-1t≤dt≤ητt.We notice that (25) is a piecewise linear nondecreasing function of dt∈[d_t,d¯t], since bτ≥0, τ=1,2,…,I, and is on a nondecreasing order. Thus, the optimal solution to the optimization problem in (3) is given by dt=d¯t. The value of the optimal objective function, χ-(t), is given as(27)χ¯t=∑i=1Ici·uit+∑m=1τ-1bm-bτp¯mt+∑m=τ+1Ibm-bτp_mt+bτdt,where τ satisfies ητ-1(t)≤d¯t≤ητ(t).Thus, (27) is the approximate analytical solution of the inner-stage optimization problem; however, it is represented as a nonlinear form, since it is nontrivial to get what τ in 1,2,…,I satisfies ητ-1(t)≤d¯t≤ητ(t).The following theorem gives a proof that (27) can be equivalently recast by series of linear constraints and thus make the inner-stage problem be aggregated as a MILP form. To aid analysis, we introduce the following variables ςτ(t), ∀τ=1,2,…,I, where bτ are permuted on a nondecreasing order:(28)ςτt=∑i=1Ici·uit+∑m=1τ-1bm-bτp¯mt+∑m=τ+1Ibm-bτp_mt+bτdt∀t=1,2,…,T.Theorem 3.
Equation (27) can be equivalently represented by a series of linear constraints.Proof of Theorem3.
From (28), ∀τ=1,2,…,I-1, ςτ+1(t)-ςτ(t) can be restated by(29)ςτ+1t-ςτt=bτ-bτ+1ητt-d¯t∀t=1,2,…,T.
From (29), we notice that when τ′ is an integer number (1≤τ′≤I) such that ητ′-1t≤d¯t≤ητ′(t), with τ increasing from 1 to τ′-1, the values of ςτ+1(t)-ςτ(t) are nonnegative; with τ increasing from τ′ to I-1, the values of ςτ+1(t)-ςτ(t) are nonpositive. Thus, the maximum value of ςτ(t) for τ=1,2,…,I is ςτ′(t) and we have χ¯(t)=ςτ′(t) from (27). In (2), the objective function is to minimize the total cost; thus (27) can be equivalently recast in terms of the following linear constraints:(30)χ¯t≥ςτt,∀t=1,2,…,T,∀τ=1,2,…,I.
Thus, the inner-stage problem in (3)-(4) is recast by the linear constraints (28) and (30).Exact Analytical Solution of k 1 l ( t ). From (18), (20), and (15), we notice that k1l(t) has a form where, in (20), k~1l(t) is to determine the minimum line power flow injected from all controllable units and, in (15), γ¯l,t is to determine the maximum line power flow injected from all wind sites; in (18), k1l(t) is to determine the maximum line power flow of the sum of k~1l(t) and γ¯l,t over [d_t,d¯t].Similar to the way we get the analytical solution of (23), based on Theorem 4 in [26], we obtain the value of the optimal objective function of the LP problem in (20), k~1l(t), as follows:(31)k~1lt=∑m=1j-1ψln^m-ψln^jp¯mt+∑m=j+1Iψln^m-ψln^jp_mt+ψln^jdt∀t=1,2,…,T,where the units i are permuted on ψln^1≤ψln^2≤⋯≤ψln^I. The integer number j, j=1,2,…,I, satisfies(32)ρjt=∑m=1jp¯mt-p_mt+∑m=1Ip_mt,ρj-1t≤dt≤ρjt.To obtain the analytical solution of (15), we firstly state (15) by its dual form with a similar structure to (23). With the Strong Duality theorem, γ¯l,t is also the optimal solution of the dual problem. We obtain the optimal solution of the dual problem of (15) (described in Theorem 4) based on [26, Theorem 4], as follows:(33)γ¯l,t=∑m=1δ-1ψln~m-ψln~δw¯tm-w_tm+∑n=1Nd^tn-dtψln~δ+∑s=1Sw_tsψln~s-ψln~δ∀t=1,2,…,T,where the wind sites are permuted on ψln~1≥ψln~2≥⋯≥ψln~S.The integer numberδ, δ=1,2,…,S, satisfies(34)ωδ,t=∑n=1Nd^tn-∑m=1δw¯tm-w_tm-∑s=1Sw_ts,ωδ,t≤dt≤ωδ-1,t.We notice from (31) and (33) that k~1l(t) and γ¯l,t are both piecewise linear functions of dt. Thus, all the values of k~1l(t) and γ¯l,t over [d_t,d¯t] can be represented by a finitely affine function of their extreme points, which separately are the values of k~1l(t) and γ¯l,t when dt=ρj(t), j=1,2,…,I, and dt=ωδ,t, δ=1,2,…,S. Similar to Theorem 3, we give a proof in [27] that the maximum values of k~1l(t) and γ¯l,t over [d_t,d¯t] are reached at dt=d_t or d¯t. Based on the definition of k1l(t), we obtain that the value of k1l(t) is reached at dt=d_t or d¯t. Furthermore, k1l(t) can be also recast as a linear form [27]. The linear representation of k2l(t) can also be stated in a similar way.Theorem 4.
λ l , t +, αl,ts, and βl,tss=1,2,…,S are the dual variables for the constraints in (15). We permute wind sites on ψln~1≥ψln~2≥⋯≥ψln~S; the dual problem is stated as follows:(35)minλl,t+,αl,ts,βl,ts∑s=1Sw¯tsαl,ts-∑s=1Sw_tsβl,ts+∑n=1Nd^tn-dtλl,t+(36)s.t.ΦtT⋮-ΦtT⋮∑s=1SΦtsαl,t1⋮αl,tSβl,t1⋮βl,tSλl,t+=Hl,t1⋮Hl,tS(37)αl,ts≥0,βl,ts≥0,∀s=1,2,…,S.
From PCA, we knowΦtT is invertible, and (36) is recast as(38)αl,t1-βl,t1⋮αl,tS-βl,tS=ψln~1⋮ψln~S-λl,t+⋮λl,t+.
The following is a proof that the optimal solution of the dual problem stated by (35)–(37) is(39)αl,t1=ψln~1-ψln~δ,⋮αl,tδ-1=ψln~δ-1-ψln~δ,αl,tδ=0,⋮αl,tS=0,λl,t+=ψln~δ,βl,t1=0,⋮βl,tδ=0,βl,tδ+1=ψln~δ-ψln~δ+1,⋮βl,tS=ψln~δ-ψln~S,where ωδ,t≤dt≤ωδ-1,t.Proof of Theorem4.
We restate (35) as (40) based on (38):(40)∑s=1Sw¯ts-w_tsαl,ts+∑s=1Sw_tsψln~s+∑n=1Nd^tn-dt-∑s=1Sw_tsλl,t+.
Ifλl,t+ is +∞ or -∞, the value in (40) is +∞; if λl,t+ decreases, the third term of (40) decreases while the first term increases. If λl,t+>ψln~δ, the increasing rate of the first term is less than the decreasing rate of the third term; thus, the whole value in (40) decreases. Otherwise, the whole value increases. Thus, the optimal solution of the dual problem is the value represented by (39).Based on the above proof,γ¯l,t can be recast by (33).


## 4. Numerical Tests

### 4.1. Motivation

In this section, we will illustrate the performance of the proposed solution approach on a large-scale system and discuss its economic efficiency and computational advantage. To do so, we firstly collect a real data set which can represent a power system with large-scale wind power output. Then, we make a design for the whole numerical tests where the test data is realistic and applied to the future facts that the wind power output will be increased as the largest power generation resource within a power system. Finally, the related test results corresponding to the computational and economic efficiency are demonstrated through the computation time, the average dispatch and start-up cost, and the unit commitment results compared with the conventional reserve method in previous publications.


### 4.2. Numerical Tests Design

#### 4.2.1. Power System Description

The test dataset is from the power system operated by Northwest China Grid (NCG). NCG has a large-scale wind resource which is mainly in Hexi Corridor within Jiuquan area of Gansu Province, and its energy has been to form into three large aggregated wind farms around this area with estimated wind capacity installation as 5,450 MW, 6,150 MW, and 1,110 MW in 2015. All the wind farms are within a 130 km “belt” surrounded by Qilian and Bei Mountains; thus their wind power outputs generally demonstrate a characteristic of high correlation. The wind penetration level is defined as the ratio of the installed wind capacity to the peak load and in the test cases it is around 35%. This penetration level is large enough to represent the future tendency that the wind power is integrated into power system and can meet the basic engineering assumption we made in Section3. In this paper, our test dataset is based on four representative wind farms and their capacity at the 2015 level. The wind data is collected from the EMS system of NCG: we selected 5-minute wind power output for a whole month and the PCA matrix is acquired for every hour. The locations of the wind sites are presented in Table 1.Table 1
Location and capacity of wind sites (MW).
Site location Bus nodes name Capacity (MW) Chang-ma Yu-men-zhen 1,110 Gan-he-kou Dong An-xi 3,500 Qiao Dong An-xi 1,950 Qiao Wan Qiao Wan 6,150  Total capacity  12,710


#### 4.2.2. Conventional Power Resource and Load Data

All the other system data including the load and conventional power resource are collected from NCG (all bus nodes are at 330 KV or above; 330 KV are the backbone of the voltage levels in NCG’s transmission grid). The load data is taken from a typical autumn day of 2010 and the units’ cost curves are obtained. The total controllable generation capacity is 67,000 MW and the system peak load is 36,000 MW. The system and network data are 221 generating units including 16 hydro plants in the grid of 126 loads, 152 nodes, and 220 transmission lines. The hydro plant is simplified as an aggregated unit which is always on line during the horizon; each hydro resource has a maximum availability for limited water and minimum output for reservoir amount constraints. Here we take hydro generation as a dispatchable resource but not a deterministic one, because it offers at least some flexibility to offset or “store” some variations in wind generation. The number of units and the capacity for each fuel type are shown in Table2.Table 2
Controllable generation mix.
Type Number of units Capacity (MW) Coal 109 43,052 Hydro 16 15,661 Other fuel 95 7,191 Import 1 1,000The detailed data for thermal units including the minimum up/down time are listed in Table3. The ramping rate of the thermal units is set to be 1/3 of the unit’s installed capacity per hour. Since we assume the hydro plant as an aggregated plant, the minimum up/down time of the hydro one are set to be the whole scheduling period, that is, 24 hours. The spinning reserve is set to be 20% of the forecasted peak load. (Based on current reliability criteria by NCG, the spinning reserve for load variation and unit outage is separately 5% and 15%.)Table 3
Thermal units data.
Type Number of units Minimum up time (hour) Minimum down time (hour) Type 1# 2 1 1 Type 2# 23 1 2 Type 3# 30 2 3 Type 4# 40 3 4 Type 5# 11 3 5 Type 6# 70 4 7 Type 7# 23 5 8 Type 8# 6 24 24


#### 4.2.3. Wind Power Forecast Data

Since the day-ahead wind forecast technique in NCG is still based on point prediction methods, we use a typical day’s wind power point forecast data to create the uncertainty set. The wind sites uncertainty bounds are set to be plus/minus a percentage value from the hourly point forecast data. Although the real uncertainty bound of probabilistic forecast is not a fixed value but a various one which is determined by point forecast value at different time period [14], we assume in this paper that the setting of all the up/down uncertainty bounds in the horizon time is based on uniform level; the reason for that is it can provide much clearer comparison of the dispatch cost variation between different probabilistic forecast levels. The system wind power uncertainty bounds are also set in a similar way, but the uncertainty level is smaller than that for single wind site’s forecast.


### 4.3. The Performance Analysis Results

#### 4.3.1. Computational Efficiency

The program for all the numerical tests is implemented in C++ on Visual Studio 2008 platform [28] and the MILP problem is solved by CPLEX 11.0 [29] on a Windows PC with an 3.40 GHz CPU and 4 GB RAM. The wind uncertainty bounds are set to be ±10% of forecasted value. The convergence tolerance is set as 0.001. The average computational time to solve the robust UC described above is around 50 seconds. A fast computation comes from the one-stage solution structure which avoids the iteration between outer and inner stage (every iteration will solve a MIP outer-stage problem); in Section 3, we obtain the analytical solution of χ-(t), k~1l(t), and γ¯l,t, stated as in (27), (31), and (33). The analytical solutions obtained can be aggregated as linear constraints into the deterministic reformulation without applying decomposition solution methods. This will reduce the enumeration time of the finite portfolios to find the optimal solution of optimization problems in (23), (20), and (15). If the convergence tolerance is reduced to 0.0001, the average computational time will be increased to 300 seconds with an average of 0.045% cost decrease in terms of the worst-case total cost; this means that a reduced convergence tolerance will cause six times more computational effects with only fractional gain in cost reduction. Therefore we set the convergence tolerance in the following subsections as 0.001.


#### 4.3.2. Economic Efficiency

Table4 demonstrates the selective unit commitment results of thermal units: we select eight unit commitment results separately from eight types of thermal units which are divided based on the same way described in Table 3 (in Table 3, the types of thermal units are divided based on their installed capacity (IC); the ICs of eight thermal units are separately smaller than 15 MW, 30 MW, 75 MW, 150 MW, 200 MW, 400 MW, 660 MW, and 1000 MW).Table 4
Selective unit commitment results of thermal units.
Type Unit commitment results Type 1# 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 Type 2# 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 Type 3# 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 Type 4# 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type 5# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 Type 6# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type 7# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type 8# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1The worst-case total cost is calculated under the different wind power uncertainty levels. Since this level is basically determined by the forecast techniques chosen by wind site operators, under same confidence probability, the level can vary for different wind sites. Table5 demonstrates the start-up cost and total cost when wind power uncertainty levels vary from 10% to 35%. With the increasing of the uncertainty bound, our finding is that, for scheduling the system against the worst-case wind realization, more controllable units would be committed and the total cost for worst-case would also increase.Table 5
The start-up costs and total costs in terms of worst case under different uncertainty level.
Uncertainty level Total cost in terms of worst case (100 M RMB) Start-up cost (100 M RMB) 10% 1.915 0.068 15% 1.921 0.068 20% 1.926 0.069 25% 1.931 0.069 30% 1.937 0.069 35% 1.950 0.070To compare the differences in cost between the proposed and conventional method (i.e., we specified a deterministic amount of reserve to cope with wind uncertainty), we introduce the penalty costs for load imbalance and line flow constraints violation under some wind power realization. The penalty price is set to be the same with 10 times the highest units’ supply price. The conventional deterministic UC problem is solved at the expected wind power generation. The comparison experiment is proceeding as follows: we first obtain the UC solution, respectively, for the proposed UC problem and the conventional one; then the ED problem is solved repeatedly for 10 wind generation realizations based on their correlation characteristics. The total dispatch cost is the sum of the start-up plus the average ED cost. The performance is compared in two aspects: the average dispatch cost and the start-up cost and the sensitivity of these costs to different uncertainty levels. Table6 reveals the start-up and the average dispatch cost for the two methods under 10% and 30% uncertainty levels. The penalty cost due to load imbalance and line flow constraints violation by the deterministic method is also shown in Table 6.Table 6
The start-up costs and real dispatch cost for proposed approach and conventional reserve method.
10% Uncertainty level 30% Uncertainty level Proposed approach start-up cost (100 M RMB) 0.068 0.069 Proposed approach real dispatch cost (100 M RMB) 1.837 1.838 Conventional reserve method start-up cost (100 M RMB) 0.068 0.068 Conventional reserve method real dispatch cost (100 M RMB) 1.839 1.846 Conventional reserve method percent of disp. cost for penalty cost 0.2% 1.0%It is founded that, from Table6, under low uncertainty level the proposed solution approach has a nearly identical total dispatch cost with that of the conventional reserve method. This is because, under low wind power uncertainty level, there is relatively low load imbalance and line flow constraint violation. Thus, if the system operator uses the conventional reserve method, the percentage of the total dispatch cost for the penalty cost is small. This result reveals that the proposed solution approach in our paper will not give a more conservative result than that of the conventional reserve method. Furthermore, when the uncertainty level increases, the feasible ED cannot be obtained by the conventional reserve method without purchasing energy outsiders, cutting wind power, or even shedding load. From the 30% uncertainty level case, we can find that although the start-up cost for the proposed approach is higher than that of the conventional reserve method, the latter one will create the penalty cost for acquiring outside power supply to meet the actual system demand realization. Thus, the total dispatch under the conventional reserve method is larger than that of the proposed approach. The economic efficiency of the proposed approach outperforms that of the conventional reserve method under a high wind power uncertainty.


## 4.1. Motivation

In this section, we will illustrate the performance of the proposed solution approach on a large-scale system and discuss its economic efficiency and computational advantage. To do so, we firstly collect a real data set which can represent a power system with large-scale wind power output. Then, we make a design for the whole numerical tests where the test data is realistic and applied to the future facts that the wind power output will be increased as the largest power generation resource within a power system. Finally, the related test results corresponding to the computational and economic efficiency are demonstrated through the computation time, the average dispatch and start-up cost, and the unit commitment results compared with the conventional reserve method in previous publications.


## 4.2. Numerical Tests Design

### 4.2.1. Power System Description

The test dataset is from the power system operated by Northwest China Grid (NCG). NCG has a large-scale wind resource which is mainly in Hexi Corridor within Jiuquan area of Gansu Province, and its energy has been to form into three large aggregated wind farms around this area with estimated wind capacity installation as 5,450 MW, 6,150 MW, and 1,110 MW in 2015. All the wind farms are within a 130 km “belt” surrounded by Qilian and Bei Mountains; thus their wind power outputs generally demonstrate a characteristic of high correlation. The wind penetration level is defined as the ratio of the installed wind capacity to the peak load and in the test cases it is around 35%. This penetration level is large enough to represent the future tendency that the wind power is integrated into power system and can meet the basic engineering assumption we made in Section3. In this paper, our test dataset is based on four representative wind farms and their capacity at the 2015 level. The wind data is collected from the EMS system of NCG: we selected 5-minute wind power output for a whole month and the PCA matrix is acquired for every hour. The locations of the wind sites are presented in Table 1.Table 1
Location and capacity of wind sites (MW).
Site location Bus nodes name Capacity (MW) Chang-ma Yu-men-zhen 1,110 Gan-he-kou Dong An-xi 3,500 Qiao Dong An-xi 1,950 Qiao Wan Qiao Wan 6,150  Total capacity  12,710


### 4.2.2. Conventional Power Resource and Load Data

All the other system data including the load and conventional power resource are collected from NCG (all bus nodes are at 330 KV or above; 330 KV are the backbone of the voltage levels in NCG’s transmission grid). The load data is taken from a typical autumn day of 2010 and the units’ cost curves are obtained. The total controllable generation capacity is 67,000 MW and the system peak load is 36,000 MW. The system and network data are 221 generating units including 16 hydro plants in the grid of 126 loads, 152 nodes, and 220 transmission lines. The hydro plant is simplified as an aggregated unit which is always on line during the horizon; each hydro resource has a maximum availability for limited water and minimum output for reservoir amount constraints. Here we take hydro generation as a dispatchable resource but not a deterministic one, because it offers at least some flexibility to offset or “store” some variations in wind generation. The number of units and the capacity for each fuel type are shown in Table2.Table 2
Controllable generation mix.
Type Number of units Capacity (MW) Coal 109 43,052 Hydro 16 15,661 Other fuel 95 7,191 Import 1 1,000The detailed data for thermal units including the minimum up/down time are listed in Table3. The ramping rate of the thermal units is set to be 1/3 of the unit’s installed capacity per hour. Since we assume the hydro plant as an aggregated plant, the minimum up/down time of the hydro one are set to be the whole scheduling period, that is, 24 hours. The spinning reserve is set to be 20% of the forecasted peak load. (Based on current reliability criteria by NCG, the spinning reserve for load variation and unit outage is separately 5% and 15%.)Table 3
Thermal units data.
Type Number of units Minimum up time (hour) Minimum down time (hour) Type 1# 2 1 1 Type 2# 23 1 2 Type 3# 30 2 3 Type 4# 40 3 4 Type 5# 11 3 5 Type 6# 70 4 7 Type 7# 23 5 8 Type 8# 6 24 24


### 4.2.3. Wind Power Forecast Data

Since the day-ahead wind forecast technique in NCG is still based on point prediction methods, we use a typical day’s wind power point forecast data to create the uncertainty set. The wind sites uncertainty bounds are set to be plus/minus a percentage value from the hourly point forecast data. Although the real uncertainty bound of probabilistic forecast is not a fixed value but a various one which is determined by point forecast value at different time period [14], we assume in this paper that the setting of all the up/down uncertainty bounds in the horizon time is based on uniform level; the reason for that is it can provide much clearer comparison of the dispatch cost variation between different probabilistic forecast levels. The system wind power uncertainty bounds are also set in a similar way, but the uncertainty level is smaller than that for single wind site’s forecast.


## 4.2.1. Power System Description

The test dataset is from the power system operated by Northwest China Grid (NCG). NCG has a large-scale wind resource which is mainly in Hexi Corridor within Jiuquan area of Gansu Province, and its energy has been to form into three large aggregated wind farms around this area with estimated wind capacity installation as 5,450 MW, 6,150 MW, and 1,110 MW in 2015. All the wind farms are within a 130 km “belt” surrounded by Qilian and Bei Mountains; thus their wind power outputs generally demonstrate a characteristic of high correlation. The wind penetration level is defined as the ratio of the installed wind capacity to the peak load and in the test cases it is around 35%. This penetration level is large enough to represent the future tendency that the wind power is integrated into power system and can meet the basic engineering assumption we made in Section3. In this paper, our test dataset is based on four representative wind farms and their capacity at the 2015 level. The wind data is collected from the EMS system of NCG: we selected 5-minute wind power output for a whole month and the PCA matrix is acquired for every hour. The locations of the wind sites are presented in Table 1.Table 1
Location and capacity of wind sites (MW).
Site location Bus nodes name Capacity (MW) Chang-ma Yu-men-zhen 1,110 Gan-he-kou Dong An-xi 3,500 Qiao Dong An-xi 1,950 Qiao Wan Qiao Wan 6,150  Total capacity  12,710


## 4.2.2. Conventional Power Resource and Load Data

All the other system data including the load and conventional power resource are collected from NCG (all bus nodes are at 330 KV or above; 330 KV are the backbone of the voltage levels in NCG’s transmission grid). The load data is taken from a typical autumn day of 2010 and the units’ cost curves are obtained. The total controllable generation capacity is 67,000 MW and the system peak load is 36,000 MW. The system and network data are 221 generating units including 16 hydro plants in the grid of 126 loads, 152 nodes, and 220 transmission lines. The hydro plant is simplified as an aggregated unit which is always on line during the horizon; each hydro resource has a maximum availability for limited water and minimum output for reservoir amount constraints. Here we take hydro generation as a dispatchable resource but not a deterministic one, because it offers at least some flexibility to offset or “store” some variations in wind generation. The number of units and the capacity for each fuel type are shown in Table2.Table 2
Controllable generation mix.
Type Number of units Capacity (MW) Coal 109 43,052 Hydro 16 15,661 Other fuel 95 7,191 Import 1 1,000The detailed data for thermal units including the minimum up/down time are listed in Table3. The ramping rate of the thermal units is set to be 1/3 of the unit’s installed capacity per hour. Since we assume the hydro plant as an aggregated plant, the minimum up/down time of the hydro one are set to be the whole scheduling period, that is, 24 hours. The spinning reserve is set to be 20% of the forecasted peak load. (Based on current reliability criteria by NCG, the spinning reserve for load variation and unit outage is separately 5% and 15%.)Table 3
Thermal units data.
Type Number of units Minimum up time (hour) Minimum down time (hour) Type 1# 2 1 1 Type 2# 23 1 2 Type 3# 30 2 3 Type 4# 40 3 4 Type 5# 11 3 5 Type 6# 70 4 7 Type 7# 23 5 8 Type 8# 6 24 24


## 4.2.3. Wind Power Forecast Data

Since the day-ahead wind forecast technique in NCG is still based on point prediction methods, we use a typical day’s wind power point forecast data to create the uncertainty set. The wind sites uncertainty bounds are set to be plus/minus a percentage value from the hourly point forecast data. Although the real uncertainty bound of probabilistic forecast is not a fixed value but a various one which is determined by point forecast value at different time period [14], we assume in this paper that the setting of all the up/down uncertainty bounds in the horizon time is based on uniform level; the reason for that is it can provide much clearer comparison of the dispatch cost variation between different probabilistic forecast levels. The system wind power uncertainty bounds are also set in a similar way, but the uncertainty level is smaller than that for single wind site’s forecast.


## 4.3. The Performance Analysis Results

### 4.3.1. Computational Efficiency

The program for all the numerical tests is implemented in C++ on Visual Studio 2008 platform [28] and the MILP problem is solved by CPLEX 11.0 [29] on a Windows PC with an 3.40 GHz CPU and 4 GB RAM. The wind uncertainty bounds are set to be ±10% of forecasted value. The convergence tolerance is set as 0.001. The average computational time to solve the robust UC described above is around 50 seconds. A fast computation comes from the one-stage solution structure which avoids the iteration between outer and inner stage (every iteration will solve a MIP outer-stage problem); in Section 3, we obtain the analytical solution of χ-(t), k~1l(t), and γ¯l,t, stated as in (27), (31), and (33). The analytical solutions obtained can be aggregated as linear constraints into the deterministic reformulation without applying decomposition solution methods. This will reduce the enumeration time of the finite portfolios to find the optimal solution of optimization problems in (23), (20), and (15). If the convergence tolerance is reduced to 0.0001, the average computational time will be increased to 300 seconds with an average of 0.045% cost decrease in terms of the worst-case total cost; this means that a reduced convergence tolerance will cause six times more computational effects with only fractional gain in cost reduction. Therefore we set the convergence tolerance in the following subsections as 0.001.


### 4.3.2. Economic Efficiency

Table4 demonstrates the selective unit commitment results of thermal units: we select eight unit commitment results separately from eight types of thermal units which are divided based on the same way described in Table 3 (in Table 3, the types of thermal units are divided based on their installed capacity (IC); the ICs of eight thermal units are separately smaller than 15 MW, 30 MW, 75 MW, 150 MW, 200 MW, 400 MW, 660 MW, and 1000 MW).Table 4
Selective unit commitment results of thermal units.
Type Unit commitment results Type 1# 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 Type 2# 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 Type 3# 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 Type 4# 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type 5# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 Type 6# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type 7# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type 8# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1The worst-case total cost is calculated under the different wind power uncertainty levels. Since this level is basically determined by the forecast techniques chosen by wind site operators, under same confidence probability, the level can vary for different wind sites. Table5 demonstrates the start-up cost and total cost when wind power uncertainty levels vary from 10% to 35%. With the increasing of the uncertainty bound, our finding is that, for scheduling the system against the worst-case wind realization, more controllable units would be committed and the total cost for worst-case would also increase.Table 5
The start-up costs and total costs in terms of worst case under different uncertainty level.
Uncertainty level Total cost in terms of worst case (100 M RMB) Start-up cost (100 M RMB) 10% 1.915 0.068 15% 1.921 0.068 20% 1.926 0.069 25% 1.931 0.069 30% 1.937 0.069 35% 1.950 0.070To compare the differences in cost between the proposed and conventional method (i.e., we specified a deterministic amount of reserve to cope with wind uncertainty), we introduce the penalty costs for load imbalance and line flow constraints violation under some wind power realization. The penalty price is set to be the same with 10 times the highest units’ supply price. The conventional deterministic UC problem is solved at the expected wind power generation. The comparison experiment is proceeding as follows: we first obtain the UC solution, respectively, for the proposed UC problem and the conventional one; then the ED problem is solved repeatedly for 10 wind generation realizations based on their correlation characteristics. The total dispatch cost is the sum of the start-up plus the average ED cost. The performance is compared in two aspects: the average dispatch cost and the start-up cost and the sensitivity of these costs to different uncertainty levels. Table6 reveals the start-up and the average dispatch cost for the two methods under 10% and 30% uncertainty levels. The penalty cost due to load imbalance and line flow constraints violation by the deterministic method is also shown in Table 6.Table 6
The start-up costs and real dispatch cost for proposed approach and conventional reserve method.
10% Uncertainty level 30% Uncertainty level Proposed approach start-up cost (100 M RMB) 0.068 0.069 Proposed approach real dispatch cost (100 M RMB) 1.837 1.838 Conventional reserve method start-up cost (100 M RMB) 0.068 0.068 Conventional reserve method real dispatch cost (100 M RMB) 1.839 1.846 Conventional reserve method percent of disp. cost for penalty cost 0.2% 1.0%It is founded that, from Table6, under low uncertainty level the proposed solution approach has a nearly identical total dispatch cost with that of the conventional reserve method. This is because, under low wind power uncertainty level, there is relatively low load imbalance and line flow constraint violation. Thus, if the system operator uses the conventional reserve method, the percentage of the total dispatch cost for the penalty cost is small. This result reveals that the proposed solution approach in our paper will not give a more conservative result than that of the conventional reserve method. Furthermore, when the uncertainty level increases, the feasible ED cannot be obtained by the conventional reserve method without purchasing energy outsiders, cutting wind power, or even shedding load. From the 30% uncertainty level case, we can find that although the start-up cost for the proposed approach is higher than that of the conventional reserve method, the latter one will create the penalty cost for acquiring outside power supply to meet the actual system demand realization. Thus, the total dispatch under the conventional reserve method is larger than that of the proposed approach. The economic efficiency of the proposed approach outperforms that of the conventional reserve method under a high wind power uncertainty.


## 4.3.1. Computational Efficiency

The program for all the numerical tests is implemented in C++ on Visual Studio 2008 platform [28] and the MILP problem is solved by CPLEX 11.0 [29] on a Windows PC with an 3.40 GHz CPU and 4 GB RAM. The wind uncertainty bounds are set to be ±10% of forecasted value. The convergence tolerance is set as 0.001. The average computational time to solve the robust UC described above is around 50 seconds. A fast computation comes from the one-stage solution structure which avoids the iteration between outer and inner stage (every iteration will solve a MIP outer-stage problem); in Section 3, we obtain the analytical solution of χ-(t), k~1l(t), and γ¯l,t, stated as in (27), (31), and (33). The analytical solutions obtained can be aggregated as linear constraints into the deterministic reformulation without applying decomposition solution methods. This will reduce the enumeration time of the finite portfolios to find the optimal solution of optimization problems in (23), (20), and (15). If the convergence tolerance is reduced to 0.0001, the average computational time will be increased to 300 seconds with an average of 0.045% cost decrease in terms of the worst-case total cost; this means that a reduced convergence tolerance will cause six times more computational effects with only fractional gain in cost reduction. Therefore we set the convergence tolerance in the following subsections as 0.001.


## 4.3.2. Economic Efficiency

Table4 demonstrates the selective unit commitment results of thermal units: we select eight unit commitment results separately from eight types of thermal units which are divided based on the same way described in Table 3 (in Table 3, the types of thermal units are divided based on their installed capacity (IC); the ICs of eight thermal units are separately smaller than 15 MW, 30 MW, 75 MW, 150 MW, 200 MW, 400 MW, 660 MW, and 1000 MW).Table 4
Selective unit commitment results of thermal units.
Type Unit commitment results Type 1# 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 Type 2# 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 Type 3# 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 Type 4# 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type 5# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 Type 6# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type 7# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type 8# 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1The worst-case total cost is calculated under the different wind power uncertainty levels. Since this level is basically determined by the forecast techniques chosen by wind site operators, under same confidence probability, the level can vary for different wind sites. Table5 demonstrates the start-up cost and total cost when wind power uncertainty levels vary from 10% to 35%. With the increasing of the uncertainty bound, our finding is that, for scheduling the system against the worst-case wind realization, more controllable units would be committed and the total cost for worst-case would also increase.Table 5
The start-up costs and total costs in terms of worst case under different uncertainty level.
Uncertainty level Total cost in terms of worst case (100 M RMB) Start-up cost (100 M RMB) 10% 1.915 0.068 15% 1.921 0.068 20% 1.926 0.069 25% 1.931 0.069 30% 1.937 0.069 35% 1.950 0.070To compare the differences in cost between the proposed and conventional method (i.e., we specified a deterministic amount of reserve to cope with wind uncertainty), we introduce the penalty costs for load imbalance and line flow constraints violation under some wind power realization. The penalty price is set to be the same with 10 times the highest units’ supply price. The conventional deterministic UC problem is solved at the expected wind power generation. The comparison experiment is proceeding as follows: we first obtain the UC solution, respectively, for the proposed UC problem and the conventional one; then the ED problem is solved repeatedly for 10 wind generation realizations based on their correlation characteristics. The total dispatch cost is the sum of the start-up plus the average ED cost. The performance is compared in two aspects: the average dispatch cost and the start-up cost and the sensitivity of these costs to different uncertainty levels. Table6 reveals the start-up and the average dispatch cost for the two methods under 10% and 30% uncertainty levels. The penalty cost due to load imbalance and line flow constraints violation by the deterministic method is also shown in Table 6.Table 6
The start-up costs and real dispatch cost for proposed approach and conventional reserve method.
10% Uncertainty level 30% Uncertainty level Proposed approach start-up cost (100 M RMB) 0.068 0.069 Proposed approach real dispatch cost (100 M RMB) 1.837 1.838 Conventional reserve method start-up cost (100 M RMB) 0.068 0.068 Conventional reserve method real dispatch cost (100 M RMB) 1.839 1.846 Conventional reserve method percent of disp. cost for penalty cost 0.2% 1.0%It is founded that, from Table6, under low uncertainty level the proposed solution approach has a nearly identical total dispatch cost with that of the conventional reserve method. This is because, under low wind power uncertainty level, there is relatively low load imbalance and line flow constraint violation. Thus, if the system operator uses the conventional reserve method, the percentage of the total dispatch cost for the penalty cost is small. This result reveals that the proposed solution approach in our paper will not give a more conservative result than that of the conventional reserve method. Furthermore, when the uncertainty level increases, the feasible ED cannot be obtained by the conventional reserve method without purchasing energy outsiders, cutting wind power, or even shedding load. From the 30% uncertainty level case, we can find that although the start-up cost for the proposed approach is higher than that of the conventional reserve method, the latter one will create the penalty cost for acquiring outside power supply to meet the actual system demand realization. Thus, the total dispatch under the conventional reserve method is larger than that of the proposed approach. The economic efficiency of the proposed approach outperforms that of the conventional reserve method under a high wind power uncertainty.


## 5. Conclusions

We cast an UC problem for a system integrated with wind resource in this paper. Given a nodal wind injection uncertainty set and the associated correlated information, the system operator can obtain a unit commitment solution which is robust against all uncertain wind power injection realizations. We propose a computationally tractable solution approach to tackle this two-stage optimization problem with mixed integer and continuous basics, aiming at fitting this framework into a real power system operation. The numerical test backs the merit of our problem statement and solution approach from economic efficiency, computational advantage, and reliability.



---

*Source: 5208290-2017-02-16.xml*