How Price-Based Frequency Regulation Impacts Stability in Power Grids: A Complex Network Perspective

.


Introduction
Climate change mitigation requires us to transform current energy systems, a major shift over the next few decades in how energy is produced and transferred in a cleaner way [1]. In modern deregulated power grids, electricity markets provide substantial benefits for energy conservation and emission reduction, especially with the help of various means of advanced metering and communications [2,3]. Particularly, with the real-time based electricity price in electricity markets, competitive activities and thus abundant interactive phenomena are introduced into physical power networks [4]. Due to such complex interactions and low efficiency of analysis in large-scale grids, however, power system engineers and operators are still prone to getting trouble in a systematic analysis of interactive phenomena, especially in system stability analysis [5], although the power engineering community has made tremendous efforts to address these problems for years.
In the fields of complex systems and complex networks, considerable attempts have also been devoted to an efficient analysis of rich interactive and collective phenomena [6,7] in power grids. Various models of power systems have increasingly been brought to higher attention recently, partially thanks to the better understanding of collective phenomena such as synchronization of power grids [7][8][9][10][11]. In particular, the Kuramoto model as a paradigm for general oscillator models has been intensively investigated recently [7], from the perspective of networked control [10], selforganized synchronization [12][13][14][15][16][17], to stability against even large perturbations [18,19]. e further investigation of the Kuramoto model with or without inertia attracts great interest from the aspects of, e.g., the optimal placement of virtual inertia for system performance [20] and for perfect synchronization [21,22], non-Gaussian frequency fluctuations described by Levy-stable laws [23]. From the viewpoint of power grid operators, it is vital yet challenging to enhance network stability in complex power grids, especially in the presence of electricity markets [4]. In particular, how the interplay between networked dynamics and real-time price regulations affects power grid stability still remains unanswered.
Here, to cope with these problems, we adopt the typical governing dynamics of power systems, the second-order Kuramoto model but additionally including generation/consumption adaption to electricity price regulation.
is is intended to maximize energy conservation and respective economic benefits and to quantify impacts of price-based dynamics on the network against even large perturbations. Our experimental results first in a low-dimensional system indeed suggest that the frequency adaption to price regulation can enhance system's stability, but the adaption can also reduce its stability. Next, we show that these findings remain valid in the real world-based China Southern Power Grid (CSG) where the existing thermal, hydro, and nuclear generation and various typical consumers are included. Furthermore, we find that, after duly decreasing the percentage of thermal generation and accordingly increasing hydrogeneration's percentage by improving economic efficiency to meet the total load demands, not only the environmental contamination is mitigated but the network also remains fairly stable.

