Computationally Efficient Market Simulation Tool for Future Grid Scenario Analysis
Shariq Riaz Gregor Verbic Archie C. Chapman





Computationally Efficient Market Simulation Tool for Future Grid Scenario Analysis
    Shariq Riaz, Graduate Student Member, IEEE,
        Gregor Verbič, Senior Member, IEEE,
        and Archie C. Chapman, Member, IEEE
Shariq Riaz, Gregor Verbič and Archie C. Chapman are with the School of Electrical and Information Engineering, The University of Sydney, Sydney, New South Wales, Australia. 	e-mails: (shariq.riaz, gregor.verbic, archie.chapman@sydney.edu.au).
Shariq Riaz is also with the Department of Electrical Engineering, University of Engineering and Technology Lahore, Lahore, Pakistan.

    December 30, 2023
==================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




This paper proposes a computationally efficient electricity market simulation tool (MST) suitable for future grid scenario analysis. The market model is based on a unit commitment (UC) problem and takes into account the uptake of emerging technologies, like demand response, battery storage, concentrated solar thermal generation, and HVDC transmission lines. To allow for a subsequent stability assessment, the MST requires an explicit representation of the number of online generation units, which affects powers system inertia and reactive power support capability. These requirements render a full-fledged UC model computationally intractable, so we propose unit clustering, a rolling horizon approach, and constraint clipping to increase the computational efficiency. To showcase the capability of the proposed tool, we use a simplified model of the Australian National Electricity Market with different penetrations of renewable generation. The results 

are verified by comparison to a more expressive and computationally-intensive binary UC, which confirm the validity of the approach for long term future grid studies. 



Electricity market, future grid, electricity market simulation tool, optimization, scenario analysis, unit commitment, stability assessment, inertia, loadability.


[A01]𝒞Set of consumers c.
[A02]𝒢Set of generators g, 𝒢^ = 𝒢^𝓈𝓎𝓃∪𝒢^ℛℰ𝒮.
[A03]𝒢^𝓈𝓎𝓃Set of synchronous generators, 𝒢^𝓈𝓎𝓃⊆𝒢.
[A04]𝒢^ℛℰ𝒮Set of renewable generators, 𝒢^ℛℰ𝒮⊆𝒢.
[A05]𝒢^𝒞𝒮𝒯Set of concentrated solar thermal generators, 𝒢^𝒞𝒮𝒯⊆𝒢^syn.
[A06]𝒢^𝓇Set of synchronous generators in region r, ⋃_𝒢^𝓇  = 𝒢^.
[A07]ℋSet of sub-horizons h.
[A08]ℒSet of power lines l, ℒ = ℒ^𝒜𝒞∪ℒ^ℋ𝒱𝒟𝒞.
[A09]ℒ^𝒜𝒞Set of AC power lines, ℒ^𝒜𝒞⊆ℒ.
[A10]ℒ^ℋ𝒱𝒟𝒞Set of HVDC power lines, ℒ^ℋ𝒱𝒟𝒞⊆ℒ.
[A11]𝒩Set of nodes n.
[A12]𝒩^𝓇Set of nodes in region r.
[A13]𝒫Set of prosumers p.
[A14]ℛSet of regions r.
[A15]𝒮Set of storage plants s.
[A16]𝒯Set of time slots t.

[D01]s_g,tNumber of online units of generator g, s_g,t∈{0,1} in BUC and s_g,t∈ℤ_+ in MST.
[D02]u_g,tInteger startup status variable of a unit of generator g, u_g,t∈{0,1} in BUC and u_g,t∈ℤ_+ in MST.
[D03]d_g,tInteger shutdown status variable of a unit of generator g, d_g,t∈{0,1} in BUC and d_g,t∈ℤ_+ in MST.
[D04]δ_n,tVoltage angle at node n.
[D05]p_l,t^Power flow on line l.


[D08]Δ p_l,tPower loss on line l.
[D09]p_g,tPower dispatch of generator g.
[D10]p_p,t^g+/-Grid/feed-in power of prosumer p.
[D11]p_s,tPower flow of storage plant s.
[D12]p_p,t^bBattery power flow of prosumer p.
[D13]e_g,tThermal energy stored in TES of generator g ∈𝒢^𝒞𝒮𝒯.
[D14]e_s,tEnergy stored in storage plant s.
[D15]e_p,t^bBattery charge state of prosumer p.
	
[P]c_g^fix/varFix/variable cost of a unit of generator g.
[P]c_g^su/sdStartup/shutdown cost of a unit of generator g.
[P]p_c,t^Load demand of consumer c.
[P]p_p,t^Load demand of prosumer p.
[P]p_n,t^rPower reserve requirement of node n.
[P]x/xMinimum/maximum limit of variable x.
[P]U_gTotal number of identical units of generator g.
[P]r^+/-_gRamp-up/down rate of a unit of generator g.
[P]τ^u/d_gMinimum up/down time of a unit of generator g.
[P]t̃Time slot offset index.
[P]ΔtTime resolution.
[P]B_lSusceptance of line l.
[P]p_g,t^RESMax. output power of renewable generator g ∈𝒢^RES.
[P]p_g,t^CSTMax. thermal power capture by generator g ∈𝒢^CST.
[P]H_gInertia of a unit of generator g.
[P]S_gMVA rating of a unit of generator g.
[P]H_n,tMinimum synchronous inertia requirement of node n.
[P]η_xEfficiency of component x.
[P]p_p,t^pvAggregated PV power of prosumer p.
[P]λFeed-in price ratio.

