Control of active power for synchronization and transient stability of power grids
Abstract
As the energy transition transforms power grids across the globe it poses several challenges regarding grid design and control. In particular, high levels of intermittent renewable generation complicate the job of continuously balancing power supply and demand, which is necessary for the grid’s stability. Although there exist several proposals to control the grid, most of them have not demonstrated to be cost efficient in terms of optimal control theory. Here, we mathematically formulate the control problem for stable operation of power grids, determining the minimal amount of control in the active power needed to achieve the constraints, and minimizing a suitable cost function at the same time. We investigate the performance of the optimal control method with respect to the uncontrolled scenario and we compare it to a simple linear control case, for two types of external disturbances. Considering case studies with two and five nodes respectively, we find that the linear control can improve the synchronization and transient stability of the power grid. However, if the synchronized angular velocity after a disturbance is allowed to differ from its initial steady state value, the linear control performs inefficiently in comparison to the optimal one.
Index terms— Optimal control, Power system control, Power system dynamics, Power system stability
I Introduction
The electrical power grid is undergoing drastic changes due to the energy transition [1, 2, 3] and sophisticated control approaches are necessary to ensure a reliable and stable operation [4]. The generation side of the grid is changing as additional renewable generators are installed to mitigate climate change, introducing fluctuations in a time scale of days [5] to subseconds [6]. In addition, the demand side is changing due to the ongoing electrification of heating and transport [7] and the introduction of demand control [8]. Regardless of these changing conditions, the grid needs to stay within strict operational boundaries to guarantee a stable electricity supply and to prevent damage to sensitive electronic devices [4].
A fundamental aspect of power system stability is the ability of interconnected synchronous machines of a power system to remain synchronized [4]. Transient stability is the ability to maintain synchronism in the face of severe transient disturbances [4], and is of great importance in preventing cascading failures [9]. Control mechanisms that balance active power and regulate frequency in the grid are key to maintaining these stability conditions. Primary controls respond within a few seconds of an event to stabilize the frequency within its permissible operating limits, after which secondary and tertiary controls restore the frequency to its nominal value [10]. Schäfer et al [11, 12] recently proposed a decentralized linear local frequency controller that regulates demand and supply on a power grid through economic incentives in order to improve the grid’s transient stability. However, it remains unclear how well and at what cost this linear control ensures that the system’s operational constraints are satisfied. In fact, some decentralized controllers for frequency response are known to be inefficient [13].
In this paper, we formulate an optimal control problem [14, 15] to keep the power grid within operational boundaries when facing known changes in demand and generation. Note that there are several predictable sources of perturbations for which a designated control reaction may be designed. For example, a recent study analysed power grid frequency fluctuations in grids of different synchronous regions and found a significant impact of trading on the frequency quality with large fluctuations every fifteen minutes [16]. Furthermore, there are predictable changes even in volatile renewable generation. For instance, in a system with large amounts of installed solar capacity, the rising or setting sun will cause similar power ramps as a recent solar eclipse caused in today’s system [17]. The control problem we formulate is wellsuited to keep the grid within operational boundaries for such deterministic perturbations.
In the following section we present a third order “fluxdecay” model which we use to describe the power grid dynamics in terms of a state vector of phase angles, frequencies and voltage amplitudes. Let be a suitable set of control variables . For a given , we quantify its cost through a cost function , and evaluate its performance with respect to various operational constraints and their tolerances . The optimal control problem for power grid synchronization is expressed mathematically as follows.
Problem:
(1) where function governs the intrinsic dynamics of the state of the grid, and is the number of nodes in its representation as a network.
Problem (1) is solved numerically using a control parametrization method [18, 19]. We illustrate the efficiency of the optimal control compared to the decentralized local frequency control for the elementary twonode system and a fivenode network motif. In particular, we compute any violations of the operational constraints and quantify the costs of control actions. Finally, we close with a conclusion and outlook.
Ii An optimal control problem for power system synchronization
Notation  Description  Units 

