Decentralized DC MicroGrid Monitoring and Optimization via Primary Control Perturbations

TL;DR

This paper presents a decentralized method for DC MicroGrid monitoring and optimization that uses local voltage measurements and primary control perturbations, eliminating the need for external communication.

Contribution

It introduces a novel decentralized estimation approach integrated within primary control loops for autonomous MicroGrid operation.

Findings

Accurate estimation of generation capacities, load demands, and line conductances.

Effective decentralized solution for economic dispatch.

Validation through simulations demonstrating practical applicability.

Abstract

We treat the emerging power systems with direct current (DC) MicroGrids, characterized with high penetration of power electronic converters. We rely on the power electronics to propose a decentralized solution for autonomous learning of and adaptation to the operating conditions of the DC Mirogrids; the goal is to eliminate the need to rely on an external communication system for such purpose. The solution works within the primary droop control loops and uses only local bus voltage measurements. Each controller is able to estimate (i) the generation capacities of power sources, (ii) the load demands, and (iii) the conductances of the distribution lines. To define a well-conditioned estimation problem, we employ decentralized strategy where the primary droop controllers temporarily switch between operating points in a coordinated manner, following amplitude-modulated training sequences. We study the use of the estimator in a decentralized solution of the Optimal Economic Dispatch problem. The evaluations confirm the usefulness of the proposed solution for autonomous MicroGrid operation.

Tables

Table 1. TABLE I: Fixed Simulation Parameters

Parameter	Value
Simulation platform	MATLAB
Reference voltage $x$ (volts)	$400$
Lower and upper voltage margins $v_{\min }$, $v_{\max }$ (volts)	$385$, $415$
Distribution network topology	cut-ring
Max. gen. capacity per DER $g$ (kW)	$1$
Max. const. conductance demand per bus $d^{\text{ca}}$ (kW)	$0.2$
Max. const. current demand per bus $d^{\text{cc}}$ (kW)	$0.2$
Max. const. power demand per bus $d^{\text{cp}}$ (kW)	$0$
Average conductance per line $y$ (S)	$1$
Sampling frequency ${\varphi }_S$ (kHz)	$50$
Sampling noise standard dev. ${\sigma }_S$ (volts/sample)	$0.1$
Transient time duration ${\tau }^{\text{transit}}$ (ms)	$2.5$
Nominal droop control params. $\tilde{x}$, $\Delta \tilde{v}$ (volts)	$400$, $15$
Max. voltage drop in $M$-phase $\Delta v$ (volts)	$15$
Total number of slots in the training epoch	$T=600$
Total number of slots in sub-phase $\alpha$	$T_{\alpha }=2N$
Total number of slots per block in sub-phase $\beta$	$L=2N$
Other $C$-phase params. ${\kappa }^{\alpha }$, ${\kappa }^{\beta }$, ${\chi }_n,n\in 𝒩$	$1$, $1$, ${\overline{v}}_n$
OED epoch duration ${\tau }^{\text{OED}}$ (s)	$300$
Backup gen./storage cost $c_{\text{source}}^{\text{extra}}/c_{\text{storage}}^{\text{extra}}$ (units/W)	$12$
DC bus signaling threshold $\xi$	$6.25\cdot {10}^{-4}$

Equations

$$[\mathbf{Y}{]}_{n,m}$$

$$i_n^{\text{cp}}\approx \frac{2d_n^{\text{cp}}}{v_n},y_n^{\text{cp}}\approx -\frac{v_n^2}{d_n^{\text{cp}}},n\in \mathcal{N}.$$

$$v_n=x_n-(y_n^{\text{va}}{)}^{-1}{i}_n,n\in \mathcal{N}.$$

$$v_{\min }<x_n\le v_{\max },y_n^{\text{va}}=\frac{g_n}{(x_n-\Delta v_n)\Delta v_n}\equiv s_n{g}_n,$$

$${\omega }_n=0,n\in \mathcal{N},$$

$${\omega }_n=v_n^2\left({\zeta }_n{y}_n^{\text{va}}+\frac{1}{x^2}{d}_n^{\text{ca}}+\sum _{m\in \mathcal{N}}{y}_{n,m}\right)-v_n\sum _{m\in \mathcal{N}}{v}_m{y}_{n,m}-v_n\left({\zeta }_n{x}_n{y}_n^{\text{va}}-\frac{1}{x}{d}_n^{\text{cc}}\right)+d_n^{\text{cp}}-(1-{\zeta }_n)p_n.$$

$$\mathbf{\theta }=[{\mathbf{g}}^{\mathsf{T}},{\mathbf{d}}^{\mathsf{T}},{\mathbf{\psi }}^{\mathsf{T}}{]}^{\mathsf{T}},$$

$${\tilde{x}}_n\equiv \tilde{x},\Delta {\tilde{v}}_n\equiv \Delta \tilde{v},{\tilde{s}}_n\equiv \tilde{s},n\in \mathcal{N}.$$

$$\mathbf{\Omega }$$

$$\displaystyle\mathbf{\Omega}=\bigg{(}\mathbf{S}\mathsf{D}(\mathbf{g})+\frac{1}{x^{2}}\mathbf{1}_{T}(\mathbf{d}^{\text{ca}})^{\mathsf{T}}\bigg{)}\odot\mathbf{V}\odot\mathbf{V}+(\mathbf{V}\mathbf{Y})\odot\mathbf{V}-\bigg{(}(\mathbf{S}\odot\mathbf{X})\mathsf{D}(\mathbf{g})-\frac{1}{x}\mathbf{1}_{T}(\mathbf{d}^{\text{cc}})^{\mathsf{T}}\bigg{)}\odot\mathbf{V}+\mathbf{1}_{T}(\mathbf{d}^{\text{cp}})^{\mathsf{T}}.$$

$$\mathbf{W}=\mathbf{V}+\mathbf{Z},$$

$$\rho (\mathsf{vec}(\mathbf{W});\mathbf{\theta })=\mathsf{N}(\mathsf{vec}(\mathbf{V}),{\sigma }^2{\mathbf{I}}_{NT}).$$

$$\mathbf{W}=\left[\begin{array}{c}\overline{\mathbf{W}}\\ {\mathbf{W}}^{\alpha }\\ {\mathbf{W}}^{\beta }\end{array}\right],{\mathbf{W}}^{\beta }=\left[\begin{array}{c}{\mathbf{W}}^{\beta ;1}\\ ⋮\\ {\mathbf{W}}^{\beta ;\overline{T}}\end{array}\right].$$

$$\mathsf{rank}({\mathbf{{\rm Y}}}_{-n})$$

$$\mathsf{rank}(\mathbf{\Gamma })$$

$$\mathsf{vec}(\overline{\mathbf{\Omega }})=\mathbf{{\rm Y}}\mathbf{\theta }=0_{\overline{T}N}.$$

$$x_n(t)=\tilde{x}+\sqrt{{\pi }_n(t)}\Delta x_n(t),s_n(t)=\tilde{s},n\in \mathcal{N},t\in {\mathcal{T}}^{\alpha /\beta },$$

$${\mathbf{X}}^{\alpha /\beta }=\tilde{x}{1}_{T^{\alpha /\beta }\times N}+{\mathbf{\Pi }}^{\alpha /\beta }\odot \Delta {\mathbf{X}}^{\alpha /\beta },$$

$${\mathbf{w}}_n^{\alpha /\beta }\approx {\tilde{v}}_n{1}_{T^{\alpha /\beta }}+({\mathbf{\Pi }}^{\alpha /\beta }\odot \Delta {\mathbf{X}}^{\alpha /\beta }){\mathbf{h}}_n+{\mathbf{z}}_n^{\alpha /\beta }.$$

$$\sqrt{{\pi }_n(t)}=\sqrt{{\pi }^{\alpha }},n\in \mathcal{N},t\in {\mathcal{T}}^{\alpha }.$$

$$\sqrt{{\pi }_n(t)}=\sqrt{{\pi }^{\beta }}({\overline{w}}_n(b)-{\chi }_n),n\in \mathcal{N},t\in {\mathcal{T}}^{\beta ;b},b\in \overline{\mathcal{T}},$$

$$(\Delta {\mathbf{X}}^{\alpha }{)}^{\mathsf{T}}{1}_{T^{\alpha }}=0_N,(\Delta {\mathbf{X}}^{\alpha }{)}^{\mathsf{T}}\Delta {\mathbf{X}}^{\alpha }={\delta }^{\alpha }{\mathbf{I}}_N,$$

$$(\Delta {\mathbf{X}}^{\beta ;b}{)}^{\mathsf{T}}{1}_L=0_N,(\Delta {\mathbf{X}}^{\beta ;b}{)}^{\mathsf{T}}\Delta {\mathbf{X}}^{\beta ;b}={\delta }^{\beta }{\mathbf{I}}_N,$$

$$\mathsf{vec}({\overline{\mathbf{W}}}_{(n)})=\frac{\sqrt{{\pi }^{\alpha }}{\delta }^{\alpha }}{\sqrt{{\pi }^{\beta }}{\delta }^{\beta }}{\mathsf{D}}^{-1}({\mathcal{X}}^{\alpha }{\mathbf{w}}_n^{\alpha })\sum _{b\in \overline{\mathcal{T}}}({\mathcal{X}}^{\beta ;b}{\mathbf{w}}_n^{\beta ;b})+\mathcal{I}\mathbf{\chi },$$

$$\rho (\mathsf{vec}({\overline{\mathbf{W}}}_{(n)});\mathbf{\theta })\approx \mathsf{N}(\mathsf{vec}(\overline{\mathbf{V}}),\mathbf{\Sigma }).$$

$$\mathbf{\Sigma }={\sigma }^2\left({\mathbf{I}}_{\overline{T}N}+\frac{{\pi }^{\alpha }({\delta }^{\alpha }{)}^2}{{\pi }^{\beta }{\delta }^{\beta }}{\mathsf{D}}^{-2}({\mathcal{X}}^{\alpha }{\mathbf{v}}_n^{\alpha })+\frac{{\pi }^{\alpha }({\delta }^{\alpha }{)}^3}{{\pi }^{\beta }({\delta }^{\beta }{)}^2}{\mathsf{D}}^{-2}({\mathcal{X}}^{\alpha }{\mathbf{v}}_n^{\alpha })\left(\sum _{b\in \overline{\mathcal{T}}}\mathsf{D}({\mathcal{X}}^{\beta ;b}{\mathbf{v}}_n^{\beta ;b})({\mathbf{I}}_N\otimes 1_{\overline{T}\times \overline{T}})\mathsf{D}({\mathcal{X}}^{\beta ;b}{\mathbf{v}}_n^{\beta ;b})\right){\mathsf{D}}^{-2}({\mathcal{X}}^{\alpha }{\mathbf{v}}_n^{\alpha })\right).$$

$$\mathbf{\vartheta }=\left[\begin{array}{c}\mathbf{\theta }\\ \mathsf{vec}(\overline{\mathbf{V}})\end{array}\right].$$

$${\hat{\mathbf{\vartheta }}}_{-n}$$

$$\mathsf{vec}(\overline{\mathbf{\Omega }})=0_{\overline{T}N},$$

$$\mathsf{vec}(\overline{\mathbf{\Omega }})\approx {\mathbf{{\rm Y}}}^{(j)}\mathbf{\theta }+{\mathbf{\Gamma }}^{(j)}(\mathsf{vec}(\overline{\mathbf{V}})-\mathsf{vec}({\overline{\mathbf{V}}}^{(j)})),$$

$${\mathbf{\theta }}_{-n}=-(({\mathbf{{\rm Y}}}_{-n}^{(j)}{)}^{\mathsf{T}}({\mathbf{\Gamma }}^{(j)}\mathbf{\Sigma }({\mathbf{\Gamma }}^{(j)}{)}^{\mathsf{T}}{)}^{-1}{\mathbf{{\rm Y}}}_{-n}^{(j)}{)}^{-1}({\mathbf{{\rm Y}}}_{-n}^{(j)}{)}^{\mathsf{T}}({\mathbf{\Gamma }}^{(j)}\mathbf{\Sigma }({\mathbf{\Gamma }}^{(j)}{)}^{\mathsf{T}}{)}^{-1}({\mathbf{\upsilon }}_n^{(j)}{g}_n+({\mathbf{\Gamma }}^{(j)}{)}^{\mathsf{T}}(\mathsf{vec}({\overline{\mathbf{W}}}_{(n)})-\mathsf{vec}({\overline{\mathbf{V}}}^{(j)}))),$$

$$\mathsf{vec}(\overline{\mathbf{V}})=\mathsf{vec}({\overline{\mathbf{W}}}_{(n)})-\mathbf{\Sigma }({\mathbf{\Gamma }}^{(j)}{)}^{\mathsf{T}}({\mathbf{\Gamma }}^{(j)}\mathbf{\Sigma }({\mathbf{\Gamma }}^{(j)}{)}^{\mathsf{T}}{)}^{-1}({\mathbf{{\rm Y}}}^{(j)}\mathbf{\theta }+({\mathbf{\Gamma }}^{(j)}{)}^{\mathsf{T}}(\mathsf{vec}({\overline{\mathbf{W}}}_{(n)})-\mathsf{vec}({\overline{\mathbf{V}}}^{(j)}))).$$

