Abstract
The stochastic nature of demand and wind generation has a considerable effect on solving the scheduling problem of a modern power system. Network constraints such as power flow equations and transmission capacities also need to be considered for a comprehensive approach to model renewable energy integration and analyze generation system flexibility. Firstly, this paper accounts for the stochastic inputs in such a way that the uncertainties are modeled as normally distributed forecast errors. The forecast errors are then superimposed on the outputs of load and wind forecasting tools. Secondly, it efficiently models the network constraints and tests an iterative algorithm and a piecewise linear approximation for representing transmission losses in mixed integer linear programming (MILP). It also integrates load shedding according to priority factors set by the system operator. Moreover, the different interactions among stochastic programming, network constraints, and prioritized load shedding are thoroughly investigated in the paper. The stochastic model is tested on a power system adopted from Jeju Island, South Korea. Results demonstrate the impact of wind speed variability and network constraints on the flexibility of the generation system. Further analysis shows the effect of loss modeling approaches on total cost, accuracy, computational time, and memory requirement.
1. Introduction
The growing concern over the increasing energy demand and the escalating price of fossil fuels have prompted efforts to develop viable alternative energy sources. Recently, the high level of penetration and integration of dispersed generation from renewable sources such as wind and solar has significantly invoked a critical issue in planning and operation of power systems. Wind energy is an intermittent resource that fluctuates depending on the wind speed and imposes difficulties in terms of both operation and planning [1, 2]. Although, wind energy output can be predicted with some degree of accuracy using a forecasting tool, there is always some uncertainty about future wind output since the weather systems are not perfectly predictable [3].
For more wind power utilization and improved demand satisfaction, it is important that the conventional fuel driven power plants need to operate with flexibility to accommodate the load and wind variability. Besides, load flow restrictions and transmission line capacities limit the extent of available wind energy that can be utilized into electrical power systems. Thus, the optimal scheduling has to be determined in order to integrate much of the wind energy without significant curtailment and at the same time respecting the network constraints. Moreover, it is essential to have precise modelling of the power loss to increase the profitability of the system operator and enhance planning schedules. A computationally fast and efficient algorithm needs to be employed for a better loss estimation. Considering all the operational and network constraints, if the electrical demand is higher than the available conventional and renewable generation at a certain time instant in the scheduling period, the planning/operational optimization tool has to deploy a load shedding strategy according to the load priority factors [4] set by the system operator. This guarantees the load and generation balance with the most cost-effective option.
For solving power system scheduling problem, deterministic and stochastic approaches are proposed by different researchers. Stochastic models incorporate a number of scenarios for wind power and load uncertainty; however, deterministic approaches assume only a single scenario for expected wind power and load. Stochastic optimization has advantage over deterministic approaches in such a way that it can help system operators to efficiently handle the uncertainty and variability from wind power and electricity demand. Among the deterministic approaches proposed in literature, Carrión and Arroyo [5] developed computationally efficient scheduling model for the unit commitment problem of thermal units using MILP. The model requires fewer binary variables and includes a precise representation of time dependent start-up cost and intertemporal constraints such as ramping limits and minimum-up/minimum-down times. The approach described in [6] uses a detailed formulation to efficiently model the power trajectories followed by a thermal unit during the start-up and shut-down processes. Large scale wind systems are included in the scheduling model proposed in [7], and the impact of wind integration in terms of cost, reliability, and environment has been studied. Furthermore, instead of MILP, other optimization problem solving approaches such as a modified subgradient (MSG) method combined with simulated annealing (SA) algorithm and a feasible modified subgradient (F-MSG) method have been proposed to solve a unit commitment problem in [8] and [9], respectively. It has been mentioned that both of these approaches resulted in a better quality of the optimal solutions. However, all the abovementioned deterministic models do not consider the stochastic behavior of load and wind power, and grid constraints have been neglected. Most importantly, for a large scale wind integration to electrical networks, possible bottlenecks and extra congestion to the transmission lines need to be determined.
Recently, García-González et al. [10] and Khodayar and Shahidehpour [11] investigated the combined optimization of a wind farm and a pumped-storage facility from the point of view of a generation company in a market environment considering price and wind uncertainties. However, these studies lack fossil fuel driven power plants such as thermal units and diesel generators in their stochastic optimization model. Similarly, the uncertain price of selling electricity to the price-taker is considered in [12]. Nevertheless, the uncertainty for electricity demand and wind generation has not been taken into account. Meibom et al. [13] established a stochastic optimization model to study the impact of large scale wind integration. Stochastic operation planning approach compared to deterministic worst-case scenario planning is discussed in [14]. Correia et al. [15] developed a scheduling method for small scale isolated power system containing wind and other dispatchable resources. In addition to load and wind uncertainties, solar power and generator outage uncertainties are addressed in [16]. However, the impact of grid constraints and transmission losses has not been assessed in the abovementioned stochastic models.
Some of optimization models formulated as either deterministic or stochastic use a DC power flow for modelling network constraints and transmission losses. In [17–20], the electrical network is represented by a DC power flow model. While line flow limits are taken into account, transmission losses are not considered here. In [20], the scheduling problem is solved in two steps. In the first step, the transmission constraints are neglected, and only one extra constraint representing the load-generation balance is added to the existing operational constraints. Then, the optimization problem is solved assuming that all the units are connected to a single bus. Afterwards, the resulting line flows are checked if the transmission line constraints are satisfied. In the second step, if the presence of line flow violations is detected in the first step, then the corresponding network constraints are incorporated and the system is solved again. In [21], transmission line constraints are neglected, and transmission losses are treated as extra demand. Similarly, network constraints and transmission losses are modeled using DC power flow and loss penalty factors in [22, 23]. Other approaches presented in [24, 25] include transmission losses as a quadratic function of bus power injections, and the scheduling problem is solved using quadratic programming. Further details of the aforementioned methods for modelling transmission losses using DC power flow are explained in Section 3.2. AC power flow based approaches in [26, 27] incorporate nonlinear power flow equations, and a more realistic optimal solution can be found using nonlinear system optimization solvers based on algorithms such as Lagrangian relaxation (LR) and genetic algorithm (GA). Nevertheless, the solution delivery time increases drastically while solving such type of stochastic scheduling problems due to large number of decision variables and constraints in every scenario and time.
As improper loss estimations may have a significant effect on the accuracy of optimal solutions and nonlinear methods can increase the computation time, in this paper the electrical network is modelled using DC power flow with precise representation of transmission losses. Two different loss modeling approaches, namely, iterative method and piecewise linear approximation, are explained in detail later in Sections 3.2.1 and 3.2.2. The former has a merit in terms of solution delivery time and the latter is better in terms of reduction in total operational cost.
Moreover, in this work, the scheduling problem is formulated in a comprehensive stochastic optimization model. The inputs are stochastic load, stochastic wind generation, thermal unit characteristics, and network parameters. Nearly all aspects in terms of operational constraints are included with the objective of minimizing the total fuel cost consumed by conventional thermal power plants. The main contributions of this paper are as follows:(i)It explicitly considers the grid constraints and tests two different approaches to model the transmission losses. The impacts of line flow constraints are neglected in [10, 11, 13–16, 21] and transmission losses are either omitted or poorly estimated in [10–23].(ii)It introduces slack variables for solution feasibility in similar manner to [13, 14]. However, in this work, by assigning different priority factors for expected lost load at each load bus, a prioritized load shedding scheme is integrated in the proposed optimization model.(iii)It accounts for the stochastic nature of electricity demand and wind generation and conducts a comprehensive analysis to study the possible interactions of grid constraints, transmission losses, and prioritized load shedding in stochastic programming. For instance, it analyses the relevancy of modeling transmission losses in a generation system subjected to wind and load uncertainty and the effect of prioritized load shedding on scheduling performance compared to that of considering transmission losses. However, this has not been discussed in either stochastic or deterministic approaches mentioned in the literature.
The rest of this paper is organized as follows. Section 2 describes the scenario generation of forecast error trajectories for load and wind power. Section 3 follows by modelling the grid constraints using DC power flow and explains different methods to model transmission losses in a scheduling problem. Section 4 deals with the network-constrained stochastic mixed integer scheduling model. Section 5 discusses the simulation and results. Finally, conclusions are given in Section 6.
2. Load and Wind Trajectory Generation
When intermittent electricity demand and wind generation are integrated into the scheduling problem, the forecast errors associated with these two inputs will significantly influence the system operator’s scheduling process. In order to account for the uncertainty arising from those forecast errors, the scheduling problem incorporates a scenario generation tool to produce several normally distributed forecast error trajectories with autocorrelation.
Considering the usual phenomena that uncertainty gets worse for increasing time horizon, the trajectories for either load or wind are generated with linearly increasing variance. It is necessary to consider multiple scenarios of wind and load forecast errors to cover the uncertainty range for the statistical parameters of different forecast trajectories such as mean values and autocorrelations. It is expected that increasing the number of scenarios should comply with a better scheduling performance. An efficient schedule should result in minimum lost load, lost reserve, and lost wind.
However, since simulating the optimization problem with a large number of scenarios increases the computational burden, a scenario reduction method using -means [28] is used. The algorithm begins by choosing arbitrary centres (i.e., predefined number of scenarios) from the forecast error trajectories. Then, Euclidean distance between the centres and each forecast trajectory is computed. Each trajectory is then assigned to the nearest centre, and each centre is recomputed as the average of all trajectories assigned to it. Then, the sum of probabilities of these trajectories is allocated to their respective centres. These three procedures (grouping, centre calculation, and probability allocation) are repeated until the distance between the centres and the other forecast trajectories is less than the predefined tolerance.
Finally, each centre that corresponds to a group of forecast error trajectories is superimposed on a given time series load and wind data. Alternatively, these representative trajectories can also be added to the outputs of load and wind forecasting tools.
3. Linearized Power Flow for Stochastic Optimization
In most of optimization problems, the network constraints are not explicitly considered. In such cases, the processing speed is fast and the problem complexity is relatively simple to solve. However, in most cases, solutions are not accurate due to poor loss approximation, and transmission lines are assumed to be operating regardless of the power level that has to be transmitted between buses. Considering its noniterative process and its reliable and distinctive solution, in this work, the network constraints are modelled using DC power flow. A simplified system formulation that is comparatively easy to obtain makes DC power flow applicable to stochastic optimization due to the involvement of volumes of computing with a number of scenarios.
3.1. DC Power Flow Constraints
Given that the admittance for all branches is provided, the DC power flow equations for net nodal power injections are expressed as And the line MW flows are given by A flat magnitude of bus voltages and high inductance to resistance ratio of transmission lines are assumed in the derivation of DC power flow.
3.2. Loss Inclusion Methods in Optimization Model
When transmitting electricity from a certain power plant through high voltage transmission lines or generator to local distribution network areas, some power is lost. Accounting for transmission losses in optimal scheduling problems is necessary to obtain improved generator power dispatch especially in some complex transmission networks with potential network congestion and notable power loss. Specific calculations are required to determine the level of the power loss.
There have been various trends for representing transmission losses in optimal scheduling problems. The system loss is preestimated prior to optimization and added to the electricity load as a constant parameter in [21]. The load-generation balance is given by
In [13], the transmission losses are represented by a fixed percentage () of the power flow via the transmission lines connected to a certain demand point (load bus) as stated in
An alternative approach proposed in [22, 23] uses loss penalty factors () in order to account for the incremental change in system losses as a result of the change in generator outputs. The transmission losses are handled in such a way that loss penalty factors are assigned to each bus directly connected to a certain group of power units. Then, the generation outputs are divided by the corresponding loss penalty factors. The expression in (5) imposes extra generation compared to the same system where no loss is considered: The -coefficients method is applied to calculate the transmission losses in [24, 25]. The total loss is assumed to be a quadratic function of the aggregate generation at each bus and expressed by (6). It is then superimposed on the system demand as shown on the right-hand side of the load-generation balance equation stated in (7). The transmission loss coefficients () for the given network need to be determined prior to optimization in order to include the loss expression inside the scheduling problem. The total loss and load-generation balance equations are given byNote that in this case is a variable, but, in the loss estimation method [21] explained above, is a known parameter.
However, modelling transmission losses using offline preestimation [13, 21] or loss penalty factors [22, 23] is not accurate and may result in an unexpected value of the generation dispatches which might incur an extra expense in the total operating cost of the system. The approach in [24, 25] involves a nonlinear expression of transmissions losses that invokes intensive computational resources when solving stochastic scheduling problems considering these nonlinear equalities. Based on DC power flow, in order to represent transmissions loss in stochastic mixed integer optimization with reasonable computational time and accuracy, an iterative algorithm and a piecewise linear approximation method are applied in this work.
Transmission losses are nonlinear functions of the generator dispatches. Applying the well-known assumptions of DC power flow, line losses are approximated as a quadratic function of line flows. However, the line flows are linear functions of phase angle differences between the sending and receiving buses. Then, the power loss on a certain line can be approximately computed as Based on (8), two different approaches for incorporating transmission losses inside the stochastic scheduling problem are discussed below.
3.2.1. Iterative Algorithm
This method encapsulates losses in the optimization problem by computing them externally using iterative algorithm and passing the loss values as parameters to the MILP solver. The procedure is summarized below.
Step 1 (input data). Load and wind data from scenario generation tool, network data, and thermal generator characteristics are fed to the optimizer.
Step 2 (initialization). The system losses are set to zero for all lines in all time instants and scenarios:
Step 3 (running the optimizer). The optimal scheduling problem is solved using commercially available software products such as CPLEX.
Step 4 (line loss calculation). The power loss in each line is calculated using (1)–(3). Then, the computed loss is divided equally and lumped as a load at sending and receiving end buses of each transmission line.
Step 5 (stopping criteria). The relative error percentage is defined in (10) as the difference between the line losses at the current iteration, , and the previous one, . According to (11), if the relative error in every dimension (, and ) is less than or equal to the specified error tolerance, the iteration halts. Otherwise, the procedure returns to Step 3. The expressions for the relative error are given by
3.2.2. Piecewise Linear Approximation Using MILP
The quadratic expression in (8) can be modelled using a piecewise linear function composed of straight line segments. For real applications, the piecewise linear function is indistinguishable from the nonlinear (quadratic) model if enough segments are used.
In addition to the variables corresponding to existing operational constraints of the stochastic scheduling problem, for every line, scenario, and time instant, binary and segment variables are required to model the piecewise linear loss function in MILP. Prior to optimization, the horizontal axis is discretized, and the break points () are identified as shown in Figure 1. Then, the binary and segment variables together with the breakpoints fully describe the linear constraints representing the piecewise linear function.
The binary variable denotes the status (active/inactive) of a segment, and the segment variable is associated with the fractional value of a segment in both coordinates. If the optimizer requires computing the loss corresponding to a certain value of the phase angle difference () that falls in between two break points, and , then the binary variable assigned to that particular segment connecting () and () is switched on. Then, the segment variable takes a real number between zero and one such that the value of the phase angle difference is precisely represented in between end points of that active linear segment.
The mathematical formulation of the equalities and inequalities associated with the linear constraints is given as follows:For a specified scenario and particular time instant, the power loss and phase angle difference for line- are computed in (12) and (13), respectively. Equation (14) ensures only one binary variable to be switched on whenever a particular piecewise segment is required, while, in (15), the fraction of the active segment is defined in between zero and one. For instance, consider a piecewise linear approximation of quadratic losses as shown in Figure 1 given that the values of and are and , respectively. Suppose the optimizer decides that the phase difference between the extremes of a certain transmission line is 0.05 , then only the second segment should be selected, and the value of the binary variable () has to be while and remain at . Once the segment is specified, the value of the segment variable is calculated as .
The main advantage of describing the piecewise linear function using (12)–(15) completely defined by binary and segment variables is that the MILP solver is always able to find a global optimum [29]. Furthermore, the piecewise linear approximation method has advantage over the iterative method because the expressions are directly embedded in the stochastic scheduling model, and it allows the MILP optimizer to consider all the constraints together at the same time. This advantage is reflected on the reduction of the total operational cost compared to the iterative method. The comparison of these loss modelling approaches by means of a numerical example is discussed in Section 5.5.
Some methods in literature use the decomposition strategy [20] to consider the line flow constraints after violations are detected. The process is iterated until all line flows in each scenario and time fall within the capacity of the corresponding transmission lines. However, since transmission loss always exits in every transmission line, the piecewise linear loss constraints (12)–(15) can not be decoupled and treated as line flow constraints in decomposition method. On the contrary, a certain transmission line with larger capacity could handle every power level that the system requires to flow via that particular line. Ignoring the line flow constraint in this case does not have effect on the optimal solution. However, omitting the line loss resulting from the same transmission line has an impact on the accuracy of the optimal solution.
Incorporating the line flow constraints using the decomposition method in a lossless stochastic scheduling problem with small number of scenarios might have some advantage in reducing the computational burden. However, when piecewise linear loss constraints are included in the stochastic scheduling problem with larger number of scenarios, the iterative procedure drastically increases the computational time since the total number of constraints is sufficiently large enough compared to the number of line flow constraints. For instance, the power system described in Section 5.1 has fourteen () transmission lines. When the stochastic scheduling problem is formulated with two scenarios () and twelve-hour period () without considering transmission losses, the number of constraints excluding the line flow constraints is . With sixteen scenarios, the number of constraints excluding the line flow constraints is grown to . The number of line flow constraints is calculated as , that is, for two scenarios and for sixteen scenarios. The proportion of the line flow constraints for the two scenario case is (=336/(15522 + 336)) of the total number of constraints. When piecewise linear losses are included in the model with ten-segment () representation, the number of constraints corresponding to binary and segment variables is computed as , that is, for sixteen scenarios. The proportion of the line flow constraints considering the piecewise losses with sixteen scenarios is (=2688/(256944 + 107520 + 2688)) of the total number of constraints. This shows that when the decomposition method is used, the reduction in the number of constraints is insignificant when a sufficient number of scenarios are used to account for uncertainty, and transmission losses are modelled in piecewise linear functions. Rather, it costs higher computational time since the problem with the rest of all operational constraints has to be resolved every time a line flow violation occurred.
In stochastic scheduling problem, in order to obtain deterministic solution for a certain scheduling period and account for updated load and wind forecasts, a repeated rolling planning process is required. In rolling planning window framework, when line flow constraints, piecewise linear loss constraints, and several other operational constraints are simultaneously considered without any decomposition as applied in this paper, the solver needs only to run twenty-four times for a planning window with one-hour nonrecourse period and a day-ahead scheduling horizon. However, for the decomposition method, suppose that a line violation does not occur after three iterations in average for each planning window and then the solver has to run (=3 24) times to deliver a final complete solution. This heavily costs extensive processing time.
4. Stochastic Scheduling Model
The scheduling model involves the determination of on/off status and power output levels of all conventional thermal plants and diesel generators in each time instant and start-up/shut-down times of each units. It also aims at resolving the magnitude of unserved load, wasted wind, and lost reserve in each time instant. The main objective is to minimize the running and penalty costs for a designated scheduling period specified by the system operator subjected to operational constraints.
The stochastic optimization problem is formulated in such a way that the uncertainty for some part of the problem is explicitly addressed using random variables. The uncertain parameters are discretized to a finite set of scenarios. The size of the problem grows larger as the number of scenarios increases. In two-stage type of stochastic problems, the decision variables are classified as nonrecourse and recourse variables. The nonrecourse variables are associated with the nonrecourse period and here-and-now first-stage decisions, and the recourse variables are associated with the recourse period and wait-and-see second-stage decisions. The decisions for the nonrecourse variables are deterministic in initial time instants of the planning window; however, decisions are not taken for the rest of recourse variables until they reside in a nonrecourse period of the following new planning window. Nevertheless, decomposing the problem for the two-stage approach causes difficulty in the formulation of the scheduling model. In this paper, the problem is solved in a single stage approach in such a way that the nonrecourse variables are imposed to exist only in a single scenario while the recourse variables are free to appear in different scenarios. For instance, for a planning window with a length of twenty-four hours, if a nonrecourse period is set in the first four initial time instants, then the nonrecourse variables take a unique value in to for the first planning window ( to ), while the value of recourse variables varies with scenarios in the recourse period, to . In the subsequent planning window ( to ), the nonrecourse period corresponds to time instants to , while the recourse period corresponds to time instants to . In such a way, a deterministic value can be obtained for all decision variables in twenty-four time instants since the they fall in a nonrecourse period of six consecutive planning windows in a rolling planning window framework.
4.1. Objective and Constraints
The objective function given in (16) includes the fuel consumption cost, shut-down cost, start-up cost, penalty cost for not fulfilling the power demand in each bus (lost load), penalty cost for not using the available wind power from each bus connected with a wind farm (lost wind), and penalty cost for not securing some portion of spinning reserve (lost reserve). The large percentage of the operating cost is associated with the fuel consumption of online thermal units which is normally considered as a quadratic function of the power generated from each unit. The piecewise linear approximation of this nonlinear cost is similar to the approach used in linearization of quadratic losses, and it is included in (16) at the end of the objective function. In (17), the power output from each online unit is expressed in terms of variables representing the piecewise linear cost segment and fraction of piecewise linear cost segment. Constraint (18) selects only one linear cost segment to be active, while in (19) the variable corresponding to fraction of piecewise linear cost segment takes a value between zero and one for active linear cost segment and zero for any inactive linear cost segment.
The equality expression in (20) assures a unique scenario and deterministic solution in nonrecourse period, and the nonrecourse variables remain the same for all scenarios in the recourse period.
The network constraints are subdivided into two options. Constraint (21) is formulated for the energy balance in iterative approach. is a parameter calculated from external algorithm and passed to the optimizer. The load-generation balance includes power generation from thermal units and wind farms, transmission loss, wasted wind, and lost demand. The alternative approach for modelling network constraints using piecewise linear approximation of the line losses is stated in (22)–(27). The formulation of modelling these losses has been already explained in detail in Section 3.2.2. In (26), the line losses are distributed equally to the sending and receiving buses. Similar to (21), the load-generation balance for the piecewise approach is given in (27). The accuracy and computational efficiency of iterative and piecewise methods are analysed in Section 5.5. Line power flows are calculated in (28). The capacity of transmission lines and bus angle boundaries are provided in (29) and (30).
Constraint (31) denotes the spinning reserve in terms of the difference between the current power output and the maximum capacity of online units. In (32), the minimum spinning reserve is defined as the sum of the largest unit and some fraction of the total load. The minimum number of units that need to be online is set by (36). Constraints (37)–(39) present the minimum and maximum power ratings of the generators.
Ramp-up/ramp-down limits critically determine the rate at which a thermal unit can increase or decrease its power generation. At a certain period, the unit may lose flexibility to deliver power within its lower/upper rating; rather it is constrained to the minimum/maximum available power output () at that time instant. The power output is constrained by ramp-up and start-up ramp rates (40), as well as by shut-down ramp rates (41) and ramp-down limits (42). For instance, suppose an online unit has a lower rating of 60 MW, upper rating of 200 MW, a ramp-up rate of 60 MW/hr, and ramp-down rate of 50 MW/hr. If the unit is generating 120 MW at current period and intended to increase generation in the following period, the maximum power that it can deliver at is only limited to 180 MW, although the upper rating of the unit is 200 MW. In a similar fashion, if the unit is intended to decrease generation in the subsequent period, the minimum possible generation at is 70 MW. However, the lower power rating of the unit is 60 MW.
A restriction on the minimum-up time of a unit is imposed in (43)–(45). The number of periods that unit needs to be initially online is computed using the formula [5]. The minimum-up time basically reflects the need to minimize thermal stresses in the generator, which can get worse if the units are forced to switch off instantly [30]. In addition, some thermal plants require to take some time to be fully operational. The minimum-down time constraints are formulated in (46)–(48). The number of periods that unit needs to be initially offline is computed using the formula [5]. The minimum-down time ensures adequate time for the required maintenance after a unit has been shut down [31]. The start-up and shut-down costs are computed in (49)-(50). These expenses are largely associated with labor and maintenance [32].
4.2. Slack Variables and Load Shedding
Slack variables are required to secure solution feasibility. In the scheduling problem, lost load (), lost reserve (), and lost wind () are defined as slack variables with a restriction to be greater than or equal to zero. In addition, in (33) to (35), these variables are limited to be lower than the corresponding load, reserve, or wind power at each time instant and scenario. The penalty associated with each slack variable is added to the running cost of thermal units in the objective function:Network constraints-iterative approach (21): Network constraints-piecewise linear approximation method (22)–(27):
Using the slack variables, prioritized load shedding is integrated in the scheduling problem. The strategy attributes different priority factors () to the load unit at each bus according to its importance determined by the system operator. This forces the optimizer to decide the location of the load to be shed and the level of demand that is forced to be unserved at the corresponding bus.
5. Simulation and Numerical Results
Simulation is carried out for a modified power system of Jeju Island [33], South Korea. The effect of network constraints on the integration of wind energy and load shedding strategy is analysed. In addition, for further installation of more wind turbines, the power system flexibility is studied. The possibility of standalone system independent of the main grid is also tested. Moreover, the loss modelling techniques discussed in Section 3.2 are compared in terms of the total cost, computational time, and memory requirement.
Prior to simulations, the length of planning window and the number of scenarios are determined considering the computational time and scheduling performance. A wide planning window usually results in better scheduling performance measured in terms of slack variables. However, the introduction of network constraints and transmission losses makes the simulation slower. Considering this, a planning window of twelve hours (half a day) with discretization of one hour is found to be the best option in terms of speed and accuracy.
The simulation comprises three software packages. MATLAB [34] is used for writing algorithms for loss computation and scenario generation tool. GPLK [35] is utilized for describing the linear mathematical programming model. CPLEX [36] is used to solve the model. The simulation is performed on a PC with Intel Core Duo (3 GHz) processor with 4 GB of RAM.
5.1. Data Source
The modified Jeju Island power system shown in Figure 2 consists of three wind farms (WF) constructed at three different locations: Hangwon (HWN), Sungsan (SSN), and Halim (HLM) currently generating 50, 30, and 20 MW, respectively. The grid from the mainland supplies 300 MW of power. Besides wind power, the system also includes one combined cycle (HLM-CC) plant, two thermal (TP) power plants (NMJ-TP and JJU-TP), and another two units composed of diesel generators (JJU-DP and NMJ-DP). The ramp-up/ramp-down limits of a unit are assumed to be equal to its lower power rating and all the units are considered to be initially online, generating their minimum level of power output. The minimum-up time and minimum-down time of units are set for one hour and two hours for smaller and larger units, respectively. Additional parameters of the thermal units for representing the quadratic cost curve and transmission line data used for developing the network constraints are presented in Tables 1 and 2, respectively. Note that the transmission line parameters are given in per unit (pu). The values are calculated using voltage and power bases of 154 kV and 100 MVA, respectively.
5.2. Load and Wind Scenarios
The wind speed timely variation is synthetically established based on specific Weibull parameters [37] estimated using the annual mean speed and standard deviation at the place where the wind farm is located. Since considering different scenarios for wind speed at every location is practically infeasible due to large number of scenarios, the sum of wind power outputs is passed to the scenario generation tool and the forecast errors are superimposed into this aggregate sum. Wind speed trajectories for three scenarios are shown in Figure 3. In the scheduling model, when the network constraints are formulated, the stochastic aggregate wind power in each time instant and scenario is distributed to the respective wind locations (buses) according to their rated power. The trajectory generation for load (Figure 4) is constructed with a growing forecast error variance of 10% after a day (24 hours).
5.3. Scheduling Outputs
The simulation is run for two months with lowest aggregate and highest aggregate demand. The real system performance in terms of lost wind, load, and reserve lies in between these two extremes. As explained in the previous section, the scheduling parameters (such as size of planning window and time discretization) are rigorously determined using a number of simulations so that the loss of accuracy due to these parameters can be minimized.
Figures 5–11 represent the scheduling outputs for a typical day of the month with highest aggregate demand. Figures 5–8 show the scheduling outputs resulted in simulating the scheduling model with nine scenarios (i.e., three scenarios for load and three scenarios for wind) using a single twelve-hour-long planning window and a one-hour-long nonrecourse period. It can be identified from the figures that, the solution for the first period is always deterministic whereas it is nondistinctive and scenario based for the rest of the planning period (i.e., recourse period), to .
However, in order to obtain a deterministic solution for the entire scheduling period ( to ), and to account for updated load and wind forecasts, a rolling window framework is used. In rolling window, new scenarios are generated based on the new load and wind forecasts, and the dispatch levels of thermal units () are reevaluated in each time period moving forward to the next time period. The final simulation outputs are denoted by the values of the nonrecourse variables in the first period of every rolling planning window. The final scheduling of thermal generation, line MW flows, and voltage angles are shown in Figures 9–11.
5.4. Impact of Network Constraints on Wind Integration
The effect of network constraints on optimal scheduling for installing extra wind power on different locations in the power system and the generation system flexibility for accommodating the wind power variability is studied in the following two conditions.
(1) The system is tested with additional 50 MW of wind power (i.e., a total of 150 MW) installed at HWN, SSN, or HLM. Therefore, three cases are stated as follows.
Case 1.
Consider the following subcases.
Case 1.1. Additional 50 MW of wind power is installed at HWN. Other thermal generations remain the same.
Case 1.2. This case is the same as Case 1.1 but additional 50 MW wind is installed at SSN.
Case 1.3. This case is the same as Case 1.1/1.2 but additional wind power is installed at HLM. These cases are compared with existing system (Case 1.4).
(2) A standalone system with additional 200 MW of wind power (i.e., a total of 300 MW) is identified.
Case 2. The system is assumed to be operating independent of the grid with 150, 90, and 60 MW of wind power generated at HWN, SSN, and HLM, respectively.
Simulations are run for April (lowest aggregated demand month) and August (highest aggregated demand month). The lower and upper bound values are associated with results for lowest and highest aggregated month, respectively. As shown in Table 3, the difference in wind speed variability at different locations has a significant impact on the optimal solution. The maximum expected lost wind occurs in Case 1.2 with an upper boundary of 6.01% of the total wind. A relatively reduced expected lost wind is noticed in Case 1.1. The boundary is tightened to %. In this case, the upper boundary of the expected lost wind and lost reserve is even slightly reduced compared to the existing system (Case 1.4). In both of these cases, much of the wind power is generated at HWN where the generation system is stressed with the same wind speed variability. However, in the other cases, much of the wind power is generated from other locations which in turn is accompanied with a different wind speed variability. The expected lost reserve in all cases is found to be % of the total spinning reserve and the expected lost load is negligible. HWN is an optimal location to install additional wind generation because of lower wind curtailment compared to that of SSN and HLM.
For the given transmission line capacities, the difference in results with and without considering network constraints is negligible. In order to demonstrate the impact of these constraints, the transmission capacity of the branch that connects HWN and SSN (line-14) is limited to 90 MW. Results are reported in Table 4. The expected lost wind in Case 1.1 is raised to % which is far worse compared to the result found when the network constraints and losses are ignored. Nearly 9–11% of the lost wind is introduced by the network constraints. Improved solution in terms of expected lost wind is noted in Case 1.2. The boundary is tightened to %. The upper boundary of the expected lost wind is even reduced in nearly 2% compared to the result reported in Table 3. However, lower expected lost load is obtained in Table 4. An increase in lost load shifts some portion of the thermal generation to secure the spinning reserves resulting in a reduced lost reserve. Moreover, some of the wind generation is utilized to cover the transmission losses which in turn results in a reduced lost wind. In all cases, the expected lost reserve is negligible and the transmission loss is bounded in between 1.45% and 1.67% of the total load. When grid constraints and losses are included in the model, a relatively reduced wind curtailment is observed at SSN (Case 1.2).
Case 2 investigates whether a standalone power system can operate independently and efficiently if the electricity supplied from the mainland grid is totally replaced by wind power. Whether the system is simulated with or without grid constraints, the expected lost wind is found to be % of the total wind as shown in Table 5. However, a significant difference in the expected lost load occurs. Without the network constraints, the lost load is limited to % of the total load and the expected lost reserve is found to be % of the total spinning reserve. On the other hand, when the constraints are included, the lost load is raised to % and the lost reserve is negligible. The transmission loss is less than 1% of the total demand in this standalone system. The expected lost load in the network-constrained system is nearly twice larger than the expected lost load when the network constraints are not included. This forces some portion of the thermal generation to be saved as spinning reserve. However, when the network constraints are not introduced, this part of the thermal generation is used to supply the load rather than securing the spinning reserve, resulting in larger lost reserve.
Note that, for the purpose of comparison, to generate similar trajectories for all the cases, the random number generator in MATLAB is set at specified stream so that it always draws the same sequence of random numbers for load and wind forecasting tool as well as for scenario generation tool.
5.5. Comparison of Loss Approximation Methods
The computational time, accuracy, memory requirements, and cost increment for different loss approximation methods compared to the original case (i.e., without considering network constraints and losses) are shown in Table 6. Since simulations take longer time (i.e., need to run 1464 times) for both of the extreme months, for purpose of comparing loss approximation methods, tests are carried out only for a typical day in one of the two extreme months in Case 1.4 with a single scenario for both wind and demand. In the iterative method, the computational time per single planning window increases by almost 35% and the total cost increment is calculated to be 11.74%. The memory requirement (measured in terms of number of variables) also increases by 30.48%. Since the losses in iterative approach are computed exactly similar to nonlinear quadratic losses and lumped as extra loads to the buses, the error is zero.
It is clearly seen that the number of segments in piecewise linear approximation has a considerable effect on computational time and accuracy. Increasing the number of segments reduces the total cost while it imposes memory burden and longer processing time. Approximations with four and six segments do not result in a reduced total cost compared to the iterative method. However, approximation with ten segments in this power system case is found to be a good tradeoff between solution delivery time and cost reduction. A reduction of 0.05% of total cost is gained compared to the iterative method. The average percentage deviation (error) from quadratic losses is about 4.9% which is a fair value, since practically transmission losses are already a small percentage of the total demand. Figure 12 shows the percentage deviation from the real quadratic losses using different number of piecewise segments. A piecewise linear loss approximation with ten segments is entirely used in all tests carried out to study the impact of network constraints in Section 5.4 and load shedding in the next subsection. Further increasing of segments yields to reduction in cost and a little improvement in accuracy but it raises the memory requirement and computational time. It should be noted that although the iterative approach is fast and converges with about five iterations, the piecewise method delivers a less expensive solution if fair number of segments are used for approximation.
5.6. Prioritized Load Shedding
To observe the impact of prioritized load shedding, different factors are applied for certain group of load buses. With decreasing order, the priorities are assigned for buses 3, (9, 4), and (2, 7, 8, 11). For purpose of demonstration, simulations are performed only for the largest aggregate demand month in Case 2. Results in Table 7 illustrate that when priority factors are applied, the less prioritized group of load buses contribute to the entire expected lost load while no load is shed in highly prioritized load buses. However, the total expected load, lost wind, and lost reserve in this specific case do not show a significant change when the system is simulated with and without priority load shedding.
Load shedding provides the balance of supply and demand. Furthermore, it guarantees a feasible solution and helps in reducing the peak power demand. In systems with load shedding capability, the installed generation capacity does not need to be designed to supply the largest demand periods, which subsequently reduces the total cost of the system significantly.
5.7. Joint Analysis of Grid Constraints, Transmission Losses, and Prioritized Load Shedding in Stochastic Scheduling
To study the possible interactions between the different aspects mentioned in the previous sections, a joint analysis is carried out for the existing system (Case 1.4) given that the transmission line capacity of line-14 is constrained to 90 MW, and results are shown in subsequent tables (Tables 8 and 9). Stochastic programming can handle the random nature of wind and demand using several trajectories constructed from scenario generation tool. The number of trajectories is reduced to a certain number of representative scenarios using -means as explained in Section 2. The number of scenarios considered affects the scheduling performance together with other factors such as unit ramp-up/ramp-down rates and minimum-up/minimum-down times, grid constraints related with transmission line capacity, and transmission losses.
When load priority is not considered and transmission losses are not modelled, the upper boundary of expected lost wind decreases from 12.64% (for one scenario) to 10.86% (for sixteen scenarios). The scheduling performance (measured in terms of lost wind) is improved as the number of scenarios increases. Likewise, when losses are considered with nonprioritized load shedding scheme, the upper boundary of the expected lost wind decreases from 12.87% (for one scenario) to 11.43% (for sixteen scenarios). Modelling the transmission losses accounts for an increase in 0.23% of lost wind for a single scenario and 0.57% for sixteen scenarios. It is observed that when the system is simulated with small number of scenarios, the effect of cumulative uncertainty of wind and demand is higher than the impact of considering the transmission losses on the scheduling performance. For instance, when losses are not considered, increasing the number of scenarios from one to four results in 1.09% reduction of the expected lost wind; however, when losses are considered for a single scenario, the expected lost wind increases from 12.64% to 12.87% (i.e., only an increment in 0.23%). On the other hand, as the number of scenarios grows, the effect of intrinsic uncertainty of wind and demand on the scheduling performance decreases compared to the impact of transmission losses. For instance, increasing the number of scenarios from nine to sixteen in ideal lossless system results in 0.12% reduction of the expected lost wind; however, when losses are not omitted and nine scenarios are used, it increases by 0.54% compared to the lossless system with the same number of scenarios.
When prioritized load shedding scheme is implemented and transmission losses are ignored, the upper boundary of expected lost wind decreases from 14.84% to 12.91% for one and sixteen scenarios, respectively. The magnitude of the lost wind increases compared to the corresponding results obtained in nonprioritized load shedding scheme. For instance, an increment in 2.08% of expected lost wind resulted for prioritized ideal lossless system simulated with nine scenarios compared to nonprioritized ideal lossless system with identical number of scenarios. When losses are considered together with prioritized load shedding, the expected lost wind increases in 0.35% and 0.63%, for one scenario and sixteen scenarios, respectively. Similar to the results observed from nonprioritized load shedding scheme, the impact of transmission losses on scheduling performance increases as the number of scenarios grows.
The expected lost load in all cases is bounded in the range . However, it shows a very slight increase in nonprioritized load shedding scheme compared to the prioritized one. The expected lost reserve is negligible, and the average transmission losses are approximately 1.55% of the total load in all cases.
In summary, it has been observed that modeling of transmission losses is relevant after the effect of wind and demand uncertainties on scheduling performance is sufficiently reduced using large number of scenarios. In this power system case, not so much is gained beyond sixteen scenarios; that is, the scheduling performance is not affected by the uncertainty any more even if the number of scenarios kept increasing. However, for small number of scenarios, the effect of transmission losses on scheduling performance is surpassed by the uncertainty in wind and load. In addition, when prioritized load shedding scheme is introduced, the magnitude of expected lost wind increases compared to the expected lost wind when load priority is not assumed. For the same number of scenarios, prioritized load shedding scheme has more significant effect on scheduling performance than the transmission losses. Note that the expected lost wind is much larger than the expected lost load and lost reserve. This is mainly related with the penalty factors associated with each slack variable in the objective function. The penalty factors are decided by the system operator according to their importance. Some system operators have a tight policy of not losing an electricity customer. In that situation, lost load is more penalized than wasted wind and lost reserve. If higher penalty factors are assigned for lost load () and lost reserve () compared to lost wind (), the optimizer tries to minimize the expected lost load and reserve while it imposes the system to lose more available wind power for the minimum possible value of the objective function.
6. Conclusion
In this paper, a comprehensive stochastic optimization model, which accounts for the random nature of load and wind generation, grid constraints, transmission losses, and thermal unit operational constraints, is presented. The loss modelling approaches are analysed and compared based on accuracy, cost minimization, and CPU processing time. Slack variables are introduced for implementing prioritized load shedding, in addition to securing solution feasibility. Moreover, the individual impacts of grid constraints, transmission losses, and prioritized load shedding and the possible interactions when jointly considered are studied under the framework of stochastic programming.
The developed stochastic scheduling model is applied to test the generation flexibility of the modified power system of Jeju Island. In the first case study that the island still maintains the grid supply, the magnitude of the wind curtailment depends on the location where much of the generation is installed. This is basically related with the difference in wind speed variability at different locations. The network constraints have also a considerable impact on increasing the wind curtailment. In some situations, the wasted wind rises up to 13.46% whereas the lost load is nearly negligible in all cases. In the standalone system, the wind curtailment is less than 4%. However, the unserved load is significantly increased, and it reaches up to 10.37% for the month with largest aggregate demand. Hence, the island’s system operator needs to install extra flexible thermal units in order to reduce the magnitude of unserved loads and wind curtailment.
The model can be adopted and applied to any size of power system network. Reactive power requirement and voltage magnitude restrictions that could have significant impact on further wind curtailment are to be considered in the future work.
Nomenclature
| : | Time step, set of time steps |
| : | Set of time steps of nonrecourse committed variables |
| : | Stochastic scenario, set of scenarios |
| : | Thermal unit, set of thermal units |
| : | Wind unit, set of wind units |
| : | Bus, set of buses |
| : | Line, set of lines |
| : | Segment of piecewise linear cost, set of linear cost segments |
| : | Segment of piecewise linear loss, set of linear loss segments |
| : | Set of thermal units connected to bus |
| : | Set of lines connected to bus |
| : | Set of wind units connected to bus . |
| : | Probability of scenario [] |
| : | Start-up cost of unit [] |
| : | Shut-down cost of unit [] |
| : | Penalty for lost load [/MW] |
| : | Penalty for lost reserve [/MW] |
| : | Penalty for lost wind [/MW] |
| : | Fraction of the total load [] |
| : | Priority factor of load at bus |
| : | Minimum-down time of unit |
| : | Minimum-up time of unit |
| : | Minimum number of units committed to be online |
| : | Number of periods unit should be initially online |
| : | Number of periods unit should be initially offline |
| : | Number of online periods of unit prior to planning |
| : | Number of offline periods of unit prior to planning |
| : | Status of unit before planning window |
| : | Ramp-up rate of unit [MW/h] |
| : | Ramp-down rate of unit [MW/h] |
| : | Upper rating of unit [MW] |
| : | Lower rating of unit [MW] |
| : | Base power [MW] |
| : | Maximum of voltage angle [rad] |
| : | Transmission capacity of line [MW] |
| : | Piecewise linear cost coordinates [MW] |
| : | Piecewise linear cost coordinates [/h] |
| : | Piecewise linear loss coordinates [rad] |
| : | Piecewise linear loss coordinates [MW] |
| : | Load [MW] |
| : | Wind [MW] |
| : | Power loss [MW] |
| : | Resistance of line [p.u] |
| : | Reactance of line [p.u] |
| : | branch incidence matrix |
| : | susceptance matrix |
| : | diagonal susceptance matrix. |
| : | Power output of unit [MW] |
| : | Voltage angle at bus [rad] |
| : | Power flow on line [MW] |
| : | Available power of unit at period [MW] |
| : | Transmission loss of line m [MW] |
| : | Piecewise linear cost segment [] |
| : | Fraction of piecewise linear cost segment [] |
| : | Piecewise linear loss segment [] |
| : | Fraction of piecewise linear loss segment [] |
| : | Lost load at bus [MW] |
| : | Lost spinning reserve in period [MW] |
| : | Lost wind [MW] |
| : | Start-up cost of unit in time step [] |
| : | Shut-down cost of unit in time step [] |
| : | Spinning reserve at period [MW]. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to sincerely thank the editor and the anonymous reviewers for their valuable comments and suggestions, which improved the quality of this paper.