Governing Dynamics.
In the context of power grids, electricity markets play a crucial role in demand side response and economic dispatch of generation. Given an appropriate governing mechanism of the real-time electricity price in the electricity market, both generators and loads are inspired to adjust their generation and consumption activities accordingly, resulting in an energy conservation and an emission reduction, as well as possible enhancement of system stability. Taking into account both physical network dynamics and market effects, all the dynamics of price regulations, node responses, and networked interactions are expressed as follows.
e time-varying price is mainly designed for frequency regulation. Supposing the price at a given time t is p(t), its control law follows [24]: where p N denotes the rated price in a steady state (without frequency deviation); K P and T l are the main parameters of the control law; and e(t) represents the average frequency deviation of the whole network, serving as the feedback signal of the overall network dynamics. In particular, e(t) is estimated as [25] where L f is the linear coefficient relating the amount of active power imbalance of the whole network to the average frequency deviation and P gi and P dj (with implicit subscript t) stand for the generated and consumed active power by the ith generator and the jth load, respectively. L f is calculated as where f N is the rated frequency (50 Hz in China and Europe and 60 Hz in the USA), M is the average inertia of all the generators, and S N � P gi is the total rated generation capacity of the whole network in the steady state. For a certain power supplier, denoted as i, given the current output P gi , the rate of expansion and/or reduction of its production is adjusted in proportion to the difference between the real-time price p(t) and its marginal cost p gi [26]. Such a response behavior yields where τ gi represents the time constant of the response of the supplier i. Qualitatively, τ gi is determined by an array of elements, including technical features, the working status for the time instant of receiving the price signal, and the means of communication of these signals [26]. Especially, technical features, primarily involving unit ramp rates, seem to have a main effect on the response time.
Intuitively, marginal cost p gi corresponds to the increment of cost with one-unit increase of the output P gi . Given the cost function of a generator, denoted by C(P gi ), p gi is estimated by the derivative of C(P gi ). In terms of thermal generation, supposing its output is P gt , the corresponding cost C(P gt ) can be simply approximated as a quadraticfunction form [27,28]: where α t , β t , and c t are the coefficients of the function. en, the marginal cost p g (P gt ) is computed by where b gt and c gt denote the constant and monomial coefficients of the linear marginal-cost function, respectively. erefore, considering the generation cost with a quadratic-function form [27,28] (see supplementary note (SN) 1), the marginal cost is determined by p gi � b gi + c gi P gi . e consumer j responds to the time-varying price p(t) in a similar manner as above, reducing its consumption once p(t) exceeds its marginal benefit under the current consumption P dj . Specifically, it is described as where τ dj and p dj correspond to its time constant of response and marginal benefit. Given a quadratic function of the consumption benefit [27,28], analogous to the marginal cost p gj , the consumer's marginal cost is computed by We use the classical power grid model [9,12,19,29,30], in coarse scale, where each oscillator i is characterized by generation P i � P gi or consumption P i � − P di , and its governing dynamics follows: 2 Complexity where θ i and ω i are the instantaneous phase and frequency of the ith oscillator, α is the viscous damping coefficient, and the line capacity K ij � K ji accounts for topological connections, with K ij > 0 if nodes i and j are connected and K ij � 0, otherwise.
Interactions between different kinds of dynamics are illustrated in the schematic diagram (Figure 1(a)). Particularly, generations P gi (3) and consumptions P di (6) are intercoupled with price regulation (1). Equations (1)-(6) can affect the node dynamics (7), but not vice versa. Rich phenomena arise due to the competition between the generated power P i and the synchronization forcing term summing neighbors' power delivery K ij sin(θ i − θ j ) along the transmission lines [12,19,30]. Functional reasons require that the limit of the transmission capacity K ij ≫ P i , which is though uneconomical. However, with the rapid growth of power demands these years, P i is expected to be extremely near but less than K ij to avoid costly transmission system expansion. Grids therein are operated at high risk of instability and even inducing blackout when suffering from large disturbances.

Two-Node Model.
As an illustration to reveal the underlying phenomena and mechanisms, we first consider a simple two-node model consisting of one generator (denoted by subscript g) and one consumer (denoted by subscript d). eir phase difference Δθ � θ g − θ d and frequency difference Δω � ω g − ω d satisfy from (7): where ΔP � P g + P d and in the stationary state. e governing dynamics of ΔP is given by equations (1)- (6). Under special assumptions on parameter values with c g � − c d � c and τ g � τ d � τ, it can be simplified by summing equations (3) and (6), which yields where C is solely determined by a perturbation and r ≡ c/τ quantifies the reactability of the system on perturbations and reflects how fast the power goes back to the normal state, i.e., the synchronization rate. For generation, r indicates the comprehensive impacts of the marginal cost and response time constant on the generator's synchronization stability: the larger c is, the more sensitive the generator becomes for power fluctuations or perturbations; the smaller τ is, the faster the generator responds for generation adjustment.
With the increase of r, effects of the electricity price on the node dynamics become more and more pronounced. r is independent of the stationary state, but its variance influences the changes in the size of the attracting basin of synchronization (see Supplementary Figure S2). e formation (9) is of course a strong simplification, yet acting as a cornerstone for the understanding of the intricate dynamics of interactions between price regulation and the physical network.
is reduced system (8) in the stationary state which has one stable fixed point, one saddle, and one stable limit cycle [13,31]. Perturbations affecting the system's dynamics can be induced by changes in generation and/or consumption, or by initial conditions of phase and frequency [32]. According to the definition of power system stability [5], the stability issue of power networks is related to physical disturbances, as load/generation variation, generator tripping, transmission line tripping, short circuit, etc. As a rule of thumb, power network dynamics emerge once it subjects to such perturbations, whereby the fluctuations of all the nodes' angular phase and speed occur.
us, imposing physical perturbations instead of directly changing initial conditions of phase or frequency seems more natural and reasonable. Although [19] attempts to imitate short-circuit disturbances, it is still considerably simplified by ignoring the changes of initial conditions of most dimensions. us, in the sequel, all the perturbations take the form of setting step changes to power consumption or generation, affecting all-dimension dynamics simultaneously.
For small disturbances, the coupling term can be linearized in the vicinity of the equilibrium point [32]. Its maximal Lyapunov exponent is expected to be negative to stabilize the dynamics of price regulations, node responses, and networked dynamics. In this case, the equilibrium point is locally stable and hence the system will return back to the normal state after a certain perturbation. For nonsmall disturbances, depending on whether located in the attracting basin of the equilibrium point, the system will run normally or will go out of power. In what follows, we focus on large perturbations and quantify the power grid stability with the interpretation as the likelihood of the system's returning back to the normal state with respect to large perturbations (see Section 4).
Initially, system (8) is in the equilibrium state. As shown in Figures 1(b) and 1(c), at time t � t 0 , it is then suddenly subjected to a physical disturbance, and the generation P g (in red) is perturbed far away from the normal state. After the sudden perturbation, the price p regulates according to the linear correlation (1), and the consumption P d is forced to adjust its difference to generation, which can be taken as the second disturbance to the networked dynamics (7). e frequency difference Δω is successively driven away from 0. Depending on whether the perturbations are located inside or outside of the attracting basin of the synchronous state, the system will either return back (in Figure 1(b)) or will start wandering periodically to the wanted synchronous state (in Figure 1(c)).

