Distributed Control and Optimization of DC Microgrids: A Port-Hamiltonian Approach

This article proposes a distributed secondary control scheme that drives a dc microgrid to an equilibrium point where the generators share optimal currents, and their voltages have a weighted average of nominal value. The scheme does not rely on the electric system topology nor its specifications; it guarantees plug-and-play design and functionality of the generators. First, the incremental model of the microgrid system with constant impedance, current, and power devices is shown to admit a port-Hamiltonian (pH) representation, and its passive output is determined. The economic dispatch problem is then solved by the Lagrange multipliers method; the Karush-Kuhn-Tucker conditions and weighted average formation of voltages are then formulated as the control objectives. We propose a control scheme that is based on the Control by Interconnection design philosophy, where the consensus-based controller is viewed as a virtual pH system to be interconnected with the physical one. We prove the regional asymptotic stability of the closed-loop system using Lyapunov and LaSalle theorems. Equilibrium analysis is also conducted based on the concepts of graph theory and economic dispatch. Finally, the effectiveness of the presented scheme for different case studies is validated with a test microgrid system, simulated in both MATLAB/Simulink and OPAL-RT environments.


I. INTRODUCTION
E LECTRIC power systems are shifting towards the use of more green technologies. To effectively integrate the renewable energy resources, energy storage systems, and electric loads into the power systems, they are interfaced with the grid via power electronic converters and are grouped in the form of microgrids (MGs) easing their control and management [1]. As a key component of modern power systems, dc microgrids have recently become more attractive [2]. They are compatible with the dc electric nature of renewable energy resources, energy storage systems, and a majority of electric loads. In addition, compared to ac MGs where control of frequency, phase, reactive power, and power quality are big challenges, control and management of dc grids are inherently simpler [2].
In dc MGs, distributed generators (DGs) and loads are connected to the grid via power converters which are either voltage-controlled (grid-forming) or current/power-controlled (grid-following). Grid-forming devices adjust the voltage of their point of common coupling (PCC) to follow a given voltage reference. The grid-following devices, on the other hand, follow some current/power references [3]. Therefore, The authors are with the Department of Electric Power Engineering, Norwegian University of Science and Technology, 7491 Trondheim, Norway (e-mail: babak.abdolmaleki@ntnu.no; gilbert.bergna@ntnu.no).
in terms of current/power, the grid-forming DGs and loads are dispatchable while the grid-following ones are nondispatchable and can be considered as constant current/power loads. In autonomous MGs, normally a cluster of dispatchable (grid-forming) generators are in charge of shaping the desired voltage level; thus, they should control the dc MG in a collaborative effort.
A common practice to dispatch the current and to adjust the voltage of grid-forming DGs in a decentralized, communication-free fashion is droop control. Despite its simple and robust functionality, this primary controller cannot guarantee desired current-sharing and voltage formation between the DGs [4]. To address these shortcomings, the droop characteristic can be corrected by a secondary controller exploiting data exchange either between the DGs and a central control unit or only between the DGs. In the former case, the controller is known as centralized secondary control which exhibits a single-point-of-failure to the system and requires a complex communication network between the DGs and the central control unit. Therefore, the distributed control techniques, using neighbor-to-neighbor inter-DG data transmissions, are preferred to the centralized ones [5].