[I]ŝ_gNumber of online units of generator g ∈𝒢^syn at start of horizon.
[I]p̂_gPower dispatch of generator g at start of horizon.
[I]û_g,tMinimum number of units of generator g ∈𝒢^syn required to remain online for time t<τ_g^u.
[I]d̂_g,tMinimum number of units of generator g ∈𝒢^syn required to remain offline for time t<τ_g^d.
[I]ê_gEnergy stored in TES of g ∈𝒢^𝒞𝒮𝒯 at start of horizon.
[I]ê_sEnergy stored in storage plant s at start of horizon.
[I]ê_p^bBattery state of charge for prosumer p at start of horizon.









§ INTRODUCTION

Power systems worldwide are moving away from domination by large-scale synchronous generation and passive consumers.
Instead, in future grids[We interpret a future grid to mean the study of national grid type structures with the transformational changes over the long-term out to 2050.] new actors, such as variable renewable energy sources (RES)[For the sake of brevity, by RES we mean “unconventional” renewables like wind and solar, but excluding conventional RES, like hydro, and dispatchable unconventional renewables, like concentrated solar thermal.], price-responsive users equipped with small-scale PV-battery systems (called prosumers), demand response (DR), and energy storage will play an increasingly important role.
Given this, in order for policy makers and power system planners to evaluate the integration of high-penetrations of these new elements into future grids, new simulation tools need to be developed. 
Specifically, there is a pressing need to understand the effects of technological change on future grids, in terms of energy balance, stability, security and reliability, over a wide range of highly-uncertain future scenarios.
This is complicated by the inherent and unavoidable uncertainty surrounding the availability, quality and cost of new technologies (e.g. battery or photo-voltaic system costs, or concentrated solar thermal (CST) generation operating characteristics) and the policy choices driving their uptake.
The recent blackout in South Australia <cit.> serves as a reminder that things can go wrong when the uptake of new technologies is not planned carefully. 

Future grid planning thus requires a major departure from conventional power system planning, where only a handful of the most critical scenarios are analyzed. To account for a wide range of possible future evolutions, scenario analysis has been proposed in many industries, e.g. in finance and economics <cit.>, and in energy <cit.>. In contradistinction to power system planning, where the aim is to find an optimal transmission and/or generation expansion plan, the aim of scenario analysis is to analyze possible evolution pathways to inform power system planning and policy making. Given the uncertainty associated with long-term projections, the focus of future grid scenario analysis is limited only to the analysis of what is technically possible, although it might also consider an explicit costing <cit.>.
In more detail, existing future grid feasibility studies have shown that the balance between demand and supply can be maintained even with high penetration of RESs by using large-scale storage, flexible generation, and diverse RES technologies <cit.>.
However, they only focus on balancing and use simplified transmission network models (either copper plate or network flow; a notable exception is the Greenpeace pan-European study <cit.> that uses a DC load flow model). This ignores network related issues, which limits these models' applicability for stability assessment. 

To the best of our knowledge, the Future Grid Research Program, funded by the Australian Commonwealth Scientific and Industrial Research Organisation (CSIRO) is the first to propose a comprehensive modeling framework for future grid scenario analysis that also includes stability assessment. The aim of the project is to explore possible future pathways for the evolution of the Australian grid out to 2050 by looking beyond simple balancing. To this end, a simulation platform has been proposed in <cit.> that consists of a market model, power flow analysis, and stability assessment, Fig. <ref>. The platform has been used, with additional improvements, to study fast stability scanning <cit.>, inertia <cit.>, modeling of prosumers for market simulation <cit.>,  impact of prosumers on voltage stability <cit.>, and power system flexibility using CST <cit.> and battery storage <cit.>.
In order to capture the inter-seasonal variations in the renewable generation, computationally intensive time-series analysis needs to be used.
A major computational bottleneck of the framework is the market simulation.



Within this context, the contribution of this paper is to propose a unified generic market simulation tool (MST) based on a unit commitment (UC) problem suitable for future grid scenario analysis, including stability assessment. The tool incorporates the following key features:

	
  * market structure agnostic modeling framework,
	
  * integration of various types and penetrations of RES and emerging demand-side technologies,
	
  * generic demand model considering the impact of prosumers,
	
  * explicit network representation, including HVDC lines, using a DC power flow model, 
	
  * explicit representation of the number of online synchronous generators,
	
  * explicit representation of system inertia and reactive power support capability of synchronous generators,
	
  * computational efficiency with sufficient accuracy.


The presented model builds on our existing research <cit.> and combines all these in a single coherent formulation.

In more detail, to reduce the computational burden, the following techniques are used building on the methods proposed in <cit.>:

	
  * unit clustering,
	
  * rolling horizon approach,
	
  * constraint clipping.
 

The computational advantages of our proposed model are shown on a simplified 14-generator model of the Australian National Energy Market (NEM) as a test grid <cit.>. Four cases for different RES penetration are run for one to seven days horizon length, and computational metrics are reported. To reflect the accuracy of the proposed MST, system inertia and voltage stability margins are used as a benchmark. In simulations, RES and load traces are taken from the National Transmission Network Developed Plan (NTNDP) data, provided by the Australian Energy Market Operator (AEMO) <cit.>.

The remainder of the paper is organized as follows: Literature review and related work are discussed in Section II, while Section
III details the MST. A detailed description
of the simulation setup is given in Section IV. In Section V results are analyzed and discussed in detail. Finally, Section VI concludes the paper.




§ RELATED WORK

In order to better explain the functional requirements of the proposed MST, we first describe the canonical UC formulation.
An interested reader can find a comprehensive literature survey in <cit.>.




 §.§ Canonical Unit Commitment Formulation

