{'fraction_non_alphanumeric': 0.04974452906584402, 'fraction_numerical': 0.01423567595359748, 'mean_word_length': 4.1694352159468435, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Energy-based Stabilization of Network Flows in Multi-machine Power Systems

This paper considers the network flow stabilization problem in power systems and adopts an output regulation viewpoint. Building upon the structure of a heterogeneous port-Hamiltonian model, we integrate network aspects and develop a systematic control design procedure. First, the passive output is selected to encode two objectives: consensus in angular velocity and constant excitation current. Second, the non-Euclidean nature of the angle variable reveals the geometry of a suitable target set, which is compact and attractive for the zero dynamics. On this set, circuit-theoretic aspects come into play, giving rise to a network potential function which relates the electrical circuit variables to the machine rotor angles. As it turns out, this energy function is convex in the edge variables, concave in the node variables and, most importantly, can be optimized via an intrinsic gradient flow, with its global minimum corresponding to angle synchronization. The third step consists of explicitly deriving the steady-state-inducing control action by further refining this sequence of control-invariant sets. Analogously to solving the so called regulator equations, we obtain an impedance-based network flow map leading to novel error coordinates and a shifted energy function. The final step amounts to decoupling the rotor current dynamics via feedback-linearziation resulting in a cascade which is used to construct an energy-based controller hierarchically.


I. INTRODUCTION
This paper aims to address two fundamental questions in power systems: how to explicitly characterize a steadystate operating point and how to achieve it via feedback. Starting from first-principles, we cast a system-theoretic treatment and arrive at the problem of stabilizing the angle configuration corresponding to an optimal network flow.
An important obstacle in stabilization of power systems dwells in the complex nonlinearity of the mechanical-toelectrical energy exchange mechanism of the rotating machines. Throughout the power systems literature, this aspect is mostly obviated while adopting reduced network and machine models of various degrees of fidelity [2], [3].
Despite the recent adoption of high order dynamical models, the analysis is still being done in a rotating frame attached to the rotor angle, a procedure difficult to extend to multi-machine networks [1], [8], [14]. In this work, we express these kind of models either in a stationary frame or in a rotating frame attached to a single independent angle variable for all machines. Thus we avoid the nonintegrability obstacle encountered in [8], while retaining the phasor interpretation, as in [13].
This work was supported in part by the ETH Zürich funds and the SNF Assistant Professor Energy Grant #160573.
Another major impediment in control design involves the choice of equilibrium with respect to which stability of the interconnected system is referred to [14]. Relating to [5], our procedure can be regarded as a reduction to a specific coupled-oscillator system, while also adopting a set-stabilization approach. Here, as in [3], the dissipation and potential functions are specified based on the network topology and in terms of the machine angles. By explicitly adressing the angle dependency in the steady-state network losses, we are able to overcome the dissipation obstacle in [3] and indeed construct a suitable energy function.
When dealing with physical systems, often complex control specifications can be formulated in terms of the level sets of an output function. For example, [4] provides a study of output synchronization of systems with relative degree one, as is the case here, while in [6] circular formations are targeted. In this spirit, we consider the consensus in angular velocity to be encoded in the passive output to be driven to zero. The resulting zero dynamics possess an invariant and attractive compact set on which the steadystate control action decomposes naturally into two angledependent components. One component balances the steadystate dissipation throughout the network, while the other acts as a gradient of an interconnection potential energy. From a power systems perspective, this function encodes the canonical network objective of inductor current minimization and capacitor voltage maximization and, as we shall see, is equivalent to that of machine angle synchronization.
The rest of the paper is organized as follows. In Section II notation is introduced, while in Section III the power system model is presented in detail. In Section IV we build up the control specifications from the system structure, while in Section V we elaborate the control strategy. Before concluding, a numerical study is performed in Section VI.

II. NOTATION
Notation will consist of the following matrices and vectors: is a matrix with all elements of the type u i on the diagonal and zeros elsewhere. In a similar fashion we use blkdiag(U i ) to lay all (square) matrices U i as blocks on the diagonal and zeros elsewhere. Denoting ⊗ as the Kronecker product, we make use of the symbol j = I ⊗ j, where I is the identity matrix in the space of dimension implied by the context. We use the symbol ∇ to denote the transpose of the Jacobian operator and to denote positive definiteness. If unspecified, we use index i to equally refer to each of the 1 up to n elements, while (·, ·) is used for the column vector concatenation. Consider the angle θ i ∈ S 1 as an object lying on the unit circle, vectorized as θ = (. . . θ i , . . .) ∈ T n , an object lying on the n-torus. Let R θi = cos θi − sin θi sin θi cos θi ∈ R 2×2 and R θ = blkdiag(R θi ) ∈ R 2n×2n denote the corresponding rotation matrices.
III. MODELLING In this section we introduce a heterogeneous multimachine power system model composed of n synchronous machines and m transmission lines which consist of inductive edges and capacitive nodes with dissipation in every element. The interconnection topology is based on a connected undirected graph defined by a signed incidence matrix E ∈ R n×m . The loads of the system are modelled as constant conductances attached to the capacitive voltage buses which are associated to each synchronous machine. To preserve clarity of presentation we omit voltage busses which are not directly attached to a generator.
Assumption 1: The three-phase system is considered to be symmetric; all resistances, inductances and capacitances have equal, positive value for each phase element.
Due to this symmetry, given a three-phase quantity z i ∈ R 3 , we shall only consider its components restricted to the two-dimensional plane orthogonal to the vector 1 ∈ R 3 , known as the αβ-frame, as follows: z i = T αβ z i , where T αβ is given by the first two rows of the power invariant Clarke

A. Synchronous machine dynamics
Consider the reference the model presented in [14], where we omit the damper windings and express the stator flux and stator current in αβ coordinates. The dynamics of the n synchronous machines are defined as Here, M = diag(M i ) ∈ R n×n denotes the rotor moments of inertia, D = diag(D i ) ∈ R n×n denotes the rotor damping coefficients, while ω ∈ R n , the rotor angular velocities. The rotor angles are represented by θ ∈ T n , the stator winding resistances by R s = diag(R si ) ⊗ I 2 ∈ R 2n×2n and the rotor winding resistances by R r = diag(R ri ) ∈ R n×n . Each generator is connected to a voltage bus v i = (v αi , v βi ) ∈ R 2 such that the stacked vector is denoted v = (. . . v i , . . .) ∈ R 2n . The machines are mechanically actuated by the rotor shaft torques u m ∈ R n and the rotor winding excitation voltages u r ∈ R n . Let λ si = (λ sα i , λ s β i ) ∈ R 2 represent the αβ component of the stator flux in machine i, such that λ s = (. . . λ si , . . .) ∈ R 2n , while λ r = (. . . λ ri , . . .) ∈ R n denote the rotor fluxes. The magnetic energy stored in all electrical machines is defined as Here, L s = diag(L si ) ⊗ I 2 ∈ R 2n×2n denotes the stator self-inductances, L m = diag(L mi ) ∈ R n×n , the mutual inductances and L r = diag(L ri ) ∈ R n×n , the rotor selfinductances. By assuming that L s L r − L 2 m 0, we have that L θ is invertible and that W e is positive. We also denote the stator current by i si = (i sα i , i s β i ), the rotor current by i ri and by τ ei , the electrical (air-gap) torque for each machine. We will find convenient the vectorized expressions for stator, rotor flux and electrical torque derived as where I r = diag(i ri ). Finally, we will also use the vectorized expressions for stator and rotor electromotive forces representing the voltages mutually induced in the stator and rotor circuits, respectively, as per Faraday's law.

B. Transmission network dynamics
The transmission system interconnecting all generators is modelled as the following linear circuit where we denote by C = diag(C i ) ⊗ I 2 ∈ R 2n×2n the bus voltage capacitances and by G = diag(G i ) ⊗ I 2 the load conductances at each node with respect to the potential of the earth. Furthermore, i t ∈ R 2m denotes the vector of transmission line currents in αβ coordinates, while L t = diag(L ti ) ⊗ I 2 ∈ R 2m×2m denotes the transmission line inductances and R t = diag(R ti ) ⊗ I 2 ∈ R 2m×2m the transmission line resistances. Here, the incidence matrix E = E ⊗ I 2 ∈ R 2n×2m describes the interconnection topology. Remark 1: (The γ-component). Notice that the α and β components of the transmission system are independent and follow identical dynamics. The main assumption for reducing to the αβ-components of the three-phase circuit is that the third component, also called γ-subsystem, is not excited and follows linear, asymptotically stable dynamics.

C. Combined system
To write system (1), (4) in port-Hamiltonian form, first define the state vector as x = (M ω, θ, λ r , λ s , Cv, L t i t ) ∈ X and the input vector as u = (u m , u r ) ∈ U, where the state space and input space are defined respectively as Further consider the following Hamiltonian function By appropriately using constant matrices J = −J , R positive semi-definite and B, equations (1), (4) can be written in port-Hamiltonian forṁ where ∇H = (ω, τ e , i r , i s , v, i t ) represents the co-energy variables. In this setting, we identify the passive output y = (ω, i r ) to serve as a starting point for control design.
IV. PROBLEM FORMULATION We begin with the observation that, for systems with well defined vector relative degree such as (7), there exists a unique feedback which renders the zero-dynamics manifold invariant, see e.g. [9]. Consider ω 0 > 0 an angular velocity to be tracked by all generators and i * r ∈ R n >0 a vector of reference rotor currents. The associated zero-dynamics manifold can be written as . The rest of this chapter develops a series of refinements of the target set Ω, together with the control action which renders it invariant, arriving in the end at a suitable control specification for the power system at large.

A. Zero dynamics analysis
We will find it convenient to denote the stator impedance, bus admittance and line impedance matrices, respectively, as as well as the transmission line Laplacian, L t = EZ −1 t E . The following result, yet intermediary, will aid us in the subsequent derivations.
Lemma 4.1: (Feedback on Ω). The unique control action u = u (x) which renders Ω invariant for (7) is given by Proof of Lemma 4.1: Using the tangency condition for a zero-level set of a function, we have that Ω is invariant if and only if the following condition holds ∀x ∈ Ω ω = 0 and d dt i r = 0 .
From (2b), (3b) we have thatλ r = d dt i r + ξ r which, together with (1a), (1c) allows us to rewrite (9) as By evaluating the expressions (3) using (9) we require that, ∀x ∈ Ω, Taking a time derivative of (2a) yieldsλ s = d dt i s + ξ s . By using this, together with (1d), in (11a), we express (11b) as which, together with (10b) yields (8a). Finally, using ω = ω 0 1 and i r = i * r in (2c) along with (10a), completes the proof. Notice that, with this feedback, the stator electro-motive force ξ s in (11a) simply becomes an Euclidean embedding of the angle θ, whenever x ∈ Ω. The zero-dynamics vector field now reads asθ with ξ s given by (11a). The next statement leverages the fact that, when ξ s is used to embed θ in Euclidean space, the zero-dynamics (13) can be written as a block-triangular system. Proposition 4.2: (Stability of zero-dynamics). There exists a unique linear mapping Π such that the set Γ defined below is globally asymptotically stable for the closed loop system (7),(8) Proof of Prop 4.2: Notice that the closed loop dynamics (7), (8), when restricted to Ω, become equivalent to (9), (13).
Using (11a) we write (13a) asξ s = ω 0 jξ s . By defining z = (i s , v, i t ) as the state of the three-phase electrical subsystem on Ω, we can equivalently write system (13) aṡ As we will see shortly, −Q −1 A is Hurwitz, and its spectrum does not intersect the imaginary axis. Hence, there exists unique solution Π ∈ R (4n+2m)×2n for the associated Sylvester equation To find Π we use the commutativity property of the Kronecker product, involving j = −j −1 . Observe that Π = (A + ω 0 jQ) −1 P solves (15) and is given by To show asymptotic stability, we define the transient component of the driven subsystem asz = z − Πξ s and, using equation (15), write its dynamics as: Choosing as Lyapunov function which is negative definite at zero, hence the transient subsystem is asymptotically stable and −Q −1 A is Hurwitz as originally assumed. Moreoverz(t) → 0 as t → ∞ here implies global asymptotic stability of Γ for (9),(13).
Remark 2: (Geometry of Γ). We have seen that the set Γ is positively invariant for the zero dynamics; notice that it is also compact. To illustrate its geometry, we express it as

B. The steady-state behaviour
Consider feedback (8), suppose that system (7) is now initialized on Γ and call such solutionsx(t) ∈ Γ, ∀t > 0. Due to control invariance of Ω one has that ω = ω 0 1 and i r = i * r for all time, but since we are also on the invariant set Γ ⊂ Ω, we have that the solutionsx(t) satisfẏ where (î s ,v,î t ) =π(θ) represent the three-phase network state variables restricted to Γ, seen as phasors, and are the electrical torques and stator electromotive forces of system (7), when restricted to Γ. We denote I * r = diag(i * ri ). Consequently, θ solely represents the state of the reduced dynamics on Γ. Notice that these dynamics can be parametrized by a constant θ dq ∈ T n , seen as the initial condition of system (18), since all solutions that satisfy (18) are of the form Consider the following set, which we call a diagonal fiber of T n through θ dq , . In this way, every operating point θ dq ∈ T n can be associated with a flow on Γ, defined as the subset We see that such a set Γ| θ dq ⊂ X corresponds to a closed curve in terms of the evolution ofx, thereby characterizing a steady-state operation of our system. By using the Π map in equation (12), we have that ξ r = 0 for the dynamics restricted to Γ. Let us now characterize how the control input (8) appears on the set Γ| θ dq . Proposition 4.3: (Steady-state input). Consider system (7) and a set-point θ dq ∈ T n . The unique control action u = u * (θ dq ) which renders Γ| θ dq invariant is given by Proof of Proposition 4.3: First, to satisfy (18c) while on Γ, we pick u * r = R r i * r . By lettingẑ = (î s ,v,î t ), denote the state of the threephase electrical subsystem on Γ, notice that the last three equations of (18) are equivalent tô Using (19b) and (22a) into (19a) allow us to express the steady-state electrical torque solely in terms of θ and i * r aŝ Finally, if we consider the commutativity property of the Kronecker product in (23), we have thatτ e = K net (θ dq )ω 0 1, which is constant on Γ| θ dq . By using it in (18b), we conclude the proof. Alternatively, this result is also obtained by using the Π map in (8) We then have that ⇧ = (A + ! 0 jQ) 1 P solves the Sylvester equation and is given by To show asymptotic stability, we define the transient component of the driven subsystem asz = z ⇧⇠ s and, using equation (14), write its dynamics as: Choosing as Lyapunov function , which is negative definite at zero, hence the transient subsystem is asymptotically stable and Q 1 A is Hurwitz as originally assumed. Moreoverz(t) ! 0 as t ! 1 here implies global asymptotic stability of for (12).

B. The steady-state locus
Remark 3: (Geometry of ). We have seen that the set is positively invariant for the zero dynamics; notice that it is also compact. To illustrate its geometry, we express it as , where the map⇡ : T n 7 ! R (4n+2m) is given bŷ ⇤ Consider feedback (8), suppose that system (7) is now initialized on and call such solutionsx(t) 2 , 8t > 0. Due to control invariance of ⌦ one has that ! = ! 0 1 and i r = i ⇤ r for all time, but since we are on the steady-state locus ⇢ ⌦ we also have that✓ for all t > 0, wherê and where ✓ represents the state of the reduced dynamics on . Define ✓ dq as the initial condition of system (19) such that all solutions that satisfy (19a) are of the form Definition 1: (Diagonal fiber). The set F| ✓dq defined below is called a diagonal fiber of T n through ✓ dq Notice that every initial condition compact subset of defined as Before stating the next result obse structure of the ⇧ map, ⇠ r = 0 for the Proof of Proposition 4.3: First, to sati we pick u ⇤ r = R r i ⇤ r . Notice that the la (19) are equivalent to the definition of According to (22b), we can interpre bined losses and loads in the three-phase state, projected as a damping factor in balance relation, and illustrated in Fig.  According to (21b), K net can be interpreted as the equivalent admittance Y net , projected via the angle θ dq as a damping factor onto the torque balance equation, as illustrated in Fig. 1. We shall shortly see that this projected damping term can be decomposed into an active part and a gradient term.

C. Definition of control specifications
The aim of this section is to show how circuit aspects play a role in selecting the set-point θ dq , allowing a further refinement of the control objective encoded by the set Γ| θ dq . Lemma 4.4: (Network potential). Consider feedback u = u * (θ dq ) from (21) and suppose that system (7) is initialized on Γ| θ dq for some θ dq ∈ T n . Consider the following energy function Then, θ dq is a critical point of S if and only if θ dq is chosen such that (21b) becomes Proof of Lemma 4.4: One can also express Y net as At this point, writting (21b) as , and splitting Y net into its symmetric and antisymmetric parts allows us to respectively decompose u * m into an active part and a gradient term. The crucial observation is that one can rewrite S via the network flow map as whereby taking the gradient with respect to θ dq yields Finally, all solutions of ∇S = 0 are critical points of (24). Whenever θ dq is chosen in such a way, one has that the asymmetric part of Y net vanishes and (21b) becomes (25).
As we will see in the next chapter, by designing controllers to actively minimize the energy function S, one can achieve, in terms of (î s ,v,î t ) variables, a saddle-point flow. In such a case, θ dq converges to a point where the energy stored in the stator impedance and in the transmission lines (edge currents) is minimized, while the energy delivered to the load (nodal voltage) is maximized. Fundamentally this is captured in the minus sign of the second term in (24) and can be seen as the min-flow max-cut duality in circuits.
Furthermore, when θ dq is not a critical point of S, the antisymmetric part of Y net would appear in (21b) to account for torque injections exchanged in-between generators and not delivered to the local loads. This is seen by left-multiplying (21b) by ω 0 1 , yielding the power balance on Γ ω 0 1 u * m = ω 0 1 Dω 0 1 +î s R sîs +v Gv +î t R tît + ω 0 1 ∇S .
Since 1 ∇S = 0 by construction, any choice of θ dq would satisfy the power balance. However, only when θ dq ∈ ∇S −1 (0), the gradient part of u * m would be zero (element by element) in (21b). In this case, the torque injections would exclusively account for local dissipation (local loads and the fair-share of line losses). When θ dq is chosen as a global minimizer of S, then, in addition, no unnecessary circulating power would be produced in the network, either because there are no loops in the transmission line graph or, conversely, because such a low energy level would not allow circulating currents through loops in the graph. Fig. 2. A plot of S(θ) for a two-generator test case, with θ 1 , θ 2 ∈ [−π, π). Notice that the minimum occurs on the zero-fiber, F| 0 ⊂ T 2 .
As seen in Figure 2, for a given θ dq , S(θ) = S(θ dq ) for all θ ∈ F| θ dq . In the aim of illustrating a global stability property, we define our control specification via the set which is closed and bounded and, depending on the graph topology, may have more than two disconnected components.

V. CONTROL DESIGN
Now we are ready to state our main control problem. Problem 1: Multi-machine synchronization problem (MMSP). Consider system (7) and a pair ω 0 > 0, i * r ∈ R n >0 . Design feedback u(x) such that the set O is globally asymptotically stable for the closed-loop system. Remark 3: (Voltage control). A general guideline for choosing the rotor current set-point i * r is to prescribe a common stator-electromotive-force amplitude reference for all generators, denoted e 0 , from which we derive i * r = (ω 0 L m ) −1 e 0 1. This baseline rotor current reference can be further adjusted to account for the voltage drop across the stator impedance and regulate instead the bus voltage amplitudes to a value closer to the set-point e 0 .

A. The hierarchical control strategy
Drawing from the distinct nature of the two inputs u r and u m , we can address the system nonlinearity hierarchically, by decoupling the rotor current and stabilizing it to a constant. This approach allows the integration in our framework of permanent-magnet synchronous machines (pmSM), which can be seen as having constant i r , and more importantly, power electronics converters which, as shown in [10], can be rendered equivalent to pmSM.
Definition 1: Multi-machine synchronizing feedback (MMSF) class. Consider system (7) and a pair ω 0 > 0, i * r ∈ R n >0 . Further consider inputs u * m (θ) from (25) and u * r from (21a). A feedback u = (u m , u r ) of the form is of multi-machine synchronizing feedback (MMSF) class if the following conditions are satisfied: (i) Rotor current regulation: The set R = {x ∈ X : i r = i * r } is globally asymptotically stable for the closed loop; (ii) Velocity and angle synchronisation: The set is globally asymptotically stable for the closed-loop dynamics restricted to the set R.
From the construction of feedback (29), we have that the stabilizing terms S r (x) and S m (x) vanish on the set O, yielding the zero-dynamics (18) with associated set of initial conditions given by the critical points of S.
Theorem 5.1: (Solvability of MMSP). Any feedback of class MMSF solves the MMSP for system (7), provided all solutions of the closed-loop system are bounded.
Proof of Thm 5.1: To prove this, we set up the conditions to apply Proposition 14 in [7], which is restated for convenience in Appendix A.
Step 1: (Hierarchy of specifications) The sets O ⊂ R ⊂ X are closed and positively invariant for the closed-loop, which is seen from Definition 1 and Lemma 4.4. Furthermore O is compact since it is also contained inside the compact Γ.
Step 2: (Hierarchical feedback) Using Definition 4 from [7] we can see that O is globally asymptotically stable relative to R, while the set R is in turn globally asymptotically stable for the closed-loop system.
Together with the hypothesis of bounded closed-loop solutions, we are able to use Proposition 14, case (b), in [7] and evaluate that the set O, in this case, is globally asymptotically stable for the closed-loop system.

B. An Energy-based solution to MMSP
Following the strategy in Theorem 5.1, we construct a solution which, on one hand enforces a cascade between the rotor current dynamics and the rest of the system via feedback-linearization and on the other, uses the input designed in (25) to account solely for the angle-dependent steady-state losses and induce a negative-gradient flow.
Proposition 5.2: (Solution to MMSP). Consider system (7) and assume that the following holds ∀θ ∈ T n D 1 4 ω 2 0 I * r (L m ⊗ e 2 )R θ Π Q 2 K −1 ΠR θ (L m ⊗ e 2 )I * r (30) where K = blkdiag(R s , G, R t ) and Q = blkdiag(L s , C, L t ). Then, the following controller is of MMSF class and solves the MMSP Proof of Proposition 5.2: Step 1: (Rotor current regulation) Using (2b) in (1c), we can express the rotor current dynamics as: By using (3b), we see that implementing (31a) assigns linear, asymptotically stable dynamics for the rotor currents which implies that the set R is rendered globally asymptotically stable in closed-loop.
Step 2: (Velocity and angle synchronization) First notice, via the decomposition of Y net , that (31b) is equivalent to (25). Consider now the following coordinate transformatioṅ where θ 0 ∈ S 1 is an auxiliary variable and where (î s ,v,î t ) ∈ R 4n+2m is defined as (î s ,v,î t ) =π(θ). Secondly, assume system (7) with input (31) is initialized on the invariant set R and consider its dynamics, written in coordinates (34) Consider now the following energy functioñ Cṽ + 1 2ĩ t L tĩt + S(θ dq ) , whose derivative along the trajectories of the closed loop system reads aṡ where Notice that (35) is negative definite at zero in (ω,ĩ s ,ṽ,ĩ t ) if and only if the following holds ∀θ dq ∈ T n D − T 12 (θ dq ) which is equivalent to (30). By LaSalle's invariance principle, we conclude that closed-loop dynamics converge to which directly implies global asymptotic stability of O for the closed-loop system (7), (31) reduced to R.
Step 3: (Boundedness) So far, we have constructed a MMSF class controller with S r (x) = ξ r and S m (x) = 0. To see that solutions of system (7) under feedback (31) are bounded, first consider the rotor current dynamics (33) and notice that i r (t) is bounded.
Next, consider the auxiliary variable ψ = R θ (L m ⊗ e 1 )i r which is a multiplication of two bounded terms, hence it is bounded for all time. Furthermore, consider the stator flux dynamics (1d) where we express, using (2a), the stator current via the auxiliary variable as i s = L −1 s (λ s − ψ). With this in mind we turn to the electrical subsystem, composed of (1d) and (4) and rewrite it aṡ which represents a linear system, asymptotically stable in open-loop, which is driven by the bounded input ψ. Hence the solution of this subsystem, is also bounded for all time.
We now have, using (2b), that λ r (t) is also bounded. What remains is to show that the ω component of the solution is bounded. This is seen by looking at the corresponding subsystem in closed-loop which is again a linear system driven by a bounded input, from which we deduce that ω(t) is bounded. Having evaluated that the entire solution x(t) is bounded for all time, we can conclude that the MMSP is solved.
Observe that the dissipation condition (30) resembles the one of Corollary 11 in [12], as well as generalizes inequality (41) in [1] to the multi-machine scenario developed here.
Remark 4: (A high-gain stabilizer) Notice that, in controller (31b), the stabilizing term S m is not active. One can construct another solution, based on feedback-linearization, which fully accounts for the electrical torque extraction (when restricted to R) and restores the negative-gradient component. In this case, the stabilizing term becomes which, compared to that of (31b), provides a high-gain error injection and does not require (30). Overall, both controllers require the implementation of either the active or the gradient part of the projected dissipation term K net (θ). Remark 5: (Projected gradient) The resulting synchronising component, namely ∇S(θ) in (31b), has the structure of a Laplacian vector field projected onto the tangent space of T n . Since it acts indirectly, via the torque balance equation (1a), the closed loop can be seen as a second-order gradient descent minimizing S over the angle variable θ.

VI. NUMERICAL EXPERIMENTS
To illustrate some practical aspects, we use a 3-machine model defined through the following set of parameters, all in corresponding S.I. units: L m,2 = 0.08 L m,3 = 0.02 L s,1 = 0.0018 L s,2 = 0.001 L s,3 = 0.0066 R s,1 = 0.166 R s,2 = 0.07 R s,3 = 0.5 The parameters of the first synchronous generator were inspired from Example 3.1 in [11], while the other two machines are slight variations of it. Furthermore, the topology encoded by E =  Figure 3, that the output regulation specifications are met, along with boundedness of solutions. In addition, the angles synchronize to zero relative differences: this is due to the fact that, for a typical choice of parameters, the syncrhonizing potential S has the global minimum corresponding to identical angles.
Notice that the angular velocity overshoot of generator 2 is induced by the long transient of its rotor current. This has been created with the purpose of showing the effects of hierarchical control. From an implementation perspective it is worth noting that the small oscillations of the control action, during mild transients, would have to be tracked by a sufficiently fast torque actuation.

VII. CONCLUSIONS
In this article, starting from a first-principle power system model, we develop a series of invariant sets in the aim of characterizing its steady-state behavior. On these sets, the steady-state control action is subsequently derived, resulting in a gradient vector field acting on the machine torque balance equation. Through the use of a steady-state network map, parametrized in terms of the angle variable, we decompose the equivalent admittance of the three-phase circuit to uncover a canonical potential function. This function incorporates the effects of bus admittance and stator and line impedances, thus representing an extension to that of the classical swing-equations. By accounting specifically for the steady-state losses, we arrive at a coordinate transformation and an associated shifted Hamiltonian which allow us to design an energy-based controller. As feedback, the controller requires machine angle measurements as well as um,1(t) D1!0 ✓dq,3(t) ✓dq,2(t) ✓dq,1(t) ir,1(t) ir,3(t) ir,2(t) ur,1(t) ur,3(t) ur,2(t) full knowledge of network parameters. For this synthesis, we propose a framework which hierarchically stabilizes the excitation current first, and then simultaneously achieves frequency consensus and an angle configuration corresponding to the optimizers of the a network potential function. The discussion on robustness of implementation and how to achieve a non-identical angle configuration, corresponding to generators exchanging power over the transmission lines, is deferred to a future publication.