A. Existing Literature and Research Gap
Distributed control of dc MGs has already been addressed in many works. A consensus-based proportional current-sharing strategy is proposed in [6] where a dynamic consensus-based estimator is additionally employed to keep the average voltage at nominal value. To reduce the communication burden and to reach faster convergence under this controller, it has respectively been modified to event-triggered and finite-time versions in [7] and [8]. In [9], a distributed optimal control scheme is proposed under which the DGs can achieve economic current-sharing. Therein, to overcome the initialization and noise robustness problems related to dynamic consensus-based estimation, a modified dynamic consensus-based average voltage observer is used which determines the voltage reference of converters so that the DG currents are shared properly. A somewhat similar control strategy to [6], but with eventtriggered communications, is proposed in [10] under which only information of the DGs' currents are communicated among them. In [11], a distributed nonlinear controller is proposed for dc MGs which, instead of droop controller, tunes the DGs voltages so their currents are shared proportionally. To bound the DG voltages within a reasonable range and to guarantee current-sharing among them to a certain degree, in [12], a containment-consensus-based controller is proposed for dc MGs. A very similar containment-based controller, but with finite-time convergence, is also proposed in [13]. In the above mentioned works, either the electric network dynamics and electrical system are not taken into account or only a simplified linear algebraic representation of the grid is considered. Consequently, the controller design and system stability may depend on the parameters of the physical system which are subject to modelling uncertainties.
One way to achieve plug-and-play (PnP) design and operation is to consider the overall system dynamics and to control the system based on energy principles. To do so, in [14], a distributed passivity-based control is proposed for buckconverter-based MGs ensuring proportional current-sharing and average voltage regulation among the DGs. Some similar versions of this controller are presented in [15]- [17] which demonstrate superior transient system performance. It should be noted that the asymptotic stochastic stability of the controller proposed in [17] has further been studied in [18] under varying loads. To reach the desired control objectives in the mentioned works in a finite-time manner, some sliding mode controllers have been developed in [19], [20]. Moreover, a few consensus-based proportional current-sharing and voltagebalancing controllers, facilitating PnP operations, have been proposed in [21], [22] for dc MGs with constant power and exponential loads where the existence and stability of the system equilibria is also studied. Following the same concept and for the sake of PnP functionality of the DGs, in [23], a distributed dynamic control strategy is proposed for voltage balancing and proportional current-sharing among parallel buck converters with the same capacity. The aforementioned works are, however, limited to buck converter-based DGs and proportional current-sharing and none of them has considered droop-controlled DGs and their economic current dispatch.

B. Contributions
Motivated by the above literature review, a distributed secondary control strategy for dc MGs with ZIP loads is proposed herein with the following noticeable features. First, in the modeling of the power system, the dynamics of transmission lines and shunt capacitors are considered, loads-and also current/power-controlled DGs-are modeled by constantimpedance-current-power (ZIP) loads, and the generators are characterized by droop-based grid-forming DGs, encompassing various types of interfacing converters. It is shown that the incremental model of the droop-based MG admits a port-Hamiltonian (pH) representation [24], and its passive output is defined. Second, drawing inspiration from the Control by Interconnection (CbI) technique of pH systems [25], a distributed consensus-based secondary controller is proposed which drives the MG to an equilibrium point where i) the DGs share optimal currents, and ii) their weighted-average voltage is the nominal voltage. The voltage weightings are directly related to coefficients of the DGs cost functions and not the electric network and loads. Third, regional asymptotic stability of the system with ZIP loads is demonstrated and it is shown that the system is globally asymptotically stable without the presence of constant power loads (CPLs). Finally, Fig. 1. A droop-based DG connected to a microgrid with ZIP load. equilibrium analysis is conducted based on the concepts of economic dispatch and graph theory. The rest of this paper is structured as follows. The MG system modeling and the control aims are formulated in Section II. Section III presents the proposed controller and the system stability and equilibrium analyses. The case studies and simulation results are given in Section IV. Finally, Section V concludes the paper and discusses future research directions.

Converter Dynamics & Internal Controllers
Throughout the paper, R n×m and R n stand for the set of n× m real matrices and n × 1 real vectors, respectively. diag{x i } indicates a diagonal matrix with x i being the corresponding diagonal arrays. col{x i } shows a column vector with the arrays x i . I is an identity matrix with appropriate dimensions. 0 and 1 are appropriate all-one and all-zero vectors or matrices. The transpose of a matrix/vector z is given by z . Given the scalar x or the vector x,x andx are their value at the equilibrium point, andx = x −x andx = x −x.