The UC problem is an umbrella term for a large class of problems in power system operation and planning whose objective is to schedule and dispatch power generation at minimum cost to meet the anticipated demand, while meeting a set of system-wide constraints. In smart grids, problems with a similar structure arise in the area of energy management, and they are sometimes also called UC <cit.>.
Before deregulation, UC was used in vertically integrated utilities for generation scheduling to minimize production costs. After deregulation, UC has been used by system operators to maximize social welfare, but the underlying optimization model is essentially the same.

Mathematically, UC is a large-scale, nonlinear, mixed-integer optimization problem under uncertainty. With some abuse of notation, the UC optimization problem can be represented in the following compact formulation <cit.>:

    minimize_𝐱_c, 𝐱_b      f_c(𝐱_c) + f_b(𝐱_b) 
    subject  to      g_c(𝐱_c) ≤𝐛
        g_b(𝐱_b) ≤𝐜
        h_c(𝐱_c) + h_b(𝐱_b)≤𝐝
       𝐱_c∈ℝ^+, 𝐱_b∈{ 0,1 }


Due to the time-couplings, the UC problem needs to be solved over a sufficiently long horizon.
The decision vector 𝐱 = {𝐱_c, 𝐱_b} for each time interval consist of continuous and binary variables. The continuous variables, 𝐱_c, include  generation dispatch levels, load levels, transmission power flows, storage levels, and transmission voltage magnitudes and phase angles. The binary variables, 𝐱_b, includes scheduling decisions for generation and storage, and logical decisions that ensure consistency of the solution.
The objective (<ref>) captures the total production cost, including fuel costs, start-up costs and shut-down costs.
The constraints include, respectively: dispatch related constraints such as energy balance, reserve requirements, transmission limits, and ramping constraints (<ref>); commitment variables, including minimum up and down, and start-up/shut-down constraints (<ref>); and constraints coupling commitment and dispatch decisions, including minimum and maximum generation capacity constraints (<ref>).  

The complexity of the problem stems from the following: (i) certain generation technologies (e.g. coal-fired steam units) require long start-up and shut-down times, which requires a sufficiently long solution horizon; (ii) generators are interconnected, which introduces couplings through the power flow constraints; (iii) on/off decisions introduce a combinatorial structure; (iv) some constraints (e.g. AC load flow constraints) and parameters (e.g. production costs) are non-convex; and (v) the increasing penetration of variable renewable generation and the emergence of demand-side technologies introduce uncertainty.
As a result, a complete UC formulation is computationally intractable, so many approximations and heuristics have been proposed to strike a balance between computational complexity and functional requirements. For example, power flow constraints can be neglected altogether (a copper plate model), can be replaced with simple network flow constraints to represent critical inter-connectors, or, instead of (non-convex) AC, a simplified (linear) DC load flow is used. 




 §.§ UC Formulations in Existing Future Grid Studies

In operational studies: the nonlinear constraints, e.g. ramping, minimum up/down time (MUDT) and thermal limits are typically linearized; startup and shutdown exponential costs are discretized, and; non-convex and non-differentiable variable cost functions are expressed as piecewise linear function <cit.>. In planning studies, due to long horizon lengths, the UC model is simplified even further. For example: combinatorial structure is avoided by aggregating all the units installed at one location <cit.>; piecewise linear cost functions and constraints are represented by one segment only; some costs (e.g. startup, shutdown and fix costs) are ignored; a deterministic UC with perfect foresight is used, and; non-critical binding constraints are omitted <cit.>[An interested reader can refer to <cit.> for a discussion on binding constraints elimination for generation planning.].
To avoid the computational complexity associated with the mixed integer formulation, a recent work <cit.> has proposed a linear relaxation of the UC formulation for flexibility studies, with an accuracy comparable to the full binary mixed integer linear formulation. 

In contrast to operation and planning studies, the computational burden of future grid scenario analysis is even bigger, due to a sheer number of scenarios that need to be analyzed, which requires further simplifications. For example, the Greenpeace study <cit.> uses an optimal power flow for generation dispatch and thus ignores UC decisions. Unlike the Greenpeace study, the Irish All Island Grid Study <cit.> and the European project e-Highway2050 <cit.> ignore load flow constraints altogether, however they do use a rolling horizon UC, with simplifications. The Irish study, for example doesn't put any restriction on the minimum number of online synchronous generators to avoid RES spillage, and the e-Highway2050 study uses a heuristics to include DR. The authors of the e-Highway2050 study, however, acknowledge the size and the complexity of the optimization framework in long term planning, and plan to develop new tools with a simplified network representation <cit.>.

In summary, a UC formulation depends on the scope of the study. Future grid studies that explicitly include stability assessment bring about some specific requirements that are routinely neglected in the existing UC formulations, as discussed next. 




§ MARKET SIMULATION TOOL



 §.§ Functional Requirements

The focus of our work is stability assessment of future grid scenarios. Thus, MST must produce dispatch decisions that accurately capture the kinetic energy stored in rotating masses (inertia), active power reserves and reactive power support capability of synchronous generators, which all depend upon the number of online units and the respective dispatch levels. 