the set of nodes where is the number of nodes in the network  –  
grid reference  
control time horizon  
set of control variables  –  
dimensional coupling matrix  
dimensional state vector  –  
variance of network angular velocities  
mean value of network angular velocities  
rotor angle at node relative to the grid reference  
angular difference between nodes and ,  
angular velocity at node  
controlled power added to node  
normalized machine voltage at node  1  
normalized mechanical power at node  
normalized voltage reference at node  
damping parameter at node  
,  machine parameters at node  , 
, ,  cost functional, constraint functional, constraint tolerance  1 
Iia Dynamics for transient stability analysis
Before focusing on optimal control, we need to discuss how to describe the intrinsic dynamics of the power grid. According to [4, p. 27], deviations in rotor angles may also lead to a progressive drop in the nodal voltages. This in turn affects the angular velocities and rotor angle values. A realistic model of the power grid should therefore take into account the effect of the rotor angles’ deviations on the voltage amplitudes. This allows us to analyse slower phenomena such as large deviations in voltage or frequency, as typically done in midterm stability studies [4, p. 34]. Therefore, in this paper we use the model in [20] which describes the power grid as a network of synchronous generators and motors governed by a set of differential equations for the rotor angle , angular velocity , and voltage at each node,
(2) 
Table I describes the notation used in this model. In particular, , and the entries of matrix describe the coupling topology of the network. Positive values of indicate net generation at node , whilst negative values of indicate net consumption. If is always positive (respectively, negative) then node is called a generator (respectively, consumer or motor) node. We refer to positive changes in the value of the control variable as incremental actions [21] since they correspond to an increase in generation or an equivalent decrease in demand. Similarly, we refer to negative changes in the value of as decremental actions [21] since they correspond to a decrease in generation or an equivalent increase in demand.
IiB Decentralized control for improved power grid stability
Decentralized controllers for largescale and complex networks such as the power grid are popular choices since they rely only on local measurements and typically have low computation and communication complexity [13]. Three decentralized controllers that appear in the literature on power grid synchronization and stability are mentioned below. They all assume that at initial time the power grid, as described by the dynamics (2), is synchronized and at a steady state with in particular.
Linear local frequency controller
The basic modelling assumption for the linear local frequency (LLF) controller in [11, 12] is that a generator or consumer node acts in response to its nodal price of electricity, which is a linear function of the angular velocity at that node. Upon linearising the response in power output at node with respect to the nodal price, the authors of [11, 12] derive a control that is directly proportional to the angular velocity ,
(3) 
where measures the willingness at node to change the active power level. This type of control effectively results in the increase of the damping parameter from to in (2). Under the assumption of constant voltages, , the authors show that by achieving an equilibrium between demand and supply through market price signals, their model can ensure synchronization and stability of the power grid over time.
Integral local frequency controller
In [22, 23] the following integral frequency control is proposed,
where . Since by (2) above, the integral frequency control can be described equivalently by,
showing that it is directly proportional to the rotor angle’s deviation from the steady state value . In the case , this control can improve transient stability and be economically efficient in a particularly defined sense (see [22, 23]). A similar statement for the general dynamics (2) with timevarying voltages is not known [24].
Bandgap local frequency controller
In [25], the following bandgap control is proposed,
where and are parameters. In this case, control is applied at node only if the angular velocity deviation is sufficiently high, , and ramps up until a maximum value which is reached when . Under the assumption of constant voltages, , the authors show that the bandgap control can improve grid stability.
IiC Operational constraints of the power grid
Let denote the dimensional controlled state variable obtained from (2) with components given by
(4) 
The dynamics of in (2) can be written compactly as
(5) 
where expressions for the components of are obtained from (2) using the assignment given in (4). Each component of the control variable corresponds to the amount of additional active power injected or withdrawn at an individual node in the network. We assume that controls are bounded: for each we have where:
(6) 
Let denote the set of all such control functions.
Synchronization
Synchronization of the power grid occurs when all of the rotor angles have the same velocity [4, p. 19]. Consequently, for every , we should have . Letting denote the vector of angular velocities, the variance of , defined by
(7) 
measures the total dispersion of from its mean , and synchronization is achieved when
(8) 
Let denote the length of the control horizon in seconds. Define the synchronization constraint loss function by
(9) 
and the total synchronization loss on by
where is a weight parameter which emphasizes the relative importance of the constraint at the final time . Note that we weight the variance as given by eq. (7) quadratically in this integral, which means small variances do not contribute much whereas large variances are strongly penalized. Other weighting schemes are also possible.
Mean angular velocity operational limits
The angular velocity quantifies the deviation of the frequency at node from the grid reference . In the United Kingdom and many other countries, the nominal frequency is 50 Hz, which gives a reference value of Hz. For reasons related to the quality of electricity supply, the frequency must respect certain operational limits. Therefore, values of the mean angular velocity should be constrained,
(10) 
In the United Kingdom, the statutory limits are 0.5Hz of the nominal value, but the operational limits are set to the stricter range of 0.2Hz from the nominal value [26]. In our setting, the operational limits on the mean angular velocity become Hz and Hz. Define the mean angular velocity constraint loss function by,
(11) 
and the total loss on for violating this constraint by,
where is a weight parameter. Note that only when is negative in eq. (11) we get a contribution.
Voltage operational limits
Since the voltages in our model are also time dependent, it is important to also take into account appropriate operational constraints on these variables. Regulations in the United Kingdom require that the steady state voltages should be kept within of the nominal voltage for systems between 1 and 132 kV, or of the nominal voltage for systems above 132 kV [27]. In our model we can take this into account with the following constraint,
(12) 
where . We define a loss function for the voltage constraint at each node by
(13) 
and the total loss on for violating this constraint by
where are weight parameters.
IiD Formulation of the optimal control problem
For the total loss is nonnegative, and is equal to zero if, equivalently, the th constraint is satisfied on . We relax this by introducing tolerance parameters , , and say that a control is feasible if it satisfies
(14) 
At an initial time , the power grid is synchronized and at a steady state, , in which various operational constraints are satisfied. If an external disturbance to the mechanical power causes the grid to become unsynchronized, we would like the control to return it close to a synchronized state within seconds with a “minimal cost” that ensures the constraint conditions (14) are satisfied. Letting denote the identity matrix, we define the following total cost for a control ,
(15) 
where denotes the transpose operator. The rate function in (15), , is the square of the Euclidean norm of the vector ,
Therefore, higher cost is assigned to controls which exert large amounts of effort over time. Moreover, adjustments in demand and generation of the same magnitude are penalized equally due to the symmetry . Note also that quadratic costs are standard in the frequency control literature [24, 23]. Using the total cost (15) and the previously defined constraint losses we formulate the optimal control problem (1).
Existence of a minimising control
We study the case in which the mechanical power is of the form,
where is an dimensional vector of steady state values at the initial time , and is a bounded piecewise continuous dimensional deterministic function that models the external disturbance at the nodes. Let denote the exact control, for all . Since the system is assumed to be at a steady state at time , the control trivially ensures the inequalities (14) are satisfied.
Each constraint function in the optimal control problem (1) represents one of the pure state constraints (8), (10) and (12) in canonical form [18, p. 70]:
(16) 
This formulation is convenient as a penalty function approach for solving the optimal control problem numerically [28, 18]. There is a corresponding maximum principle which provides necessary conditions for a solution to the optimal control problem [14], but does not prove that there exists such a solution. For completeness we provide the following result.
Theorem 1.
In physical terms, Theorem 1 tells us that there is a feasible control strategy that minimizes the cost. Condition (17) is a technical assumption needed to apply a general existence theorem from optimal control theory (see [15, p. 483] for a discussion on this). It is also a reasonable one since we aim to control the power grid within a bounded range of operational values. When all operational constraints are satisfied at all times exactly, and (17) is easily verified from the pure state constraints (8), (10) and (12).
Proof.
First note that the dynamics for given explicitly by (2) are finitely generated in the following sense [15, p. 478]:
(18) 
where , , is a dimensional vectorvalued function. In particular, for we have where is the dimensional vector with 1 at the th coordinate and zeros elsewhere. The optimal control problem also has the following additional properties:

Each function () in (18) is measurable in time , continuous in the state and satisfies,
(19) where is a constant and is the Euclidean norm. For condition (19) is clearly satisfied since . For the growth condition (19) follows from the dynamics (2), boundedness of the mechanical power , and compactness of .

The set of control values described by (6) is compact and convex.

The mapping is continuous, nonnegative and convex.
If there exists a compact set such that the following stronger version of (17) holds:
then the proof is completed by applying Theorem 23.11 in [15]. Otherwise, we proceed by noticing that may be written as,
(20) 
where, in analogy to (18), is finitely generated, measurable in and continuous in . If then (20) clearly holds. If we use the fact that is continuously differentiable then obtain (20) by elementary calculus and using .
Let be the dimensional vector with components for , and for has dynamics:
Notice that and the dynamics for are also finitely generated,
where , , is a dimensional vectorvalued function. We can now verify that there exist compact sets and such that,
and the proof is completed by applying Theorem 23.11 in [15]. ∎
Description of approximation method
The control parametrization method [18, 19] is a numerical procedure based on the maximum principle [14, 15, 28] which gives an approximate solution to the optimal control problem. Let denote a sequence of time points used to partition the control horizon into subintervals, where is an integer. Let denote the subset of controls of the form
(21) 
Each component , , of a control is parametrized by an dimensional vector . The optimal control problem is then approximated by a constrained nonlinear optimization problem over the bounded dimensional parameter space defining the control vector . Standard optimization packages can be used to solve this problem, and we used the Sequential Least Squares Programming (SLSQP) routine in Python. For the numerical simulations we use equidistant partitioning points , , with .
Iii Simulations of controlled networks
For the numerical simulations, we concentrate on the following two case studies. We model a twonode system and the fivenode network motif, shown in Figure (a)a, consisting of two generators at nodes and connected to three motors at the remaining nodes. In Appendix A we list parameter values for the control problem in Table III, and parameter values for the twonode and fivenode systems in Tables IV and V respectively. We consider two types of disturbances which only affect node 1: a persistent one and a temporary one on the interval . For the twonode system, node 1 is perturbed by a disturbance of . For the fivenode system a disturbance of is used as shown in Figure (b)b. The simulated disturbances reflect a sudden increase in demand, or equivalent loss in generation, that can disappear after a few seconds (temporary case) or persist for much longer (persistent case).
Besides the exact and optimized controls, we consider a linear frequency control which is defined by
(22) 
Note that in a synchronized system the same level of control is applied to all nodes. Starting from , which is the uncontrolled case, we simulate the controlled system and check whether the inequalities (14) are satisfied, incrementing the value of by until this becomes true. Table II below shows the values of we obtained.
Temporary Disturbance  Persistent Disturbance  