II. MICROGRID MODELING AND CONTROL OBJECTIVES A. Electric Network, Generators, and ZIP Load Models
Let N e , E e , and G e , with the cardinalities n N e , n E e , and n G e , be the sets of buses, transmission lines, and grid-forming (voltage-controlled) generators, respectively. Suppose that the transmission lines are modeled by serial resistor-inductor pairs, the buses are modeled by shunt capacitors and ZIP loads, and each generator is modeled by a controllable voltage source which is connected to the grid via a transmission line (See Fig. 1).
The described electric network can be modeled as two graphs M e and M G e where the buses and transmission lines play the roles of their nodes and edges, respectively. Consider e are its node set, edge set, and incidence matrix, respectively. Similarly, the graph M G e = (N e , G e , B G e ) can be defined with the same node set but different edge set G e = {1, · · · , n G e } and incidence e . An incidence matrix describes the network graph topology by determining the connections between the bus voltages and line currents. For the electric network, one should first consider an arbitrary current-flow direction for every line (edge); if current of jth line enters node k then b kj = 1, if it leaves node k then b kj = −1, otherwise b kj = 0. Similarly, if ith DG injects current to bus k via an output connector, then b Ge ki = 1; otherwise, b Ge ki = 0. Note that in this work, the generators are assumed to only inject current to the loads and network and not to absorb it, i.e., b Ge ki = −1 is not considered.
According to Fig. 1 and based on the system incidence matrices, the dynamics of the droop-based microgrid system are as follows.
where L Ge i , R Ge i , and I Ge i , ∀i ∈ G e are inductance, resistance, and current of ith generator transmission line; L Ee j , R Ee j , and I Ee j , ∀j ∈ E e are jth line inductance, resistance, and current; C Ne k , I L k , and V Ne k , ∀k ∈ N e are capacitance, load current, and voltage at bus k; G cte k ≥ 0, I cte k , and P cte k are respectively constant conductance, current, and power values of the ZIP load at bus k; V nom , V i , and V ref i are nominal voltage and ith DG voltage and its reference value, respectively; R D i and u i are respectively the droop coefficient and correction term (input) of ith generator.
There are various types of internal current and/or voltage controllers for converters which are normally designed to be very fast. Hence, in secondary control and optimization design and studies, the following assumptions are usually required.
Assumption 1: The grid-forming (voltage-controlled) generators can be modeled as controllable voltage sources so that Therefore, considering the well-known droop equation the grid-forming units are characterized by the algebraic The grid-following (current-/power-controlled) converters are considered as negative constant current/power loads in the ZIP load model. Let us define the following global matrices and vectors.
one can write the system in the following form.
where B G e and B e are the incidence matrices defined in the preamble of this subsection; are the skew-symmetric and symmetric component of F, respectively.
Assumption 3: The system Σ (2) has a unique equilibrium pointx. Moreover,ū (resp.ȳ) are the equilibrium control (resp. equilibrium output) of (2) at the equilibrium point where The incremental model of the system Σ forx = x −x and u = u −ū can then be written as the PH system below.
. Proposition 1: With the storage function H(x) = 0.5x Qx, and the passive outputỹ with respect toũ, the systemΣ (3) is passive in the following domain.
Proof: Since J = −J , the derivative of the storage function along the trajectories of (3) iṡ On the other hand, the matrixR(x) is positive definite for all x ∈ D. Therefore, the systemΣ (3) is passive with the given storage function [26].

B. Economic Dispatch and Near-Nominal Voltage Formation
Let C i (I Ge i ) = α i (I Ge i ) 2 + β i (I Ge i ) + γ i be ith generator cost function, where α i , β i , and γ i are its parameters. If I demand is the total current demand in the power network, then the economic current dispatch problem can be written as the following optimization problem.
This optimization problem can be solved by Lagrangian method with the following Lagrangian function [27].
where λ is dual variable or Lagrange multiplier. The primal problem is convex; hence, if Slater's condition is satisfied, then the Karush-Kuhn-Tucker (KKT) conditions provide necessary and sufficient conditions for primal-dual optimality of the points as follows [27].
Primal feasibility: ∂L/∂λ = 0, This implies that considering a feasible equality constraint in the problem, the KKT optimality conditions are boiled down to the stationary condition [27] lim t→∞ λ i = λ j = λ opt (6) where λ i = ∂C i /∂I Ge i = 2α i I Ge i + β i is the incremental cost (Lagrange multiplier) of ith DG, and λ opt is its optimal value. This condition is known as equal incremental costs (EIC) criteria [28]. Due to the fact that current sharing in power networks depends on the bus-voltage differences and not the absolute values of voltages, theoretically speaking, the above mentioned optimality condition can be satisfied in many voltage levels; i.e., λ opt can have various values depending on the weighted average of voltages. However, in practice the voltages must be as close as possible to the network's nominal voltage. Therefore, the controller should also guarantee a nearnominal voltage formation which can be formulated as where w i > 0, ∀i ∈ G e are voltage weightings which are defined later. Remark 1: A special choice of the cost function parameters is α i = 0.5/I rated i , β i = γ i = 0 which turns (6) into the equal current ratios criteria (I Ge i /I rated ) underlining the proportional current-sharing, studied in the literature (See e.g., [11]- [23]).