For the sake of illustration, consider a generation plant consisting of three identical (synchronous) thermal units, with the following characteristics: (i) constant terminal voltage of 1pu; (ii) minimum technical limit P_min = 0.4pu; (iii) power factor of 0.8; (iv) maximum excitation limit E_fd^max = 1.5pu; and (v) normalized inertia constant H = 5. We further assume that in the over-excited region, the excitation limit is the binding constraint, as shown in Fig. <ref>. Observe that the maximum reactive power capability depends on the active power generated, and varies between Q_n at P_max = 1pu and Q_max at P_min.
We consider three cases defined by the total active power generation of the plant: (i) 0.8pu, (ii) 1.2pu, and (iii) 1.6pu. 
The three scenarios correspond to the rows  in Fig. <ref>, which shows the active power dispatch level P, reactive power support capability Q, online active power reserves R, and generator inertia H. 
The three columns show feasible solutions for three different UC formulations: all three units are aggregated into one equivalent unit (AGG), standard binary UC (BUC) when each unit is modeled individually, and the proposed market simulation tool (MST). A detailed comparison of the three formulations is given in Section V.

 

Although the results are self-explanatory, a few things are worth emphasizing. In case (i), aggregating the units into one equivalent unit (AGG) results in the unit being shut down due to the minimum technical limit. The individual unit representation (BUC), on the other hand, does allow the dispatch of one or two units, but with significantly different operational characteristics. In cases (ii) and (iii), the total inertia in the AGG formulation is much higher, which has important implications for frequency stability. A similar observation can be made for the reactive power support capability, which affects voltage stability. Also, dispatching power from all three units results in a significantly higher active power reserve. And last, a higher reactive power generation due to a lower P reduces the internal machine angle, which improves transient stability.

In conclusion, a faithful representation of the number of online synchronous machines is of vital importance for stability assessment. An individual unit representation, however, is computationally expensive, so the computational burden should be reduced, as discussed in the following section. Next, an explicit network representation is required. An AC load flow formulation, however, is nonlinear (and non-convex), which results in an intractable mixed-integer nonlinear problem. Therefore, we use a DC load flow representation with a sufficiently small voltage angle difference on transmission lines. Our experience shows that an angle difference of 30 results in a manageable small number of infeasible operating conditions that can be dealt with separately.




 §.§ Computational Speedup
 
The MST is based on the UC formulation using constant fixed, startup, shutdown and production costs. To improve its computational efficiency, the dimensionality of the optimization problem is reduced employing: (i) unit clustering <cit.> to reduce the number of variables needed to represent a multi-unit generation plant; (ii) a rolling horizon approach <cit.> to reduce the time dimension; and (iii) constraint clipping to remove most non-binding constraints.  



  §.§.§ Unit Clustering

Linearized UC models are computationally efficient for horizons of up to a few days, which makes them extremely useful for operational studies. For planning studies, however, where horizon lengths can be up to a year, or more, these models are still computationally too expensive. Our work builds on the clustering approach proposed in <cit.>, where identical units at each generation plant are aggregated by replacing binary variables with fewer integer variables. The status of online units, startup/shutdown decisions and dispatched power are tracked by three integer variables and one continuous variable per plant per period, as opposed to three binary and one continuous variable per unit per period. Further clustering proposed in <cit.> is not possible in our formulation because of the explicit network representation required in the MST. 



  §.§.§ Rolling Horizon
 
Solving the UC as one block, especially for long horizons, is computationally too expensive. This can be overcome by breaking the problem into several smaller intervals called sub-horizons <cit.>. To ensure accuracy and consistency of the solution, a proper overlap between sub-horizons is maintained and the terminating state of the previous sub-horizon is used as the initial condition of the next sub-horizon. The minimum sub-horizon length depends on the time constants associated with the decision variables. While these might be in the order of hours for thermal power plants, they can be significantly longer for energy storage. Large-scale hydro dams, for example, require horizon lengths of several weeks, or even months. In our research, however, the sub-horizon length is up to a few days to cater for thermal energy storage (TES) of CST plants and battery storage. The optimization of hydro dams is not explicitly considered, however it can be taken into account heuristically, if needed.



  §.§.§ Constraint Clipping

The size of the problem can be reduced by removing non-binding constraints, which doesn't affect the feasible region. For instance, an MUDT constraint on a unit with an MUDT less than the time interval is redundant[This is especially the case when the time resolution is coarse. In our studies, the time step is one hour. In operational studies, where the resolution can be as short as five minutes, constraint clipping is less useful.]. Similarly, a ramp constraint for flexible units is redundant if the time step is sufficiently long. With a higher RES penetration, in particular, where backup generation is provided by fast-ramping gas turbines, this technique can significantly reduce the size of the optimization problem, and hence improves the computational performance due to a larger number of units with higher ramp rates and smaller MUDTs. It should be noted that optimization pre-solvers might not able to automatically remove these constraints. 




 §.§ MST UC Formulation

























  §.§.§ Objective function

The objective of the proposed MST is to minimize total generation cost for all sub-horizons h:

    Ωminimize ∑_t∈𝒯^∑_g∈𝒢^( c_g^fix s_g,t
    +c_g^su u_g,t + c_g^sd  d_g,t +c_g^var  p_g,t),
	
where Ω = {s_g,t,u_g,t,d_g,t,p_g,t, p_s,t, p_l,t} are the decision variables of the problem, and c_g^fix, c_g^su, c_g^sd, and c_g^var are fixed, startup, shutdown and variable cost, respectively.
As typically done in planning studies <cit.>, <cit.>, the costs are assumed constant to reduce the computational complexity. The framework, however, also admits a piece-wise linear approximation proposed in <cit.>.



  §.§.§ System Constraints

System Constraints[All the constraints must be satisfied in all time slots t, however, for sake of notational brevity, this is not explicitly mentioned.] include power balance constraints, power reserve and minimum synchronous inertia requirements.