2 Nodes  15  280 
5 Nodes  30  27 
As discussed in Section IIA, the linear frequency control effectively increases the damping parameter in the intrinsic grid dynamics (2). Comparing the values of in Table II with those for in Tables IV and V, we see that those for are relatively large. Moreover, since the range of feasible angular velocities is , we see that the range of implied values for the linear frequency control for is already , which is very large in comparison to the other active power parameters so that consumers could be required to generate power instead of consuming it.
Iiia Twonode system
We show in Figures 2, 3 and 4 respectively, trajectories for the controlled active power, average voltage, and the angular velocity mean and standard deviation in the twonode system.
Uncontrolled system
The individual voltages and angular velocities both have periodic orbits that lie outside of the operational constraints. For the voltages this is still reflected when averaged across the nodes, as Figures (a)a and (b)b show. Figures (a)a and (b)b show that the system loses its synchronism and, even when the disturbance disappears, the unsynchronized behaviour persists for the entire control horizon. For the angular velocities, however, the oscillations largely cancel out when averaged across the nodes so this behaviour is not reflected in the top panels of Figures (a)a and (b)b.
Linear local frequency control
Trajectories for the linear control applied to the temporary disturbance correspond to the coefficient value , whereas those for the persistent one correspond to . Figure 2 shows that the control initially adds power to node 1, which is nearly equal in magnitude to the disturbance, and eventually the level of additional power at this node decreases whilst the amount added to the node 2 increases. For both disturbances, Figure 4 shows that the angular velocity remains close to its initial value. For the temporary disturbance, the system is already close to being synchronized shortly after the disturbance ends, whereas for the persistent one synchronization improves only slowly over time.
Optimized control
The optimized control initially adds about half the amount of power to node 1 as the linear one whilst noticeably reducing the amount of power at node 2. This results in a slower decrease in as it approaches the lower operational limit (cf. Figure 4). In the temporary case, the system is returned close to its initial steady state values after the disturbance has ended. In the persistent case, once reaches the lower operational limit it is kept there while the disturbance persists. Figure 3 shows that the optimized control takes a similar approach to ensuring the voltage stays within its limits. Notice that in Figure 2 the optimized and LLF controls are qualitatively similar and of the same sign at node 1, but there is a sustained negative level of control at node 2 during the time that there is a positive level for the LLF control.
IiiB Fivenode system
We show in Figures 5, 6 and 7 respectively, trajectories for the controlled active power, average voltage, and the angular velocity mean and standard deviation in the fivenode system.
Uncontrolled system
Linear local frequency control
The synchronization results are similar to the twonode case, with angular velocities kept close to their initial values. The results show that larger networks with higher levels of inertia can provide a more distributed response, thereby requiring lower values of the coefficient . Indeed, for persistent disturbance we used , whereas for the twonode system we used . Consequently, the combined response at the unperturbed nodes leads to lower control costs.
Optimized control
The trajectories are quite similar to those for the twonode system. In particular, as the optimized control changes the net active power, the mean angular velocity follows its natural direction of descent or ascent within the operational limits until a particular level. The angular velocity is then kept at this level whilst the disturbance persists. Furthermore, the combined action at the unperturbed nodes is generally of the opposite type to that taken at the perturbed node. That is, when there is an increase (respectively, decrease) in there is typically a decrease (respectively, increase) in at the same time. Somewhat surprising is the extent to which the LLF and optimized controls agree, at least qualitatively. Differences between the two mainly concern the initial response to the disturbance:

the LLF control causes a combined increase in active power at the unperturbed nodes;

the optimized control causes a combined decrease in active power at the unperturbed nodes.
Another interesting observation is that the controlled power and limiting angular velocities can all be different for the LLF, uncontrolled and optimized cases, yet the limiting voltages across these cases appear to be equal in Figure 6.
IiiC Comparison of control costs
In Figure 8 we show the cost for the linear, optimized and exact controls associated with the trajectories displayed above. While the optimized control satisfies the constraints with smallest cost, we note that its costs are much lower than the linear frequency control’s, partly because it is not forced to keep the system at its initial steady state value. If solutions that are sufficiently close to the initial steady state values are more desirable, this can be taken into account in a variety of ways, including:

Adjusting the constraint costs , , so that deviations from the initial state are penalized explicitly, perhaps increasingly as the remaining time for control, , decreases.

Penalising deviations from the initial state by generalising (15) to
(23) where and,

for an dimensional square matrix and an dimensional vector , is the quadratic form ;

, and are appropriately chosen diagonal matrices with positive definite and , positive semidefinite.

Iv Conclusion and outlook
The linear local frequency control acts as a primary response service to keep the grid frequency close to its nominal value. In order to guarantee stable operation, its control coefficient generally needs to be large. Using the analogy in [11, 12], this means generators and consumers must have a high willingness to adjust their output in response to changes in the reference angular velocity. In terms of power output, this leads to response profiles that almost exactly counteract the disturbance, at least in the initial response phase.
Our results suggest that more efficient controllers distribute the controlled response amongst all nodes in the power grid. Moreover, this response need not be homogeneous throughout the network, but could simultaneously involve incremental actions (net increase in power) at some nodes and decremental ones (net decrease in power) at others.
Trajectories associated with the optimized control show that as it changes the net active power, the mean angular velocity follows its natural direction of descent, or ascent as appropriate, within the operational limits until a point is reached, possibly at the boundary, at which the power grid is synchronized and active power is balanced within the network. For planned events, or events that will occur with very high probability at an anticipated future time, an optimal centralized control such as this one can be used as an initial response, and can be coupled with a simple distributed or decentralized control, such as the linear local frequency control, to smooth out additional unknown perturbations.
Designing decentralized and distributed controllers that are economically efficient is subject of ongoing work (see [13, 23, 29], for instance). Decentralized stochastic control (see [30, 31]), which generalizes our methodology by incorporating uncertainties and different information structures amongst multiple controllers, is likely to become an important theoretical tool for understanding how these controllers work. Its relevance to the stability and control of power systems has been illustrated recently in [32]. While the numerical results presented here to illustrate how the method works were obtained for two simple case studies, namely a 2node and a 5node network, our method can be applied to more realistic and larger networks, for which the problem discussed in our paper is of utmost relevance: controlling the system not by arbitrary means but in the most cost efficient way, such that the operational constraints of the power grid keep on being satisfied.
Appendix A Tables of parameter values
Parameter  Value  Units 