III. CONTROLLER DESIGN, CLOSED-LOOP SYSTEM EQUILIBRIUM, AND STABILITY ANALYSIS
In this section, a distributed controller is proposed for the droop-based MG system to satisfy the control objectives in (6) and (7). The proposed controller relies on both the local and neighborhood measurements of the generators; hence, the generators need to exchange information through a communication network as described next.

A. Communication Network Model
A communication network between the generators can be modeled as a undirected graph with generators and communication links being its nodes and edges, respectively. Consider are its node set, edge set, and adjacency matrix, respectively. If nodes i and j exchange data, then they are neighbors, (j, i) ∈ E c , and a ij = a ji > 0; otherwise, nodes i and j are not neighbors, (j, i) / ∈ E c , and a ij = a ji = 0. Let N i = {j|(j, i) ∈ E c } and d i = j∈Ni a ij be neighbor set and in-degree of node i, respectively. The Laplacian matrix of M c is then

B. The Distributed Consensus-Based Control System
The consensus algorithm [29] is an effective technique to perform a distributed solution of the KKT condition in optimization problems (the control objective (6)) [5]. Accordingly, we choose the distributed consensus-based integral controlleṙ where u c i is the data shared between the DGs; x i is the controller state; k I i > 0 is the integral gain; a ij is the communication weight between DGs i and j, defined in the previous subsection. Let us define x c = col{x c i } ∈ R n G e , u c = col{u c i } ∈ R n G e , and k I = diag{k i } ∈ R n G e ×n G e . With the Hamiltonian H c (x c ) = 0.5x c k −1 I x c , this controller can then be represented as the PH system below.
where L is the Laplacian matrix of the communication network. Now if Σ c (8b) has a feasible equilibrium point, then the incremental model of this linear system can be written as Therefore, one can writeḢ c (x c ) =ỹ cũc . Hence, with the storage function H c (x c ), the control systemΣ c (8c) is also passive (lossless) with the inputũ c and outputỹ c .