Power balance:
Power generated at node n must be equal to the node power demand plus the net power flow on transmission lines connected to the node:

    ∑_g∈𝒢_n^p_g,t =  
    ∑_c ∈𝒞_np_c,t^ + ∑_p ∈𝒫_n p_p,t^g+ - ∑_p ∈𝒫_n p_p,t^g- + ∑_s ∈𝒮_np_s,t + ∑_l ∈ℒ_n(p_l,t +Δ p_l,t),

where 𝒢_n, 𝒞_n, 𝒫_n, 𝒮_n, ℒ_n represent respectively the set of generators, consumers, prosumers[Price-responsive users equipped with small-scale PV-battery systems.], utility storage plants and lines connected to node n. 

Power reserves:
To cater for uncertainties, active power reserves provided by synchronous generation g ∈𝒢^syn are maintained in each region r:

    ∑_g ∈{ (𝒢^syn-𝒢^CST) ∩𝒢^r} (p_g s_g,t - p_g,t) + 
    ∑_g ∈{𝒢^CST∩𝒢^r}min(p_g s_g,t - p_g,t,e_g,t-p_g,t)  ≥∑_n ∈𝒩_r^p_n,t^r.

For synchronous generators other than concentrated solar thermal (CST), reserves are defined as the difference between the online capacity and the current operating point. For CST, reserves can either be limited by their online capacity or energy level of their thermal energy system (TES).
Variable s_g,t in (<ref>) represents the total number of online units at each generation plant, and 𝒢^r and 𝒩_r represent the sets of generators and nodes in region r, respectively.

Minimum synchronous inertia requirement: 
To ensure frequency stability, a minimum level of inertia provided by synchronous generation must be maintained at all times (more details are available in <cit.>) in each region r:  

    ∑_g ∈{𝒢^syn∩𝒢^r}^ s_g,tH_g S_g ≥∑_n ∈𝒩_r^H_n,t.




  §.§.§ Network constraints

Network constraints include DC power flow constraints and thermal line limits for AC lines, and active power limits for HVDC lines.  

Line power constraints: 
A DC load flow model is used for computational simplicity for AC transmission lines[A sufficiently small (∼30) voltage angle difference over a transmission line is used to reduce the number of nonconvergent AC power flow cases.]:

    p_l,t^x,y = B_l(δ_x,t - δ_y,t),    l ∈ℒ^𝒜𝒞,

where the variables δ_x,t and δ_y,t represent voltage angles at nodes x ∈𝒩 and y ∈𝒩, respectively.

Thermal line limits:
Power flows on all transmission lines are limited by the respective thermal limits of line l:

    |  p_l,t|≤p_l,

where p_l represents the thermal limit of line l.



  §.§.§ Generation constraints

Generation constraints include physical limits of individual generation units.
For the binary unit commitment (BUC), we adopted a UC formulation requiring three binary variables per time slot (on/off status, startup, shutdown) to model an individual unit. In the MST, identical units of a plant are clustered into one individual unit <cit.>. 
This requires three integer variables (on/of status, startup, and shutdown) per generation plant per time slot as opposed to three binary variables per generation unit per time slot in the BUC, as discussed in Section III.B of A Computationally Efficient Market Model for Future Grid Scenario Studies.

Generation limits: 
Dispatch levels of a synchronous generator g are limited by the respective stable operating limits:

    s_g,tp_g≤ p_g,t≤ s_g,tp_g,     g ∈𝒢^syn.

The power of RES[For the sake of brevity, by RES we mean “unconventional” renewables like wind and solar, but excluding conventional RES, like hydro, and dispatchable unconventional renewables, like concentrated solar thermal.] generation is limited by the availability of the corresponding renewable resource (wind or sun):

    s_g,tp_g≤ p_g,t≤ s_g,tp_g,t^RES,     g ∈𝒢^RES.


Unit on/off constraints: 
A unit can only be turned on if and only if it is in off state and vice versa:

    u_g,t-d_g,t=s_g,t-s_g,t-1,     t ≠ 1,  g ∈𝒢^syn.
	
In a rolling horizon approach, consistency between adjacent time slots is ensured by:

    u_g,t-d_g,t=s_g,t - ŝ_g,     t =1,  g ∈𝒢^syn,

where ŝ_g is the initial number of online units of generator g. Equations (<ref>) and (<ref>) also implicitly determine the upper bound of u_g,t and d_g,t in terms of changes in  s_g,t.

Number of online units:
Unlike the BUC, the MST requires an explicit upper bound on status variables:

    s_g,t≤U_g.
	


Ramp-up and ramp-down limits:
Ramp rates of synchronous generation should be kept within the respective ramp-up (<ref>), (<ref>) and ramp-down limits (<ref>), (<ref>):

    p_g,t - p_g,t-1≤ s_g,tr^+_g,     t ≠ 1, g ∈{𝒢^syn | r^+_g < p_g}, 
        p_g,t - p̂_g≤ s_g,tr^+_g,     t =1, g ∈{𝒢^syn | r^+_g < p_g},  
       p̂_g - p_g,t≤ s_g,t-1r^-_g,      t ≠ 1, g ∈{𝒢^syn | r^-_g < p_g},  
       p̂_g - p_g,t≤ŝ_gr^-_g,     t =1, g ∈{𝒢^syn | r^-_g < p_g}.

In the MST, a ramp limit of a power plant is defined as a product of the ramp limit of an individual unit and the number of online units in a power plant s_g,t. If s_g,t is binary, these ramp constraints are mathematically identical to ramp constraints of the BUC.   
If a ramp rate multiplied by the length of the time resolution Δt is less than the rated power, the rate limit has no effect on the dispatch, so the corresponding constraint can be eliminated.
Constraints explicitly defined for t=1 are used to join two adjacent sub-horizons in the rolling-horizon approach.