1  
1  
1  
1  
1  
Parameter  Value  Units 

Parameter  Value  Units 

References
 [1] G. Boyle, Renewable Energy: Power for a Sustainable Future, 2^{nd} ed. Oxford, UK: Oxford University Press, 2004.
 [2] F. Ueckerdt, R. Brecha, and G. Luderer, “Analyzing major challenges of wind and solar variability in power systems,” Renew. Energ., vol. 81, pp. 1–10, Sept. 2015. DOI: 10.1016/j.renene.2015.03.002.
 [3] J. A. Turner, “A realizable renewable energy future,” Science, vol. 285, no. 5428, pp. 687–689, Jul. 1999. DOI: 10.1126/science.285.5428.687.
 [4] P. Kundur, Power System Stability and Control, N. J. Balu and M. G. Lauby, eds. New York, NY, USA: McGrawHill, 1994.
 [5] D. Heide, L. Von Bremen, M. Greiner, C. Hoffmann, M. Speckmann, and S. Bofinger, “Seasonal optimal mix of wind and solar power in a future, highly renewable Europe,” Renew. Energ., vol. 35, no. 11, pp. 2483–2489, Nov. 2010. DOI: 10.1016/j.renene.2010.03.012.
 [6] P. Milan, M. Wächter, and J. Peinke, “Turbulent character of wind energy,” Phys. Rev. Lett., vol. 110, no. 13, Mar. 2013, art. no. 138701. DOI: 10.1103/PhysRevLett.110.138701.
 [7] K. Dennis, K. Colburn, and J. Lazar, “Environmentally beneficial electrification: The dawn of emissions efficiency,” The Electricity Journal, vol. 29, no. 6, pp. 52–58, Jul. 2016. DOI: 10.1016/j.tej.2016.07.007.
 [8] P. Palensky and D. Dietrich, “Demand Side Management: demand response, intelligent energy systems, and smart Loads,” IEEE T. Ind. Inform., vol. 7, no. 3, pp. 381–388, Aug. 2011. 10.1109/TII.2011.2158841.
 [9] B. Schäfer, D. Witthaut, M. Timme and V. Latora, “Dynamically induced cascading failures in supply networks,” arXiv preprint, Jul. 2017. arXiv:1707.08018.
 [10] D. Kirschen and G. Strbac, Fundamentals of Power System Economics. Chicester, UK: John Wiley & Sons, 2004.
 [11] B. Schäfer, M. Matthiae, M. Timme and D. Witthaut, “Decentral smart grid control,” New J. Phys., vol. 17, art. no. 015002, Jan. 2015.
 [12] B. Schäfer, C. Grabow, S. Auer, J. Kurths, D. Witthaut and M. Timme, “Taming instabilities in power grid networks by decentralized control,” Eur. Phys. JSpec. Top., vol. 225, no. 3, pp. 569–582, May 2016. DOI: 10.1140/epjst/e201550136y.
 [13] D. K. Molzahn, F. Dörfler, H. Sandberg, S. H. Low, S. Chakrabarti, R. Baldick and J. Lavaei, “A survey of distributed optimization and control algorithms for electric power systems,” IEEE T. Smart Grid, vol. 8, no. 6, pp. 2941–2962, Nov. 2017. DOI:10.1109/TSG.2017.2720471.
 [14] M. R. Hestenes, Calculus of Variations and Optimal Control Theory. New York, NY, USA: John Wiley & Sons, 1966.
 [15] F. Clarke, Functional Analysis, Calculus of Variations and Optimal Control (Graduate Texts in Mathematics Vol. 264). London, UK: Springer, 2013.
 [16] B. Schäfer, C. Beck, K. Aihara, D. Witthaut and M. Timme, “NonGaussian power grid frequency fluctuations characterized by Lévystable laws and superstatistics,” Nature Energy, to be published, 2018. DOI: 10.1038/s415600170058z.
 [17] R.G. Harrison and E. Hanna, “The solar eclipse: a natural meteorological experiment,” Phil. Trans. R. Soc. A., vol. 374, no. 2077, Aug. 2016. DOI: 10.1098/rsta.2015.0225.
 [18] K. L. Teo and C.J. Goh, A Unified Computational Approach to Optimal Control Problems. Essex, UK: Longman Scientific & Technical, 1991.
 [19] V. Rehbockt, K. L. Teo, L. S. Jennings and H. W. J. Lee, “A survey of the control parametrization and control parametrization enhancing methods for constrained optimal control problems,” in Progress in Optimization (Applied Optimization Vol. 30). Boston, MA, USA: Springer US, 1999, pp. 247–275. DOI:10.1007/9781461332855_13.
 [20] K. Schmietendorf, J. Peinke, R. Friedrich and O. Kamps, “Selforganized synchronization and voltage stability in networks of synchronous machines,” Eur. Phys. JSpec. Top., vol. 223, no. 12, pp. 2577–2592, Oct. 2014. DOI: 10.1140/epjst/e2014022098.
 [21] D. Z. Szabó and R. Martyr, “Real option valuation of a decremental regulation service provided by electricity storage,” Phil. Trans. R. Soc. A., vol. 375, no. 2100, Aug. 2017. DOI: 10.1098/rsta.2016.0300.
 [22] C. Zhao, E. Mallada, and F. Dörfler, “Distributed frequency control for stability and economic dispatch in power networks,” in Proc. Am. Control Conf., Chicago, IL, USA, 2015, pp. 2359–2364. DOI: 10.1109/ACC.2015.7171085.
 [23] E. Weitenberg, Y. Jiang, C. Zhao, E. Mallada, C. De Persis and F. Dörfler, “Robust decentralized secondary frequency control in power systems: merits and tradeoffs,” arXiv preprint, Nov. 2017. arXiv:1711.07332.
 [24] S. Trip, M. Bürger and C. De Persis, “An internal model approach to (optimal) frequency regulation in power grids with timevarying voltages,” Automatica, vol. 64, pp. 240–253, Feb. 2016. DOI: 10.1016/j.automatica.2015.11.021.
 [25] S. Auer, C. Roos, J. Heitzig, F. Hellmann and J. Kurths, “The contribution of different electric vehicle control strategies to dynamical grid stability,” arXiv preprint, Aug. 2017. arXiv:1708.03531.
 [26] Electricity Transmission System Operations. National Grid UK, 2017. [Online]. Available: http://www2.nationalgrid.com/uk/industryinformation/electricitytransmissionsystemoperations/. Accessed on: July, 2017.
 [27] The Electricity Supply Regulations 1988. The National Archives, 2017. http://www.legislation.gov.uk/uksi/1988/1057/made/data.xht?wrap=true. Accessed on: July, 2017.
 [28] R. F. Hartl, S. P. Sethi and R. G. Vickson, “A survey of the maximum principles for optimal control problems with state constraints,” SIAM Rev., vol. 37, no. 2, pp. 181–218, Jun. 1995. DOi: 10.1137/1037043.
 [29] T. Stegink, C. De Persis and A. van der Schaft, “A unifying energybased approach to stability of power grids with market dynamics,” IEEE T. Automat. Contr., vol. 62, no. 6, pp. 2612–2622, Jun. 2017. DOI: 10.1109/TAC.2016.2613901.
 [30] A. Mahajan and M. Mannan, “Decentralized stochastic control,” Ann. Oper. Res., vol. 241, no. 1–2, pp. 109–126, Jun. 2016. DOI: 10.1007/s1047901416520.
 [31] C. D. Charalambous and N. U. Ahmed, “Centralized versus decentralized optimization of distributed stochastic differential decision systems with different information structurespart I: a general theory,” T. Automat. Contr., vol. 62, no. 3, pp. 1194–1209, Mar. 2017. DOI: 10.1109/TAC.2016.2575818.
 [32] R. Singh, P. R. Kumar and L. Xie, “Decentralized control via dynamic stochastic prices: the independent system operator problem,” arXiv preprint, Jun. 2016. arXiv:1605.08926v2.
 [33] K. Schmietendorf, J. Peinke and O. Kamps, “On the stability and quality of power grids subjected to intermittent feedin,” arXiv preprint, Nov. 2016. arXiv:1611.08235.