Full text

I Introduction

Since their inception, MicroGrids (MGs) have evolved substantially, particularly in the domain of low voltages (LV), leading to variety of use cases and topologies [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]: from small clusters of distributed energy resources (DERs) serving houses or buildings, to large meshes of small MGs covering large areas, such as neighborhoods, industrial complexes and remote villages. As a result, the future smart grid (SG) is envisioned as a mesh of interconnected autonomous MG systems. It is also within the field of MGs where direct current (DC) power networks have experienced a renaissance due to the seamless integration with DC renewable generation, DC energy storage systems and DC smart loads [2, 3, 4]. Hence, LV DC MGs are considered as a solution for residential and industrial use cases.

A distinctive characteristic of DC MGs is the use of programmable DC/DC and AC/DC power electronic converters (PECs) to connect the DERs to the DC distribution system. PECs are digital signal processors (DSPs) that allow for software implementation of advanced control systems [2, 3]. Leveraging on the advanced features of PECs the control system design also shifted from simple strategies, suitable for small systems [12, 13, 14], to modular hierarchical architectures where several interacting control layers dynamically respond to state variations on different time scales and pursue various complementary objectives [3, 4, 15, 16, 17, 18, 19, 20, 21, 22, 23]. Specifically, the MG control plane is organized into dual-layer architecture, comprising primary and upper control layer [3, 15]. The primary control is decentralized and deals with high frequency dynamic compensation and state regulation [3]. The upper control layer deals with slow, global changes in the MG by providing updated primary control references and is implemented in distributed/centralized fashion [15, 16, 17, 18, 19, 20, 21, 22, 23]. An exemplary upper layer application is the Optimal Economic Dispatch (OED), which aims to compute the optimal dispatch policies that minimize the total generation cost while keeping the load balanced [17].

The standard design assumption is that the feedback of the upper control layer is closed via an external communication system, usually via off-the-shelf wireless technologies [3, 17, 21]. However, this approach was challenged recently due to several issues [2, 3]. First, the distributed power systems, particularly MGs, are dynamic and ad-hoc in nature, thus the installation of communication hardware may prove impractical and cost inefficient. Second, the external communication system reduces the resilience of the overall MG system, as it becomes a factor in the system reliability/availability. Finally, there is a growing concern about the cyber-security of power systems that exploit external communications, as the related security threats and attacks might severely compromise their stability and operation, leading to blackouts, equipment damage, data theft and investment losses [24, 25, 26, 27].

A straightforward solution would be to remove the upper layer completely and run the DC MG only with primary control without any further coordination. However, the approach is not suitable for advanced MG topologies, as it can not foster optimal and sustainable regulation. The DC bus signaling has been introduced as an enhancement of the above idea [12, 13, 14]. It uses the variations of the steady state bus voltage as an implicit coordination signal that tells the DERs how to behave in specific conditions. The idea is motivated by the fact that DC systems are inherently tolerant to steady state voltage variations, allowing for voltage ripples of up to $10\%$ [2, 3, 5]. Each PEC monitors the local voltage and if it the crosses predefined threshold, the PEC takes predefined actions. This approach has reliability, availability and security advantages over traditional networked design and requires only software modifications of the PECs. However, it is configuration-dependent, performing well in environments with predictable loads, but not in large, dynamic and general-purpose MGs. Moreover, the range of upper layer applications that can be supported is limited. Another alternative to wireless communications is to use conventional powerline communications (PLC) [28]. This way, some of the security concerns can be alleviated as now an attacker would need physical access to the MG. Nevertheless, PLCs are still essentially an external communication system coupled to the control of the MG, as they require installation of dedicated modems.

Motivated by the shortcomings of the above approaches, we propose a decentralized dual-layer control architecture for autonomous DC MGs in which each primary controller locally acquires the information required for the operation of the upper layer and determines the updated primary control references without the support of external communication enabler. To support the majority of applications, the upper control layer requires information about: i) the generation capacities of the dispatchable DERs, ii) the demands of the loads, and iii) the conductance matrix of the distribution network [17, 19]. This information can be inferred from local voltage observations, since the bus voltages are functionally related to the MG parameters through a non-linear model. To extract these parameters, the PECs deliberately move the MG through a sequence of sub-optimal states via coordinated and amplitude-modulated perturbations of the primary control parameters, referred to as training sequences. This way, the PECs obtain sequences of local bus voltage measurements from which the required information can be uniquely estimated, provided that the training sequences satisfy sufficiency criteria. To this end, we formulate a constrained Maximum Likelihood (ML) estimation problem that estimates the MG parameters jointly with the state of the DC MG. To solve the non-convex optimization problem, we develop an iterative algorithm and compare its performance against the Cramer-Rao Lower Bound (CRLB). We illustrate the practical potential of the method by applying it in decentralized OED (DOED) and we show how to minimize the operational cost by optimizing the design of the training sequences. The proposed solution does not rely on any additional communication hardware, as it exploits the signal processing capabilities of the PECs and its locally available voltage measurements, such that it can be implemented only in software.

The rest of the paper is organized as follows. Section II gives an overview of the main contributions. Section III introduces the system model. Section IV presents the training protocol and formulates the decentralized system identification problem. Section V is the pivotal section of the paper, presenting our take to the problem formulated in Section IV. Section VI introduces the periodic DOED protocol. Section VII presents the results and Section VIII concludes the paper.

Notation: Column vectors and matrices are denoted by lowercase and uppercase bold letters, e.g., $\mathbf{a}\in\mathbb{R}^{N\times 1}$ and $\mathbf{A}\in\mathbb{R}^{N\times M}$. $\mathbf{a}_{-n}\in\mathbb{R}^{(N-1)\times 1}$ is obtained from $\mathbf{a}$ by removing the element at position $n$. Similarly, $\mathbf{A}_{-m}\in\mathbb{R}^{N\times(M-1)}$ is obtained from $\mathbf{A}$ by removing the $m$-th column $\mathbf{a}_{m}$. $(\cdot)^{\mathsf{T}}$, $(\cdot)^{\dagger}$, $\mathsf{vec}(\cdot)$, $\mathsf{dim}(\cdot)$, $\mathsf{rank}(\cdot)$, $\mathsf{trace}(\cdot)$ and $\|\cdot\|_{l}$ denote the transpose, the pseudo-inverse, the vectorization, the dimension, the rank, the trace and the $l$-norm of the argument. $\otimes$ denotes the Kroneker product while $\odot$ and $\oslash$ denote the Hadamard (element-vise) product and division of vectors/matrices of adequate dimensions. The vectors $\mathbf{1}_{N}$, $\mathbf{0}_{N}$ and $\mathbf{e}_{n},~{}n\in\mathcal{N}$, denote the all-one, all-zero and the principal coordinate vector, $\mathbf{1}_{N\times M}$, $\mathbf{0}_{N\times M}$ denote the $N\times M$ all-one and all-zero matrices, and $\mathbf{I}_{N}$ is the $N\times N$ identity matrix. $\mathsf{D}(\mathbf{a})$ denotes diagonal matrix with the entries of $\mathbf{a}$ on the main diagonal. We frequently use the identity $\mathsf{vec}(\mathsf{D}(\mathbf{a}))=\mathbf{O}_{N}\mathbf{a}$ where the $N^{2}\times N$ matrix $\mathbf{O}_{N}=\sum_{n=1}^{N}\mathbf{e}_{n}\otimes(\mathbf{e}_{n}\mathbf{e}_{n}^{\mathsf{T}})$.

II Overview of Contributions

The proposed solution is illustrated in Fig. 1. We consider a generic DC MG model with multiple buses, described in Sections III and IV. We assume that the MG does not have access to reliable external communication resources. The physical state of DC MGs is characterized by the steady state bus voltages. We introduce a parameter vector $\boldsymbol{\theta}$ that collects all system variables whose values are determined by exogenous influences; this includes the generation capacities of the DERs, the load demands and the distribution network topology, i.e., the conductance matrix, see Section IV-A. Using the power balance equation, we represent the bus voltages thorough a non-linear and implicit model, parametrized by $\boldsymbol{\theta}$, see Section III-B. Evidently, $\boldsymbol{\theta}$ varies with time; to respond to its variations on different time scales, the DC MG is governed by a hierarchical control system, comprising primary and upper control layer. The primary control is decentralized: several controllers regulate the bus voltages, using only local feedbacks without exchanging any information with peer controllers. They are very fast and capable of responding to high frequency variations in $\boldsymbol{\theta}$. Popular primary controller in DC MGs is the Voltage Source Converter (VSC) with voltage droop control, which is reminiscent to the widespread frequency droop control in AC systems, but defined over the DC voltage; it is therefore standard practice to refer to it simply as droop controller [2, 3]. The upper control layer, on the other hand, responds to less frequent changes in $\boldsymbol{\theta}$ that affect the global behavior of the system; examples include changes of the load/generation profile, faults, attacks, etc. Its main role is to adapt the system to the new conditions by computing updated optimal control references for the primary controllers; all upper layer control applications require full/partial knowledge of $\boldsymbol{\theta}$ to determine the control references that adequately reflect the new conditions [15, 16, 17, 18, 19, 20, 21, 22, 23].

Unlike conventional centralized networked control solutions, where the upper control layer is supported by an external communication enabler, we propose a decentralized control architecture that relies solely on the DSP capabilities of the PECs: namely, in our solution the upper control layer is implemented locally within each PEC, and uses only the locally available state measurements, as depicted in Fig. 1. The solution comprises two main functional blocks, i.e., the monitoring and optimization, executed sequentially.

Monitoring. This functional block exploits the fact that the steady state bus voltages are functionally related with $\boldsymbol{\theta}$ through the power balance equation; hence, each controller can compute a local estimate of $\boldsymbol{\theta}$. The key challenge is that it is impossible to infer $\boldsymbol{\theta}$ by using only local measurements of a single realization of the state, as the system is not observable and the estimation is ill-conditioned. To address this, the monitoring block comprises two procedures: (1) coordinated decentralized training [29, 30] via primary control perturbations, see Section IV, and (2) Joint System Identification and State Estimation (J-SISE), see Section V. During training, the controllers perturb the values of the local droop control parameters, for a limited period of time, following predetermined training sequences. This generates a sequence of different realizations of the state. The controllers collect the local measurements of the state sequence and modulate them into the perturbation signals, see Section V-C. In other words, the relation between the primary control perturbation signals and the induced state deviations is interpreted as the input-output relation of an implicit communication channel [31, 32, 33, 34, 35], through which the controllers exchange their local observations. Hence, the training sequences are used both for generating multiple states and communicating the local state observations. If the training sequences satisfy sufficiency criteria, see Section V-B, each controller is able to compute unique estimate $\hat{\boldsymbol{\theta}}$ using the steady state voltage measurements acquired during training and the J-SISE algorithm, see Section V-D. The J-SISE is formulated as non-convex, constrained ML optimization problem in classical estimation framework which we solve via iterative algorithm based on partially linearized constraints and evaluate its performance using the CRLB, see Sections V-E and VII-B.

Optimization. The local estimates $\hat{\boldsymbol{\theta}}$ are used as inputs to an energy management application which computes updated primary control references, see Fig. 1. Any application for which $\boldsymbol{\theta}$ is sufficient can be applied. We focus on DOED with linear generation cost model, since a simple, decentralized closed form solution is available in this case [17, 35]. To this end, we design periodic protocol, detailed in Section VI, where the controllers first perform training and obtain $\hat{\boldsymbol{\theta}}$ via J-SISE, then re-dispatch. Finally, we show how to minimize the operational cost of the protocol by calibrating the training parameters, see Section VII-C.

We conclude by highlighting the benefits of the proposed solution. First and foremost, it promotes the principle of self-sustainability in SG as it reuses the DSP features of the available power electronics and obviates critical reliance on external communication system. Further, the optimization block is not limited only to OED, as the knowledge of $\boldsymbol{\theta}$ allows each controller to solve locally a great deal of energy management optimizations (even if they do not have decentralized formulation) such as Optimal Power Flow (OPF), Unit Commitment (UC) and security-related applications, such as Fault Detection and Diagnosis (FDD) [21, 19]. This flexibility strengthens the autonomous operation of the DC MG. Finally, the developed framework can be adapted for arbitrary DC MG systems, as discussed in Section VII-B.

III System Model

The terminology and the notation system applied to the model is standardly used in power engineering literature [3, 19]. Section IV introduces compact, matrix notation of the power balance equation which is easier to manipulate later on; this can be also seen as a standalone contribution, as this is the first work that introduces such compact notation for droop-controlled DC MG.

III-A General Multiple-Bus DC MicroGrid

III-A1 Buses and Distribution Network