Minimum up and down times: 
Steam generators must remain on for a period of time  τ_g^u once turned on (minimum up time):

    s_g,t≥∑_t̃=τ_g^u-1^0 u_g,t-t̃,     t ≥τ_g^u,   g ∈{𝒢^syn | τ_g^u > Δt}, 
        s_g,t≥∑_t̃=t-1^0u_g,t-t̃  + û_g,t,     t < τ_g^u,  g ∈{𝒢^syn | τ_g^u > Δt}.
		
Similarly, they must not be turned on for a period of time  τ_g^d once turned off (minimum down time):

    s_g,t≤U_g - ∑_t̃=τ_g^d-1^0 d_g,t-t̃,     t ≥τ_g^d,  g ∈{𝒢^syn | τ_g^d > Δt}, 
        s_g,t≤U_g - ∑_t̃=t-1^0d_g,t-t̃ - d̂_g,t,     t < τ_g^d,  g ∈{𝒢^syn | τ_g^d > Δt}.
	
Similar to the rate limits, if the minimum up and down times are smaller than the time resolution Δt, the corresponding constraints can be eliminated. 
Due to integer nature of discrete variables in the MST, the definition of the MUDT constraints in the RH approach requires the number of online units for the last τ^u/d time interval to establish the relationship between the adjacent sub-horizons. If the τ_g^u/d is smaller than time resolution Δt, then these constraints can be eliminated. 



  §.§.§ CST constraints:

CST constraints include TES energy balance and storage limits.                

TES state of charge (SOC) 
determines the TES energy balance subject to the accumulated energy in the previous time slot, thermal losses, thermal power provided by the solar farm and electrical power dispatched from the CST plant:

    e_g,t=η_ge_g,t-1+p_g,t^CST-p_g,t,	    t ≠ 1,  g ∈𝒢^CST, 
        e_g,t=η_gê_g+p_g,t^CST-p_g,t,      t=1,  g ∈𝒢^CST,
	
where, p_g,t^CST is the thermal power collected by the solar field of generator g ∈𝒢^CST. 

TES limits:	Energy stored is limited by the capacity of a storage tank:

    e_g≤e_g,t≤e_g,	   g ∈𝒢^CST.




  §.§.§ Utility storage constraints

Utility-scale storage constraints include energy balance, storage capacity limits and power flow constraints. The formulation is generic and can capture a wide range of storage technologies.  

Utility storage SOC limits determine the energy balance of storage plant  s:

    e_s,t=η_se_s,t-1+p_s,t,	    t ≠ 1, 
        e_s,t=η_sê_s+p_s,t,	    t=1.
	

Utility storage capacity limits:
Energy stored is limited by the capacity of storage plant s:

    e_s≤e_s,t≤e_s.
	

Charge/discharge rates limit the charge and discharge powers of storage plant s:

    p_s^- ≤p_s,t≤p_s^+,

where p_s^- and p_s^+ represent the maximum power discharge and charge rates of a storage plant, respectively.



  §.§.§ Prosumer sub-problem

The prosumer sub-problem captures the aggregated effect of prosumers. It is modeled using a bi-level framework in which the upper-level unit commitment problem described above minimizes the total generation cost, and the lower-level problem maximizes prosumers' self-consumption. The coupling is through the prosumers' demand, not through the electricity price, which renders the proposed model market structure agnostic. As such, it implicitly assumes a mechanism for demand response aggregation. The Karush-Kuhn-Tucker optimality conditions of the lower-level problem are added as the constraints to the upper-level problem, which reduces the problem to a single mixed integer linear program.

The model makes the following assumptions: (i) the loads are modeled as price anticipators; (ii) the demand model representing an aggregator consists of a large population of prosumers connected to an unconstrained distribution network who collectively maximize self-consumption; (iii) aggregators do not alter the underlying power consumption of the prosumers; and (iv) prosumers have smart meters equipped with home energy management systems for scheduling of the PV-battery systems, and, a communication infrastructure is assumed that allows a two-way communication between the grid, the aggregator and the prosumers. More details can be found in <cit.>.

Prosumer Objective function: 
Prosumers aim to minimize electricity expenditure:

    p_p^g+/–, p_p^bminimize∑_t∈𝒯^ p_p,t^g+ - λ p_p,t^g-,

where λ is the applicable feed-in price ratio. In our research, we assumed λ = 0, which corresponds to maximization of self-consumption.

The prosumer sub-problem is subject to the following constraints:

Prosumer power balance: 
Electrical consumption of prosumer p, consisting of grid feed-in power, p_p,t^g-, underlying consumption, p_p,t^, and battery charging power, p_p,t^b, is equal to the power taken from the grid, p_p,t^g+, plus the power generated by the PV system, p_p,t^pv:

    p_p,t^g+ + p_p,t^pv =  p_p,t^g- +  p_p,t^ + p_p,t^b.


Battery charge/discharge limits: 
Battery power should not exceed the charge/discharge limits:

    p_p^b-≤p_p,t^b≤p_p^b+,
	
where p_b^- and p_b^+ represent the maximum power discharge and charge rates of the prosumer's battery, respectively.		

Battery storage capacity limits: 
Energy stored in a battery of prosumer p should always be less than its capacity:

    e_p^b≤e_p,t^b≤e_p^b.


Battery SOC limits:
Battery SOC is the sum of the  power inflow and the SOC in the previous period:

    e_p,t^b = η_p^be_p,t^b +  p_p,t^b,      t ≠ 1, 
       e_p,t^b = η_p^bê_p^b  +  p_p,t^b,     t=1,
		