Complexity
To quantify the effects of electricity markets on power system stability, without varying the location of equilibrium points, we calculate the system's antidisturbance capability (ADC) with respect to the synchronization rate r. As the stationary state of systems (1)-(6) is τ-independent but c-dependent, r is varied via the reaction time τ. ADC vs. r is shown in Figure 2(a). e projections of two trajectories consisting of one returning back to the normal state and the other oscillating periodically in the phase space (Δθ, Δω) are shown in the inset. Simulations are conducted with T � 400 perturbations for P g randomly selected within the range of [− 20, 20] at each r. System stability is approximately estimated according to equation (11) with the standard error e < 0.03. Intuitively, strengthening electricity markets might enhance system's stability.
While given C/r constant, various transitions could occur depending crucially on the value of r. For high values of r, the system exhibits robustness against perturbations and has high volumes (ADC) of basins of attraction. If starting from large values of r, ADC increases smoothly and then gradually becomes globally stable with ADC � 1 for further increases of r. For small values of r, synchronization rate decreases and this induces the shrinking of the attraction basin, as shown in Figure 2(b). Interestingly, if starting from a rather small value of r, ADC first decreases and approaches a minimum value, and only then it increases smoothly (for detailed interpretation, see SN 2). is can also be reflected by the emerging red region in the low-left corner of Figure 2(b).
is also indicates that enhancing market effects may undermine power grid stability under certain circumstances, which needs to be avoided for both power network planning and operation.
ese results are obtained from the lowdimensional system, the simple two-node model of a strong simplification. However, can similar results be observed in high-dimensional systems?   (8), trajectories of the price p, generation P g , and consumption P d are shown in (b) (resp. in (c)) given a perturbation suffered at time t 0 located inside (resp. outside) the attracting basin of the synchronous state. Parameters are given as K P � 50, T l � 20,  (3) and (6) yet with different characteristics (see Tables 1 and 2), except for nuclear plants which are independent of price regulations. As stated in Section 4, each node stands for an aggregated plant or consumer, corresponding to a series of actual generators or loads. Under this circumstance, perturbations of P g on individual loads represent sudden tripping of small-scale generators in practice or the total fluctuation of power generation of the corresponding aggregated plant. In what follows, we especially address these two aspects: effects of electricity markets on the grid stability and on environmental pollution.