C. Control by Interconnection of the Incremental Systems
Now that the incremental model of both physical and control systems are represented as PH systems, one can couple them through the following subsystem.
where b = col{b i } and b c = col{b c i } are constant vectors in R n G e ; r and w are square matrices belonging to R n G e ×n G e . Assumption 4: The systems Σ (2) and Σ c (8b) have feasible equilibrium points which are coupled through the subsystem Σ I (9a) as follows.
If Assumption 4 holds, then the incremental model of Σ I (9a) can be written as the following lossy interconnection subsystem [25].
Therefore, one hasỹ The configuration used for control by interconnection [25] of the systems is shown in Fig. 2. Proposition 2: Consider the PH system Σ (2) coupled with the controller Σ c (8b) through the interconnection subsystem Σ I (9a) (See Fig. 2). If Assumption 4 holds and the matrix R D + r is positive-definite, then the equilibrium point of the closed-loop system is asymptotically stable in the region Fig. 2. Block (circuit) diagram of the control by interconnection scheme for both non-incremental (a) and incremental (b) system models.
Proof: Consider the total storage function H t (x t ) = H(x) + H(x c ) for the incremental model of the closed loopsystem which has a minimum at the equilibrium point. Taking its derivative and using (3), (8c), and (9d), one haṡ According to (4), [G cte − G P (q Ne )] is positive-definite for allx t ∈ S, if the closed-loop system has a feasible equilibrium point (Assumption 4 holds). Moreover, the matrices R Ee and R Ge are also positive-definite. Therefore, if R D + r is positive-semi-definite, then T(x t ) is positive-definite anḋ H t ≤ 0, ∀x t ∈ S which proves that the equilibrium point is stable in S [30], with Lyapunov function H t . On the other hand, positive-definiteness of H t ensures that ∃ζ > 0 such that the level set Ω ζ = {x t ∈ S : H t ≤ ζ} is bounded. Since its requirements are all satisfied, LaSalle's theorem can be applied [30]. According to LaSalle's theorem, every solution starting in Ω ζ converges to the largest invariant set, say M, in E = {x t ∈ Ω ζ :Ḣ t = 0}; i.e.,x t ∈ M ⊆ E as t → ∞. Since T(x t ) is positive-definite ∀x t ∈ S and H(x) is quadratic, according to (10), one can write E = {x t ∈ Ω ζ :x = 0, ∇H(x) = 0} implying thatẋ = 0. Therefore, by using (3), (8c), and (9c) it is easy to observe that the motion in this invariant set is governed byẋ c = 0, ∀x t ∈ E. In other words, the largest invariant set in E is the equilibrium point; i.e., M = {x t ∈ Ω ζ :x = 0,x c = 0}. Therefore, LaSalle's theorem implies asymptotic stability of the equilibrium point in S.

D. Equilibrium (Steady State) Analysis
Proposition 3: Let Assumption 4 hold. Then, if the communication network is connected, the KKT optimality condition in (6) and the near-nominal voltage formation in (7) with w i = α −1 i are simultaneously achieved at the equilibrium point of the closed-loop system.
Proof: According to (8b), at equilibrium point one has Lū c = 0, where we used the fact that k I is positive-definite. If the communication network is connected, then L has a simple zero eigenvalue [29] and thereforeū c = u opt 1 is the unique solution of Lū c = 0, where u opt is the consensus value. Thus, according to (9b) one can write Let us define λ = col{λ i } ∈ R n G e , β = col{β i } ∈ R n G e and α = diag{α i } ∈ R n G e ×n G e . The KKT condition (6) can then be written asλ = 2αȳ + β = λ opt 1.

E. Implementation of the Proposed Controller
The proposed controller is presented in matrix form so far. However, in what follows, to better understand its practical implementation, it is formulated in a non-matrix format in terms of the required measurements, parameters, and communication data. If w = 0.5α −1 , r = −R D + k P 2αL2α, b c = β, and b = −k P 2αLβ, then considering (2) and defining z λ = col{z λ i } = −Lλ, and z c = col{z c i } = −Lx c , the control system (8b) coupled with the interconnection subsystem (9a) can be written as which can be written in the following scalar format. Fig. 3 depicts a schematic diagram of the proposed controller described in (12). One can see that except for x c j and λ j , received from the neighboring DGs, the other parameters and variables are locally available for each DG.