where ê_p^b represents the initial SOC and is used to establish the connection between adjacent sub-horizons.





§ SIMULATION SETUP

The case studies provided in this section compare the computational efficiency of the proposed MST with alternative formulations. For detailed studies on the impact of different technologies on future grids, an interested reader can refer to our previous work <cit.>.




 §.§ Test System

We use a modified 14-generator IEEE test system that was initially proposed in <cit.> as a test bed for small-signal analysis. The system is loosely based on the Australian National Electricity Market (NEM), the interconnection on the Australian eastern seaboard. The network is stringy, with large transmission distances and loads concentrated in a few load centres. Generation, demand and the transmission network were modified to meet future load requirements. The modified model consists of 79 buses grouped into four regions, 101 units installed at 14 generation plants and 810 transmission lines.




 §.§ Test Cases

To expose the limitations of the different UC formulations, we have selected a typical week with sufficiently varying operating conditions.
Four diverse test cases with different RES penetrations are considered. 
First, RES0 considers only conventional generation, including hydro, black coal, brown coal, combined cycle gas and open cycle gas. The generation mix consists of 2.31 hydro, 39.35 of coal and 5.16 of gas, with the peak load of 36.5. To cater for demand and generation variations, 10 reserves are maintained at all times. The generators are assumed to bid at their respective short run marginal costs, based on regional fuel prices <cit.>.

Cases RES30, RES50, RES75 consider, respectively, 30, 50 and 75 annual energy RES penetration, supplied by wind, PV and CST. Normalized power traces for PV, CST and wind farms (WFs) for the 16-zones of the NEM are taken from the AEMO's planning document <cit.>. The locations of RESs are loosely based on the AEMO's 100% RES study <cit.>. 




 §.§ Modeling Assumptions

Power traces of all PV modules and wind turbines at one plant are aggregated and represented by a single generator. This is a reasonable assumption given that PV and WF don't provide active power reserves, and are not limited by ramp rates, MUDT, and startup and shutdown costs, which renders the information on the number of online units unnecessary.

Also worth mentioning is that RES can be modeled as negative demand, which can lead to an infeasible solution. Modeling RES (wind and solar PV) as negative demand is namely identical to preventing RES from spilling energy. Given the high RES penetration in future grids, we model RES explicitly as individual generators. 
Unlike solar PV and wind, CST requires a different modeling approach. Given that CST is synchronous generation it also contributes to spinning reserves and system inertia. Therefore, the number of online units in a CST plant needs to be modeled explicitly.

An optimality gap of 1% was used for all test cases. Simulation were run on Dell OPTIPLEX 9020 desktop computer with Intel(R) Core(TM) i7-4770 CPU with 3.40 clock speed and 16B RAM. 




§ RESULTS AND DISCUSSION

To showcase the computational efficiency of the proposed MST, we first benchmark its performance for different horizon lengths against the BUC formulation employing three binary variables per unit per time slot and the AGG formulation where identical units at each plant are aggregated into a single unit, which requires three binary variables per plant per time slot. 
We pay particular attention to the techniques used for computational speedup, namely unit clustering, rolling horizon, and constraint clipping. Last, we compare the results of the proposed MST with BUC and AGG formulations for voltage and frequency stability studies.




 §.§ Binary Unit Commitment (BUC)
 
We first run the BUC for horizon lengths varying from one to seven days, Fig. <ref> (top).

As expected, with the increase in the horizon length, the solution time increases exponentially. For a seven-day horizon, the solution time is as high as 25000 (7). Observe how the computational burden is highly dependent on the RES penetration. The variability of the RES results in an increased cycling of the conventional thermal fleet, which increases the number of on/off decisions and, consequently the computational burden. In addition to that, a higher RES penetration involves an increased operation of CST. This poses an additional computational burden due to the decision variables associated with TES that span several time slots.
In summary, the computational burden of the BUC renders it inappropriate for scenario analysis involving extended horizons.




 §.§ Aggregated Formulation (AGG)

Aggregating identical units at a power plant into a single unit results in a smaller number of binary variables, which should in principle reduce the computational complexity. 
Fig. <ref> confirms that this is mostly true, however, for RES50-HL7 the computation time is higher than in the BUC formulation. The reason for that is that, in this particular case, the BUC formulation has a tighter relaxation than the AGG formulation and, consequently, a smaller root node gap. Compared to the MST formulation, with a similar number of variables than the AGG formulation, the MST has considerably shorter computation time due to a smaller root node gap. 
 
In terms of accuracy, the AGG formulation works well for balancing studies <cit.>. On the other hand, the number of online synchronous generators in the dispatch differs significantly from the BUC, which negatively affects the accuracy of voltage and frequency stability analysis, as shown later. Due to a large number of online units in a particular scenario, a direct comparison of dispatch levels and reserves from each generator is difficult. Therefore, we compare the total number of online synchronous generators, which serves as a proxy to the available system inertia. Fig. <ref> shows the number of online generators of four different RES penetration levels for a horizon length of seven days. For most of the hours there is a significant difference between the number of online units obtained from the BUC and the AGG formulation.     


In conclusion, despite its computational advantages, the AGG formulation is not appropriate for stability studies due to large variations in the number of online synchronous units in the dispatch results. In addition to that, the computational time is comparable to the BUC in some cases.


We now evaluate the effectiveness of the techniques for the computational speedup.



  §.§.§ Unit Clustering