China Southern Power Grid (CSG
(i) Let us now first revisit generation/consumption dynamics and settings of physical disturbances referred above, so as to quantify the effects of electricity markets on the grid stability ADC as a function of   the synchronization rate r. Synchrony in the grid is achieved when the maximal absolute deviation of the frequencies ω i vanishes in time. e networked model (7) exhibits the essential feature of the system, the coexisting stationary state, as soon as it comprises of the two-node model [12]. After perturbation, individual nodes behave differently. And after some transition, actually only nodes around the birthplace of the perturbation are observably driven away from the synchronous state, while others remain close to it, as depicted in Figure 3(a), where the perturbation is imposed on node 0. Following the process outlined in the previous lowdimensional model, concretely, simulations are conducted via perturbing power generation (P g ) of thermal and hydro plants within the range of [− 20, 20]. Our results suggest that, as depicted in Figure 3(b), variation of ADC with respect to r illustrates qualitatively a similar trend (as in Figure 2(a)): sinking first and then rising subsequently. For each node i, we estimate single-node stability ADC i which expresses the probability of all nodes returning back to the normal state after the ith node has been subjected to small and even large perturbations. Projection of ADC i on each node i with its color strength proportional to 1 − ADC i is shown in Figure 4. Clearly, the low degree nodes (in dark red) are likely to undermine power grid stability [33,34].
(ii) Next we seek to decarbonize the power plants.
Would it be possible to decrease the thermal generation and to increase the hydrogeneration, while simultaneously keeping the power grid still stable and robust? On the assumption that the stationary state of price p remains as a constant, in order to achieve this goal, one might improve the economic efficiency of hydro and thermal generators, e.g., decreasing its marginal-cost coefficient c g by an equipment upgrade and technological improvement.
Assume that both hydro and thermal generation units in CSG are updated with more advanced technology, corresponding to the parameter shift of c g , which are artificially decreased by 30% and 10%, respectively. is induces the reallocation of power flows. As shown in the insets of Figure 4, our results suggest that after improving the economic efficiency, the thermal generation over the total generation is decreased from 55% to 47%, whereas the hydrogeneration's percentage is increased from 39% to 48%. e price remains almost constant as expected. is of course further affects the system's performance, especially its stability against large perturbations. erefore, we re-estimate the single-node stability ADC i ′ for all nodes listed in Supplementary Tables S1 and S2. We find that few nodes' stability is even enhanced (see green-colored values in Supplementary Table S1), while some nodes' stability is decreased (in red). But the whole power grid still remains fairly stable. erefore, on the whole, both tasks of decarbonization of power production and maintenance of system capability are accomplished in this way. However, the results also reveal that the single-node stability of some individual nodes is somewhat decreased. It is crucial to identify overloaded nodes to avoid cascading failures [35,36]. In practice, the update process of equipment on these nodes is suggested to be avoided or deferred, so that a better trade-off between the decarbonization and system stabilization can be achieved.

Discussion
We have investigated effects of real-time price-based frequency regulation on power grid stability against even large perturbations. Specifically, in the coexistence regime of the two-node model and in a linear-stable regime of the markets (see SN 2), we have found that the system stability can be undermined via increasing the control of the electricity price mechanism.
e results remain valid in the rather large China Southern Power Grid with more practical scenarios, composing of three types of generators (thermal, hydro, and nuclear plants) and of consumers (residential, commercial, and industrial loads). Of course, power grid stability analysis is rather complicated and can be influenced by various aspects, e.g., parameter settings, ranges of perturbations, and network topologies. erefore, we have validated the results in the two-node model with an even larger range of perturbations (see SN 2), in the northern European power grid (see SN 3) and additionally, in the China Southern Power Grid with a different power allocation with a single type of generator and of consumer (see SN 4). All these simulations exhibit the local minimum phenomenon of system antidisturbance capability with a continuous increase of the synchronization rate.
Generators react accordingly with different characteristics and have different environmental impacts. Given the specific scenario in the China Southern Power Grid, we have demonstrated a simple but feasible approach, involving updates of hydro and thermal generation units, to achieve the dual target of decarbonization and stability maintenance. Our experiments confirm the possibility of reducing the thermal production's percentage and increasing the percentage of hydroproduction, and meanwhile the grid remains fairly stable afterwards.
What remains to be investigated next? Primarily, renewable energy generation, including wind power and solar energy, and more realistic models of consumers instead of synchronous motors deserve further research efforts [37]. Additionally, several efforts should be taken into account, e.g., more detailed network structures and the relationship between the overall decarbonization benefits and the system's basin of attraction.

Methods and Materials
4.1. Stability Measure. We take the following general approach to quantify system antidisturbance capability (ADC), which is intended to evaluate how stable the synchronous state is against large physical perturbations for generations and/or consumption P, especially defined as ADC � χ(P)ρ(P)dP, (10) where χ(P) is an indicator function with χ(P) � 1, if P belongs to the attracting basin of the equilibrium point and χ(P) � 0, otherwise; ρ is a density function normalized as ρ(P)dP � 1. Its magnitude suggests how likely the system comes back to the stable fixed point after suffering from large perturbations. Numerical Monte Carlo methods are employed to estimate the likelihood. In particular, initially the system is synchronized and we draw T random initial conditions according to the density function ρ(P), solve the system by numerical integration, and then count the number of initial conditions driving the system back to the synchronous state. is yields with standard error e � ����������������� (ADC(1 − ADC))/T. Note that perturbations are imposed physically on power production or demand of individual nodes instead of the changing phase and frequency as before in [19].

Topology of China Southern Power Grid. China Southern
Power Grid (CSG) is one of the two large-scale power grids in China, consisting of Yunnan, Guizhou, Guangxi, Guangdong, and Hainan provincial subnetworks. As shown in Figure 4, the main topology of the bulk alternating-current (AC) electric network of CSG is extracted for this case study. is real-world network preserves the primary structure and characteristics of the whole system. On the basis of the planning data of the Dispatching Center of CSG in 2014, it is simplified by the following equivalences. Similar generators in one location are equivalently represented by a large generator [38]; to avoid trivial and complicated network modeling, all the low-voltage distribution networks are aggregated into corresponding high-voltage substations. e final reduced network includes 57 generators and 111 loads, connected by 194 edges. In particular, the generators are mainly comprised of hydro, thermal, and nuclear plants, whose deployments are based on actual situations. It is because these three forms of generations still occupy a major proportion of the total power generation in southern China, though an increasing penetration of wind power and photovoltaic generation has been witnessed in the last decade. e loads also involve industrial, commercial, and residential categories, and, from the perspective of high-voltage grids, each load node contains all the three.

Parameters Setting.
Within the three categories of generation in CSG, thermal and hydro generators are likely to respond to the real-time price, while the outputs of nuclear ones are not affected by price at any time, due to safety consideration and technical restriction. In terms of thermal plants, owing to the inertia of thermodynamic systems and technical requirements of mild output, their ramp rates are relatively limited; thus, their τ gi is fairly large. On the contrary, hydroplants are permitted to operate in a more flexible way, which means that a smaller τ gi is approachable. Taking into account the economic cost and environmental expense in China [39], the monomial coefficient c gi of thermal generation is supposed to be about 1.5 ∼ 3 times larger than that of the hydrocounterpart. Given specific c gi , b gi is determined by the rated price p N and output P giN at the stationary status. Concrete ranges of τ gi , b gi , and c gi are summarized in Table 1 [24,39,40].
In the light of practical experiences, concrete ranges of parameters of the 3 typical kinds of consumers are listed in Table 2 supplementary table S2). e two pie charts show the percentage of the three kinds of generator plants before and after improving the economic efficiency of hydropower plants by decreasing their c g from 1 to 0.7 and decreasing thermal c g from 1 to 0.9. Parameters are the same as in Figure 3 with τ g � 3.0 for thermal plants and τ g � 0.2 for hydroplants.
of the production process should be guaranteed first, and so as to ensure the quality of its products, it is not that willing to change its amount of power consumption even given a certain degree of marginal benefit. However, residents care more about the varying price, and they are most likely to adjust their consumption once receiving the latest price information, aiming at gaining larger benefits. As for shopping malls or supermarkets, sometimes (during the slack period) they tend to react to the varying price timely, while sometimes (during the peak period) not. erefore, from a comprehensive point of view, residential consumers are deemed to be most sensitive to price variation, commercial ones are relatively mediate, whereas industrial users are least sensitive.

Data Availability
e China Southern Power Grid and other data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.

Authors' Contributions
P.J. and L.Z. designed the study and prepared the manuscript. L.Z. prepared the data. W.L. and P.J. carried out the analysis. All authors discussed the results and contributed to editing the manuscript. J.K., W.L., and C.L. supervised the study.