Distributed Frequency Regulation for Heterogeneous Microgrids via Steady State Optimal Control

Distrib uted Frequenc y Regulation for Heterogeneous Microgrids via Steady State Optimal Control Lukas K ¨ olsch, Manuel Dupuis, Kirtan Bhatt, Stefan Krebs, and S ¨ oren Hohmann Institute of Contr ol Systems, Karlsruhe Institute of T echnology (KIT) , Karlsruhe, Germany lukas.koelsch@kit.edu, manuel.dupuis@student.kit.edu, kirtan.bhatt@student.kit.edu, stefan.krebs@kit.edu, soeren.hohmann@kit.edu Abstract —In this paper , we present a model-based frequency controller for microgrids with nonzero line resistances based on a port-Hamiltonian formulation of the micr ogrid model and real-time dynamic pricing. The controller is applicable for con ventional generation with synchronous machines as well as for power electronics interfaced sources and it is robust against power fluctuations from uncontrollable loads or volatile regen- erative sources. The price-based formulation allows additional requir ements such as active power sharing to be met. The capability and effectiveness of our procedure is demonstrated by means of an 18-node exemplary grid. Index T erms —fr equency regulation, steady state optimal con- trol, microgrid, port-Hamiltonian systems, distributed control I . I N T RO D U C T I O N A. State of Resear ch The energy transition motiv ates a worldwide trend towards renew able energy generation which should substitute the con- ventional po wer plants in the future. A key aspect of rene wable sources is their distributed and volatile nature compared to the centralized and well predictable character of conv entional power plants [1], [2]. This also results in a necessary change in the control schemes of the power network. So far , frequency control, i.e. the regulation of the imbalance between power generation and demand, has been the task of the transmission system operator using a hierarchy of primary , secondary and tertiary frequency control layers: In the first layer , frequency deviations and thus power imbalances are prevented from further increasing, in the second layer , the nominal state is restored, and in the third layer , an economic optimization is carried out. Both secondary and tertiary layer are each gov erned by a central controller . Howe ver , the distributed nature of renew able energy gen- eration encourages the application of a distributed frequency control scheme between multiple agents which are able to handle the control task in parallel [3]. For this reason, steady state optimal control by real-time dynamic pricing poses an ad- vantageous control concept especially for lar ge scale networks, since it enables communication of network imbalances via a price signal, see [5] for a survey on current research directions regarding frequency regulation. This kind of controller features a distributed architecture for frequency restoration based on neighbor-to-neighbor communication and local measurements This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—project number 360464149. which is able to reach a desired economic optimum at steady state and thus provides a unifying approach incorporating all three control layers [6]. A common assumption made in pre vious publications on dynamic pricing methods for frequency regulation, e.g. [7], is that all power lines are lossless, which is in fact an incorrect assumption especially for medium and low voltage grids. Practically , if controllers like the one proposed in [7] are used, a synchronous frequency ω is achie ved which deviates from the nominal frequency ω n [6]. There are publications like [2] which take nonzero line resistances into account by proposing a generalized droop control method. Howe ver , these methods rely on a fixed R/X ratio for the designed generalized droop control concept, which makes it necessary to compensate for this simplifying assumption by an additional control layer . B. Main Contributions A gradual transition from today’ s power network with large, centralized power plants to wards a future network with a lar ge number of small, decentralized generation units is essential. For this purpose, we present a model-based, distributed, steady state optimal frequency controller that is applicable for heterogeneous and lossy microgrids with both types of network connectors, i.e. synchronous generators for con ventional power plants as well as inv erters for rene wable sources. The controller design is based on [7], [8] and our previous work [6] by inte grating an in verter model based on [9] in the underlying microgrid model. Both system and controller are represented as a port-Hamiltonian system, which results in a closed-loop system that is again port-Hamiltonian. Due to the port-Hamiltonian structure, stability of the closed-loop system can finally be characterized using a (shifted) passivity property . The remainder of this paper is structured as follows. In section II, we deri ve a port-Hamiltonian model of a hetero- geneous microgrid consisting of a mixture of con ventional generation with synchronous machines, renewable generation via po wer electronics interfaced sources, and uncontrollable consumers or producers. In section III, we formulate a price- based, distrib uted frequency controller which is robust against power demand fluctuations. In section IV, we demonstrate the performance of the controller under hea vy load changes by means of an 18-node test network and in section V, we summarize our results and provide an outlook on future research directions. I I . M I C RO G R I D M O D E L Microgrids consist of con ventional and regenerativ e gener- ators as well as consumers that are all physically connected together via a lossy and meshed electrical network. Accordingly , the microgrid is modeled as a directed graph consisting of three different types of nodes: 1) Synchr onous gener ator nodes which are connected to synchronous generators of con ventional po wer plants. 2) In verter nodes which are connected to power electronics interfaced sources. 3) Load nodes which are characterized by a giv en and uncontrollable activ e and reactiv e power demand. T o set up a dynamic model of the microgrid in port- Hamiltonian formulation in the course of this section, at first all notational conv entions as well as symbols used are outlined (section II-A) and the model assumptions and sim- plifications used are listed (section II-B), before submodels for synchronous generators (section II-C), in verters (section II-D) and load nodes (section II-E) are deriv ed. Interconnection with lossy lines (section II-F) finally results in an ov erall model in port-Hamiltonian representation (section II-G), which forms the basis (“plant model”) for the controller design. A. Notational Pr eliminaries V ector a = col i { a i } = col { a 1 , a 2 , . . . } is a column vector of elements a i , i = 1 , 2 , . . . and matrix A = diag i { a i } = diag { a 1 , a 2 , . . . } is a (block-)diagonal matrix of elements a i , i = 1 , 2 , . . . . The ( n × n ) -identity matrix and ( n × n ) -zero matrix are denoted by I n and 0 n , respecti vely . Steady state (i.e. equilibrium) variables are marked with an ov erline. The microgrid is modeled by a directed graph G p = ( V , E p ) with V = V G ∪ V I ∪ V ` being the set of n G = |V G | generator nodes , n I = |V I | in verter nodes, and n ` = |V ` | load nodes. The physical interconnection of the nodes is represented by the incidence matrix D p ∈ R n × m p with n = n G + n I + n ` and m p = |E p | . Incidence matrix D p can be subdi vided as follows D p =   D p G D p I D p`   , (1) where submatrices D p G , D p I and D p` correspond to the generator , inv erter, and load nodes, respectiv ely . W e note j ∈ N i if node j is a neighbor of node i , i.e. j is adjacent to i in the undirected graph. Positiv e semidefiniteness of a matrix is denoted by  0 , whereas nonnegati vity of a scalar is denoted by ≥ 0 . Subscript p denotes the plant v ariables, i.e. the v ariables of the microgrid model, whereas subscript c denotes the variables of the contr oller . A list of symbols used for parameters and state variables of the microgrid model is giv en in T able I. T ABLE I L I ST O F M I C RO G RI D P A R A M ET E R S A N D S T AT E V A RI A B LE S A i positiv e generator, inv erter and load damping constant B ii negati ve of self-susceptance B ij negati ve of susceptance of line ( i, j ) C DC capacitance in DC circuit of the inv erter D p incidence matrix of microgrid G DC conductance in DC circuit of the inv erter G ij negati ve of conductance of line ( i, j ) i DC current source in DC circuit of inv erter i αβ output current of inv erter L i deviation of angular momentum from nominal value M i ω n M i moment of inertia p i sending-end activ e power flow p g,i activ e power generation p `,i activ e power demand q i sending-end reactiv e power flow q `,i reactiv e power demand T e electrical torque at generator T m mechanical torque at generator u DC input voltage of the inverter U i magnitude of transient internal voltage U f ,i magnitude of excitation voltage X d,i d-axis synchronous reactance X 0 d,i d-axis transient reactance θ i bus voltage phase angle ϑ ij bus voltage angle dif ference θ i − θ j Φ overall transmission losses τ U,i open-circuit transient time constant of synchronous machine ω i deviation of bus frequency from nominal value ω n ω I virtual frequency of inverter B. Modeling Assumptions In accordance with [7]–[9], we make the follo wing model- ing assumptions for the microgrid model and the controller: (A1) The grid is operating around the nominal frequency ω n = 2 π · 50 Hz . (A2) The grid is a balanced three-phased system and the lines are represented by its one-phase π -equi valent circuits. (A3) Subtransient dynamics of the synchronous generators is neglected. (A4) The matching controller of the in verters presented in section II-D has fast dynamics compared to the price- based frequency controller . Howe ver , we make the follo wing less restrictiv e assumptions: (A5) Power lines are lossy , i.e. ha ve nonzero resistances. (A6) Loads do not hav e to be constant. (A7) Excitation voltages of the generators do not have to be constant. C. Dynamic Model of Generator Nodes For generator node i ∈ V G , the third-order generator model (”flux-decay” model), described in local dq-frame, appropri- ately represents the transient dynamic behavior [10]: ˙ θ i = ω i , (2) M i ˙ ω i = − A i ω i + p g ,i − p `,i − p i , (3) τ U,i ˙ U i = U f ,i − U i +  X d,i − X 0 d,i  I d,i . (4) According to [11], the stator d-axis current I d,i can be formu- lated as I d,i = U j cos ( θ i − θ j ) − U i X d,i (5) with transient internal voltage U i and terminal voltage U j . In power system literature, it is a common assumption that the stator resistance can be neglected and thus lossless reactiv e power flow q i = U 2 i X d − U i U j X d cos ( θ i − θ j ) (6) can be used to describe the generator dynamics. W ith (6), the identity q i U i = − U i X d,i + U j X d,i cos ( θ i − θ j ) = I d,i (7) holds. Substituting (4) with (7), this yields the generator model as in [8]: ˙ θ i = ω i , i ∈ V G , (8) M i ˙ ω i = − A i ω i + p g ,i − p `,i − p i , i ∈ V G , (9) τ U,i ˙ U i = U f ,i − U i +  X d,i − X 0 d,i  U − 1 i q i , i ∈ V G . (10) W ithout loss of generality , we assume that the generated power p g ,i is controllable, while the po wer demand p `,i is uncontrollable. Thus p g ,i serves as control input and p `,i is a disturbance input. D. Dynamic Model of In verter Nodes The in verters are regulated by v DC ∼ ω matching control [9] to mimic the dynamic behavior of a synchronous generator . This is done by exploiting the structural similarities between kinetic energy of the rotor of a synchronous generator and electric energy stored in the DC-side capacitor of an inv erter . The relev ant equation describing the used 3-phase DC/A C in verter dynamics in αβ -frame is C DC ˙ u DC = − G DC u DC + i DC − 1 2 i > αβ m αβ , (11) with m αβ being an A C power electronics modulation signal generated by the controller to match the behavior of the synchronous generator presented abov e [12]. For this purpose, the internal model, also in αβ -frame, equals [12] ˙ θ = ω, (12) M ˙ ω = − Aω + T m − T e , (13) where the expression for the electrical torque T e can be expressed by T e = − L m i f i > αβ  − sin ( θ ) cos ( θ )  , (14) with the stator-to-rotor mutual inductance L m and the exci- tation current i f . Hence, the internal generator model to be matched by the in verter is ˙ θ = ω, (15) M ˙ ω = − Aω + T m + L m i f i > αβ  − sin ( θ ) cos ( θ )  . (16) For this purpose, the modulation signal is chosen to m αβ = µ ·  − sin ( θ I ) cos ( θ I )  , (17) with θ I as a “virtual” rotor angle and a constant gain µ > 0 . The “virtual” frequency thus results in ω I = ˙ θ I . Furthermore, a linking between rotational speed, or fre- quency respectiv ely , and power consumption needs to be created by construction. This is achie ved by appyling a propor- tional controller for ω I based on local measurement of u DC : ω I = η u DC , (18) again with a constant gain η > 0 . By setting µ to µ = − 2 η L m i f (19) and using (17), the last term in (11) can be reformulated to 1 2 i > αβ m αβ = − η L m i f i > αβ  − sin ( θ I ) cos ( θ I )  = η T e, I , (20) with the ”virtual” electrical torque T e, I . By inserting (20) in (11), substituting u DC w .r .t (18) and dividing by η , the in verter model can be formulated as ˙ θ I = ω I , (21) C DC η 2 ˙ ω I = − G DC η 2 ω I + i DC η − T e, I . (22) T o further highlight the resemblance between the modulated in verter equations and the synchronous generator model, as in [12], the coef ficients in (22) can be interpreted as virtual inertia M ∗ I = C DC η 2 and virtual damping A ∗ I = G DC η 2 . Moreov er, the DC current i DC is chosen according to [9] as i DC = η A ∗ I ω n + η · p g ω I . (23) Furthermore, the virtual electrical torque can be described as T e, I = p I ω I , with p I being the po wer input of the inv erter . Under the assumption that no power is dissipated in the in verter i ∈ V I , i.e. p I ,i = p `,i + p i , and use of (18) and (23), this allows a reformulation of (22) as M ∗ I ,i ˙ ω I ,i = − A ∗ I ,i ( ω I ,i − ω n ) + 1 ω I ( p g ,i − p `,i − p i ) . (24) If (A1) holds, multiplying (24) with ω n yields ˙ θ I ,i = ω I ,i , i ∈ V I , (25) M I ,i ˙ ω I ,i = − A I ,i ω I ,i + p g ,i − p `,i , i ∈ V I (26) with M I = M ∗ I ω n , A I = A ∗ I ω n , and ω I expressing the de- viation of frequency from its nominal v alue ω n . In particular , the structural similarities between (26) and the swing equation (9) can be clearly seen now . E. Dynamic Model of Load Nodes The loads are modeled by an activ e power consumption which consists of both a frequency-dependent part with load damping coefficients A i ≥ 0 and a frequency-independent part p ` , as well as frequency-independent reactiv e consumption q ` [11]: ˙ θ i = ω i , i ∈ V ` , (27) 0 = − A i ω i − p `,i − p i , i ∈ V ` , (28) 0 = − q `,i − q i , i ∈ V ` . (29) F . T ransmission Lines Generator , in verter , and load nodes are interconnected via transmission lines, which are modeled by the lossy A C po wer flow equations [10] p i = X j ∈N i B ij U i U j sin( ϑ ij ) + G ii U 2 i + X j ∈N i G ij U i U j cos( ϑ ij ) , i ∈ V , (30) q i = − X j ∈N i B ij U i U j cos( ϑ ij ) + B ii U 2 i + X j ∈N i G ij U i U j sin( ϑ ij ) , i ∈ V (31) with Y = G + j B being the admittance matrix and ϑ ij = θ i − θ j being the voltage angle deviation between two adjacent nodes. Note that by definition of the admittance matrix, G ij < 0 and B ij > 0 if nodes i and j are connected via a resistive- inductiv e line [10]. G. Overall Model The equations for generator (8)–(10), in verter (25)–(26) and load nodes (27)–(29) can be summarized in a compact notation as follows: ˙ θ i = ω i , i ∈ V , (32) ˙ L i = − A i ω i + p g ,i − p `,i − p i , i ∈ V G ∪ V I , (33) τ d,i ˙ U i = U f ,i − U i − X d,i − X 0 d,i U i · q i , i ∈ V G , (34) 0 = − A i ω i − p `,i − p i , i ∈ V ` , (35) 0 = − q `,i − q i , i ∈ V ` . (36) The interconnection of these node dynamics with the power flow equations (30)–(31) leads to the overall model of the microgrid, which is presented in port-Hamiltonian form. First, the “plant” state vector x p of the microgrid is defined as x p = col { ϑ , L G , L I , U g , ω ` , U ` } (37) with voltage angle deviations ϑ = D p θ and L G , L I being the vectors of angular momentum deviations L i = M i · ω i of generator and in verter nodes, respectiv ely . The state vector is used to set up the plant Hamiltonian H p ( x p ) = 1 2 X i ∈V G M − 1 i L 2 i + U 2 i X d,i − X 0 d,i ! + 1 2 X i ∈V I M − 1 i L 2 i − 1 2 X i ∈V B ii U 2 i − X ( i,j ) ∈E B ij U i U j cos( θ i − θ j ) + 1 2 X i ∈V ` ω 2 `,i , (38) which describes the total energy stored in the system. The first row of (38) represents the shifted kinetic energy of the rotors and the magnetic energy of the generator circuits, the second row represents the “virtual” kinetic ener gy at inv erter nodes, the third row represents the magnetic energy of transmission lines and the fourth row represents the local deviations of load nodes from nominal frequency . Using the Hamiltonian H p ( x p ) and its gradient ∇ H p ( x p ) , equations (30)–(36) can be written as follows:         ˙ ϑ ˙ L G ˙ L I ˙ U g 0 0         =                 0 D > p G D > p I 0 D > p` 0 − D p G 0 0 0 0 0 − D p I 0 0 0 0 0 0 0 0 0 0 0 − D p` 0 0 0 0 0 0 0 0 0 0 0         | {z } J p −         0 0 0 0 0 0 0 A G 0 0 0 0 0 0 A I 0 0 0 0 0 0 R g 0 0 0 0 0 0 A ` 0 0 0 0 0 0 b U `         | {z } R p         ∇ H p −         0 ϕ G ϕ I % G ϕ ` % `         | {z } r p +         0 0 0 0 0 I 0 0 0 − b I G 0 I 0 0 − b I I 0 0 ˆ τ U 0 0 0 0 0 0 − b I ` 0 0 0 − I 0               p G p I U f q ` p `       , (39) with A G = diag i { A i } , i ∈ V G , (40) A I = diag i { A i } , i ∈ V I , (41) A ` = diag i { A i } , i ∈ V ` , (42) R G = diag i  X di − X 0 di τ U,i  , i ∈ V G , (43) b U ` = diag i { U i } , i ∈ V ` , (44) ϕ G = col i  G ii U 2 i + X j ∈N i G ij U i U j cos( ϑ ij )  ,i ∈ V G , (45) ϕ I = col i  G ii U 2 i + X j ∈N i G ij U i U j cos( ϑ ij )  ,i ∈ V I , (46) ϕ ` = col i  G ii U 2 i + X j ∈N i G ij U i U j cos( ϑ ij )  ,i ∈ V ` , (47) % G = col i  R g ,i X j ∈N i G ij U i U j sin( ϑ ij )  , i ∈ V G , (48) % ` = col i  X j ∈N i G ij U i U j sin( ϑ ij )  , i ∈ V ` , (49) ˆ τ U = diag i { 1 /τ U,i } , i ∈ V G , (50) b I G =  I n g × n g 0 n g × n i 0 n g × n `  , (51) b I I =  0 n i × n g I n i × n i 0 n i × n `  , (52) b I ` =  0 n ` × n g 0 n ` × n i I n ` × n `  . (53) Note that J p = − J > p and R p  0 . Hence, this is a port- Hamiltonian descriptor system [13] with a nonlinear dissipa- tiv e relation due to r p 6 = 0 [14]. I I I . P R I C E - B A S E D C O N T RO L L E R A. Contr ol Objective In the following controller design, the control variables p g = col { p G , p I } are to be regulated in such a way that the steady-state frequency de viation ω from the nominal frequency ω n is zero at each node while, at the same time, the steady- state generation p g is optimal with respect to an objecti ve function being to be defined. In [15], it is shown that a necessary condition for ω = 0 is that the ov erall resistive losses Φ = X i ∈V G ii U 2 i + 2 · X ( i,j ) ∈E G ij U i U j cos( ϑ ij ) (54) are equal to the net sum of generation and load, i.e. Φ ! = X i ∈V G p G ,i + X i ∈V I p I ,i − X i ∈V p `,i . (55) The abov e condition serves as a fundamental constraint for any equilibrium that the closed-loop system is supposed to attain. As shown e.g. in [8], for a giv en p ` , the allocation of acti ve power injections p g is a solution of (55) if and only if there exists a ν ∈ R m c such that D c ν = b I > G p G + b I > I p I − p ` − ϕ , (56) with D c being an arbitrary incidence matrix of a commu- nication graph G c = ( V , E c ) with m c = |E c | edges and ϕ = col { ϕ G , ϕ I , ϕ ` } . This alternativ e formulation by means of (55) will result in a distributed controller with control variables p g ,i being only dependent on v ariables of node i or on variables that are adjacent with respect to the communication graph. The aim is now to design a controller in such a way that the closed-loop equilibrium, i.e. the steady state, is a solution to the optimization problem min p G , p I , ν C ( p G , p I ) sub ject to (56) (OP) with C ( p G , p I ) being an arbitrary , strictly con vex cost func- tion. Note that (OP) is a con vex optimization problem since (56) is linear-af fine with respect to optimization variables ( p G , p I , ν ) . B. Distributed Contr ol Algorithm The primal-dual gradient method for conv ex optimization problems [15]–[17] is used to deri ve a controller that solves (OP) in steady state. T o simplify the notation we define b I g = " b I G b I I # (57) and by letting [6, Proposition 1] apply , we get the distributed controller τ g ˙ p g = −∇ C ( p g ) + b I g λ + u c , (58) τ λ ˙ λ = D c ν − b I > g p g + p ` + ϕ , (59) τ ν ˙ ν = − D > c λ . (60) Diagonal matrices τ g , τ λ , τ ν > 0 are used to adjust the con- ver gence behavior of the respectiv e variable: The smaller the τ v alue, the faster the con vergence and the lar ger the transient amplitudes. u c is an additional controller input which is later chosen in such a way that a po wer-preserving interconnection of plant and controller is achiev ed [14], resulting in a closed- loop system which is again port-Hamiltonian. By defining the controller state x c = col { τ g p g , τ λ λ , τ ν ν } and the controller Hamiltonian H c ( x c ) = 1 2 x > c τ − 1 c x c (61) with τ c = diag { τ g , τ λ , τ ν } , (62) controller equations (58)–(60) ha ve the port-Hamiltonian rep- resentation ˙ x c =   0 b I g 0 − b I > g 0 D c 0 − D > c 0   | {z } J c ∇ H c −   ∇ C − ϕ 0   | {z } r c +   u c p ` 0   . (63) This representation no w provides a straightforward way to set up and analyze the closed-loop system. C. Closed-Loop System W ith the ne w composite Hamiltonian H ( x p , x c ) = H p ( x p ) + H c ( x c ) and by choosing u c = − col { ω G , ω I } as control input [6]–[8], the interconnection of plant and controller results in the closed-loop descriptor system E ˙ x = ( J − R ) ∇ H − r + F u (64) with E = diag { I 3 n G +2 n I + n + m + m c , 0 2 n ` } , (65) J =                0 b I g 0 0 − e I > G − e I > I 0 0 0 − b I > g 0 D c 0 0 0 0 0 0 0 − D > c 0 0 0 0 0 0 0 0 0 0 0 D > p G D > p I 0 D > p` 0 e I G 0 0 − D p G 0 0 0 0 0 e I I 0 0 − D p I 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 − D p` 0 0 0 0 0 0 0 0 0 0 0 0 0 0                , (66) R = diag { 0 n G + n I + n + m c , R p } , (67) r = col { r c , r p } , (68) F =               0 0 0 0 0 I 0 0 0 0 0 0 0 0 − b I G 0 0 − b I I ˆ τ U 0 0 0 0 − b I ` 0 − I 0               , (69) u = col { U f , q ` , p ` } (70) e I G =  I n G × n G 0 n G × n I  (71) e I I =  I n I × n G 0 n I × n I  . (72) Due to J = − J > and R  0 the system is again port- Hamiltonian. Denote each equilibrium x of (64) by x . This equilibrium has two salient properties which are presented in the following two propositions: Proposition 1. At each equilibrium of (64) , the fr equency deviation ω i , i ∈ V , fr om nominal fr equency ω n is zer o. Pr oof. Let 1 > =  1 · · · 1  be the all-ones ro w vector . With ˙ ϑ = 0 at steady state, the first ro w of (39) equals 0 = D > p ω . Since D p is the incidence matrix of a connected graph, this implies that each row of vector ω has the same v alue, i.e. ω = ω · 1 , and thus each node of the microgrid is synchronized to a common frequency ω . Since 1 > ϕ = Φ and ˙ λ = 0 , left-multiplying (59) with 1 > yields 0 = 0 − X i ∈V G ∪V I p g ,i + X i ∈V p `,i + Φ , (73) i.e. condition (55) is fulfilled at steady state. Moreov er, a comparison of (30) and (54) shows that P i ∈V p i = Φ . No w left-multiplying (33) with 1 > equals 0 = − X i ∈V G ∪V I A i ω i = − ω · X i ∈V G ∪V I A i (74) and since A i > 0 , it follows that ω must be zero. Proposition 2. At each equilibrium x of (64) , the marginal prices ar e equal, i.e. ∇ C ( p g ,i ) = ∇ C ( p g ,j ) for i, j ∈ V G ∪ V I . Pr oof. W ith ˙ ν = 0 at steady state, (60) equals 0 = − D > c λ . Since D c is the incidence matrix of a connected graph, this implies that each ro w of λ has the same value, i.e. λ = λ · 1 , Moreov er, with ω = 0 from Proposition 1, (58) leads to λ = ∇ C ( p g ) at steady state and hence all marginal prices are equal. Proposition 2 shows that the closed-loop system fulfills the well-known economic dispatch criterion [5] at steady state. Note that in this context, λ can be interpreted as a price signal . The stability of the closed-loop equilibrium can be in- vestigated by e xploiting the port-Hamiltonian structure (64) with its (shifted) passivity property: W ith dissipation vector R ( x ) = R ∇ H ( x ) + r , equation (64) reads as follows: E ˙ x = J ∇ H ( x ) − R ( x ) + F u , (75) with each equilibrium x fulfilling 0 = J ∇ H ( x ) − R ( x ) + F u (76) for a constant input vector u . Since H ( x ) is a con ve x and nonnegati ve function, the shifted Hamiltonian [14] H ( x ) := H ( x ) − ( x − x ) > ∇ H ( x ) − H ( x ) (77) is positive definite with minimum H ( x ) = 0 . Thus the shifted closed-loop dynamics, i.e. (75) minus (76), can be expressed in terms of H ( x ) as follo ws: E ˙ x = J ∇ H ( x ) − [ R ( x ) − R ( x )] + F [ u − u ] . (78) As a result, stability of x is giv en if the shifted passivity property [14] [ ∇ H ( x ) − ∇ H ( x )] > [ R ( x ) − R ( x )] ≥ 0 (79) is satisfied. Note that for r p = 0 , i.e. lossless microgrids, (79) is always fulfilled due to strict con vexity of C ( p G , p I ) . I V . S I M U L A T I O N A. Case Study The price-based steady state optimal controller presented in the previous section is now demonstrated by means of an 18-node exemplary microgrid with base voltage of 10 kV , n G = n I = 7 and n ` = 4 , see Fig. 1. Generator nodes are represented by black nodes, in verter nodes are represented by gray nodes, and load nodes are represented by white nodes. All parameters of the microgrid can be found in T ables II to V and are given in p.u., except τ U,i , which is gi ven in seconds. The numericals values for the parameters of generator nodes, Fig. 1. Network topology of exemplary microgrid. T ABLE II P A R A M ET E R S O F G E N ER ATO R N O D E S i 1 2 3 4 5 6 7 A i 1.6 1.22 1.38 1.42 1.4 1.3 1.3 B ii -2.67 -6.97 -4.0 -2.1 -3.5 -5.5 -7.2 M i 5.2 3.98 4.49 4.22 4.4 4.5 5.15 X d,i 0.02 0.03 0.03 0.025 0.02 0.024 0.03 X 0 d,i 0.004 0.006 0.005 0.005 0.003 0.0044 0.0068 τ U,i 6.45 7.68 7.5 6.5 6.9 7.2 6.88 load nodes and transmission lines are based on those provided in [6], [15] and the parameter values of in verter nodes base upon [9]. W ithout loss of generality , yet for sake of simplicity , we assume constant R/X ratios γ , i.e. G ij = − γ · B ij for each line ( i, j ) , and τ c = 0 . 01 · I . Moreov er , we choose D c to be identical to the plant incidence matrix D p after it has been pointed out in [6] that the specific choice of D c has little influence on the con vergence speed to the desired equilibrium. The simulations were carried out in W olfram Mathematica 12.0. B. Cost function and input signals The cost function is chosen to C ( p g ) = 1 2 X i ∈V G ∪ V I 1 w i · p 2 g ,i , (80) with weighting factors ω 1 = 1 , w 2 = 1 . 1 , w 3 = 1 . 2 and so on. Bearing in mind Proposition 2, this specific choice of C ( p g ) T ABLE III P A R A M ET E R S O F I N V ER TE R N O D E S i 8 9 10 11 12 13 14 A i 1.5 1.7 1.55 1.6 1.4 1.65 1.25 B ii -6.2 -7.1 -4.5 -4.2 -4.5 -6.05 -7.1 M i 4 3.85 6 5.55 4.1 3.9 4.32 T ABLE IV P A R A M ET E R S O F L OA D N O DE S i 15 16 17 18 A i 1.45 1.35 1.5 1.7 B ii -2.05 -2.2 -1.5 -2.1 T ABLE V P A R A M ET E R S O F T R A NS M I S SI O N L I N E S B 1 , 2 1.27 B 1 , 14 1.4 B 2 , 3 1.4 B 2 , 14 2.25 B 2 , 15 2.05 B 3 , 4 1.1 B 4 , 5 1.0 B 5 , 6 2.5 B 6 , 7 3.0 B 7 , 8 2.7 B 7 , 9 1.5 B 7 , 17 3.0 B 8 , 9 3.5 B 9 , 10 1.5 B 9 , 18 2.1 B 10 , 11 3.0 B 11 , 12 1.2 B 12 , 13 3.3 B 13 , 14 1.25 B 14 , 15 2.2 0 25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 − 0 . 2 − 0 . 1 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 t in s p `,i in p.u. i = 15 i = 16 i = 17 i = 18 Fig. 2. Stepwise increase at load nodes. as a weighted sum of squares leads to activ e po wer sharing [18] in steady state, i.e. a proportional share p g ,i /w i = p g ,j / w j = const. for all i, j ∈ V G ∪ V I . The initial values of input vector u and state vector x are chosen such that the closed-loop system starts in synchronous mode with ω ( t = 0) = 0 . At regular intervals of 100 s , a step of + 0 . 5 p.u. occurs sequentially at each load node, as shown in Fig. 2. γ is set to one. C. Results Fig. 3 shows the node frequencies for each i ∈ V . Starting from synchronous mode with a frequency of 50 Hz at each node, a de viation of the local frequency in the range of about − 0 . 45 Hz to +0 . 1 Hz occurs immediately after the load jumps, before being resynchronized again and being regulated to 50 Hz by the controller . The con vergence speed of the individual frequencies to the common frequency of 50 Hz is independent of which node the load jump occurred at. Fig. 4 shows the corresponding activ e power generation p g at generator and in verter nodes. After each load step, the controllers automatically increase p g to compensate for the additional demand. Remarkably , the individual power injections p g ,i are equidistant from each other at steady state, regardless of the total generation, thus acti ve po wer sharing is evident. The decay time of both frequency deviation and power regulation is about 40 s . It can be further accelerated by choosing smaller entries within matrix τ c of the controller . V . C O N C L U S I O N A N D F U T U R E W O R K In this paper we presented a model-based, steady state op- timal controller for heterogeneous microgrids. The underlying microgrid model for the controller can consist of a mixture of con ventional synchronous generators, power electronics 0 25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 49 . 6 49 . 7 49 . 8 49 . 9 50 50 . 1 t in s f i in Hz i = 1 i = 2 i = 3 i = 4 i = 5 i = 6 i = 7 i = 8 i = 9 i = 10 i = 11 i = 12 i = 13 i = 14 i = 15 i = 16 i = 17 i = 18 Fig. 3. Frequency regulation after step increase at load nodes. 0 25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 1 1 . 2 1 . 4 1 . 6 1 . 8 2 2 . 2 2 . 4 2 . 6 t in s p g ,i in p.u. i = 1 i = 2 i = 3 i = 4 i = 5 i = 6 i = 7 i = 8 i = 9 i = 10 i = 11 i = 12 i = 13 i = 14 Fig. 4. Active power generation at generator and inverter nodes. interfaced sources and uncontrollable loads. In contrast to state-of-the-art approaches, the controller ensures asymptotic stability of the equilibrium at nominal frequenc y of 50 Hz e ven with nonzero line resistances. The controller also provides an automatic solution to an optimization problem with a user- definable cost function. As shown in a simulation example, activ e power sharing can thus be achiev ed, for instance. Howe ver , other optimization problems can also be addressed, e.g. minimal total po wer input or minimal total line losses. The closed-loop dynamics can be formulated as a port-Hamiltonian system and thus asymptotic stability of the o verall system can be shown using a (shifted) passivity property . In future research, an integrated voltage regulation will be incorporated into the existing controller scheme. Furthermore, the presented controller will be applied to a benchmark system with significantly lar ger amount of nodes and real-world generation and load profiles to further illustrate the feasibility for large-scale systems. R E F E R E N C E S [1] G. Kariniotakis, L. Martini, C. Caerts, H. Brunner, and N. Retiere. Challenges, innov ative architectures and control strategies for future networks: the web-of-cells, fractal grids and other concepts. CIRED - Open Access Pr oceedings Journal , 2017(1):2149–2152, oct June 2017. [2] H. Bevrani and S. Shokoohi. An intelligent droop control for simulta- neous voltage and frequency regulation in islanded microgrids. IEEE T rans. Smart Grid , 4(3):1505–1513, 2013. [3] S. Patnaik, J. J. P . C. Rodrigues, and A. Gawanmeh. Cyber-Physical Systems for Next-Generation Networks . IGI Global, 2018. [4] F . D ¨ orfler , J. W . Simpson-Porco, and F . Bullo. Breaking the hierarchy: Distributed control and economic optimality in microgrids. IEEE T rans. Contr ol Netw . Syst. , 3(3):241–253, 2016. [5] F . D ¨ orfler , S. Bolognani, J. W . Simpson-Porco, and S. Grammatico. Distributed control and optimization for autonomous power grids. Pr oc. Eur opean Control Conference (ECC) , pages 2436–2453, 2019. [6] L. K ¨ olsch, K. Bhatt, S. Krebs, and S. Hohmann. Steady-state optimal frequency control for lossy power grids with distributed communica- tion. In IEEE International Conference on Electrical, Control and Instrumentation Engineering , 2019. to appear . Preprint available at https://arxiv .org/abs/1909.07859. [7] T . W . Stegink, C. De Persis, and A. J. van der Schaft. Stabilization of structure-preserving power networks with market dynamics. IF A C- P apersOnLine , 50:6737–6742, 2017. [8] T . Stegink, C. De Persis, and A. J. van der Schaft. A unifying energy- based approach to stability of power grids with market dynamics. IEEE T rans. Autom. Contr ol , 62(6):2612–2622, 2017. [9] P . Monshizadeh, C. De Persis, T . Stegink, N. Monshizadeh, and A. J. van der Schaft. Stability and frequency re gulation of in verters with capaciti ve inertia. In Pr oc. IEEE Conf. Decision Contr ol , pages 5696–5701, 2017. [10] J. Machowski, J. W . Bialek, and J. R. Bumby . P ower system dynamics: Stability and contr ol . W iley , 2012. [11] I. Boldea. Synchr onous generators . The electric generators handbook. CRC/T aylor & Francis, Boca Raton, FL, 2006. [12] T . Jouini, C. Arghir , and F . D ¨ orfler . Grid-friendly matching of syn- chronous machines by tapping into the dc storage. IF A C-P apersOnLine , 49(22):192–197, 2016. [13] C. Beattie, V . Mehrmann, H. Xu, and H. Zwart. Port-Hamiltonian descriptor systems. arXiv e-prints , 1705.09081, 2017. [14] A. J. van der Schaft. L2-Gain and P assivity T echniques in Nonlinear Contr ol . Springer International Publishing, Cham, 2017. [15] S. Trip, M. B ¨ urger , and C. De Persis. An internal model approach to (optimal) frequency regulation in power grids with time-v arying voltages. Automatica , 64:240–253, 2016. [16] A. Jokic. Price-based optimal contr ol of electrical power systems . PhD thesis, Eindhoven University of T echnology , 2007. [17] A. Jokic, M. Lazar, and P . P . J. van den Bosch. On constrained steady- state regulation: Dynamic kkt controllers. IEEE Tr ans. Autom. Contr ol , 54(9):2250–2254, 2009. [18] J. Schiffer . Stability and P ower Sharing in Micr ogrids . PhD thesis, TU Berlin, 2015.