In unit clustering, binary variables associated with the generation unit constraints are replaced with a smaller number of integer variables, which allows aggregating several identical units into one equivalent unit, but with the number of online units retained. This results in a significant reduction in the number of variables and, consequently, in the computational speedup. Compared to the BUC, the number of variables in the MST with this technique alone reduces from 24649 to 5990 for RES75 with a horizon length of seven days. Therefore, the solution time for RES75-HL7 reduces from 25000 in the BUC to 450 in MST with unit clustering alone.



  §.§.§ Rolling Horizon Approach

A rolling horizon approach splits the UC problem into shorter horizons. Given the exponential relationship between the computational burden and the horizon length, as discussed in Section <ref>, solving the problem in a number of smaller chunks instead of in one block results in a significant computational speedup. The accuracy and the consistency of the solution are maintained by having an appropriate overlap between the adjacent horizons. However, the overlap depends on the time constants of the problem. Long term storage, for example, might require longer solution horizons. The solution times for different RES penetrations are shown in Table <ref>. Observe that in the RES75 case, the effect of rolling horizon is much more pronounced, which confirms the validity of the approach for studies with high RES penetration.













  §.§.§ Constraint Clipping

Eliminating non binding constraints can speedup the computation even further. Table <ref> shows the number of constraints for different scenarios with and without constraint clipping. Observe that the number of redundant constraints is higher in scenarios with a higher RES penetration. The reason is that a higher RES penetration requires more flexible gas generation with ramp rates shorter than the time resolution (one hour in our case). Note that the benefit of constraint clipping with a shorter time resolution will be smaller.





 §.§ MST Computation Time and Accuracy
  
The proposed MST outperforms the BUC and AGG in terms of the computational time by several orders of magnitude, as shown in Fig. <ref> (bottom). The difference is more pronounced at higher RES penetration levels. For RES75, the MST is more than 500 times faster than the BUC. In terms of the accuracy, the MST results are almost indistinguishable from the BUC results, as evident from Fig. <ref> that shows the number of online synchronous units for different RES penetration levels. Minor differences in the results stem from the nature of the optimization problem. Due to its mixed-integer structure, the problem is non-convex and has therefore several local optima. Given that the BUC and the MST are mathematically not equivalent, the respective solutions might not be exactly the same. The results are nevertheless very close, which confirms the validity of the approach for the purpose of scenario analysis. The loadability and inertia results presented later further support this conclusion.
 



 §.§ Stability Assessment

To showcase the applicability of the MST for stability assessment, we analyze system inertia and loadability that serve as a proxy to frequency and voltage stability, respectively. More detailed stability studies are covered in our previous work, including small-signal stability <cit.>, frequency stability <cit.>, and voltage stability <cit.>.



  §.§.§ System inertia

Fig. <ref> (bottom) shows the system inertia for the BUC, AGG and the proposed MST, respectively, for RES0. Given that the inertia is the dominant factor in the frequency response of a system after a major disturbance, the minuscule difference between the BUC and the MST observed in Fig. <ref> validates the suitability of the MST for frequency stability assessment. The inertia captured by the AGG, on the other hand, is either over or under estimated and so does not provide a reliable basis for frequency stability assessment.



  §.§.§ Loadability Analysis

The dispatch results from the MST are used to calculate power flows, which are then used in loadability analysis[The loadability analysis is performed by uniformly increasing the load in the system until the load flow fails to converge. The loadability margin is calculated as the difference between the base system load and the load in the last convergent load flow iteration.]. Fig. <ref> (top) shows loadability margins for the RES0 scenario for different UC formulations. Observe that the BUC and the MST produce very similar results. The AGG formulation, on the other hand, gives significantly different results. From hours 95 to 150, in particular, the AGG results show that the system is unstable most of the time, which is in direct contradiction to the accurate BUC formulation. 
Compared to the inertia analysis, the differences between the formulations are much more pronounced.
Unlike voltage, frequency is a system variable, which means that it is uniform across the system. In addition to that, inertia only depends on the number of online units but not on their dispatch levels.
Voltage stability, on the other hand, is highly sensitive both to the number of online units and their dispatch levels, which affects the available reactive power support capability, as illustrated in Fig. <ref>. 
Close to the voltage stability limit, the system becomes highly nonlinear, so even small variations in dispatch results can significantly change the power flows and, consequently, voltage stability of the system. One can argue that in comparison to BUC  the proposed MST result in the more conservative loadability margin, although this is not always the case (around hour 85, the MST is less conservative).






§ CONCLUSION

This paper has proposed a computationally efficient electricity market simulation tool based on a UC problem suitable for future grid scenario analysis. The proposed UC formulation includes an explicit network representation and accounts for the uptake of emerging demand side technologies in a unified generic framework while allowing for a subsequent stability assessment. We have shown that unit aggregation, used in conventional planning-type UC formulations to achieve computational speedup, fails to properly capture the system inertia and reactive power support capability, which is crucial for stability assessment. To address this shortcoming, we have proposed a UC formulation that models the number of online generation units explicitly and is amenable to a computationally expensive time-series analysis required in future grid scenario analysis. To achieve further speedup, we use a rolling horizon approach and constraint clipping.

The effectiveness of the computational speedup techniques depends on the problem structure and the technologies involved so the results cannot be readily generalized. The computational speedup varies between 20 to more than 500 times, for a zero and 75% RES penetration, respectively, which can be explained by a more frequent cycling of the conventional thermal units in the high-RES case. The simulation results have shown that the computational speedup doesn't jeopardize the accuracy. Both the number of online units that serves as a proxy for the system inertia and the loadability results are in close agreement with more detailed UC formulations, which confirms the validity of the approach for long term future grid studies, where one is more interested in finding weak points in the system rather than in a detailed analysis of an individual operating condition.






IEEEtran