IV. CASE STUDIES AND RESULTS
To show the effectiveness of the proposed controller, it is tested on a 48-Volt meshed dc MG, powered by six DGs. The DGs with odd (resp. even) numbers are interfaced to the grid via buck (resp. boost) converters, which are depicted in Fig. 4 by circles (resp. squares). The electrical and control specifications of the MG shown in Fig. 4 are given in Table I. Remark 2: According to Assumption 1, to design the sec-  ondary controller, the converters are modeled by an equivalent zero-order model as in (1e); thus, the converter dynamics and its internal voltage controller are hidden in Fig. 1 under the dashed blue box. However, in the simulations, Linear Quadratic Regulator (LQR) controller technique is used for the voltage V i to track its reference V ref i [31]. Fig. 5 depicts the converter dynamics and the internal voltage controller. The resistance R i , inductance L i , and capacitance C i of all the converters are 0.1Ω, 2.64mH, and 2.2mF , respectively; the input voltage to the converters V in i of the DGs 1 to 6 are 80, 25, 100, 20, 80, 25 V , respectively; I i , ζ i , and V i are the states of the system, m i is the duty cycle given to the PWM generator to produce the switching signal g i with frequency of 5kHz. To design proper feedback gain matrix K i ∈ R 3×3 , the linearized second-order average model of converters augmented with a voltage-tracker integrator, is used where the output current of the converter capacitor I Ge i is considered as an external disturbance, along the lines of [31].
A. Controller Performance: Activation and Load Change Fig. 6 depicts the performance of the MG under the proposed controller in different stages. Before t = 5s, the MG is operated without the proposed secondary control. Therefore, the DGs voltages are settled away from their nominal voltages so that their average value is deviated from the nominal value 48V . Moreover, the incremental costs of the DGs have different values which underlines the KKT condition is not satisfied. After activating the controller at t = 5s, the DGs reach a consensus on their incremental costs and at the same time they form their voltages around the nominal value with a weighted average of nominal voltage. It should be noted that before t = 14s, only constant impedance and constant current loads are energized. To emphasis the resiliency of the controller, at t = 14s, the constant power loads at all the buses are activated. One can see that the DGs reach an agreement on a new optimal incremental cost higher than the previous one, which returns to the previous value after deactivating the constant power loads at t = 19s. It should be emphasized that, over the load change transitions, the average voltage remains unchanged and only transient voltage drifts from the nominal voltage are observed.

B. Controller Performance: Plug-and-Play Ability
To show the DGs plug-and-play ability under the proposed controller, the 4th DG is disconnected from the grid at t = 24s and it is connected back to the grid at t = 29s. To do so, a corresponding circuit breaker is opened at t = 24s to disconnect the DG physically and the communication links related to the DG are all interrupted. Moreover, before closing the breaker at t = 29s, all the communication links are restored and both sides of the breaker are voltage-synchronized for seamless connection of the DG. According to Fig. 6, after disconnecting 4th DG from the grid, other DGs inject more current so they reach consensus on a new optimal incremental cost. Furthermore, one can see that the average voltage of the remaining five DGs still operate at the nominal value while the fourth DG voltage drops to the voltage of the bus number 4. It is also shown that after connecting it back to the grid, the DG immediately participates in the current sharing and voltage formation tasks as before.

C. Real-Time Results From OPAL-RT
To verify the real-time effectiveness of the proposed controller, the previous system is built and loaded to an OPAL-RT OP5600 real-time simulator, shown in Fig. 7. It should be pointed out that, therein, the detailed switching model of the Buck and Boost converters with switching frequency of 5kHz are employed. The selected IGBTs and Diodes have internal resistance of 1mΩ and forward voltage of 0.8V . The other (passive) components of the converters and their inner voltage controllers are exactly the same as described in the preamble of this Section (Remark 2). Fig. 8 indicates alignment of the real-time system responses with the simulation results in Section IV-A. Due to the input limitation of the oscilloscope only the results for the DGs 1, 2, 5, and 6 are given. After activating the controller, the incremental costs reach a consensus and the voltages reach a formation around the nominal value so that their weighted average settles at the nominal value. The results for the load increase scenario further approves the effectiveness of the proposed control in reaching the current-sharing and voltageformation control goals, under severe load changes.

V. CONCLUSIONS
A distributed secondary control technique is proposed for dc MGs with ZIP loads which drives the MG to a point where the KKT optimality condition is satisfied for all the DGs and their weighted average voltage is the nominal value. The closedloop system (the MG engaged with the proposed controller) is formulated in a port-Hamiltonian representation which is shown to be asymptotically stable by using Lyapunov and LaSalle theorems. It is also shown that the system is globally asymptotically stable without the constant power loads. The effectiveness of the proposed controller for different case studies is verified by adapting it to a test system through both non-real-time and real-time simulations. It should be noted that for the theoretical analyses each DG is modeled by an equivalent zero order model as a controllable voltage source, while, in MATLAB/Simulink simulation and OPAL-RT model the average model and detailed switching model are used, respectively. All in all, the theoretical analyses and case studies demonstrate effectiveness of the proposed controller in achieving the desired control goals.