A DC MG is a collection of DERs and loads, connected to low voltage DC distribution system, see Fig. 2. The distribution system consist of $N\geq 1$ buses, indexed in the set $\mathcal{N}=\left\{1,...,N\right\}$. Each bus $n$ in steady state is characterized by a bus voltage $v_{n}$, and all DERs and loads connected to bus $n$ measure the same voltage $v_{n}$. The distribution line connecting buses $n$ and $m,~{}n\neq m$ has a line conductance denoted by $y_{n,m},~{}y_{n,m}\equiv y_{m,n}\geq 0$ [3]. The topology of the distribution system is specified via the symmetric $N\times N$ conductance matrix $\mathbf{Y}$ with elements:

[TABLE]

III-A2 Distributed Energy Resources

We model each DER as separate bus, i.e., we assume that each bus hosts at most one DER; hence, the total number of DERs is $N$ and they are indexed in the set $\mathcal{N}$. This modeling choice simplifies the notation without losing generality; in fact, if DERs $n$ and $m$ are connected to the same physical point, i.e., the same bus, by definition $y_{n,m}=\infty$. The $n-$th DER has current $i_{n}$ and power output $p_{n}=v_{n}i_{n}$. We assume that the DERs in the MG are small-scale power sources such as renewables (RESs) or distributed generators (DGs) based on traditional fossil fuel. Each DER $n$ has an instantaneous generation capacity $g_{n}\geq 0$, and the output power $p_{n}$ should satisfy $0\leq p_{n}\leq g_{n}$.

III-A3 Loads

The $n-$th bus hosts a collection of loads, represented through an aggregate model as a mixture of three components (also known as ZIP load model [36]): 1) constant conductance $y_{n}^{\text{ca}}={x^{-2}}{d_{n}^{\text{ca}}}$, 2) constant current $i_{n}^{\text{cc}}={x^{-1}}{d_{n}^{\text{cc}}}$, and 3) constant power component $d_{n}^{\text{cp}}$, see Fig. 2. The quantities $d_{n}^{\text{ca}}$, $d_{n}^{\text{cc}}$ and $d_{n}^{\text{cp}}$ are the instantaneous power demands of the components at a rated voltage $x$. For a given $d_{n}^{\text{cp}}$, the constant power component in steady state is approximated with an equivalent positive current source in parallel with negative conductance and the electrical parameters are [3]:

[TABLE]

III-A4 Primary Control

The DERs use PECs to interface the buses; the bus voltage $v_{n}$ and/or current $i_{n}$, i.e., power $p_{n}$ are locally controlled through decentralized primary controller, which is a software program executed by the PEC [3]. Two primary control schemes, i.e., modes are commonly used, see Fig. 3: 1) a closed loop Voltage Source Converter (VSC), and 2) an open loop Current Source Converter (CSC). VSC regulates the bus voltage and current of the DER as the loads/generation in the system change in order to keep the bus voltage within predefined margins and foster fair power sharing. It contains fast inner and slow outer control loops. An inner control loop consists of a cascade of voltage and current loops with control bandwidth of the order of several tens of kHz, equal to the sampling frequency $\phi_{S}$ of the converter. Its role is to maintain the output bus voltage $v_{n}$ to specific reference value, dictated by the outer control loop. The outer control loop is closed via filtered current feedback, and is slower than the inner control loop by an order of magnitude. The current feedback generates the reference value for the inner voltage loop, via the following steady state control law:

[TABLE]

This is known as decentralized droop control for DC MGs [3, 15] with two controllable parameters: the reference voltage $x_{n}$ and the virtual conductance $y_{n}^{\text{va}}$. Their values are set (i) to keep the bus voltage, as closely as possible to the rated voltage $x$, within predefined margins $v_{\max}\leq v_{n}\leq v_{\min}$ for any $n\in\mathcal{N}$, and (ii) to enable fair power sharing among DERs based on their instantaneous generation capacities [3]. Fig. 4 depicts a widespread droop control law that meets the above conditions, with droop control parameters set as follows:

[TABLE]

where $s_{n}\equiv((x_{n}-\Delta v_{n})\Delta v_{n})^{-1}$ is the droop slope (in volts*-2*). The configuration $y_{n}^{\text{va}}=s_{n}g_{n}$ enables proportional power sharing among the DERs. When the DER operates close to its capacity, the maximal voltage drop is $\Delta v_{n}$, $0<\Delta v_{n}\leq x_{n}-v_{\min}$. In steady state, the droop-controlled VSC units are modeled as voltage sources in series with virtual conductance, see Fig. 2.

The other primary control mode CSC does not have outer control loop and inner voltage loops, see Fig. 3. The reference for the inner current loop is generated via a separate algorithm that gets as an input fixed power reference [3]. Hence, a CSC acts as a constant power component, neither participating in voltage regulation nor power sharing. It is modeled as a negative current source and parallel conductance, as in (4) but with opposite sign. It is architecturally equivalent to a negative constant power load, see Fig. 2.

The subsets of DERs operating in VSC/CSC, denoted respectively by $\mathcal{N}^{\text{V}}/\mathcal{N}^{\text{C}}$, are determined dynamically by the upper layer application, see Section VI for an example. To support this dynamic operation, each converter is assumed to have dual mode, and is capable to switch between VSC and CSC control mode seamlessly [3, 37], see also Fig. 3.

III-B Steady State Equations

A DC MG is governed by Ohm’s and Kirchhoff’s laws, resulting in a system of $N$ steady state power balance equations for $N$ buses:

[TABLE]

with $\omega_{n}$ given by:

[TABLE]

The binary variable $\zeta_{n}$ in (8) is $1/0$ if DER $n$ is configured in VSC/CSC control mode, respectively. The system of equations is quadratic in the bus voltages, such that, in general, a closed form solution for $v_{n},~{}n\in\mathcal{N}$ is not possible. The non-linear nature of the power balance equations stems from the presence of constant power components [38], both constant power loads and CSCs. Hence, in the case when $d_{n}^{\text{cp}}=0$ for all $n$ and $\mathcal{N}^{\text{C}}=\emptyset$, the system (7) becomes linear in the bus voltages. Another special case with closed-form solution is the Single-Bus DC MG which we have studied separately [34] due to its practical importance.

IV Problem Formulation and Training Epoch

The DC MG is not connected to an external communication system and the PECs only have the local voltage/current measurements to work with. To learn (i) the generation capacities of remote DERs, (ii) the power demands of the loads and, (iii) the conductances of the distribution lines, the controllers need to solve a decentralized system identification problem, formulated below.

Before we begin, we list the main assumptions:

1. ($A_{1}$)

The primary controllers are fully synchronized to a common time reference. 2. ($A_{2}$)

No prior knowledge on the generation capacities, load demands or the conductance matrix is used. 3. ($A_{3}$)

The rate of load/generation/topology variations is an order of magnitude smaller than the frequency of the primary controllers.

IV-A Parameter Vector

Let $\mathbf{g}=[g_{1},\ldots,g_{N}]^{\mathsf{T}}$ be a $N\times 1$ vector that collects the instantaneous generation capacities of all DERs in the MG. Similarly, the instantaneous load demands are collected in separate $N\times 1$ vectors: $\mathbf{d}^{\text{ca}}=[d_{1}^{\text{ca}},\ldots,d_{N}^{\text{ca}}]^{\mathsf{T}}$, $\mathbf{d}^{\text{cc}}=[d_{1}^{\text{cc}},\ldots,d_{N}^{\text{cc}}]^{\mathsf{T}}$ and $\mathbf{d}^{\text{cp}}=[d_{1}^{\text{cp}},\ldots,d_{N}^{\text{cp}}]^{\mathsf{T}}$. The $3N\times 1$ load demand vector is defined as $\mathbf{d}=[(\mathbf{d}^{\text{ca}})^{\mathsf{T}},\;(\mathbf{d}^{\text{cc}})^{\mathsf{T}},\;(\mathbf{d}^{\text{cp}})^{\mathsf{T}}]_{n\in\mathcal{N}}^{\mathsf{T}}$. Further, we observe that $\mathbf{Y}$ is fully specified by its supra(infra)-diagonal elements, see (3). We organize these elements in a vector $\boldsymbol{\psi}=[\ldots,y_{n,m},\ldots]^{\mathsf{T}}$, $n,m\in\mathcal{N}$, $m>n$, with dimension $\mathsf{dim}(\boldsymbol{\psi})=\frac{1}{2}N(N-1)\times 1$. Using $\boldsymbol{\psi}$, we can write $\mathbf{Y}$ as the weighted Laplacian $\mathbf{Y}=\mathbf{A}\mathsf{D}(\boldsymbol{\psi})\mathbf{A}^{\mathsf{T}}$, where $\mathbf{A}\in\left\{-1,0,1\right\}^{N\times\mathsf{dim}(\boldsymbol{\psi})}$ is the oriented incidence matrix [39].

The deterministic parameter vector $\boldsymbol{\theta}$ is defined as:

[TABLE]

with dimension $\mathsf{dim}(\boldsymbol{\theta})=\frac{1}{2}N(N+7)\times 1$. From the discussion in Section III-B, the steady state bus voltage $v_{n}$ depends on $\boldsymbol{\theta}$, see eq. (6), (8). This suggests that an arbitrary controller can infer the parameter vector $\boldsymbol{\theta}$ locally, using local measurements of the steady state bus voltage (see also [40] and references therein for similar approaches). However, it is impossible to determine $\boldsymbol{\theta}$ uniquely in classical, non-Bayesian estimation framework, using only a single observation of the local steady state bus voltage. To address this issue, the following subsection introduces a technique based on decentralized training via primary control perturbations.

IV-B Training Protocol and Training Sequences

We introduce a dedicated training epoch of predefined duration, in which (i) all controllers switch to VSC mode using a droop control law of the form (6), (ii) perturb their local droop control parameters, causing deviations of the bus voltages, and (iii) measure the local bus voltage response, collecting sequences of steady state bus voltage measurements. The training epoch design uses the assumptions $(A_{1})$ and $(A_{2})$. Specifically, the time axis during the training epoch is divided into $T$ time slots, see Fig. 5, and all controllers are synchronized to this structure. We index each slot with $t\in\mathcal{T}=\left\{1,\ldots,T\right\}$. The slot duration $\tau$ complies with the control bandwidth of the primary control loops, allowing the bus to reach a steady state after a transient time $\tau^{\text{transit}}\ll\tau$, yielding $\phi_{S}(\tau-\tau^{\text{transit}})$ voltage samples per slot for each controller, see Fig. 5. The system constant $\tau^{\text{transit}}$, usually several milliseconds [3], is determined by the sampling frequency $\phi_{S}$ and the line capacitors. Following $(A_{2})$, $\boldsymbol{\theta}$ can be assumed to remain constant during the training epoch.

We use $\tilde{\cdot}$ to denote the unperturbed, i.e., nominal droop control parameters during the training epoch; we use the law (6) with equal reference voltages and droop slopes:

[TABLE]

In slot $t$, all controllers simultaneously perturb the reference voltages and droop slopes, according to perturbation signals $x_{n}(t)\neq\tilde{x}$, $s_{n}(t)\neq\tilde{s}$, $n\in\mathcal{N}$; they are organized in $T\times N$ training matrices $\mathbf{X}$, $\mathbf{S}$, defined as $[\mathbf{X}]_{t,n}=x_{n}(t)$ and $[\mathbf{S}]_{t,n}=s_{n}(t)$, $n\in\mathcal{N}$, $t\in\mathcal{T}$. The columns $\mathbf{x}_{n}$/$\mathbf{s}_{n}$ of $\mathbf{X}$/$\mathbf{S}$, correspond to the training sequence injected by controller $n$.

IV-C Steady State Bus Voltages and Measurement Vectors

The steady state bus voltage $\tilde{v}_{n}$ corresponds to the nominal, unperturbed, droop parameters $\tilde{x}_{n}$, $\tilde{s}_{n}$. The steady state bus voltage response in the $t-$th slot is $v_{n}(t)\neq\tilde{v}_{n},~{}n\in\mathcal{N}$. The $T\times N$ steady state bus voltage matrix $\mathbf{V}$ is defined as $[\mathbf{V}]_{t,n}=v_{n}(t)$, $n\in\mathcal{N}$, $t\in\mathcal{T}$. The following proposition characterizes $\mathbf{V}$ in terms of $\mathbf{X}$, $\mathbf{S}$ and $\boldsymbol{\theta}$:

Proposition 1**.**

The steady state of DC MG during the training epoch is characterized by the implicit power balance equation:

[TABLE]

where $\mathbf{\Omega}:[v_{\min},\;v_{\max}]^{T\times N}\times\mathbb{X}\times\mathbb{S}\times\mathbb{R}^{\mathsf{dim}(\boldsymbol{\theta})}\mapsto\mathbf{0}_{T\times N}$ is defined as $[\mathbf{\Omega}]_{t,n}=\omega_{n}(t),~{}n\in\mathcal{N},~{}t\in\mathcal{T}$, and given by:

[TABLE]

The subsets $\mathbb{X}\subset\mathbb{R}^{T\times N}$ and $\mathbb{S}\subset\mathbb{R}^{T\times N}$ comprise all training matrices $\mathbf{X}$ and $\mathbf{S}$ that keep $\mathbf{V}$ within $[v_{\min},\;v_{\max}]^{T\times N}$.

Proof.

See Appendix A. ∎

The power balance equation (11) reflects the requirement to keep the system balanced and stable, i.e., in a valid (albeit suboptimal) operating point, in each slot during training. It also gives an implicit relation between $\mathbf{V}$ and $\boldsymbol{\theta}$, since (11) cannot be solved in closed form for $\mathbf{V}$.

The $n-$th controller measures the $n$-th column $\mathbf{v}_{n}$ of $\mathbf{V}$ during the training epoch. The noisy measurement obtained by controller $n$ in slot $t$ is an average of multiple voltage samples collected during the steady state period of the slot, and can be written as $w_{n}(t)=v_{n}(t)+z_{n}(t)$ with $z_{n}(t)$ denoting the additive noise. The $T\times N$ bus-voltage measurements matrix $\mathbf{W}$, with $[\mathbf{W}]_{t,n}=w_{n}(t)$, $n\in\mathcal{N}$, $t\in\mathcal{T}$, is given as:

[TABLE]

where $\mathbf{Z}$ represents the noise and $\mathsf{vec}(\mathbf{Z})$ is a zero-mean, white Gaussian random vector with standard deviation $\sigma$ [41], such that the probability density function (pdf) of $\mathsf{vec}(\mathbf{W})$ is:

[TABLE]

The decentralized system identification problem for DC MGs is about devising an efficient and unbiased estimator of the local parameter vector $\boldsymbol{\theta}_{-n}$, denoted with $\hat{\boldsymbol{\theta}}_{-n}$, using only local bus voltage measurements $\mathbf{w}_{n}$, for any $n\in\mathcal{N}$.

IV-D Relaxing Assumptions $(A_{1})-(A_{3})$

We briefly discuss the implications that arise when assumptions $(A_{1})-(A_{3})$ are no longer valid; addressing these implications is out of the paper’s scope. We start with $(A_{1})$, as the strongest assumption. Maintaining precise synchronization among the controllers on the level of slot and training epoch can be easily achieved if the PECs are equipped with GPS modules. Alternatively, one can use common decentralized network synchronization approaches, typically used in sensor networks [42]. Since the method operates in a time scale in the order of milliseconds, it should be significantly easier to maintain (at least coarse) synchronization for long periods of time. Finally, if synchronization is not possible, and the controllers inject perturbation signals without any prior coordination, then the formulation of the problem should be modified accordingly to account for asynchronous training. For instance, the parameter vector should be extended to include binary variables that capture the activity patterns of the controllers and the start times of individual training sequences, as well as their end times in case of variable training sequence durations.

Assumption $(A_{2})$ simply casts our problem in classical estimation framework. In practice, prior knowledge is always available to some extent; in fact, $\boldsymbol{\theta}$ can be assumed to evolve over time following a stochastic process, paving the way for formulating the identification problem in sophisticated Bayesian filtering/prediction framework [43]. Nevertheless, the analysis of the non-Bayesian case naturally comes first.

We use assumption $(A_{3})$ to postulate that $\boldsymbol{\theta}$ remains fixed during training, which is not true in general. In practice, $\boldsymbol{\theta}$ might change at any time due to load/generation variation or a system fault. To incorporate this notion we should reformulate the problem accordingly. One way is to first relax assumption $(A_{2})$ and model the dynamic evolution of $\boldsymbol{\theta}$ via stochastic process, where relaxing assumption $(A_{3})$ arises naturally. We can avoid relaxing $(A_{2})$ and still use the classical framework as presented in the paper, but with modified definition of the parameter vector. For instance, let us assume that $\boldsymbol{\theta}$ has changed no more than $J\geq 0$ times during training; then, the parameter vector should comprise $J+1$ different values for $\boldsymbol{\theta}$ as defined in (9), in addition to the time instances when the changes have occurred. Such formulations in the literature are known as model change detection, see [44].

V Decentralized Generation, Demand and Topology Estimation

V-A Preliminaries and Notation

In the case when the controllers do not not have any knowledge of the steady state bus voltages at remote buses, the system is not observable; hence $\boldsymbol{\theta}_{-n}$ cannot be uniquely identified in classical, non-Bayesian sense (see Appendix B).

Motivated by the ideas in [31], we propose a decentralized solution that splits the slots into two consecutive training phases: (i) measurement phase, denoted as ${M}$-phase, and (ii) communication phase, denoted as ${C}$-phase. The slots in the ${C}$-phase are used to disseminate the local steady state voltage measurements acquired in the ${M}$-phase to remote controllers via amplitude modulation of the reference voltage perturbation signals. Each controller then uses a sequential-type of demodulator to process the local bus voltage measurements acquired in the $C$-phase and acquire full knowledge of the portion of $\mathbf{W}$ that corresponds to the ${M}$-phase. If the training matrices in the ${M}$-phase satisfy predefined conditions, elaborated in subsection V-B, then knowing only the ${M}$-phase portion of $\mathbf{W}$ is sufficient to uniquely estimate the parameter vector locally.

The temporal organization of the proposed training protocol is depicted in Fig. 6, see also Fig. 7. The $C$-phase is further split into sub-phases $\alpha$ (channel estimation sub-phase) and $\beta$ (modulation and demodulation sub-phase). The $M$-phase contains the first $\overline{T}$ slots, indexed in $\overline{\mathcal{T}}=\left\{1,\ldots,\overline{T}\right\}$, the sub-phase $\alpha$ takes the subsequent $T^{\alpha}$ slots indexed in $\mathcal{T}^{\alpha}=\left\{\overline{T}+1,\ldots,\overline{T}+T^{\alpha}\right\}$ and the sub-phase $\beta$ comprises the remaining $T^{\beta}=T-\overline{T}-T^{\alpha}$ slots indexed in $\mathcal{T}^{\beta}=\left\{\overline{T}+T^{\alpha}+1,\ldots,T\right\}$. The sub-phase $\beta$ is further split into $\overline{T}$ blocks, one for each slot in the $M$-phase, see Fig. 6; hence, the blocks are indexed in $\overline{\mathcal{T}}$. Each block is formed by $L$ consecutive time slots, such that $L\overline{T}=T^{\beta}$. We write $\mathcal{T}^{\beta}=\cup_{b\in\overline{\mathcal{T}}}\mathcal{T}^{\beta;b}$ where $\mathcal{T}^{\beta;b}=\left\{\overline{T}+T^{\alpha}+(b-1)L+1,\ldots,\overline{T}+T^{\alpha}+bL\right\},~{}b\in\overline{\mathcal{T}}$ is the set indexing the slots in block $b$. As elaborated in subsection V-C, in block $b$, the controllers disseminate the measurements obtained in slot $b$ in the $M$-phase, see Fig. 7. We introduce notation corresponding to (sub-)phase-wise and block-wise partition of the matrices $\mathbf{W}$, $\mathbf{X}$, $\mathbf{S}$, $\mathbf{V}$ and $\mathbf{\Omega}$. Take the measurement matrix $\mathbf{W}$ as an example (analogous notation applies to $\mathbf{X}$, $\mathbf{S}$, $\mathbf{V}$ and $\mathbf{\Omega}$); it can be partitioned as, see Fig. 7:

[TABLE]

The $\overline{T}\times N$ matrix $\overline{\mathbf{W}}$, with $[\overline{\mathbf{W}}]_{t,n}=w_{n}(t),~{}n\in\mathcal{N},~{}t\in\overline{\mathcal{T}}$, contains the steady state bus voltage measurements from the $M$-phase; $\mathbf{W}^{\alpha}$, $\mathbf{W}^{\beta}$ as well as each of the matrices $\mathbf{W}^{\beta;b},~{}b\in\overline{\mathcal{T}}$ are defined analogously. $\overline{\mathbf{w}}_{n}$ denotes the $n-$th column of $\overline{\mathbf{W}}$; analogous notation applies to the other matrices.

V-B Sufficient Excitation

The purpose of the $C$-phase is to enable each controller to learn $\overline{\mathbf{W}}$, which is sufficient to generate locally a unique estimate of $\boldsymbol{\theta}_{-n}$ for any $n\in\mathcal{N}$, if and only if the Jacobians of $\mathsf{vec}(\overline{\mathbf{\Omega}})$ w.r.t. $\boldsymbol{\theta}_{-n}$ and $\mathsf{vec}(\overline{\mathbf{V}})$, denoted with $\mathbf{\Upsilon}_{-n}$ and $\mathbf{\Gamma}$, respectively, satisfy the rank conditions:

[TABLE]

for any $n\in\mathcal{N}$. The sufficient excitation conditions provide practical guidelines for designing the training matrices $\overline{\mathbf{X}}$ and $\overline{\mathbf{S}}$; this is further discussed in subsection V-F.

We note that the vectorization of $\overline{\mathbf{\Omega}}$ is linear in $\boldsymbol{\theta}$:

[TABLE]

In fact, it can be shown that it is always linear in $\mathbf{d}$ and $\boldsymbol{\psi}$; however, the linearity in $\mathbf{g}$ is a direct corollary of the virtual resistance configuration (6) for proportional power sharing based on the instantaneous generation capacities. This result is useful for finding good initial estimates of $\boldsymbol{\theta}_{-n}$ which will be used to initialize the iterative algorithm.

V-C Training Phases and Sub-phases

In the $M$-phase, the $n-$th controller obtains $\overline{\mathbf{w}}_{n}$, the $n$-th column of $\overline{\mathbf{W}}$. Learning the remaining columns $\overline{\mathbf{w}}_{n},~{}n\neq m$ and obtaining local copy of $\overline{\mathbf{W}}$, denoted with $\overline{\mathbf{W}}_{(n)}$, is done in the $C$-phase where controller $n$ disseminates $\overline{\mathbf{w}}_{n}$ to remote controllers by modulating the amplitudes of the reference voltage deviations and, in the same time, demodulates $\overline{\mathbf{w}}_{m},~{}m\neq n$ from the locally available measurements $\mathbf{w}_{n}^{\alpha/\beta}$ via sequential demodulator, see Fig. 7.

In the $C$-phase, we adopt the following perturbation signals:

[TABLE]

where $\Delta x_{n}(t)\in[-1,+1]$ is the reference voltage perturbation and $\sqrt{\pi_{n}(t)}>0$ is the perturbation amplitude; hence, the droop slopes in the $C$-phase are kept fixed to the nominal value and the communication channel is established via the reference voltage perturbation signals. The $C$-phase training matrices $\mathbf{X}^{\alpha}$ and $\mathbf{X}^{\beta}$ can then be written as follows:

[TABLE]

where $\Delta\mathbf{X}$ and $\mathbf{\Pi}$ are the reference voltage perturbation and perturbation amplitude matrices, defined as $[\Delta\mathbf{X}]_{t,n}=\Delta x_{n}(t)$ and $[\mathbf{\Pi}]_{t,n}=\sqrt{\pi_{n}(t)}$, $n\in\mathcal{N},~{}{t}\in\mathcal{T}^{\alpha/\beta}$, respectively. To facilitate the design of the demodulator, we make the following small signal assumption: the reference voltage deviation amplitudes in the $C$-phase are relatively small w.r.t. the nominal reference voltage, i.e., $\pi_{n}(t)\ll{\tilde{x}}_{n}$, $n\in\mathcal{N}$, $t\in\mathcal{T}^{\alpha/\beta}$. Using Taylor’s series expansion, the signal collected by controller $n$ in the $C$-phase can be written as:

[TABLE]

The model above defines the input-output relation of a real, linear, synchronous communication channel with channel vector given by the gradient $\mathbf{h}_{n}$ (evaluated at the nominal droop values) which contains the real coefficients of the equivalent linear channels that controller $n$ sees to the other controllers; in localized and strongly connected MGs, the entries in $\mathbf{h}_{n}$ do not differ significantly (see also [33]), i.e., the channel (20) experiences strong all-to-all property.

We use the linear model to design sequential transceiver that operates as follows. First, in sub-phase $\alpha$, controller $k$ estimates $\mathbf{h}_{k}$; for this purpose, we fix the perturbation amplitudes to be all known and equal constants:

[TABLE]

Then, in sub-phase $\beta$ the controllers disseminate the information acquired in the $M$-phase via the following linear amplitude modulation (without any additional error protection):

[TABLE]

where $\pi^{\beta}$ and $\chi_{n}$ are known positive constants. Clearly, $\pi_{n}(t)$ remains fixed in block $b\in\overline{\mathcal{T}}$, carrying the information about $\overline{w}_{n}(b)$ by embedding it into the amplitude of the perturbation signal $\Delta\mathbf{x}_{n}^{\beta;b}$. The controllers operate in full duplex transmission mode, simultaneously broadcasting and receiving one voltage measurement per block to/from all other controllers.111The scheme suits well channels with strong all-to-all property, i.e., channels where the gains in $\mathbf{h}_{k}$ do not differ significantly; this is the case for small and localized MGs. As the system grows in size and scope, the all-to-all property ceases to be valid and one should consider applying more sophisticated digital modulation/demodulation and scheduling schemes, including error protection coding; see [32, 33] for alternatives.

To guarantee the uniqueness of the local copies $\overline{\mathbf{W}}_{(n)}$, we restrict the columns of the reference voltage perturbation matrices $\Delta\mathbf{X}^{\alpha}$ and $\Delta\mathbf{X}^{\beta;b}$ to be zero mean and orthogonal:

[TABLE]

where $\delta^{\alpha}=\|\Delta\mathbf{x}_{n}^{\alpha}\|_{2}^{2}\leq T^{\alpha}$, $\delta^{\beta}=\|\Delta\mathbf{x}_{n}^{\beta;b}\|_{2}^{2}\leq L$, for every $n\in\mathcal{N}$, $b\in\overline{\mathcal{T}}$. We note that the above assumptions are a bit restrictive. Given the perturbation signals (21) and (22) in sub-phases $\alpha$ and $\beta$, the sufficient conditions for uniqueness of $\overline{\mathbf{W}}_{(n)}$ for any $n\in\mathcal{N}$ are $\mathsf{rank}(\Delta\mathbf{X}^{\alpha})=\mathsf{rank}(\Delta\mathbf{X}^{\beta;b})=N$ for any $b\in\overline{\mathcal{T}}$; however, we use (23), (24) for convenience, namely, to obtain compact expression for ${\overline{\mathbf{W}}}_{(n)}$ without loosing generality. Replacing (21) and (22) in (20) and using assumptions (23), (24), we derive $\overline{\mathbf{W}}_{(n)}$:

Proposition 2**.**

The local estimators of $\mathsf{vec}(\overline{\mathbf{W}})$ are given by:

[TABLE]

for any $n\in\mathcal{N}$; for notational brevity, we used $\boldsymbol{\mathcal{X}}^{\alpha}=(\Delta\mathbf{X}^{\alpha})^{\mathsf{T}}\otimes\mathbf{1}_{\overline{T}}$, $\boldsymbol{\mathcal{X}}^{\beta;b}=(\Delta\mathbf{X}^{\beta;b})^{\mathsf{T}}\otimes\mathbf{e}_{b}$, $\boldsymbol{\mathcal{I}}=\mathbf{I}_{N}\otimes\mathbf{1}_{\overline{T}}$ and $\boldsymbol{\chi}=[\chi_{1},\ldots,\chi_{N}]^{\mathsf{T}}$.

Proof.

See Appendix C. ∎

By the end of the training epoch, the $n-$th controller has a local copy of the $M$-phase measurement matrix $\mathsf{vec}(\overline{\mathbf{W}}_{(n)})$; if the sufficient excitation conditions (15), (16) hold, then $\mathsf{vec}(\overline{\mathbf{W}}_{(n)})$ is sufficient to estimate $\boldsymbol{\theta}_{-n}$. Formulating an ML estimation problem using $\mathsf{vec}(\overline{\mathbf{W}}_{(n)})$ requires knowledge of the pdf $\rho(\mathsf{vec}(\overline{\mathbf{W}}_{(n)});\boldsymbol{\theta})$; however, obtaining the closed from expression is tedious since (25) involves ratios of non-zero Gaussian random variables. Therefore, we derive Gaussian approximation for the pdf of $\mathsf{vec}(\overline{\mathbf{W}}_{(n)})$ based on first-order perturbation-theoretic approach (see supplementary material). Using Neumann series expansion, we get:

[TABLE]

The covariance matrix $\mathbf{\Sigma}$ can be computed via the first-order approximation (see Appendix D) and is given by:

[TABLE]

The approximation is valid for sub-phase $\alpha$ signals satisfying $\mathbf{0}_{T^{\alpha}}<\mathbf{w}_{n}^{\alpha}<2\mathbf{v}_{n}^{\alpha}$. In practice, this is expected to be satisfied as the probability that $\mathbf{w}_{n}^{\alpha}$ is negative or larger than $2\mathbf{v}_{n}^{\alpha}$ is negligible. In light of this, one can easily verify that the Gaussian approximation converges to the true distribution of $\mathsf{vec}(\overline{\mathbf{W}}_{(n)})$ in the limit $\mathbf{w}_{n}^{\alpha}\rightarrow\mathbf{v}_{n}^{\alpha}$. Expression (27) also captures the effect of the $C$-phase and the transmission schemes we adopted there on the uncertainty in the local copies $\overline{\mathbf{W}}_{(n)}$; specifically, the initial uncertainty in $\overline{\mathbf{W}}$, represented with the first term in (27), increases due to (1) measurement noise in sub-phase $\beta$ (second term) and, (2) the uncertainty in the channel estimates induced in sub-phase $\alpha$ (third term).

V-D Joint System Identification and State Estimation

By the end of the training epoch, the $n-$th controller has $\overline{\mathbf{W}}_{(n)}$ and the $C$-phase measurement vectors $\mathbf{w}_{n}^{\alpha}$ and $\mathbf{w}_{n}^{\beta}$. The reference voltage training matrix $\mathbf{X}^{\alpha}$ is deterministic, so $\mathbf{w}_{n}^{\alpha}$ can still be useful when formulating the estimation problem. On the other hand, the training matrix $\mathbf{X}^{\beta}$ in sub-phase $\beta$ is modulated with $M$-phase measurements; since controller $n$ knows only the noisy copy $\overline{\mathbf{W}}_{(n)}$, it is impossible to reconstruct $\mathbf{X}^{\beta}$ perfectly which makes $\mathbf{w}_{n}^{\beta}$ of no further use. The optimal ML that uses all available information should be defined over an augmented vector, comprising $\mathsf{vec}(\overline{\mathbf{W}}_{(n)})$ and $\mathbf{w}_{n}^{\alpha}$. Including $\mathbf{w}_{n}^{\alpha}$ increases the dimensionality of the problem, but the numerical investigations indicate that it does not yield any practically significant performance gain. We therefore omit $\mathbf{w}_{n}^{\alpha}$ from the ML for clarity of exposition.

The relation between the steady state bus voltages and the parameter vector is defined implicitly in Proposition 1; therefore, we define a joint system identification and state estimation (J-SISE) problem via constrained ML estimation [43, 45]. We introduce the joint parameter/state vector:

[TABLE]

We define $\hat{\boldsymbol{\vartheta}}_{-n},~{}n\in\mathcal{N}$ as the globally optimal solution to:

[TABLE]

formulated w.r.t. the true distribution of $\mathsf{vec}(\overline{\mathbf{W}}_{(n)})$. The problem (29) is neither convex nor concave due to the quadratic nature of the constrains that contain bilinear terms in the decision variables. Since $\mathsf{vec}({\mathbf{\Omega}})$ is sufficiently differentiable in $\boldsymbol{\vartheta}_{-n}$, the constrained optimization problem (29) can be restated as an unconstrained one using the Lagrange method of multipliers [46]. Using the Gaussian approximation (26) and applying the Karush-Kuhn-Tucker (KKT) conditions yields a non-linear system of equations. Using the result of the following proposition, we propose Algorithm 1 based on partially linearized constraints to solve the system iteratively [45]. Specifically, denote $\boldsymbol{\vartheta}^{(j)}$ in the $j$-th iteration and let:

[TABLE]

be the linear approximation of $\mathsf{vec}(\overline{\mathbf{\Omega}})$ around $\boldsymbol{\vartheta}^{(j)}$. The Jacobians ${\mathbf{\Upsilon}}^{(j)}$, ${\mathbf{\Gamma}}^{(j)}$ are evaluated in $\boldsymbol{\vartheta}^{(j)}$. We obtain the following result:

Proposition 3**.**

If the sufficient excitation conditions (15), (16) are satisfied in $\boldsymbol{\vartheta}^{(j)}$, the global solution to (29) after substituting the power balance constraint with (30) is given by:

[TABLE]

Proof.

See Appendix E. ∎

The algorithm starts with an initial guess $\boldsymbol{\vartheta}_{-n}^{(0)}$. Then, we apply Proposition 3 iteratively; the solutions (31), (32) in each iteration serve as an input for the next iteration until convergence. In order to apply Algorithm 1, controller $n$ should know the covariance matrix $\mathbf{\Sigma}$ up to a scaling factor, i.e., knowledge of the noise variance $\sigma^{2}$ is not necessary. To ensure fast convergence, we propose the following initialization: once $\overline{\mathbf{W}}_{(n)}$ is locally available, a reasonable initial estimate of the state $\mathsf{vec}(\overline{\mathbf{V}}^{(0)})$ can be obtained via eq. (13):

[TABLE]

Then, we evaluate ${\mathbf{\Upsilon}}$ in $\mathsf{vec}(\overline{\mathbf{V}}^{(0)})$ and solve (17) for $\boldsymbol{\theta}_{-n}$:

[TABLE]

where ${\boldsymbol{\upsilon}}_{n}^{(0)}$ is the $n$-th column of ${\mathbf{\Upsilon}}^{(0)}$. It can be easily verified that $\boldsymbol{\vartheta}_{-n}^{(0)}$ satisfies the KKT conditions and is a stationary point of the objective in (29). Section VII shows that (34) is unbiased but not efficient estimator of $\boldsymbol{\theta}_{-n}$. In this regard, Algorithm 1 serves to refine the initial estimate $\boldsymbol{\vartheta}_{-n}^{(0)}$ and further reduce its covariance.

V-E Performance

The Mean Squared Error matrix of the unbiased estimator of ${\boldsymbol{\vartheta}}_{-n}$ is defined as:

[TABLE]

$\text{MSE}(\hat{\boldsymbol{\theta}}_{-n})$ and $\text{MSE}(\mathsf{vec}(\hat{\overline{\mathbf{V}}}))$ are defined analogously. In stead of deriving the MSE matrix directly, we use the CRLB inequality to bound it and derive an approximate lower bound using the Gaussian approximation (26). Referring to the optimization problem (29), a straightforward way to bound $\text{MSE}(\hat{\boldsymbol{\vartheta}}_{-n})$ is to use the constrained CRLB [47]. Let $\mathbf{O}$ denote the $\mathsf{dim}(\boldsymbol{\vartheta}_{-n})\times\mathsf{dim}(\boldsymbol{\theta}_{-n})$ matrix whose columns form the orthonormal basis for the null space of the Jacobian $[{\mathbf{\Upsilon}}_{-n},\;{\mathbf{\mathbf{\Gamma}}}]$. Then, $\text{MSE}(\hat{\boldsymbol{\vartheta}}_{-n})$ can be bounded as follows [47]:

[TABLE]

where $\mathbf{0}$ denote all-zero matrices of adequate dimensions. $\mathbf{O}$ is is computed numerically, as it is unavailable in closed form. To bound the MSE matrices of $\hat{\boldsymbol{\theta}}_{-k}$ and $\mathsf{vec}(\hat{\overline{\mathbf{V}}})$ separately, we need to perform numerical block inversion of the right-hand side of (36); the following proposition gives alternative and simpler closed form expressions for these bounds:

Proposition 4**.**

The MSE matrices $\textbf{MSE}(\hat{\boldsymbol{\theta}}_{-n})$ and $\textbf{MSE}(\mathsf{vec}(\hat{\overline{\mathbf{V}}}))$ can be bounded from below as follows:

[TABLE]

where $\boldsymbol{\mathcal{J}}$ denotes the Fisher Information Matrix of $\boldsymbol{\theta}_{-n}$.

Proof.

See Appendix F. ∎

The expressions (37) and (38) can be verified to be asymptotically tight; it can be shown that if Algorithm 1 converges to the global optimum, the MSE matrix is of the same analytical form as (37) and (38), but evaluated at $\hat{\boldsymbol{\vartheta}}_{-n}$. Conversely, expressions (37) and (38) prove the asymptotic efficiency of Algorithm 1.

V-F Discussion

We take a closer look on few crucial aspects that set the applicability boundaries of the proposed method. We consider the sufficient excitation conditions, outlined in subsection V-B; they provide guidelines for designing the training sequences and they determine the overall duration of the training epoch. A straightforward way to guarantee (15), (16) is to ensure that $\overline{T}\geq N^{-1}\mathsf{dim}(\boldsymbol{\theta}_{-n})$ and $\mathsf{rank}(\overline{\mathbf{X}})=N$ and/or $\mathsf{rank}(\overline{\mathbf{S}})=N$. The minimal duration of the $C$-phase is determined by the conditions for uniqueness of ${\overline{\mathbf{W}}}_{(n)}$, such that the total duration of the training epoch (in slots) $T=\overline{T}(1+L)+T^{\alpha}$ in a system with $N$ DERs is lower bounded as:

[TABLE]

The lower bound on $T$ can be attained by random training sequences. Alternatively, when using deterministic codes such as orthogonal Walsh-Hadamard sequences, meeting the rank conditions and, possibly additional conditions such as (23), (24) might require more time slots than $T_{\min}$.

The frequency of the training epoch should match the requirements of the upper layer application. If the application runs periodically, then the training epoch should be invoked in each period, preferably at the beginning, while in event-triggered applications, the training epoch should be invoked whenever the application is triggered. While the frequencies should be equal, the total duration (in seconds) $T\tau$ is expected to constitute only a fraction $0<\gamma<1$ of the average time $\tau^{\text{u.app}}$ between two consecutive application runs. Then, we have the following upper bound on the slot duration:

[TABLE]

where $\tau_{\max}$ is obtained by fixed $\gamma=1$ and $T=T_{\min}$.

Further, since the proposed method is developed in classical estimation framework, each controller requires perfect knowledge of the training matrices $\overline{\mathbf{X}}$, $\overline{\mathbf{S}}$, $\Delta{\mathbf{X}}^{\alpha}$ and $\Delta{\mathbf{X}}^{\beta}$. This means that the training matrices should be designed a priori, delivered to the controllers and kept fixed afterward (via hard-coding for instance). Relaxing this condition requires adequate modifications of the problem formulation, which is out of the paper’s scope. For instance, if no prior knowledge is available, we have no choice but to model the training matrices are deterministic unknowns and modify the definition of the parameter vector to include them.

The method can only identify buses that host at least one DER whose primary controller engages in decentralized training. In other words, buses that host only loads are unidentifiable. However, we can still apply the method in MGs with (potentially many) load buses; in this case, the method identifies the Kron-reduced conductance matrix, which is obtained by isolating the DER buses in the original network and applying block inversion on the original conductance matrix. Analyzing the structure of the Kron-reduced conductance matrix, the controllers might be able to deduce some information on the original conductance matrix, see [48].

VI Decentralized OED via Training

We illustrate the practical potential of the proposed system identification method by applying it in decentralized OED (DOED) as the most common upper layer application in power systems. In OED, each DER $n\in\mathcal{N}$ is assigned a monotonic and convex cost function $c_{n}(p_{n})$ that determines the cost of the output power $p_{n}$ of DER $n$. The aim of the OED is to find the optimal local output powers, referred to as optimal dispatch policies $p_{n}^{*},~{}n\in\mathcal{N}$ that minimize the total cost $\sum_{n\in\mathcal{N}}c_{n}(p_{n})$ such that the total load demand $d^{\star}=\mathbf{1}_{3N}^{\mathsf{T}}\mathbf{d}$ is balanced and the box constraints on the output powers are satisfied:

[TABLE]

where $c(\mathbf{p})=\sum_{n\in\mathcal{N}}c_{n}(p_{n})$, $\mathbf{p}=[p_{1},\ldots,p_{N}]^{\mathsf{T}}$ and $\mathbf{v}=[v_{1},\ldots,v_{N}]^{\mathsf{T}}$. Distributed MGs with small-scale DERs typically use linear cost functions [17]. Hence, we adopt $c_{n}(p_{n})=a_{n}p_{n}$ where $a_{n}$ is the constant marginal cost of the $n-$th DER per unit of injected/stored power. Without loss of generality, the costs are ordered as $a_{n}\leq a_{n+1},~{}n\in\mathcal{N}$, which divides the DERs in several ordered cost groups based on the marginal costs. The optimal solution to (41) is the following decentralized program:

[TABLE]

for any $n\in\mathcal{N}$ (see also [17, 35]). Specifically, the total load demand is first filled with the capacities of the DERs from the cheapest cost groups, until the third condition in (45) is met. Then, the DERs from the cost group that meets this condition share the remaining net load demand proportionally to their local capacities while the DERs from the remaining, most expensive cost groups do not inject power. The DERs that satisfy the first condition in (45) are operated at a constant power (at capacity) and their local controllers are configured in CSC mode (forming the subset $\mathcal{N}^{\text{C}}$), whereas the DERs that satisfy the third condition have flexible power outputs and their local controllers are configured in VSC mode, tuned for proportional power sharing (forming $\mathcal{N}^{\text{V}}$).

Knowing $\boldsymbol{\theta}$, specifically $\mathbf{g}$ and $d^{\star}$, is sufficient for implementing the decentralized program (45). We design a OED protocol in which the controllers utilize decentralized training and Algorithm 1 to acquire the information necessary to execute (45). Fig. 8 illustrates the temporal organization of the protocol. The OED typically runs periodically, every $5-30$ minutes depending on the average rate of change of $\mathbf{g}$ and/or $\mathbf{d}$ [3, 17]. Therefore, we (i) divide the time axis into periodic OED epochs, each of of duration $\tau^{\text{OED}}$, and (ii) assume that $\boldsymbol{\theta}$ changes independently at the beginning of and OED epoch and remains fixed throughout the epoch [17]. In each epoch, the DERs locally run the program (45) using up-to-date information about the generation capacities and load demands. To obtain this information, a fraction of the total duration $\tau^{\text{OED}}$ of the OED epoch is allocated for decentralized training, see Fig. 8. The OED epoch is split into a training epoch of duration $T\tau$ and an optimal operation epoch of total duration $\tau^{\text{OED}}-T\tau$. In the training epoch, the DER controllers perform decentralized training and estimation as described in Sections IV and V. At the end of the training epoch, controller $n$ obtains $\hat{\boldsymbol{\theta}}_{-n}$, used at the beginning of the optimal operation epoch to determine the local dispatch policy $\hat{p}_{n}^{*}$, i.e., to determine which condition in (45) is satisfied locally. Hence, each DER individually decides its primary control configuration via (45) using $\hat{\boldsymbol{\theta}}_{-n}$ and configures the local controller accordingly, forming the subsets $\hat{\mathcal{N}}^{\text{V}}\subset\mathcal{N}$ and $\hat{\mathcal{N}}^{\text{C}}\subset\mathcal{N}$. We use $\hat{\cdot}$ to denote that (45) is solved using $\hat{\boldsymbol{\theta}}_{-n}$.

Implicit in the derivation of the decentralized program (45) is the assumption that the MG is balanced $d^{\star}\leq\sum_{m\in{\mathcal{N}}^{\text{C}}\cup{\mathcal{N}}^{\text{V}}}g_{m}$. However, the stochastic renewable generation might sometimes violate the balance condition. Moreover, due to estimation errors in $\hat{\boldsymbol{\theta}}_{-n}$, the resulting dispatch policies $\hat{{p}}_{n}^{*}$ will in general differ from ${p}_{n}^{*}$, attainable only when $\boldsymbol{\theta}$ is known perfectly; hence, $\hat{\mathcal{N}}^{\text{C/V}}\neq\mathcal{N}^{\text{C/V}}$ in general. This leads to slightly suboptimal MG operation, but it might also violate the balance condition even when $\mathcal{N}^{\text{C/V}}$ satisfy it. This results in loss of voltage regulation as the bus voltage quickly (i) drops towards the lower margin $v_{\min}$ when the net load demand is positive $d^{\star}>\sum_{m\in\hat{\mathcal{N}}^{\text{C}}\cup\hat{\mathcal{N}}^{\text{V}}}g_{m}$ or (ii) rises towards the upper margin $v_{\max}$ when the net-load demand is negative $d^{\star}<\sum_{m\in\hat{\mathcal{N}}^{\text{C}}}g_{m}$. Clearly, additional generation/storage capacity is necessary to balance the remaining demand. We employ a solution based on classical DC bus signaling, where a backup source/storage is activated if the bus voltage crosses certain thresholds [12, 13]. The marginal costs of the backups are denoted with $c_{\text{source}}^{\text{extra}}/c_{\text{storage}}^{\text{extra}}$ per unit generated/stored power; these values are always larger than the largest marginal cost among the DERs in $\mathcal{N}$, i.e., $c_{\text{source}}^{\text{extra}}/c_{\text{storage}}^{\text{extra}}>c_{N}$. In normal operating conditions, the MG is balanced, the backups are not active, and the bus voltage is regulated by the DERs in $\hat{\mathcal{N}}^{\text{V}}$, using the droop control law (6) with parameters:

[TABLE]

dimensioned to maintain the bus voltages in a tight region around the rated voltage $x$, i.e., in the interval $[(1-\xi)x,(1+\xi)x]$ with $\xi$ being a small positive number. If the bus voltage drops below $(1-\xi)x$, it signals power deficit and the backup source is activated and configured in droop-controlled VSC mode, using (6) with parameters set as:

[TABLE]

maintaining the bus voltages in $[v_{\min},(1-\xi)x]$. Conversely, if the voltage rises above $(1+\xi)x$, it signals power surplus and the storage is activated and also configured in droop-controlled VSC mode, using (6) with parameters set as:

[TABLE]

maintaining the bus voltages in $[(1+\xi)x,v_{\max}]$. Fig. 9 summarizes the complete operational dynamics of the proposed system on a single $v-i$ diagram. Note that installing backup generation/storage is standard practice when dimensioning standalone systems [3, 12, 13]. In grid-connected systems, the grid can be used as backup, effectively acting as ideal voltage source with infinite generation/storage capacity [3].

VII Evaluation

VII-A General Simulation Description and Design Parameters

Table I summarizes the numerical values of the simulation parameters that remain fixed in all simulation studies; the values of the remaining parameters are provided in the captions of the respective plots. We consider a line, i.e., cut-ring distribution network topology, where all buses are connected to two other buses except for buses $n=1$ and $n=N$ that are connected to a single bus each. As it is a regular practice for any power system, the MG is dimensioned to operate over a range of load demands. For simplicity, we use $d_{n}^{\text{c}\cdot}\leq d^{\text{c}\cdot}$ for any $n\in\mathcal{N}$ (“$\cdot$” stands for either “a”, “c” or “p”); similarly, $g_{n}\leq g$ for any $n\in\mathcal{N}$, see Table I.

The measurement noise variance $\sigma^{2}$ after averaging $\phi_{S}(\tau-\tau^{\text{transit}})$ samples per slot, see Fig. 5, can be computed as:

[TABLE]

where $\sigma_{S}^{2}$ is the noise variance of the PECs’ ADCs [41].

The number of slots $\overline{T}$ in the $M$-phase for fixed $T^{\alpha}=2N$ and $L=2N$, see Table I, is determined from the total number of slots $T=(1+L)\overline{T}+T^{\alpha}$ which is also fixed:

[TABLE]

The perturbation signals are set as (see also Fig. 10):

[TABLE]

The binary sequences $\Delta x_{n}(t)\in\left\{-1,+1\right\},~{}t\in\overline{\mathcal{T}}$ are formed by tossing a fair coin for any $n\in\mathcal{N}$. This is done a priori, i.e., $N$ binary Bernoulli sequences of length $\overline{T}$ are generated, confirmed to satisfy (15), (16) and stored. The droop slope perturbation laws (52) ensure that the bus voltages will not drop below $x-\Delta v\geq v_{\min}$ or rise above $x+\Delta v\leq v_{\min}$ as long as $\sqrt{\pi}<\Delta v$, see Fig. 10. The reference voltage training sequences in sub-phase $\alpha$ and block $b$ in sub-phase $\beta$ have fixed length of $2N$ slots and are set as:

[TABLE]

Hence, $\delta^{\alpha}=\delta^{\beta}=2$. We also fix $\sqrt{\pi^{\alpha}}=\kappa^{\alpha}\sqrt{{\pi}}$ and $\sqrt{\pi^{\beta}}=\kappa^{\beta}\sqrt{\pi}$, where $0<\kappa^{\alpha},\kappa^{\beta}\leq 1$ are set to keep the reference voltage deviation amplitudes in the $C$-phase relatively small, ensuring that the model (20) is valid for any $\sqrt{\pi}\in(0,\Delta v)$.

The performance of J-SISE w.r.t. the MSE and the performance of DOED w.r.t. the cost, are determined by the configuration of the training epoch, which in turn is determined by variety of factors such as slot duration, number of slots, nominal droop control parameters, training matrices and deviation amplitudes. With all specifications listed above and in Table I, most of these factors are kept fixed in our evaluations and the design parameters of the training epoch are the slot duration $\tau$ and the reference voltage deviation amplitude $\sqrt{\pi}$. Next, we evaluate the performance of the J-SISE in terms of the design parameters and show how to find their optimal values w.r.t. DOED.

VII-B J-SISE Performance

First, we investigate the performance, the scalability and the convergence properties of Algorithm 1 w.r.t. $\boldsymbol{\theta}_{-n}$ from the perspective of controller $n=1$ and compare it against CRLB. We fix the generation capacities of all DERs to have equal values, i.e., $g_{n}=g,~{}n\in\mathcal{N}$ and we do the same with the load components $d_{n}^{\text{ca}}=d_{n}^{\text{ca}},~{}d_{n}^{\text{cc}}=d_{n}^{\text{cc}},~{}d_{n}^{\text{cp}}=d_{n}^{\text{cp}}$ and the line conductances $y_{n,m}=y$ for all $n,m\in\mathcal{N}$. We use the Relative Root Mean Squared Error (RRMSE) metric, derived from the MSE matrix as follows:

[TABLE]

To evaluate the MSE matrix, we use statistical average of individual MSE matrices, obtained for $1000$ different realizations of the noise matrix $\mathbf{Z}$. “$\cdot$” in the above definition stands for either the full vector $\boldsymbol{\theta}_{-n}$ or its constituent vectors, i.e., $\mathbf{g}_{-n}$, $\mathbf{d}$ or $\boldsymbol{\psi}$; in either case, the RRMSE is interpreted as the standard deviation of the estimation error per component of the vector that is used as argument. Note that, when applied to a constituent vector of $\boldsymbol{\theta}_{-n}$, we plug the diagonal block of the MSE matrix corresponding to that particular constituent vector. To compute the corresponding lower bound on the RRMSE, we use the CRLB matrix in (54) instead of the MSE matrix.

We focus particularly on the RRMSE as function of $\sqrt{\pi}$, since RRMSE decreases linearly with ${\tau}$ in the log-domain, see eq. (49)). Fig. 11 depicts the performance of J-SISE for each of the constituent vectors of $\boldsymbol{\theta}_{-n}$, i.e., $\mathbf{g}_{-n}$, $\mathbf{d}$ and $\boldsymbol{\psi}$ against the corresponding lower bounds, for $N=6$ DERs. We have evaluated the lower bounds using both, the constrained CRLB (36) and expression (37) from Proposition 4, and they both yield numerically identical results. Empty markers correspond to the initial estimate that initializes Algorithm 1, obtained via (34), while filled markers correspond to $\hat{\boldsymbol{\theta}}_{-n}$ after Algorithm 1 converges. As expected, J-SISE is efficient and attains the CRLB as $\sqrt{\pi}$ increases, except for values very close to $\Delta v$; here, the RRMSE hits a turning point, after which it increases sharply as a result of the fact that when $\sqrt{\pi}\rightarrow\Delta v$, the droop slope $s_{n}(t)$ grows arbitrarily large and the virtual resistance $y_{n}^{\text{va}}\rightarrow 0$. Hence, the controller starts to behave as an ideal voltage source with infinite capacity, pushing the bus voltages to a fixed value $x-\Delta v$ and making the MG insusceptible to reference voltage perturbations.

We further observe that the generation capacities, Fig. 11LABEL:sub@results1a, and the line conductances, Fig. 11LABEL:sub@results1c, can be identified with very high precision (less than $1\%$ of the true value). In contrast, the RRMSE of the load demands of individual components, Fig. 11LABEL:sub@results1b, is several orders of magnitude higher. We conclude that, identifying the individual components of the loads with satisfactory performance might require excessive (even prohibitive) training epoch durations to suppress the noise. However, in many upper layer applications, detailed knowledge on the individual load component demands is not necessary and knowing only the total bus demand $d_{n}^{\star}=d_{n}^{\text{ca}}+d_{n}^{\text{cc}}+d_{n}^{\text{cp}}$ is sufficient [17, 35]; in such case, an estimate of the total load demand vector $\mathbf{d}^{\star}=[d_{1}^{\star},\ldots,d_{N}^{\star}]^{\mathsf{T}}$, comprising the total demands at each bus, can be obtained from $\hat{\mathbf{d}}$ via $\hat{\mathbf{d}}^{\star}=[\mathbf{I}_{N},\mathbf{I}_{N},\mathbf{I}_{N}]\hat{\mathbf{d}}$. Fig. 11LABEL:sub@results1b shows that $\hat{\mathbf{d}}^{\star}$ can be identified with a precision comparable to the one achieved for the generation capacities and line conductances.

The improvement of $\hat{\boldsymbol{\theta}}_{-n}$ w.r.t. $\boldsymbol{\theta}_{-n}^{(0)}$ given with (34), is also evident, clearly showing that the initial estimate is not efficient. The numerical results (not shown here due to space limitations) show that the average of $\boldsymbol{\theta}_{-n}^{(0)}-\hat{\boldsymbol{\theta}}_{-n}$ converges to zero vector asymptotically. We conclude that the initial estimate $\boldsymbol{\theta}_{-n}^{(0)}$ is indeed unbiased estimator of $\boldsymbol{\theta}_{-n}$ and can be still used in practice even though it is not efficient, particularly, when $\sqrt{\pi}$ is of the same order as/smaller than $\sigma$ or for small $N$. In the first case, Algorithm 1 does not converge, see Fig. 11, and $\boldsymbol{\theta}_{-n}^{(0)}$ remains as the only reasonable choice. The second case can be more clearly observed in Fig. 12 that investigates the performance of the framework for increasing number of buses; we see that for small number of buses (e.g. $N=2$), the RRMSE of the initial estimate approaches the CRLB; in such case, the gain from applying Algorithm 1 is marginal, and $\boldsymbol{\theta}_{-n}^{(0)}$ is sufficient for all practical purposes.

From Fig. 12, we also observe that, the performance of J-SISE tends to deteriorate as the number of buses increases, which is expected due to the increase of $\mathsf{dim}(\boldsymbol{\theta}_{-n})$. A straightforward way to improve the performance of Algorithm 1 and make the estimation error arbitrarily small for large $N$, is to increase $\tau$. However, note that (29) treats the vector $\boldsymbol{\psi}$ as full vector, when in fact it may be sparse, containing many zero entries. This might prove to be problematic as the size of the MG scales, i.e., as the number of buses increases since larger distribution systems are significantly sparser [39], so estimating $\boldsymbol{\psi}$ as if it is full vector might lead to performance degradation [49]. So, an appropriate way to improve the performance when $N$ is large (which is out of the scope of this work) is to modify (29) by adding sparsity constraint on $\boldsymbol{\psi}$ and apply a common relaxation method [49].

Finally, we comment on the convergence speed of Algorithm 1; in all tested cases, that is for $N\leq 12$, Algorithm 1 converges already after $10$ iterations. This remarkable result can be mainly attributed to the fact that the initial estimates $\boldsymbol{\chi}^{(0)}$, $\boldsymbol{\theta}_{-n}^{(0)}$, given with eq. (33), (34), respectively, form a stationary point of the optimization problem (29) (see subsection V-D). The additional fact that they are also (asymptotically) unbiased, implies that $\boldsymbol{\theta}_{-n}^{(0)}$ must lie in a neighborhood around $\hat{\boldsymbol{\theta}}_{-n}$, possibly being an inflection point from which it can easily converge to the global optimum only after several iterations.

VII-C Optimizing the Cost Trade-off in DOED

The results presented in the previous subsection do not consider (i) the effect that the estimation error has on the upper layer applications, and (ii) the effect that the power dissipation during training has on the overall performance of the MG. In other words, improving the performance of J-SISE, which is desirable from the perspective of the upper layer application, comes at the “price” of increased power dissipation during training, either by using large perturbation amplitudes or long slot durations, which in turn compromises the performance of the upper layer control application. This leads to a fundamental trade-off between the performance of J-SISE, which is determined by the configuration of the training epoch, and the performance of the application. Our goal is to (i) show how to characterize this trade-off via utility function that jointly captures the performances of J-SISE and the upper layer application, and (ii) provide guidelines on how to design optimal training epochs, namely, how to choose $\tau$ and $\sqrt{\pi}$ such that the utility function is optimized.

As a case study, we take the DOED protocol, described in subsection VI, noting that the approach described below can be applied to any upper layer application. The performance of specific DOED policy vectors $\mathbf{p}$ is assessed via the cost $c(\mathbf{p})=\mathbf{a}^{\mathsf{T}}\mathbf{p}$. The cost of the optimal policy $\mathbf{p}^{*}$ is $c^{*}=c(\mathbf{p}^{*})+c^{\text{extra}}=\mathbf{a}^{\mathsf{T}}\mathbf{p}^{*}+c^{\text{extra}}$ with $c^{\text{extra}}$ denoting any extra cost entailed by activating backups in case the MG is unbalanced. $c^{*}$ is in fact the minimal cost, attainable only when $\mathbf{g}$ and ${d}^{\star}$ are perfectly known to each controller. However, when running the DOED protocol using the estimated parameter vector, see subsection VI, the cost of the resulting dispatch policy vector $\hat{\mathbf{p}}^{*}$ should also account for (i) the fact that $\hat{\mathbf{p}}^{*}\neq\mathbf{p}^{*}$, i.e., the DOED policy $\hat{\mathbf{p}}^{*}$ is, in general, suboptimal, (ii) the fact that $\hat{\mathbf{p}}^{*}$ is valid only in the optimal operation epoch within the OED epoch, and (iii) the power dissipation incurred in the training epoch. We denote this cost with $\hat{c}^{*}$ and we write:

[TABLE]

where the $T\times N$ matrix $\mathbf{P}$ is defined as $[\mathbf{P}]_{t,n}=p_{n}(t),~{}n\in\mathcal{N},~{}t\in\mathcal{T}$ and $p_{n}(t)$ is the output power of DER $n$ in slot $t$. The first term corresponds to the cost of training, whereas the second gives the actual cost of $\hat{\mathbf{p}}^{*}$. We define the Relative Cost Increase (RCI) $\hat{\mu}$, relative to the optimal cost $c^{*}$:

[TABLE]

The RCI can be interpreted as a measure of the additional monetary charge that the the community served by the MG will be subjected to when operating autonomously using the proposed DOED protocol, without any access to external communication enabler. We observe that $\hat{\mu}$ is a random variable whose pdf is parametrized w.r.t fixed $\boldsymbol{\theta}$. In practice, it is desirable to optimize the performance of the upper layer application over the range of $\boldsymbol{\theta}$, which the MG is foreseen to operate in. Therefore, we choose the average RCI, denoted by ${\mu}$ and computed as an average of $\hat{\mu}$ over $\boldsymbol{\theta}$, to be the utility function for the DOED. The aim is to find the optimal training epoch configuration parameters, namely, $\tau$ and $\sqrt{\pi}$, that minimize the average RCI:

[TABLE]

Computing $\mu$ in closed form is far from trivial; therefore, we resort to Monte-Carlo simulation, run the DOED protocol for $100000$ different values of $\boldsymbol{\theta}$ and use the statistical average of the individual RCIs as an estimate of $\mu$. In each trial, $\boldsymbol{\theta}$ is generated independently from the uniform distribution, i.e., $\mathbf{g}\in\mathsf{Unif}[\mathbf{0}_{N},g\mathbf{1}_{N}]$, $\mathbf{d}^{\text{c}\cdot}\in\mathsf{Unif}[\mathbf{0}_{N},d^{\text{c}\cdot}\mathbf{1}_{N}]$, where $g$, $d^{\text{c}\cdot}$ are given in Table I; note that we keep the line conductances fixed to $y$ as the topology changes very infrequently compared to the generation and the load.

Rewriting $p_{n}(t)=\tilde{p}_{n}+\Delta p_{n}(t)$, where $\tilde{p}_{n}$ is the output power of DER $n$ corresponding to the nominal droop parameters, the time average of the power dissipation $\sum_{t\in\mathcal{T}}\Delta p_{n}(t)\approx 0$. We conclude that with (56) and linear OED cost function, it is difficult to asses the impact of power dissipation during training. Therefore, we introduce a quadratically-modified RCI (QRCI), denoted with $\hat{\eta}$:

[TABLE]

where the $T\times N$ matrix $\mathbf{Q}$ is defined as $[\mathbf{Q}]_{t,n}=(p_{n}(t)-\tilde{p}_{n})^{2},~{}n\in\mathcal{N},~{}t\in\mathcal{T}$ and $0<q\leq\frac{\tau}{\tau^{\text{OED}}c^{*}}=q_{\max}$. In similar way as $\mu$, we define the average QRCI, denoted with $\eta$ and restate the optimization problem (57) with $\eta$ as utility function.

The results are presented in Fig. 13. We observe that within the investigated domain, the average RCI, see Fig. 13a is a convex function of $\tau$ and $\sqrt{\pi}$. Specifically, for fixed $\sqrt{\pi}$, $\mu$ decrease as $\tau$ increases due to the effect of noise suppression, see (49). In this regime, the duration of the training epoch is still very short relative to $\tau^{\text{OED}}$, such that the first term in (55) is negligible and the RCI is dominated by the second term which decreases towards $c^{*}$ as the estimation error is reduced. However, $\mu$ hits a turning point when $\tau$ and, consequently, the duration of the training epoch become long enough such that the first term in (55) starts to dominate over the second; after this, it makes no sense to keep increasing $\tau$ as $\mu$ will also increase. Conversely, for fixed $\tau$, $\mu$ decreases as $\sqrt{\pi}$ increases until it hits the turning point after which it starts to increase quickly; evidently, this is happening when we get very close to $\Delta v$. As discussed in the previous subsection, the performance of J-SISE starts to deteriorate when $\sqrt{\pi}\rightarrow\Delta v$, pushing the second term in (55) away from its lower bound $c^{*}$. Hence, within the domain of interest, the average RCI for an MG, specified in Table I and the caption of Fig. 13a, is minimized when $\sqrt{\pi}\approx 8.8$ volts and $\tau\approx 13$ milliseconds. The minimized average RCI is $\mu^{*}\approx 0.008$; in other words, the average increase of the cost is less than $1\%$ of the optimal cost $c^{*}$. This increase, besides being completely tolerable by the OED [3], it is also comparable to the additional operating cost charges imposed by mobile operators when employing wireless cellular solution not including the cost of installing dedicated communication hardware [3, 17].

Similarly as the average RCI, the average QRCI, see Fig. 13b, is also a convex function of $\tau$ and $\sqrt{\pi}$ within the investigated domain with behaviour governed by the same reasoning we used on the average RCI. However, the minimum this time moves closer to the down-left corner due to the second term in (58). Specifically, $\eta$ is minimized when $\sqrt{\pi}\approx 4$ volts and $\tau\approx 5.1$ milliseconds with average RCI $\mu\approx 0.015$, i.e., still around $1\%$ of $c^{*}$.

VIII Concluding Remarks

We introduced autonomous system identification solution, based on temporary primary control perturbations and iterative ML-based algorithm for DC MGs and without access to an external communication system. The method is implemented in a decentralized manner within the primary droop controllers of the PECs and enables the controllers to learn i) the generation capacities of power sources, ii) the load demands, and iii) distribution network topology using only local bus voltage measurements. The key enabling tool is the decentralized training where the primary controllers inject small, amplitude-modulated training sequences that complete the rank of the estimation problem and enable regaining full system observability. We evaluated the performance of the ML-based algorithm, showing that we can achieve high reliability in DC MGs of small to moderate size ($N\leq 12$). Then, we showcased the potential of the solution in fully decentralized OED where the controllers perform training periodically and reconfigure according to the locally estimated information. Last but not least, we illustrated an elaborate methodology for designing training epochs that optimize the operational cost of an autonomous DC MG.

Although we focused on DC MGs and we used several assumptions that simplified the developments, the same design principles introduced in this paper can be applied to any cyber-physical system with dual-layer control architecture that does not not have access to external communication resources, under broader circumstances. Such investigations are part of our on-going and future work.

Appendix A Proof of Proposition I

The power balance condition in each slot states that:

[TABLE]

Recall that during training all DERs are in droop-controlled VSC mode configured for proportional power sharing. Hence $\zeta_{n}=1$ for any $n\in\mathcal{N}$. In such case (59) can be rewritten as:

[TABLE]

Let $\boldsymbol{\omega}^{t}$ be defined as $[\boldsymbol{\omega}^{t}]_{n}=\omega_{n}(t)$; we get:

[TABLE]

where $\mathbf{v}^{t}$, $\mathbf{x}^{t}$ and $\mathbf{v}^{t}$, defined as $[\mathbf{v}^{t}]_{n}=v_{n}(t)$, $[\mathbf{x}^{t}]_{n}=x_{n}(t)$ and $[\mathbf{s}^{t}]_{n}=s_{n}(t)$, represent the $t$-th rows of $\mathbf{V}$, $\mathbf{X}$ and $\mathbf{S}$, respectively. Stacking $(\boldsymbol{\omega}^{t})^{\mathsf{T}}$ vertically for each $t\in\mathcal{T}$, we get the power balance matrix $\mathbf{\Omega}$:

[TABLE]

yielding the compact form (12) which completes the derivation.

Appendix B $\boldsymbol{\theta}_{-n}$ is not identifiable when the system is not observable

We consider the following situation: controller $n$ knows only $\mathbf{w}_{n}$ and knows $\mathbf{X}$ and $\mathbf{S}$ completely. In other words, the controllers do not exchange any local steady state voltage measurements as in the proposed solution, i.e., the $C$-phase training matrices are completely deterministic and known. Hence, all other columns $\mathbf{w}_{m},m\neq n$ are not observable. Since the power balance equation concerning the observable voltages $\boldsymbol{\omega}_{n}=\mathbf{0}_{T}$ also includes and depends on $\mathbf{v}_{m},m\neq n$ (as a result of the fact that the buses are connected through $\mathbf{Y}$) and if classical, non-Bayesian framework is employed (without exploiting any prior knowledge), $\mathbf{v}_{m},m\neq n$ should be treated as unknown parameters in the same way as the generation capacities, load demands and line conductances. Therefore, the parameter vector $\boldsymbol{\theta}$ should be redefined as:

[TABLE]

with

[TABLE]

The sufficient excitation conditions in this case should be restated in term of $\mathbf{\omega}_{n}$ since only $\mathbf{v}_{n}$ is observable; we get:

[TABLE]

where $\mathbf{\Upsilon}_{-n}$ and $\mathbf{\Gamma}$ are the Jacobians of $\mathbf{\omega}_{n}$ w.r.t. $\boldsymbol{\theta}_{-n}$ and $\mathsf{vec}(\mathbf{V})$, respectively. It becomes immediately evident that $\mathbf{\Upsilon}_{-n}$ is a fat matrix, i.e., $\mathsf{dim}(\mathbf{\Upsilon}_{-n})=T\times\mathsf{dim}(\boldsymbol{\theta}_{-n})$ with column rank at most $T<\mathsf{dim}(\boldsymbol{\theta}_{-n})$; hence, the first sufficient excitation condition is not satisfied and $\boldsymbol{\theta}_{-n}$ cannot be uniquely identified.

Equivalently, one can look at the same problem from the perspective of the constrained ML optimization. Namely, the joint parameter/state vector now is:

[TABLE]

The constrained ML optimization problem should be formulated over $\mathbf{\omega}_{n}$ since only $\mathbf{v}_{n}$ is observable:

[TABLE]

Clearly, the number of linearly independent equality constraints is at most $T<\mathsf{dim}(\boldsymbol{\theta}_{-n})$, yielding an ill-conditioned optimization problem that does not converge to any meaningful solution.

Appendix C Proof of Proposition II

Controller $n$ derives the channel estimator $\hat{\mathbf{h}}_{n}$ using the measurement vector from sub-phase $\alpha$, i.e., $\mathbf{w}_{n}^{\alpha}$. Replacing $\sqrt{\pi_{n}(t)}=\sqrt{\pi^{\alpha}},~{}n\in\mathcal{N},~{}t\in\mathcal{T}^{\alpha}$ in the linear model

[TABLE]

we get:

[TABLE]

Using the above, $\hat{\mathbf{h}}_{n}$ is obtained by solving the linear least squares problem:

[TABLE]

Using $\sqrt{\pi_{n}(t)}=\sqrt{\pi^{\beta}}(\overline{w}_{n}(b)-\chi_{n}),~{}n\in\mathcal{N},~{}t\in\mathcal{T}^{\beta;b},~{}b\in\overline{\mathcal{T}}$, (69) can be rewritten as:

[TABLE]

where we used the commutative property of the product $\mathsf{D}(\overline{\mathbf{w}}^{b}-\boldsymbol{\chi})\mathbf{h}_{n}$. Note that $\overline{\mathbf{w}}^{b}$ is the $b$-th row of the $M$-phase measurement matrix $\overline{\mathbf{W}}$ and contains the data transmitted by the controllers in block $b$. Using the the channel estimate, controller $n$ obtains a local copy of $\overline{\mathbf{w}}^{b}$, denoted with $\overline{\mathbf{w}}_{(n)}^{b}$ by solving the following linear least squares problem:

[TABLE]

Note that $\overline{\mathbf{W}}_{(n)}=\sum_{b\in\overline{\mathcal{T}}}\mathbf{e}_{b}(\overline{\mathbf{w}}_{(n)}^{b})^{\mathsf{T}}$; so we get:

[TABLE]

. Vectorizing the above, we obtain:

[TABLE]

which completes the derivation.

Appendix D Derivation of the Gaussian approximation of $\rho(\mathsf{vec}(\overline{\mathbf{W}});\boldsymbol{\theta})$

Let $\mathbf{w}_{n}^{\alpha}=\mathbf{v}_{n}^{\alpha}+\Delta\mathbf{w}_{n}^{\alpha}$ where $\Delta\mathbf{w}_{n}^{\alpha}\sim\mathsf{N}(\mathbf{0}_{T^{\alpha}},\sigma^{2}\mathbf{I}_{T^{\alpha}})$. Similarly, $\mathbf{w}_{n}^{\beta;b}=\mathbf{v}_{n}^{\beta;b}+\Delta\mathbf{w}_{n}^{\beta;b}$ where $\Delta\mathbf{w}_{n}^{\beta;b}\sim\mathsf{N}(\mathbf{0}_{L},\sigma^{2}\mathbf{I}_{L})$ for $b\in\overline{\mathcal{T}}$. Then, (25) can be written as:

[TABLE]

where $(a)$ follows from the Neumann expansion valid on the subset $\mathbf{0}_{T^{\alpha}}<\mathbf{w}_{n}^{\alpha}<2\mathbf{v}_{n}^{\alpha}$. From (85), we see that $\mathsf{vec}(\overline{\mathbf{W}}_{(n)})$ can be approximated with Gaussian random vector with mean:

[TABLE]

and covariance matrix:

[TABLE]

Appendix E Proof of Proposition III

The Lagrange method of multipliers casts the original constrained ML problem into an unconstrained as follows:

[TABLE]

where we used the Gaussian approximation for the pdf $\rho(\mathsf{vec}(\overline{\mathbf{W}}_{(n)});\boldsymbol{\theta})$. $\boldsymbol{\lambda}$ is $\overline{T}N\times 1$ vector of multipliers. Applying the KKT conditions to (91) after replacing the power balance constraint with its first order approximation, we get the following system of equations:

[TABLE]

The above system is linear in $\boldsymbol{\vartheta}_{-n}$ and can be solved efficiently; the derivation of the solution follows similar steps as in (LABEL:31). Multiplying (92) with $\mathbf{\Gamma}^{(j)}\mathbf{\Sigma}$ yields:

[TABLE]

which is substituted in (94) to yield:

[TABLE]

Solving for $\boldsymbol{\lambda}$ gives:

[TABLE]

Multiplying (97) with $(\mathbf{\Upsilon}_{-n}^{(j)})^{\mathsf{T}}$ on both sides, gives:

[TABLE]

which, after replacing $\mathbf{\Upsilon}^{(j)}\boldsymbol{\theta}=\mathbf{\Upsilon}_{-n}^{(j)}\boldsymbol{\theta}_{-n}+\boldsymbol{\upsilon}_{n}^{(j)}g_{n}$ and solving for $\boldsymbol{\theta}_{-n}$ gives (31). Finally, replacing (97) in (92) and solving for $\mathsf{vec}(\overline{\mathbf{V}})$ produces (32), completing the proof.

Appendix F Proof of Proposition IV

Recall that the implicit function theorem governs the existence of an explicit solution of the system of power balance equations $\omega_{n}=0,~{}n\in\mathcal{N}$ of the following form:

[TABLE]

Hence, again by the implicit function theorem, the solution of the $M$-phase power balance equation $\overline{\mathbf{\Omega}}=\mathbf{0}_{\overline{T}\times N}$ exists and can be written in the following form:

[TABLE]

where the $\overline{T}\times N$ matrix $\overline{\mathbf{F}}$ is defined as $[\overline{\mathbf{F}}]_{b,n}=f_{n}(b)~{}n\in\mathcal{N},~{}b\in\overline{\mathcal{T}}$. If $\overline{\mathbf{F}}$ is available in closed form, the $M$-phase measurement matrix (i.e., its vectorization) can be written explicitly in terms of $\boldsymbol{\theta}$ as:

[TABLE]

Using the above, we derive the CRLB. In particular, the MSE matrix of $\hat{\boldsymbol{\theta}}_{-n}$ can be bounded from below as:

[TABLE]

where $\boldsymbol{\mathcal{J}}(\boldsymbol{\theta}_{-n})$ is the Fisher Information Matrix (FIM) defined as:

[TABLE]

Using the Gaussian approximation for the pdf of $\mathsf{vec}(\overline{\mathbf{W}})$, the FIM can be approximated with the following Grammian:

[TABLE]

Applying the implicit function theorem, we obtain the following expression for the Jacobian $\nabla_{\boldsymbol{\theta}_{-n}}\mathsf{vec}(\overline{\mathbf{F}}(\boldsymbol{\theta}))$:

[TABLE]

Substituting the above in (104) gives expression (37). To bound the MSE matrix of $\mathsf{vec}(\hat{\overline{\mathbf{V}}})$, we use (100), i.e., the fact that $\mathsf{vec}(\overline{\mathbf{V}})$ is a transformed version of $\boldsymbol{\theta}$ and apply the corresponding CRLB formula, i.e.:

[TABLE]

completing the proof.