# Combined Transmission and Distribution State-Estimation for Future Electric
Grids

Amritanshu Pandey, Shimiao Li, and Larry Pileggi All authors are affiliated
with Electrical and Computer Engineering Department at Carnegie Mellon
University.

###### Abstract

Proliferation of grid resources on the distribution network along with the
inability to forecast them accurately will render the existing methodology of
grid operation and control untenable in the future. Instead, a more
distributed yet coordinated approach for grid operation and control will
emerge that models and analyzes the grid with a larger footprint and deeper
hierarchy to unify control of disparate T&D grid resources under a common
framework. Such approach will require AC state-estimation (ACSE) of joint T&D
networks. Today, no practical method for realizing combined T&D ACSE exists.
This paper addresses that gap from circuit-theoretic perspective through
realizing a combined T&D ACSE solution methodology that is fast, convex and
robust against bad-data. To address daunting challenges of problem size
(million+ variables) and data-privacy, the approach is distributed both in
memory and computing resources. To ensure timely convergence, the approach
constructs a distributed circuit model for combined T&D networks and utilizes
node-tearing techniques for efficient parallelism. To demonstrate the efficacy
of the approach, combined T&D ACSE algorithm is run on large test networks
that comprise of multiple T&D feeders. The results reflect the accuracy of the
estimates in terms of root mean-square error and algorithm scalability in
terms of wall-clock time.

###### Index Terms:

circuit-theoretic, combined transmission and distribution networks,
distributed computing, globally convergent, least-absolute value, state-
estimation.

##  I Introduction

Today, the operation of the continental power grids is highly segregated with
control of city-wide distribution grids isolated from control of nation-wide
transmission grids. Such an approach is inherently sub-optimal, resulting in
excessive use of carbon-intensive grid resources and marginalized reliability
during natural disasters and other stressed scenarios. Tomorrow’s grid,
therefore, will rely on an approach that enables joint operation and control
of transmission and distribution (T&D) resources. However, the present state
of the art methods are inadequate for joint operation and control of T&D grids
since they were developed for grid sizes that span up to tens of thousands of
transmission nodes whereas a real-time combined T&D network would be
represented by potentially hundreds of millions of nodes. The latter requires
new models and methods that support re-envisioning the paradigm of grid
operation to enable concurrent operation of both transmission and distribution
grids in the most efficient, reliable, and resilient manner with inclusion of
all grid resources.

Amongst the foremost and most fundamental analyses for joint T&D grid
operation is the combined T&D AC state-estimation (ACSE) problem that would
provide situational awareness for all T&D resources and parameterize the T&D
network models for other analyses. No solution methodology exists today that
can estimate the combined T&D grid states robustly due to both the sheer size
of the problem in terms of number of variables, as well as other practical
requirements. Combined T&D ACSE must accept sensor data from transmission and
distribution grids as input and produce an estimate of grid states in near
real-time. Disparate ACSE methods exist today for transmission [1], [2], [3]
and distribution grids [4] that are not generically extensible to estimation
of combined T&D network states. Furthermore, no standard method exists that
can synchronously access and incorporate measurements from a wide array of
sensors in the T&D grids. This is partly due to the privacy needs of various
T&D stakeholders and partly due to the lack of a physics-based approach that
can incorporate any type of sensor measurement into the grid model. The
difficulty of the problem is further exacerbated by other challenges, such as
inaccurate mapping of distribution grid topology and cyber-security concerns
for devices ranging from critical grid assets to low-privilege access IoT
devices.

To enable an envisioned scenario of joint control and operation of all T&D
resources, in this paper, we propose a new circuit-theoretic combined T&D ACSE
approach. By building it on a strong theoretic foundation, this approach
enables distributed evaluation of the ACSE while supporting any existing or
future T&D grid components (including measurements) without loss of
generality. We begin by distributed circuit-theoretic construction of the ACSE
problem by independently modeling different T&D networks (with measurements
included) as positive-sequence or three-phase equivalent circuits (see Section
III). To couple these independent sub-network together, we develop a novel
supervisory layer (see Section IV) to result in a combined T&D ACSE model.

Building upon the measurement-dependent combined T&D circuit model, Section V
formulates the circuit-theoretic ACSE problem. The formulation has two key
properties: 1) convexity and 2) implicit robustness against bad data. To
formulate a convex ACSE problem, Substitution Theorem [5] and feature mapping
are applied when constructing the circuit models for measurement devices. The
resulting circuit models are linear and impose affine network constraints. To
ensure implicit robustness against bad-data, we build upon our prior work
[6]-[5] to utilize an inherent property of grid bad-data: that it is sparsely
populated. Under the hold of this property, we minimize the weighted least
absolute value (WLAV) of measurement error/noise to obtain a sparse residual
vector that localizes suspicious bad-data.

To make it scale well, as the number of T&D grid states can be upwards of tens
of millions, we develop a distributed ACSE computing strategy (regarding both
memory and computing resources) away from traditional centralized approaches.
To ensure timely solution of the complete problem, we apply theory of
diakoptics (or node tearing) to breakdown the overall problem into sub-
problems and we solve the overall problem without relaxation with Gauss-Seidel
Newton (GSN) method. The efficacy and scalability of this approach is
demonstrated in Section VII), by testing the distributed combined T&D ACSE
approach on large T&D networks. To demonstrate the robustness of WLAV
estimation under presence of bad-data, we compare the results from WLAV
algorithm with a baseline of weighted least square (WLS) objective function.

##  II Prior Work

Combined T&D ACSE is a subject that has not been studied extensively
(exceptions include [7]) even though traditional ACSE for transmission
networks is a workhorse of day-to-day grid operation. The unique facets that
make combined T&D ACSE far more challenging than the traditional transmission
or distribution ACSE are:

  1. 1.

dependency on a network model that is more complex and orders-of-magnitude
larger than the traditional transmission or distribution network models, and;

  2. 2.

a need for a distributed estimation algorithm that is fast and resilient
against bad or compromised measurement data, which are increasing in numbers
due to growing attack surface.

Hence, driven by the needs of these characteristics, the related prior work
for combined T&D ACSE can be categorized into two broad categories: i)
approaches for combined modeling of T&D networks and ii) approaches for
distributed ACSE using these models.

###  II-A Combined T&D Modeling

The proposed approach for combined T&D ACSE is dependent on the underlying
combined T&D models. There is active research on development of such combined
T&D models [8], [9], [10], [11] [12]; however, most of the research focus has
been on building models for solving planning problems rather than state
estimation that requires optimization and inclusion of measurement data. These
models, when combined, represent transmission and distribution grids with a
mix of positive sequence and three-phase network components. Two modeling
approaches stand out for such models: i) co-simulation approach [8], [9], [10]
and ii) co-modeling or integrated modeling approach [11] [12].

In the co-simulation approach, separate tools model the transmission and
distribution components with disparate formulations. Combined simulations are
then run by communicating the coupling information between the tools. To
achieve convergence, separate tools develop independent heuristics with a
myopic view of its model rather than considering the complete T&D context.
Such an approach poses challenges for fast and timely convergence of the
combined simulation as it is dependent on the success of each individual sub-
simulation and the corresponding methodology.

Alternatively, co-modeling or integrated modeling approach enables a common
framework for developing combined T&D models. Literature for such frameworks
is limited due to lack of common formulations for modeling both positive
sequence transmission and three-phase distribution networks ([11], [13], [14]
and [15]). One such approach by the authors overcomes this issue by modeling
both transmission and distribution grid elements as circuit elements [12].
Preliminary results based on these equivalent circuit models have demonstrated
robust convergence for large T&D networks for power flow problems [12]. Here,
we will extend the use of these circuit-theoretic models for combined T&D ACSE
problem with optimization capabilities.

The combined T&D ACSE approach must incorporate measurement data from vast
array of T&D sensors. Although vast literature exists for incorporating
standard metering devices on transmission networks (remote terminal units
(RTUs) and phasor measurement units (PMUs) [3], [16] and [17]), far fewer have
discussed the incorporation of measurements on the distribution grid ([18],
[19] and [4]) especially those consisting of new-age devices such as smart
meters, EV charging etc. ([20] and [21]). In Section IV, we will develop a
generic methodology to include measurements from wide array of sensing devices
from both transmission and distribution grids into the ACSE framework.

###  II-B Distributed State-Estimation Algorithms

Even with the lack of substantive research in combined T&D ACSE, there is
plethora of active research in the field of distributed SE [22], [23], [24],
[25]) and distributed optimization for power grids [26], [27], [28], [29] that
target some facets of the described problem. However, most of this research is
not directly applicable to solving combined T&D ACSE due to its inability to
efficiently scale to large combined T&D networks with both positive sequence
and three-phase components. A successful combined T&D ACSE algorithm must be
fast and robust. It should be significantly faster than the future T&D market
dispatch interval and it should provide certain convergence guarantees. Due to
the implicit non-linearities in the traditional formulations [3], [1], [2],
these conditions cannot be satisfied.

Figure 1: Contour plots for the solution space of weighted-least squares SE
formulation.

We consider the ubiquitous WLSE algorithm with power mismatch based
measurement functions [1], [2] to elaborate further. Due to the non-convex
solution space (see the contour plot in Figure 1), the algorithm has multiple
peaks and valleys that represent the various local optimal solutions for the
ACSE problem (e.g. in the figure. green dot represents a local solution
whereas blue dot represents the global optimum). Depending on the choice of
initial conditions and the corresponding convergence trajectory, the WLS
formulation may converge to a local optimum or a global optimum. In case of
convergence to a local optimum, we may observe high residuals with estimated
states far from the true network states. Even worse it could converge to a
saddle point (shown in yellow) again resulting in high differences between
estimated and true network states.

Some existing literature [30], [31] addresses these underlying challenges. For
instance, many present provably convergent algorithms [32], [31]; albeit
almost all of them are only based on phasor measurement units (PMU) [32], [17]
making them impractical in a realistic T&D grid scenarios. There are methods
that incorporate network topology to converge to a more physical solution, but
they are known to face numerical difficulties [30] with added convergence
challenges. There is also recent works on hybrid approach [16] that include
both RTU and PMU measurements into its framework; however, these are still
formulated as a non-convex optimization problem, thereby facing the same
challenges as those described above. More recently, the authors have proposed
centralized circuit-theoretic WLSE and WLAV approaches for transmission grid
ACSE [5]. These tackle some challenges with traditional transmission ACSE
formulation but does not directly extend to distribution grids or distributed
computing frameworks.

##  III Equivalent Circuit Modeling for Combined T&D Networks

Figure 2: Notional figure of combined T&D network. Figure 3: Notional figure
of combined T&D network with measurement model: measurement components have
been replaced by circuit models so that the coupled T&D network is equivalent
to interconnected sub-circuits.

The combined T&D ACSE framework developed within this paper is based on a
combined T&D model with its data distributed across multiple memory units and
ownership across multiple stakeholders. We construct the model following a
circuit-theoretic approach [12] that represents transmission network
components as positive sequence circuits and distribution network components
as three-phase circuits. To couple the various sub-circuit models
($S\subseteq\textbf{S}$), we develop a supervisory layer (preferably on cloud)
that communicates between various sub-network entities through coupling
circuits to enable a robust distributed optimization framework. We apply
Substitution Theorem to map measurement data into circuit models via replacing
any measured element in the combined T&D grid (see Figure 2) with a mix of
independent voltage and current sources (see Figure 3). To model measurement
noise, we augment the sources with slack injections to satisfy Kirchhoff’s
current laws. The aggregated measurement-dependent circuit includes complete
topology of the T&D network, including the zero-injection nodes which
constitute a majority of the distribution system. This constricts the solution
space of the estimator and drives it to a more physically meaningful solution.

###  III-A Measurement Models for Combined T&D Networks

The electric grid consists of a variety of sensing equipment measuring various
system states at different sampling rates. Most of these measurements are
communicated back to grid T&D control rooms owned by entities such as local
grid operators, utilities and data aggregators. This section describes the
mapping of measurement data into measurement circuit models, which are stored
at the local entity level. Even though variety of sensing devices exist in the
T&D grid (e.g., RTUs, smart meters, IoTs etc.), all measured parameters of
interest from these devices can be categorized into one of the following four
categories without loss of generality:

  * •

Root mean square (RMS) power injection and flow measurements

  * •

RMS voltage nodal measurements

  * •

Phasor voltage nodal measurements

  * •

Phasor current injection and flow measurements

For the transmission network each measurement parameter represents a positive
sequence quantity whereas for distribution network, it represents a three-
phase quantity. Starting from fundamental block of these four categories, we
can construct the measurement circuit model for any grid measuring device. We
describe two examples in detail: (i) transmission grid PMU and (ii)
distribution grid smart meter. Following the same approach, other circuit
models for other measurement devices can be constructed.

####  III-A1 Phasor measurement units on the transmission network

Figure 4: Real part of PMU injection circuit model (left); Real part of a
smart-meter circuit model for wye connected load. (right).

Phasor measurement units (PMUs) are highly accurate measurement devices with
clock synchronization through a common GPS signal. These devices can capture
voltage and current waveforms concurrently with other such PMU devices. This
allows a PMU device (identified by $\rho$ in subscript of its parameters) to
measure voltage and current magnitudes as well as the corresponding angle
differences with respect to a system reference. Injection and flow
measurements constitute two critical measurements by a PMU device. To further
demonstrate mapping of a PMU device to an equivalent circuit model, we make
use of injection measurements as an example.

For injection measurements, a PMU can provide real-time meter readings of
positive sequence $(P)$ voltage and injection current phasors:
$(I_{\rho}^{R,P},I_{\rho}^{I,P},V_{\rho}^{R,P},V_{\rho}^{I,P})$ located at a
set of network nodes $N_{\rho}\subseteq N$. By substitution theorem, we map
these PMU measurements into an equivalent circuit model (shown in Figure 4
(left)). The model is represented by independent current sources with values
corresponding to measured real and imaginary currents. To further include
measurement noise such that the KCL is always satisfied, we add additional
slack current sources $\eta$ to the model such that:

$\displaystyle I_{i}^{R,P}=I_{\rho,i}^{R,P}+{\eta_{\rho,i}^{R,P}},$
$\displaystyle\quad\forall i\in N_{\rho}$ (1) $\displaystyle
I_{i}^{I,P}=I_{\rho,i}^{I,P}+{\eta_{\rho,i}^{I,P}},$
$\displaystyle\quad\forall i\in N_{\rho}$ (2)

To provide robustness against bad data, this paper builds a WLAV method that
minimizes least absolute value of measurement error, i.e., solving a solution
that satisfies $\min||\eta||_{1}$. To convert this L1-norm term into an
equivalent but differentiable form, we split the noise term $\eta$ into 2 non-
negative variables ${\eta^{+},\eta^{-}}$, such that
$\eta=\eta^{+}-\eta^{-},||\eta||_{1}=\eta^{+}+\eta^{-},$ and
$\eta^{+},\eta^{-}\succeq 0$. Thus the model in (1)-(2) is equivalently
represented as:

$\displaystyle
I_{i}^{R,P}=I_{\rho,i}^{R,P}+{\eta_{\rho}^{R,P}}^{+}-{\eta_{\rho}^{R,P}}^{-},$
$\displaystyle\quad\forall i\in N_{\rho}$ (3) $\displaystyle
I_{i}^{I,P}=I_{\rho,i}^{I,P}+{\eta_{\rho}^{I,P}}^{+}-{\eta_{\rho}^{I,P}}^{-},$
$\displaystyle\quad\forall i\in N_{\rho}$ (4)

Models in the rest of this paper are also represented with $\eta^{+},\eta^{-}$
variables to fit in the WLAV-based problem definition.

Even though the measured current values can completely parameterize the states
of the measured node, PMUs also provide redundant measurements in terms of
voltage phasors $(V_{\rho}^{R,P},V_{\rho}^{I,P})$. We include these
measurements into the formulation for improving estimation quality as
additional constraints (see eqns. (14), (14f), (14g)). Following the same
methodology, positive sequence and three-phase models for PMU flow
measurements can also be constructed.

####  III-A2 Smart meter for three-phase wye connected load

Smart meters or AMI infrastructure are becoming the most ubiquitous sensing
devices on the distribution grid. Although they are rarely used for grid
situational awareness today, in the near future, this is likely to change.
These smart meters (identified by $r$ subscript in its parameters) measure
three-phase complex voltage magnitudes along with real and reactive power
quantities; however, unlike PMU or $\mu$PMUs devices they cannot measure
relative angles for the voltage or current quantities. Depending on the
connection type, these meters can measure single, triplex or three-phase
voltages and powers. As an example, we show construction of an equivalent
circuit model for a smart meter that is wired across a wye connected three-
phase load. The meter for each phase $(\delta\in{A,B,C})$ measures nodal
voltage ($|V_{r}^{\delta}|$) and power quantities ($P_{r}^{\delta}$ and
$Q_{r}^{\delta}$) on the network nodes $N_{r}\subseteq N$. However, these
measured parameters are noisy and therefore the true parameters for each phase
$\delta$ are represented by:

$\displaystyle P_{i}^{\delta}=P_{r,i}^{\delta}+\Delta P_{r,i}^{\delta},$
$\displaystyle\quad\forall i\in N_{r}$ (5) $\displaystyle
Q_{i}^{\delta}=Q_{r,i}^{\delta}+\Delta Q_{r,i}^{\delta},$
$\displaystyle\quad\forall i\in N_{r}$ (6)
$\displaystyle|V_{i}^{\delta}|=|V_{r,i}^{\delta}|+\Delta|V_{r,i}^{\delta}|,$
$\displaystyle\quad\forall i\in N_{r}$ (7)

To represent the measured parameters as an equivalent circuit, we perform a
feature transformation:

$\displaystyle
G_{r,i}^{\delta}=\frac{P_{r,i}^{\delta}}{{|V_{r,i}^{\delta}|}^{2}},$
$\displaystyle\quad\forall i\in N_{r}$ (8) $\displaystyle
B_{r,i}^{\delta}=\frac{Q_{r,i}^{\delta}}{{|V_{r,i}^{\delta}|}^{2}},$
$\displaystyle\quad\forall i\in N_{r}$ (9)

and to capture noise elements in the model to satisfy KCL constraints, we add
slack injections ($\eta_{r}^{\delta}$). As in the case of PMU model, we
separate the noise terms into positive and negative components to facilitate
WLAV estimation:

$\displaystyle
I_{r,i}^{R,\delta}=G_{r,i}^{\delta}V_{i}^{R,\delta}-B_{r,i}^{\delta}V_{i}^{I,\delta}+{\eta_{r,i}^{R,\delta}}^{+}-{\eta_{r,i}^{R,\delta}}^{-},$
$\displaystyle\quad\forall i\in N_{r}$ (10) $\displaystyle
I_{r,i}^{I,\delta}=G_{r,i}^{\delta}V_{i}^{I,\delta}+B_{r,i}^{\delta}V_{i}^{R,\delta}+{\eta_{r,i}^{I,\delta}}^{+}-{\eta_{r,i}^{I,\delta}}^{-},$
$\displaystyle\quad\forall i\in N_{r}$ (11)

To construct a measurement device model, (10)-(11) are mapped into an
equivalent circuit which then replaces the wye-connected load that it senses
(see Figure 3). Figure 4 (right) represents the real circuit for the smart
meter measuring three-phase wye connected load.

The described modeling approach for smart meter does a feature transformation
that reduces three measured parameters $(|V_{r,i}|,P_{r,i},Q_{r,i})$ to two
$(G_{r,i},B_{r,i})$. This could result in loss of critical information.
Therefore, we can further minimize the L1-norm of square of voltage magnitude
noise $(|V_{r,i}|^{2}-{V^{R}_{i}}^{2}-{V^{I}_{i}}^{2}|)$ in the objective to
preserve this information; however, doing so nullifies the convexity of the
overall formulation.

####  III-A3 Assumptions on model observability

Global observability of T&D network is a necessary condition for running a
combined T&D ACSE. To satisfy this condition, the developed T&D model must
satisfy the following:
The system matrix $(Y_{S})$ of each sub-network $(S\subseteq\textbf{S})$ that
represents nodal balance equations at all studied nodes must be full-ranked
(i.e. an inverse ${(Y_{S})}^{-1}$ exists).

More generally stated, this implies whether the set of flow and injection
measurements in the sub-network are sufficient to describe all the states in
the grid. Transmission networks not only have sufficient measurements to
describe the state of the network, but they generally have redundant
measurements to further improve the quality of estimates [3]. We include these
redundant measurements into the framework (see (14)) as constraints. However,
ensuring observability for distribution networks remains an active research
area [33], [34]. Here we assume that we have access to measurements at all non
zero-injection nodes in the studied distribution grids. These can include high
accuracy substation-based intelligent electronic devices (IEDs) and $\mu$PMU
devices, as well as smart meters or newer cheap IoT devices. Nonetheless,
there is a potential to explore techniques such as pruning of networks
(topology simplification or reduction), use of forecast as measurements, use
of low quality AMI data and other methods such matrix completion [35] for
improving combined T&D model observability. However, this is not covered
within the scope of this paper alongside discussion of techniques for time-
synchronization of this data.

##  IV Supervisory Layer for Combined T&D ACSE

Section III discussed the construction of measurement-dependent circuit models
for distinct T&D utilities. In this section, we describe how these T&D
disjoint models are coupled to create a unified combined T&D circuit model. To
couple various distributed T&D models, we construct a supervisory layer that
is presumably hosted on the cloud and communicates information about boundary
nodes between various distinct utilities. In [12], we developed a similar
supervisory layer for combined T&D power flow. Here, the approach is extended
to enable online distributed optimization in presence of measurements. To
capture knowledge of boundary nodes, supervisory layer consists of a set of
coupling circuits, which are an aggregation of controlled current and voltage
sources. These circuits couple (or unify) the otherwise partitioned circuit
models.

Figure 5: Primal coupling circuit between T&D network models.

Figure 5 depicts the equivalent circuit for the coupling circuit between
positive sequence T and three-phase D circuits. The set of nodes on left of
the figure represent the nodes for the coupling circuit on the transmission
side for bus $m$, whereas the set of nodes on the right of the figure
represent the nodes on the distribution end for three-phase bus $m$. As shown
in Figure 5, the controlled current sources $(I_{R}^{P}$ and $I_{I}^{P})$ and
controlled voltage sources
$(v_{R}^{A},v_{R}^{B},v_{R}^{C},v_{I}^{A},v_{I}^{B},v_{I}^{C})$ represent the
coupling variables between the transmission and distribution sub-circuits.
Specifically, the controlled current sources $(I_{R}^{P}$ and $I_{I}^{P})$
represent the positive sequence currents to be delivered to the distribution
sub-circuit, whereas controlled voltage sources represent the nodal voltages
for the distribution grid $(v^{ABC})$ that mirrors the positive sequence
transmission voltages $(V^{P})$ at the point of interconnection (POI).
Importantly, this port enables connection of any number of transmission and
distribution networks without loss of generality and together they include the
set of nodes ${N}_{c}^{P}$ and ${N}_{c}^{ABC}$ to represent all nodes for
coupling circuits at transmission and distribution side, respectively. The
theory of symmetrical components is used to derive the expressions for
coupling currents and voltages required to model the port. The positive
sequence model of the power grid assumes balanced operation of the grid, and
therefore, ignores the zero and negative sequence components of voltages and
currents. Hence, by ignoring the zero and negative sequence currents consumed
by the distribution feeder at POI, we first model the current coupling between
T&D (as depicted in Figure 5) through the relationship characterized by
controlled current sources, see (12) where $\alpha$ denotes the value of
$\frac{2\pi}{3}$ radians:

$\begin{bmatrix}I_{R}^{P}\\\
I_{I}^{P}\end{bmatrix}=\frac{1}{3}\begin{bmatrix}1&0&\alpha&0&\alpha^{2}&0\\\
0&1&0&\alpha&0&\alpha^{2}\end{bmatrix}\begin{bmatrix}i_{R}^{A}\\\ i_{I}^{A}\\\
i_{R}^{B}\\\ i_{I}^{B}\\\ i_{R}^{C}\\\ i_{I}^{C}\end{bmatrix}$ (12)

Similarly, the voltage coupling between transmission and distribution (also
depicted in Figure 5) is modeled by formulating the distribution end voltages
as a function of transmission POI voltages, as (13) shows:

$\begin{bmatrix}v_{R}^{A}\\\ v_{I}^{A}\\\ v_{R}^{B}\\\ v_{I}^{B}\\\
v_{R}^{C}\\\ v_{I}^{C}\end{bmatrix}=\begin{bmatrix}1&0\\\ 0&1\\\
\alpha^{2}&0\\\ 0&\alpha^{2}\\\ \alpha&0\\\ 0&\alpha\\\
\end{bmatrix}\begin{bmatrix}V_{R}^{P}\\\ V_{I}^{P}\end{bmatrix}$ (13)

The coupling circuits above describe the relationship between primal states of
T&D boundary nodes. However, for combined T&D ACSE (i.e. distributed
optimization problem), we derive similar relationships for dual optimization
variables as well (See Section V). In this approach, we ensure that minimum
data has to communicated from each utility to the supervisory layer. So far it
only includes states corresponding to the boundary nodes. Additionally, if
zero injection nodes are chosen as boundary nodes, then various utilities do
not have to share any measurement data with the supervisory layer but rather
only the solutions of voltages and currents at the boundary node. This helps
satisfy the privacy needs of the various utilities. Next we describe the
distributed combined T&D ACSE algorithm that actuates upon this unified
circuit model.

##  V Distributed Combined T&D weighted-LAV ACSE Formulation

The general idea behind the proposed combined T&D ACSE approach is to minimize
the norm of slack injections in measurement devices subject to the Kirchhoff’s
Current Laws (KCL) at each node in the unified circuit model. Depending on the
chosen norm, the approach can be described as weighted least square estimation
(WLSE) or weighted least absolute value (WLAV) estimation problem. While WLSE
algorithm is more prominent due to its equivalence to maximum likelihood
estimation (MLE) under presence of Gaussian noise, it is not robust against
bad-data. As many distribution grid measurements are likely to be erroneous, a
robust estimator like the WLAV is preferred, which is statistically equivalent
to MLE when noise is modeled as Laplacian distribution.

The solution approach inputs measurement and network data in the form of a
distributed T&D circuit model. It outputs robust estimates for the T&D grid
states. The network constraints are affine due to the use of circuit-theoretic
modeling approach (see Section III and [6] for details). The resultant
optimization problem is convex with a linear objective, affine equality
constraints and convex inequality constraints. The convexity of the
formulation in addition to the natural weak coupling between T&D networks
allows efficient use of distributed algorithms for solution of the problem.
This helps ensure timely convergence of the combined T&D ACSE algorithm with
provable guarantees of convergence; necessary conditions to be of practical
use. The formulation is mathematically described in (14). Note that inequality
constraints exist in addition to KCL at all nodes due to transformation of
non-differentiable L1 terms of the WLAV objective to mathematically equivalent
continuous twice-differential terms through separation of each noise term into
two variables $(+/-)$. Having this form allows use of the distributed primal-
dual interior point (PDIP) algorithm (14) with the Newton’s method.

$\displaystyle\min_{\eta,V}\quad$ $\displaystyle\sum_{i\in
N_{\rho}}w_{\rho,i}\left({\eta_{\rho,i}^{R,\delta}}^{+}+{\eta_{\rho,i}^{R,\delta}}^{-}+{\eta_{\rho,i}^{I,\delta}}^{+}+{\eta_{\rho,i}^{I,\delta}}^{-}\right)$
$\displaystyle+\sum_{i\in
N_{\rho}}w_{\rho,i}\left({\eta_{\rho,i}^{VR,\delta}}^{+}+{\eta_{\rho,i}^{VR,\delta}}^{-}+{\eta_{\rho,i}^{VI,\delta}}^{+}+{\eta_{\rho,i}^{VI,\delta}}^{-}\right)$
$\displaystyle+\sum_{i\in
N_{r}}w_{r,i}\left({\eta_{r,i}^{R,\delta}}^{+}+{\eta_{r,i}^{R,\delta}}^{-}+{\eta_{r,i}^{I,\delta}}^{+}+{\eta_{r,i}^{I,\delta}}^{-}\right)$
$\displaystyle+\sum_{ij\in
N_{\rho}}w_{\rho,ij}\left({\eta_{\rho,ij}^{R,\delta}}^{+}+{\eta_{\rho,ij}^{R,\delta}}^{-}+{\eta_{\rho,ij}^{I,\delta}}^{+}+{\eta_{\rho,ij}^{I,\delta}}^{-}\right)$
$\displaystyle+\sum_{ij\in
N_{r}}w_{r,ij}\left({\eta_{r,ij}^{R,\delta}}^{+}+{\eta_{r,ij}^{R,\delta}}^{-}+{\eta_{r,ij}^{I,\delta}}^{+}+{\eta_{r,ij}^{I,\delta}}^{-}\right)$
(14a) subject to: Nodal constraints for current injection phasor measurements
$\displaystyle
V^{T}Y_{i}^{R}+I_{\rho,i}^{R,\delta}+{\eta_{\rho,i}^{R,\delta}}^{+}-{\eta_{\rho,i}^{R,\delta}}^{-}=0,$
$\displaystyle\quad\forall i\in N_{\rho}$ (14b) $\displaystyle
V^{T}Y_{i}^{I}+I_{\rho,i}^{I,\delta}+{\eta_{\rho,i}^{I,\delta}}^{+}-{\eta_{\rho,i}^{I,\delta}}^{-}=0,$
$\displaystyle\quad\forall i\in N_{\rho}$ (14c) Nodal constraints for RMS
power injection measurements $\displaystyle
V^{T}Y_{i}^{R}+G_{r,i}^{\delta}V_{i}^{R,\delta}-B_{r,i}^{\delta}V_{i}^{I,\delta}+{\eta_{r,i}^{R,\delta}}^{+}-{\eta_{r,i}^{R,\delta}}^{-}=0,$
$\displaystyle\quad\forall i\in N_{r}$ (14d) $\displaystyle
V^{T}Y_{i}^{I}+G_{r,i}^{\delta}V_{i}^{I,\delta}+B_{r,i}^{\delta}V_{i}^{R,\delta}+{\eta_{r,i}^{I,\delta}}^{+}-{\eta_{r,i}^{I,\delta}}^{-}=0,$
$\displaystyle\quad\forall i\in N_{r}$ (14e) Constraints for phasor voltages
$\displaystyle
V_{i}^{R,\delta}-V_{\rho,i}^{R,\delta}+{\eta_{r,i}^{VR,\delta}}^{+}-{\eta_{r,i}^{VR,\delta}}^{-}=0,$
$\displaystyle\quad\forall i\in N_{\rho}$ (14f) $\displaystyle
V_{i}^{I,\delta}-V_{\rho,i}^{I,\delta}+{\eta_{r,i}^{VI,\delta}}^{+}-{\eta_{r,i}^{VI,\delta}}^{-}=0,$
$\displaystyle\quad\forall i\in N_{\rho}$ (14g) Nodal constraints for zero-
injection (ZI) nodes $\displaystyle V^{T}Y_{i}^{R}=0,$
$\displaystyle\quad\forall i\in N_{ZI}$ (14h) $\displaystyle
V^{T}Y_{i}^{I}=0,$ $\displaystyle\quad\forall i\in N_{ZI}$ (14i) Nodal
constraints for RMS flow measurements $\displaystyle
I_{ij}^{R,\delta}+G_{r,ij}^{\delta}V_{i}^{R,\delta}-B_{r,ij}^{\delta}V_{i}^{I,\delta}+{\eta_{r,ij}^{R,\delta}}^{+}-{\eta_{r,ij}^{R,\delta}}^{-}=0$
(14j) $\displaystyle
I_{ij}^{I,\delta}+G_{r,ij}^{\delta}V_{i}^{I,\delta}+B_{r,ij}^{\delta}V_{i}^{R,\delta}+{\eta_{r,ij}^{I,\delta}}^{+}-{\eta_{r,ij}^{I,fl,\delta}}^{-}=0$
(14k) $\displaystyle ij\in E_{r}\subseteq{\\{(i,j)\quad|\quad i\in N_{r},j\in
N},j\neq i\\}$ Nodal constraints for phasor flow measurements $\displaystyle
I_{\rho,ij}^{R,\delta}-I_{ij}^{R,\delta}+{\eta_{\rho,i}^{R,\delta}}^{+}-{\eta_{\rho,i}^{R,\delta}}^{-}=0$
(14l) $\displaystyle
I_{\rho,ij}^{I,\delta}-I_{ij}^{I,\delta}+{\eta_{\rho,i}^{I,\delta}}^{+}-{\eta_{\rho,i}^{I,\delta}}^{-}=0$
(14m) $\displaystyle ij\in E_{\rho}\subseteq{\\{(i,j)\quad|\quad i\in
N_{\rho},j\in N},j\neq i\\}$ Positivity of all noise terms
$\displaystyle\eta^{+/-}\succeq 0$ (14n)

In (14), $N_{\rho}$ is a subset of network nodes $N$ to which phasor
measurement devices are connected $(N_{\rho}\subseteq N)$ and $N_{r}$ is a
subset of network nodes $N$ to which RMS measurement devices are connected
$({N_{r}\subseteq N})$, $E_{\rho}$ and $E_{r}$ denotes the set of edges
(branches) connected to the phasor and RMS measurement devices. $w$ is the
weight associated to trustworthiness of each measurement and is a function of
device type and tolerance. ${\delta}$ is the set of phases that for three-
phase networks represent phases: $\\{A,B,C\\}$ and for transmission networks
represent the positive sequence component: $\\{P\\}$. To solve the problem, we
implement a distributed PDIP algorithm augmented with circuit-theoretic
limiting methods [15] and [36] to ensure timely convergence. We describe the
algorithm in Algorithm 1. In contrast to general purpose optimization tools
where the entire vectors of primal and dual optimization variables are limited
in the PDIP approach [37], [38], we individually damp the primal and dual
variables based on their physical characteristics. For instance, grid-physics
indicates that a Newton step for the voltage variable of 2 pu exhibits non-
physical behavior and therefore we significantly damp it. Similarly for
ensuring dual feasibility, we apply limiting on the $\mu$ dual variables (see
[39] for details). Following such methodology in the past, we have
demonstrated robust convergence of power flow and AC-optimal power flow
problems [15], [39]. Here we apply these heuristics to the circuit based
distributed WLAV-ACSE algorithm.

Figure 6: Contour plots of the solution space for WLAV-ACSE with and without
bad-data.

The circuit-theoretic distributed WLAV-ACSE formulation in (14) has
significant advantages over the traditional PQV-based WLSE/WLAV formulations.
Traditional PQV formulations are highly non-convex (also see Figure 1) and are
unable to scale to robust distributed frameworks. The non-convexity also makes
it difficult to obtain timely solution for large combined T&D ACSE problems.
In contrast, the developed circuit-theoretic formulation is convex with only
one optimum solution (see Figure 6) and therefore can provide fast convergence
guarantees for the grid-physics constrained optimization problem. In addition,
the proposed WLAV approach is robust against bad-data samples and can converge
to a solution closer to true state even in presence of many bad-data samples
(see right of Figure 6 and [6]).

##  VI Distributed Computing Framework for Combined T&D WLAV ACSE Problem

To solve the problem described in (14), we develop a distributed primal dual
interior point approach that can operate on scattered data spread across
multiple utilities. This distributed approach is mandated by practical
realities of T&D grid problems:

  * •

Stakeholders representing separate T&D utilities may not intend to share all
internal network data with each other.

  * •

Solving an iteration matrix with hundreds of millions of nodes may not be
viable under a single machine single memory computing scheme.

The inherent weak coupling between the T&D sub-networks (i.e., most instances
of T&D coupling assumes radial connection) renders distributed solution
algorithm apt for this problem while satisfying the necessary provisions of
speed, robustness, and resilience. More specifically, it enables domain-based
decomposition ([40], [41] and [42]) of the overall problem through node-
tearing such that the solution of the overall problem is obtained via
concurrent runs of decomposed sub-problems.

To describe the solution methodology, we first symbolically derive the
perturbed first-order optimality conditions (KKT) in (15) for the problem
described in (14). The solution to the perturbed KKT equations in (15), for a
convex problem described in (14), provides the WLAV estimate $X^{*}$ of the
combined T&D ACSE problem as well as the global optimum of the underlying
optimization problem. However, for combined T&D ACSE the set of equations can
be very large and data scattered across many utilities. Therefore, instead of
solving the perturbed KKT condition centrally, we decompose them into many
sub-problems through node-tearing. To do so, we first represent it in the
bordered-block diagonal form as shown in Figure 7.

$\displaystyle\nabla_{\lambda}\mathcal{L}=g\left(X\right)=0$ (15a)
$\displaystyle\nabla_{X}\mathcal{L}=\nabla_{X}f_{obj}(X)+\nabla_{X}^{T}g\left(X\right)\mathit{\lambda}-\mathit{\mu}=0$
(15b) $\displaystyle-\mathit{\mu}\
\odot\left(X-\underline{X}\right)+\epsilon=0$ (15c)
$\displaystyle\mathit{\mu}\geq\mathbf{0}$ (15d) $\displaystyle
X\geq\underline{X}$ (15e)

where, $X\in\\{V^{R,I},\eta^{R,I}\\}$ represent the set of system states,
$f_{obj}$ denotes the objective function, $g(X)$ denotes the set of equality
constraints, $X\geq\underline{X}$ represents the lower bound constraints for
the $\eta^{+/-}$ variables, $\lambda$ is the set of Lagrangian multipliers of
equality constraints, and $\mu$ denotes the Lagrangian multipliers for
inequality constraints.

Figure 7: Combined T&D ACSE solution matrix in BBD structure, where
$H_{xx}=\nabla_{X}^{2}f_{obj},A=\nabla_{X}g(X)$.

To apply the distributed approach to (15), we divide the solution variables
$X$ into $|\textbf{S}|$ sets of
$X^{S}_{int}\in\\{V_{int,S}^{R,I},\eta_{S}^{R,I}\\}$ and
$X^{S}_{ext}\in\\{V_{ext,S}^{R,I}\\}$, where $S\in\textbf{S}$ and S is the set
of all the subsystems. Similarly, the set of dual variables are separated into
distinct sets as well, with $\lambda$ separated into
$\lambda^{S}_{int}\in\\{\lambda_{int,S}^{R,I}\\}$,
$\lambda^{S}_{ext}\in\\{\lambda_{ext,S}^{R,I}\\}$, and $\mu$ are separated
into $|\textbf{S}|$ sets of $\mu_{S}$. The overall problem (15) is decomposed
into $|\textbf{S}|$ sub-problems:

$\displaystyle\forall S\in\textbf{S}:$
$\displaystyle\nabla_{\lambda}^{S}\mathcal{L}=g_{S}\left(X^{S}_{int},(X^{S}_{ext})^{k}\right)=0$
(16a) $\displaystyle\nabla_{X}^{S}\mathcal{L}=\nabla_{X}f_{obj}(X^{S}_{int})+$
$\displaystyle\quad\nabla_{X}^{T}g_{S}\left(X^{S}_{int},(X^{S}_{ext})^{k}\right)(\lambda_{int}^{S},(\lambda^{S}_{ext})^{k})-\mu^{S}_{int}=0$
(16b)
$\displaystyle-\mathit{\mu^{S}_{int}}\odot\left(X^{S}_{int}-\underline{X^{S}_{int}}\right)+\epsilon=0$
(16c) $\displaystyle\mathit{\mu_{int}^{S}}\geq 0$ (16d) $\displaystyle
X\geq\underline{X}$ (16e)

Input: T&D network models, measurements, weights, interconnection nodes

Output: Transmission and distribution grid estimates

1 Read T&D input files

2Domain-based Decomposition: Obtain a set of subsystems $\\{S\\}$ by domain-
based decomposition.

3Initialize epoch $k=0$.

4for  _epoch k=1 to N_ do

5 for  _each sub-system S in parallel_ do

6 if  _k=1_ then

7 then initialize coupling port variables $X_{ext}^{S,k},\lambda_{ext}^{S,k}$

8

9 Solve state estimation sub-problem $S$ in (16) with PDIP algorithm (using
Newton’s method).

10 Transfer solutions of $S$ to the supervisory layer.

11

12 Supervisory Layer: Collect limited information for all sub-systems S and
calculate/update $X_{ext}^{S}$ based on (12, 13, 17) to be used in
$(k+1)^{th}$ epoch.

13 Check convergence and break if converged

Algorithm 1 WLAV Combined T&D ACSE

In the described estimation algorithm (also see Algorithm 1), each individual
entity solves the independent sub-problem $S$ iteratively through Newton’s
method to complete one epoch. After each iteration the external variables
$X_{ext}$ and $\lambda_{ext}^{S}$ are communicated to the supervisory layer
that solves a set of coupling equations (12)-(13) to describe the Gauss step:

$f_{ext}(X^{k},\lambda^{k})=0$ (17)

For the $(k+1)^{th}$ epoch, external variables $X_{ext}$ calculated from the
solution of $k^{th}$ epoch are used. Convergence is achieved when the norm of
change in external variables between two subsequent epochs $(\left\lVert
X_{ext}^{k+1}-X_{ext}^{k}\right\rVert)$ is less than a predefined tolerance
$\epsilon$. The algorithm satisfies a key requirement (18) for the given
domain-based decomposition, which is that the dimension of internal states
within any sub-circuit is much larger than the dimension of external states
within that sub-circuit. This condition is naturally satisfied in T&D
circuits.

$\displaystyle\text{dim}(X^{S}_{int})\gg\text{dim}(X^{S}_{ext})\quad\forall
S\in\\{1,....,|\textbf{S}|\\}$ (18)

Another nice-to-have property for such a framework is for individual entities
to keep measurement data $(z_{s})$ private without having to share it with
external entities. This condition can be satisfied by applying the domain-
based decomposition with zero-injection nodes as the boundary between
different sub-networks.

##  VII Experiments

###  VII-A Problem Setup and Creation of Synthetic Measurements

This section evaluates the performance of the developed combined T&D ACSE
algorithm. To generate synthetic measurements for the ACSE experiments, we
first ran a set of combined T&D power flow simulations following the
methodology in [12]. The output of these represented synthetic true states for
the simulated T&D configuration. We generated synthetic measurements by adding
Gaussian noise to these true states. With access to synthetic measurements,
for transmission networks, we constructed positive sequence measurement models
at all injection nodes following methodology in Section III. Additionally, we
also constructed line flow measurement models from $P_{line}^{P},Q_{line}^{P}$
measurements for some randomly selected flow measurements. For distribution
network, we constructed three-phase measurement models for all load nodes
(from $|V|^{\delta},P_{inj}^{\delta},Q_{inj}^{\delta}$). For the second group
of experiments with bad-data, we created measurement samples with synthetic
bad-data by adding large deviations to a few randomly selected measurements in
each of the sub-networks.

For all experiments, we evaluate the performance of the ACSE algorithm through
use of root mean square error $RMSE$ metric, which is defined as follows:

$RMSE=\sqrt{||x_{est}-x_{true}||_{2}^{2}/\text{dim}(x)}$

where $x_{est}$ is a vector of estimated states in the T&D network and
$x_{true}$ is a vector of true states obtained from the combined T&D power
flow simulation. A lower value of RMSE suggests that the grid estimate is
closer to its true magnitude. To run the distributed ACSE framework, we made
use of Google cloud infrastructure with multiple compute instances. The
details are: #cores: 16, Total RAM: 128 GB and CPUs: Intel Cascade Lake.

###  VII-B Experiment 1: Combined T&D ACSE (WLAV and WLSE) with Gaussian Noise

For the first experiment, we ran the WLAV- and WLSE- ACSE algorithms on two
combined T&D networks. The first combined T&D network consisted of a 10k node
synthetic Texas transmission network (ACTIVSg10k) coupled with 7 realistic
distribution feeders (see taxonomy feeders in [43]). The second network
consisted of a larger transmission network with 25k nodes (ACTIVSg25k) with an
additional meshed urban distribution feeder. The results from the ACSE runs
are documented in Table I. As expected, with no presence of bad-data the
performance of both WLSE and WLAV estimation algorithms is similar, with WLAV
algorithm providing marginally better estimates.

###  VII-C Experiment 2: Combined T&D ACSE (WLAV and WLSE) with Bad-data

In the second experiment, we added bad-data to a few randomly selected
measurements ($<$15 in each sub-network) to the same T&D networks that were
evaluated in Experiment 1. We repeated the experiment in presence of bad-data
for both combined T&D networks and documented the results in Table I. The
results show that the RMSE performance metric for the WLAV algorithm is lower
than the WLSE algorithm. This implies that the WLAV estimates are closer to
the true grid states than the WLSE estimates, hence validating the property of
the WLAV algorithm to be robust against bad-data.

TABLE I: RMSE metric comparison on combined T&D networks

T&D CASE |  | Circuit-theoretic
---
WLAV SE
| Circuit-theoretic
---
WLS SE
| ACTIVSg10k + 7 D-nets1
---
0.00465 | 0.00536
| ACTIVSg10k + 7 D-nets1
---
(with bad-data)
0.00556 | 0.01509
| ACTIVSg25k + 8 D-nets2
---
0.00848 | 0.00893
| ACTIVSg25k + 8 D-nets2
---
(with bad-data)
0.00881 | 0.01335

* D-net: Distribution network.
1 Taxonomy feeders GC-12.47-1 and R1-12.47-3 represent the set of distribution
feeders.
2 In addition to the taxonomy feeders, a larger urban meshed distribution
system Network_Model_Case1 is used.
* Both circuit-theoretic WLAV and WLS SEs are proposed within this paper. The difference is that WLAV-SE minimizes $||\eta||_{1}$ and is converted to (14) with $\eta^{+/-}$ as noise variables, while WLS SE minimizes $||\eta||_{2}^{2}$.

To further validate the efficacy of the combined T&D SE algorithms in presence
of bad-data, we plot the true (obtained from combined T&D power flow) as well
as estimated voltages at the various (7) point of interconnections (POIs)
between T&D sub-networks in Figure 8. For the estimated values, we plot
voltages obtained from both WLSE and WLAV ACSE algorithms. We see that
although both WLAV and WLSE provide good estimates for the POI voltages, WLAV
algorithm provides more accurate estimates as is inferred from the distance
between true and estimated values. This validates that the WLAV algorithm
provides more accurate estimates under presence of bad-data as compared to
WLSE algorithm.

Figure 8: Voltage estimates at point of interconnection for combined T&D model
ACTIVSg10k + 7 D-nets.

###  VII-D Experiment3: Scalability

Figure 9: Combined T&D WLAV ACSE algorithm scalability.

Finally, to demonstrate the scalability of the proposed algorithm, we ran
combined T&D ACSE simulations on a set of networks with solution vector
dimension varying from a few thousands to a million+ and plotted the
simulation time (in sec.) on the vertical axis in Figure 9. We observed that
the speed of the simulation is highly dependent on the size of the largest
sub-network. For instance, the time for the million+ dimension system is found
to be dependent on the run-time of the 70k+ node transmission network.
Nonetheless, even with the prototype implementation of the tool, we
demonstrated that the large T&D networks can be solved within 15 min.

##  VIII Conclusion

This paper developed a novel circuit-theoretic approach for combined T&D ACSE
problem. The approach is distributed both in memory and compute resources to
support the practical needs of the problem with its data stored across
multiple utilities and sheer number of the solution variables. The approach is
able to obtain a truer or more physical solution by capturing the nodal
relationships for all T&D nodes in the network. Also, by formulating the
distributed estimation problem with weighted least-absolute value objective,
the approach can implicitly reject bad-data. In the results, we demonstrate
these features of the approach through runs on a set of large T&D networks
both with and without bad-data.

### Reproducibility

All data for the test cases is available in github repo:
git@gitlab.com:amritanshup/sugar_combinedtd_acse_testcases.git

## References

  * [1] Ali Abur and A.. Exposito  “Power System State Estimation: Theory and Implementation”  In _New York: Marcel Dekker_ , 2004
  * [2] A. Monticelli  “Electric power system state estimation”  In _Proceedings of the IEEE_ , 2000
  * [3] F.. Schweppe and J. Wildes  “Power System Static-State Estimation, Part I: Exact Model”  In _IEEE Transactions on Power Apparatus and Systems_ PAS-89.1, 1970, pp. 120–125
  * [4] Mesut E. Baran and Arthur W. Kelley  “State Estimation for Real-Time Monitoring of Distribution Systems”  In _IEEE Transactions on Power Systems_ , 1994  DOI: 10.1109/59.336098
  * [5] S. Li, A. Pandey, S. Kar and L. Pileggi  “A Circuit-Theoretic Approach to State Estimation”  In _2020 IEEE PES Innovative Smart Grid Technologies_ , 2020, pp. 1–5
  * [6] Shimiao Li, Amritanshu Pandey and Larry Pileggi  “A WLAV-based robust hybrid state estimation using circuit-theoretic approach”  In _2021 IEEE PES General Meeting (To Appear)_ , 2021  URL: https://arxiv.org/pdf/2011.06021.pdf
  * [7] H.B. Sun and B.M. Zhang  “Global state estimation for whole transmission and distribution networks”  In _Electric Power Systems Research_ 74.2, 2005, pp. 187–195
  * [8] Karanjit Kalsi et al.  “Integrated Transmission and Distribution Control”, 2013  DOI: 10.2172/1117088
  * [9] K. Anderson, J. Du, A. Narayan and A. Gamal  “GridSpice: A distributed simulation platform for the smart grid”  In _2013 Workshop Modeling Simulation Cyber-Phys. Energy Syst. Proc._ , 2013, pp. 1–5
  * [10] Bryan Palmintier et al.  “Final Technical Report: Integrated Distribution-Transmission Analysis for Very High Penetration Solar PV”, 2016  DOI: 10.2172/1239539
  * [11] H. Sun et al.  “Master–Slave-Splitting Based Distributed Global Power Flow Method for Integrated Transmission and Distribution Analysis”  In _IEEE Transactions on Smart Grid_ 6.3, 2015, pp. 1484–1492
  * [12] A. Pandey and L. Pileggi  “Steady-State Simulation for Combined Transmission and Distribution Systems”  In _IEEE Transactions on Smart Grid_ 11.2, 2020, pp. 1124–1135
  * [13] G.. Chang, S.. Chu and H.. Wang  “An Improved Backward/Forward Sweep Load Flow Algorithm for Radial Distribution Systems”  In _IEEE Transactions on Power Systems_ 22.2, 2007, pp. 882–884
  * [14] P… Garcia et al.  “Three-phase power flow calculations using the current injection method”  In _IEEE Transactions on Power Systems_ 15.2, 2000, pp. 508–514
  * [15] A. Pandey et al.  “Robust Power Flow and Three-Phase Power Flow Analyses”  In _IEEE Transactions on Power Systems_ 34.1, 2019, pp. 616–626
  * [16] S. Zhang, X. Luo and Y. Li  “Performance Evaluation of ISO-NE’s Hybrid State Estimation”  In _2019 IEEE Power Energy Society General Meeting (PESGM)_ , 2019, pp. 1–5
  * [17] A.. Phadke, J.. Thorp and K.. Karimi  “State Estimation with Phasor Measurements”  In _IEEE Transactions on Power Systems_ 1.1, 1986, pp. 233–238  DOI: 10.1109/TPWRS.1986.4334878
  * [18] Anggoro Primadianto and Chan Nan Lu  “A Review on Distribution System State Estimation”  In _IEEE Transactions on Power Systems_ , 2017
  * [19] Mesut E. Baran and Arthur W. Kelley  “A branch-current-based state estimation method for distribution systems”  In _IEEE Transactions on Power Systems_ , 1995
  * [20] A. Alimardani et al.  “Distribution System State Estimation Based on Nonsynchronized Smart Meters”  In _IEEE Transactions on Smart Grid_ 6.6, 2015, pp. 2919–2928
  * [21] Jouni Peppanen et al.  “Leveraging AMI Data for Distribution System Model Calibration and Situational Awareness”  In _IEEE Transactions on Smart Grid_ , 2015
  * [22] L. Xie, D. Choi, S. Kar and H.. Poor  “Fully Distributed State Estimation for Wide-Area Monitoring Systems”  In _IEEE Transactions on Smart Grid_ 3.3, 2012, pp. 1154–1169
  * [23] A. Minot, Y.. Lu and N. Li  “A Distributed Gauss-Newton Method for Power System State Estimation”  In _IEEE Transactions on Power Systems_ 31.5, 2016, pp. 3804–3815
  * [24] Ye Guo, Wenchuan Wu, Boming Zhang and Hongbin Sun  “A distributed state estimation method for power systems incorporating linear and nonlinear models”  In _International Journal of Electrical Power & Energy Systems_ 64, 2015, pp. 608–616  DOI: https://doi.org/10.1016/j.ijepes.2014.07.053
  * [25] R. Ebrahimian and R. Baldick  “State estimation distributed processing [for power systems]”  In _IEEE Transactions on Power Systems_ 15.4, 2000, pp. 1240–1246  DOI: 10.1109/59.898096
  * [26] D.. Molzahn et al.  “A Survey of Distributed Optimization and Control Algorithms for Electric Power Systems”  In _IEEE Transactions on Smart Grid_ 8.6, 2017, pp. 2941–2962
  * [27] Emiliano Dall’Anese, Hao Zhu and Georgios B. Giannakis  “Distributed optimal power flow for smart microgrids”  In _IEEE Transactions on Smart Grid_ , 2013  DOI: 10.1109/TSG.2013.2248175
  * [28] B. Zhang, A… Lam, A.. Domínguez-García and D. Tse  “An Optimal and Distributed Method for Voltage Regulation in Power Distribution Systems”  In _IEEE Transactions on Power Systems_ 30.4, 2015, pp. 1714–1726
  * [29] A. Kargarian et al.  “Toward Distributed/Decentralized DC Optimal Power Flow Implementation in Future Electric Power Systems”  In _IEEE Transactions on Smart Grid_ 9.4, 2018, pp. 2574–2594  DOI: 10.1109/TSG.2016.2614904
  * [30] A. Gjelsvik, S. Aam and L. Holten  “Hachtel’s Augmented Matrix Method - A Rapid Method Improving Numerical Stability in Power System Static State Estimation”  In _IEEE Transactions on Power Apparatus and Systems_ PAS-104.11, 1985, pp. 2987–2993  DOI: 10.1109/TPAS.1985.318939
  * [31] D.. Haughton and G.. Heydt  “A Linear State Estimation Formulation for Smart Distribution Systems”  In _IEEE Transactions on Power Systems_ 28.2, 2013, pp. 1187–1195
  * [32] Murat Göl and Ali Abur  “LAV based robust state estimation for systems measured by PMUs”  In _IEEE Transactions on Smart Grid_ , 2014
  * [33] Siddharth Bhela, Vassilis Kekatos and Sriharsha Veeramachaneni  “Enhancing observability in distribution grids using smart meter data”  In _IEEE Transactions on Smart Grid_ , 2018  arXiv:1612.06669
  * [34] Seiun Park, Deepjyoti Deka and Michael Chcrtkov  “Exact topology and parameter estimation in distribution grids with minimal observability”  In _20th Power Systems Computation Conference, PSCC 2018_ , 2018  arXiv:1710.10727
  * [35] P.. Donti et al.  “Matrix Completion for Low-Observability Voltage Estimation”  In _IEEE Transactions on Smart Grid_ 11.3, 2020, pp. 2520–2530
  * [36] A. Pandey et al.  “Robust Convergence of Power Flow Using TX Stepping Method with Equivalent Circuit Formulation”  In _2018 Power Systems Computation Conference (PSCC)_ , 2018, pp. 1–7  DOI: 10.23919/PSCC.2018.8442680
  * [37] A. Wachter and L.. Biegler  “On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming”  In _Mathematical Programming_ 106.1, 2006, pp. 25–57
  * [38] Thomas F. Coleman and Yuying Li  “An interior trust region approach for nonlinear minimization subject to bounds”  In _SIAM Journal on Optimization_ , 1996  DOI: 10.1137/0806023
  * [39] M. Jereminov, A. Pandey and L. Pileggi  “Equivalent circuit formulation for solving ac optimal power flow”  In _IEEE Transactions on Power Systems_ 34.3, 2019
  * [40] Gabriel Kron  “A Set of Principles to Interconnect the Solutions of Physical Systems”  In _Journal of Applied Physics_ 24.8, 1953, pp. 965–980  DOI: 10.1063/1.1721447
  * [41] Arthur Richard Newton and Alberto L. Sangiovanni-Vincentelli  “Relaxation-Based Electrical Simulation”  In _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , 1984
  * [42] M. Vlach  “LU decomposition and forward-backward substitution of recursive bordered block diagonal matrices”  In _IEE Proceedings G: Electronics Circuits and Systems_ , 1985
  * [43] Kevin P Schneider et al.  “Modern Grid Initiative Distribution Taxonomy Final Report”, 2008  DOI: 10.2172/1040684


# Joint Spectrum Partitioning and Power Allocation for Energy Efficient Semi-
Integrated Sensing and Communications

Ammar Mohamed Abouelmaati, Sylvester Aboagye, Member, IEEE, and Hina Tabassum,
Senior Member, IEEE This work was supported by NSERC Discovery Grant funded by
the Natural Sciences and Engineering Research Council of Canada (NSERC). A. M.
Abouelmaati and H. Tabassum are with the Department of Electrical Engineering
and Computer Science at York University, Toronto, Canada; emails:
<EMAIL_ADDRESS>hinat@yorku.ca).S. Aboagye was with the Department of
Electrical Engineering and Computer Science at York University, Toronto,
Canada and is currently with the School of Engineering, University of Guelph,
Guelph, Canada. e-mail: (saboagye@uoguelph.ca).

###### Abstract

With spectrum resources becoming congested and the emergence of sensing-
enabled wireless applications, conventional resource allocation methods need a
revamp to support communications-only, sensing-only, and integrated sensing
and communication (ISaC) services together. In this letter, we propose two
joint spectrum partitioning (SP) and power allocation (PA) schemes to maximize
the aggregate sensing and communication performance as well as corresponding
energy efficiency (EE) of a semi-ISaC system that supports all three services
in a unified manner. The proposed framework captures the priority of the
distinct services, impact of target clutters, power budget and bandwidth
constraints, and sensing and communication quality-of-service (QoS)
requirements. We reveal that the former problem is jointly convex and the
latter is a non-convex problem that can be solved optimally by exploiting
fractional and parametric programming techniques. Numerical results verify the
effectiveness of proposed schemes and extract novel insights related to the
impact of the priority and QoS requirements of distinct services on the
performance of semi-ISaC networks.

###### Index Terms:

Semi-integrated sensing and communication, resource allocation, mutual
information, energy efficiency.

##  I Introduction

Integrated Sensing and Communication (ISaC) enables simultaneous communication
and sensing on network resources. However, using the entire spectrum for ISaC
may be impractical due to existing allocations for various applications (e.g.,
L-band for air traffic control, S-band for weather observation). Also, not all
users constantly require ISaC services. Hence, there is a need for a resource
allocation framework supporting communication-only, sensing-only, and ISaC
services referred to as semi-ISaC network.

Recently, there have been few studies that focused on resource allocation in
ISaC systems. In [1], an alternating optimization-based power and spectrum
partitioning (SP) scheme was developed to optimize probability of detection
(PD), Cramer-Rao bound (CRB), and the posterior CRB under constraints on the
communication data rate. A power allocation (PA) scheme to maximize PD while
satisfying minimum data rate was proposed in [2]. The authors in [3]
considered communication-centric energy efficiency (EE) maximization of an
ISaC system by optimizing beamforming vector under CRB constraint. In [4], the
authors examined SP in an ISaC network to maximize the sum rate. However, no
sensing performance metric was considered. Also, [5] focused on PA to maximize
signal-to-interference-plus-noise ratio (SINR) for a radar and communication
coexistence scenario and [6] optimized beamforming vectors to enhance PD with
constraints on the data rate requirements.

Unlike an ISAC system, the semi-ISAC system offers additional degree of
freedom to network operators to dynamically prioritize users based on their
required services (i.e., communication only, sensing only, or ISAC), minimum
quality-of-service (QoS) requirements, and available network resources.
Recently, in [7], outage probability analysis was performed for a single-
input-single-output (SISO) semi-ISaC system and it was shown to outperform an
ISaC system in terms of outage probability, ergodic rate, and radar estimation
information rate.

To date, none of the existing studies explored the joint optimization of SP
and PA for ISaC systems, in general, and semi-ISaC systems, in particular.

In this letter, we develop a unified sensing and communication optimization
framework that jointly optimizes PA and SP for semi-ISAC networks under base-
station (BS) power budget, total bandwidth, and users’ QoS requirements while
considering interference from clutters and users’ priorities. We consider two
objectives of practical relevance, i.e., maximizing (i) the aggregate radar
mutual information (MI) and data rate, and (ii) EE. We establish the joint
convexity of the former problem in spectrum and power allocation variables.
The latter problem is a non-convex problem, thus we transform it into a
concave-convex fractional problem and propose a Dinkelbach-based algorithm
which converges to the global optimal solution. Simulation results affirm the
optimized semi-ISaC framework’s efficacy over existing benchmarks. Also, we
illustrate how varying user requirements and priorities impact resource
allocation choices.

Figure 1: Considered semi-ISaC system with sensing-only, ISaC, and
communications-only services.

##  II System Model and Assumptions

We consider the downlink (DL) of an orthogonal multiple access-based semi-ISaC
system as illustrated in Fig. 1. The BS, operating in full-duplex mode,
partitions its available spectrum into three portions for three distinct
services, namely, $\tau_{1}$ for sensing-only where BS senses the surrounding
environment between it and user 1, $\tau_{2}$ for ISaC where BS senses the
surrounding environment between it and user 2 while communicating with that
user, and $\tau_{3}$ for communications-only with user 3. Note that the range
resolution of the ISaC system is considered sufficiently accurate to avoid the
interference between two targets (e.g., radar target and ISaC user) [7].
Typical sensing applications include high-accuracy localization and tracking,
augmented human sensing, simultaneous imaging and mapping, and space
surveillance [8, 9].

We define the proportion bandwidths $\tau_{s}$ such that the total bandwidth
can be expressed as $W={\tau_{1}W}+{\tau_{2}W}+{\tau_{3}W},$ where
$\tau_{s}\in[0,1]$ for $s=1,2,3,$ and $\sum_{s=1}^{3}\tau_{s}=1$, as depicted
in Fig. 1. Without loss of generality, we assume $J$ clutter scatterers
observed at the radar receiver due to the illumination of the radar
transmitter. Similar to [10, 7, 6], we assume that the number of clutters and
channel state information can be estimated perfectly.

The path-loss model (PLM) on the communication channels is given as
$\mathcal{L}_{c}(d_{s})=\frac{G_{\mathrm{tx}}d_{s}^{-\alpha_{c}}c^{2}}{(4\pi
f_{c})^{2}}$[11], where $G_{\mathrm{tx}}$ is the transmit antenna gain,
$s\in\\{2,3\\}$ is the index of the ISaC user and the communications-only
user, $d_{s}$ denotes the link distance, $\alpha_{c}$ is the path-loss
exponent (PLE), $f_{c}$ is the carrier frequency, and $c=3\times 10^{8}$ m/s
is the speed of light. On the other hand, the two-way PLM for the radar echoes
is given as
$\mathcal{L}_{r}(d_{1})=\frac{G_{\mathrm{tx}}d_{1}^{-2\alpha_{r}}\sigma_{\rm
RCS}\lambda^{2}}{(4\pi)^{3}}$, where $\alpha_{r}$ is the PLE, $\lambda$ is the
wavelength of the carrier, and $\sigma_{\rm RCS}$ is the target radar cross
section (RCS). By replacing $d_{1}$ with $d_{j}$, the PLM can be customized
for the clutter. In what follows, we describe the received signal models for
sensing, ISaC, and communication scenarios.

####  II-1 Sensing-Only

The BS sends a radar signal to the target. The radar echo signal at the BS in
the presence of $J$ clutter components can be expressed as [5, 7]:

$y_{1}(t)={h_{1,d}h_{1,u}\sqrt{P_{1}\mathcal{L}_{r}(d_{1})}z_{1}(t-\widetilde{\tau_{1}})}+I_{\mathrm{c}_{1}}+n_{1}(t),$
(1)

where
$I_{\mathrm{c}_{1}}={\sum_{j=1}^{J}h_{j,d}h_{j,u}\sqrt{P_{1}\mathcal{L}_{r}(d_{j})}z_{1}(t-\widetilde{\tau_{j}})}$
is the clutter interference, $P_{1}$ is the allocated transmit power,
$z_{1}(t)$ is the radar signal, $n_{1}(t)$ is the additive White Gaussian
noise (AWGN) with variance $\sigma_{1}^{2}=k_{B}T_{\rm temp}\tau_{1}W$, where
$k_{B}$ is the Boltzmann constant, $T_{\rm temp}$ is the absolute temperature,
$\widetilde{\tau_{1}}$ and $\widetilde{\tau_{j}}$ are the time delay of the
echo coming from user 1 and the $j$-th clutter, respectively. It is assumed
that the time delays can be perfectly compensated through synchronization [7].
The DL and uplink (UL) small-scale fading channel power gains between the
radar target and BS are given by $|h_{1,d}|^{2}$ and $|h_{1,u}|^{2}$,
respectively. The corresponding signal-to-clutter-plus-noise ratio (SCNR) is
thus given by:

$\gamma_{1}=\frac{{P_{1}\mathcal{L}_{r}(d_{1})\left|g_{1}\right|^{2}}}{{\sum_{j=1}^{J}P_{1}\mathcal{L}_{r}(d_{j})\left|\zeta_{j}\right|^{2}}+{\sigma_{1}^{2}}},\\\
$ (2)

where $|g_{1}|^{2}=|h_{1,d}|^{2}|h_{1,u}|^{2}$ is the cascaded channel power
between the desired target and the BS, and
$|\zeta_{j}|^{2}=|h_{j,d}|^{2}|h_{j,u}|^{2}$ is the cascaded channel power of
the $j$-th interfering clutter.

####  II-2 ISaC

The framework for perceptive mobile networks–a joint communication and sensing
enabled mobile network–is adopted here [9], where a single transmitted signal
is used for communication and sensing, simultaneously. Specifically, the ISaC
user utilizes the received signal for communication while downlink active
sensing–where the BS uses reflected communication signals from its own
transmitted signal for sensing–happens at the BS. Full-duplex technology
allows the BS to perform sensing and communication over the same time-
frequency resources [12]. Similar to [11, 13], we assume perfect self-
interference cancellation. The received signal at the ISaC user is given by:

$\displaystyle
y_{2}^{d}(t)={h_{2,d}\sqrt{P_{2}\mathcal{L}_{c}(d_{2})}z_{2}(t)}+n_{2}(t),$
(3)

and the corresponding SNR is expressed as:

$\gamma_{2}^{d}=\frac{{P_{2}\mathcal{L}_{c}(d_{2})|h_{2,d}|^{2}}}{{\sigma_{2}^{2}}},$
(4)

where $P_{2}$ is the transmit power of the ISaC signal $z_{2}{(t)}$,
$n_{2}(t)$ is AWGN with variance $\sigma_{2}^{2}=k_{B}T_{\rm temp}\tau_{2}W$,
and $|h_{2,d}|^{2}$ is the channel power. The reflected echo at the BS in the
presence of $J$ clutter components can be expressed as:

$y_{2}^{u}(t)={h_{2,d}h_{2,u}\sqrt{P_{2}\mathcal{L}_{r}(d_{2})}z_{2}(t-\widetilde{\tau_{2}})}+I_{\mathrm{c}_{2}}+n_{2}(t),$
(5)

where
$I_{\mathrm{c}_{2}}={\sum_{j=1}^{J}h_{j,d}h_{j,u}\sqrt{P_{2}\mathcal{L}_{r}(d_{j})}z_{2}(t-\widetilde{\tau_{j}})}$,
and the corresponding SCNR is given by:

$\gamma_{2}^{u}=\frac{{P_{2}\mathcal{L}_{r}(d_{2})|g_{2}|^{2}}}{\displaystyle\sum_{j=1}^{J}{P_{2}\mathcal{L}_{r}(d_{j})\left|\zeta_{j}\right|^{2}}+{\sigma_{2}^{2}}},$
(6)

where $|g_{2}|^{2}=|h_{2,d}|^{2}|h_{2,u}|^{2}$ is the cascaded channel power
gain.

####  II-3 Communication Signal

The received signal at the communications-only user is given by:

$y_{3}(t)={h_{3}\sqrt{P_{3}\mathcal{L}_{c}(d_{3})}z_{3}(t)}+n_{3}(t),$ (7)

and the corresponding SNR is
$\gamma_{3}=\frac{P_{3}\mathcal{L}_{c}(d_{3})\left|h_{3}\right|^{2}}{\sigma_{3}^{2}},$
where $z_{3}{(t)}$ is the transmit signal, $|h_{3}|^{2}$ is the small-scale
fading channel power gain, $n_{3}(t)$ is AWGN with variance
$\sigma_{3}^{2}=k_{B}T_{\rm temp}\tau_{3}W$, and $P_{3}$ is its allocated
transmit power.

##  III MI, Data Rate, and EE Optimization

###  III-A MI and Data Rate Optimization

For the radar sensing and ISaC transmission, we use MI as a metric that
measures how much information about the channel is conveyed to the radar
receiver [14, 15]. Specifically, after receiving $y(t)$, the priori
uncertainty of the target/ISaC user decreases since there is some information
about the target/ISaC user contained in the cascaded channel gain $g$. The MI
between $y(t)$ and $g$ is given as $I(y;g)=H(y)-H(y|g)=H(y)-H(n_{o})$ [14,
15], where $n_{o}$ is the noise and interference at the receiver. Assuming
both $H(y)$ and $H(n_{o})$ are Gaussian, their differential entropy is given
as $H(y)=\frac{1}{2}\log_{2}\left(2\pi eE\\{|y|^{2}\\}\right)$ and
$H(n_{o})=\frac{1}{2}\log_{2}\left(2\pi eE\\{|n_{o}|^{2}\\}\right)$ [15, 16].
The MI can be given as:

$\displaystyle
I=H(y)-H(n_{o})=\frac{1}{2}\log_{2}\left(E\left\\{\frac{|y|^{2}}{|n_{o}|^{2}}\right\\}\right)$
(8)

For radar sensing, the MI at the BS according to (1) can be given as
$I_{1}=\tau_{1}W\log_{2}(1+\gamma_{1})$. Based on (3) and (5), the MI at the
ISaC user (i.e., communication performance) and at the BS (i.e., sensing
performance) can be given as $I_{2}^{d}=\tau_{2}W\log_{2}(1+\gamma_{2}^{d})$
and $I_{2}^{u}=\tau_{2}W\log_{2}(1+\gamma_{2}^{u})$, respectively. Thus,
$I_{2}=I_{2}^{d}+I_{2}^{u}$.  For communications, according to (7), the data
rate can be determined by $I_{3}=\tau_{3}W\log_{2}(1+\gamma_{3})$.

The unified optimization framework to maximize the weighted sum of the MI and
data rate via the joint optimization of SP and PA can be mathematically
formulated as follows:

$\centering\begin{array}[]{ll}\mathop{\max}\limits_{{\bf
P},{\boldsymbol{\tau}}}\displaystyle\sum_{s=1}^{3}\Gamma_{s}\tau_{s}\mathrm{log}_{2}(1+\gamma_{s})&{\rm
s.t.}\\\ \mathbf{C}_{1}:\tau_{1}W\log_{2}(1+\gamma_{1})\geq R_{\rm
r},&\mathbf{C}_{5}:\sum\limits_{s=1}^{3}\tau_{s}=1,\\\
\mathbf{C}_{2}:\tau_{2}W\log_{2}(1+\gamma_{2}^{d})\geq R_{\rm
c},&\mathbf{C}_{6}:\sum\limits_{s=1}^{3}P_{s}\leq P_{\max},\\\
\mathbf{C}_{3}:\tau_{2}W\log_{2}(1+\gamma_{2}^{u})\geq R_{\rm
r},&\mathbf{C}_{7}:\tau_{s},P_{s}>0,\,\forall s,\\\
\mathbf{C}_{4}:\tau_{3}W\log_{2}(1+\gamma_{3})\geq R_{\rm c},&\\\
\end{array}\@add@centering$ (9)

where $\Gamma_{s}\geq 0$ is the priority of the $s$-th user, $R_{r}$ and
$R_{c}$ are the minimum sensing and communication QoS requirements, and
$P_{\max}$ is the maximum allowable transmit power. The constraints
$\mathbf{C}_{1}$, $\mathbf{C}_{2}$, $\mathbf{C}_{3}$, and $\mathbf{C}_{4}$ are
the QoS requirements for radar target, ISaC user’s downlink communication,
ISaC user’s uplink sensing, and communication user, respectively. Constraint
$\mathbf{C}_{5}$ ensures that the sum of the bandwidth coefficients must be 1.
Constraint $\mathbf{C}_{6}$ is the transmit power budget. Constraint
$\mathbf{C}_{7}$ ensures that the SP and PA variables are non-negative.

It is shown in the sequel that the optimization problem in (9) is convex. Let
${\mathcal{M}}(x,q)=Wx\log_{2}(1+\frac{a_{s}q}{b_{s}q+c_{s}x})$, dom
${\mathcal{M}}=\\{x,q|0<x\leq 1,q>0\\}$, where $a_{s}$, $b_{s}$, and $c_{s}$
are non-negative constants. ${\mathcal{M}}(x,q)$ represents the MI expression
with clutter interference for the radar target ($s=1$) and ISaC user ($s=2$).
For $s=1$, $a_{1}=G_{\mathrm{tx}}d_{1}^{-2\alpha_{r}}\sigma_{\rm
RCS}\lambda^{2}\left|g_{1}\right|^{2}$,
$b_{1}=\sum_{j=1}^{J}G_{\mathrm{tx}}d_{j}^{-2\alpha_{r}}\sigma_{\rm
RCS}\lambda^{2}\left|\zeta_{j}\right|^{2}$ and $c_{1}=(4\pi)^{3}k_{B}T_{\rm
temp}W$. For $s=2$, $a_{2}=G_{\mathrm{tx}}d_{2}^{-2\alpha_{r}}\sigma_{\rm
RCS}\lambda^{2}\left|g_{2}\right|^{2}$,
$b_{2}=\sum_{j=1}^{J}G_{\mathrm{tx}}d_{j}^{-2\alpha_{r}}\sigma_{\rm
RCS}\lambda^{2}\left|\zeta_{j}\right|^{2}$, and
$c_{2}=(4\pi)^{3}k_{B}T_{\mathrm{temp}}W$.

Also, ${\mathcal{K}}(x,q)=Wx\log_{2}(1+\frac{dq}{ex})$ dom
${\mathcal{K}}=\\{x|0<x\leq 1,q>0\\}$, where $d,e>0$, denotes the data rate of
the communication user ($s=3$). $d$ is
$d=G_{\mathrm{tx}}d_{3}^{-\alpha_{c}}c^{2}\left|h_{3}\right|^{2}$, and $e$ is
given by $e=(4\pi f_{c})^{2}k_{B}T_{\mathrm{temp}}W$.

We show that the multivariable functions ${\mathcal{M}}$ and ${\mathcal{K}}$
are concave by analyzing their Hessian matrices. The Hessian matrix for
${\mathcal{M}}(x,q)$ is given by (10).

$H_{\mathcal{M}}=\left[\begin{matrix}-\frac{W(a^{2}xq^{2}c^{2}+2axq^{2}bc^{2}+2aq^{3}b^{2}c+2a^{2}q^{3}bc)}{\ln\left(2\right)\left(xc+aq+qb\right)^{2}\left(xc+qb\right)^{2}}\\\
\\\
\frac{W(2axq^{2}b^{2}c+2a^{2}xq^{2}bc+a^{2}x^{2}qc^{2}+2ax^{2}qbc^{2})}{\ln\left(2\right)\left(xc+aq+qb\right)^{2}\left(xc+qb\right)^{2}}\end{matrix}\quad\begin{matrix}\frac{acW\left(2x^{2}qbc+ax^{2}qc+2xq^{2}b^{2}+2axq^{2}b\right)}{\ln\left(2\right)\left(xc+aq+qb\right)^{2}\left(xc+qb\right)^{2}}\\\
\\\
-\frac{ax^{2}cW\left(2qb^{2}+2aqb+axc+2xbc\right)}{\ln\left(2\right)\left(xc+aq+qb\right)^{2}\left(xc+qb\right)^{2}}\end{matrix}\right].$
(10)

Clearly, the determinant of the first leading principle $D_{1}$, denoted
$|D_{1}|$, is less than zero. The determinant of the second leading principle
$|D_{2}|=|H_{\mathcal{M}}|=0$. This concludes that the MI function is concave.

Similarly for ${\mathcal{K}}(x,q)$, the Hessian matrix is given by

$H_{\mathcal{K}}=\begin{bmatrix}-\frac{a^{2}q^{2}}{\ln\left(2\right)x\left(xb+aq\right)^{2}}&\frac{a^{2}q}{\ln\left(2\right)\left(xb+aq\right)^{2}}\\\
\frac{a^{2}q}{\ln\left(2\right)\left(xb+aq\right)^{2}}&-\frac{a^{2}x}{\ln\left(2\right)\left(xb+aq\right)^{2}}\end{bmatrix}.$
(11)

For this matrix, it can be shown that $|D_{1}|<0$ and
$|D_{2}|=|H_{\mathcal{K}}|=0$. This concludes that the data rate expression is
concave. Since the sum of concave functions is concave and the constraint sets
are either concave or affine, we can conclude that the optimization problem in
(9) is a convex problem which can be solved using standard convex optimization
solvers.

###  III-B Energy Efficiency Optimization

EE in semi-ISaC systems can be defined as the ratio of the sum of the data
rate and MI to the total power consumed. Mathematically, it can be written as:

$EE=\frac{{\mathcal{A}}({\bf P},{\boldsymbol{\tau}})}{{\mathcal{B}}({\bf
P})}=\frac{\sum_{s=1}^{3}\Gamma_{s}\tau_{s}\mathrm{log}_{2}(1+\gamma_{s})}{P_{1}+P_{2}+P_{3}+\omega}$
(12)

where $\omega$ denotes the circuit power consumption. The unified joint SP and
PA to maximize EE in semi-ISaC networks can be formulated as:

$\centering\begin{array}[]{l}\mathop{\max}\limits_{{\bf
P},{\boldsymbol{\tau}}}\,\,EE\\\
{\rm{s}}{\rm{.t}}{\rm{.}}\,\,\,\mathbf{C}_{1},\mathbf{C}_{2},\mathbf{C}_{3},\mathbf{C}_{4},\mathbf{C}_{5},\mathbf{C}_{6},\mathbf{C}_{7}.\end{array}\@add@centering$
(13)

Algorithm 1 Optimal SP and PA in Semi-ISaC Algorithm

1: set $j=0$, $\delta>0,\eta)j=0,F(\eta)>0$;

2: while $F(\eta_{j})>\delta$ do

3: Solve (14) to obtain $P_{j}$ and $\tau_{j}$;

4: Update $j\leftarrow j+1$;

5: Determine
$\eta_{j}\leftarrow{{\mathcal{A}}({P_{j}},{{\tau_{j}}})}/{{\mathcal{B}}({P_{j}},{{\tau_{j}}})}$;

6: Set $F(\eta_{j})$ to objective value of (14)

7: end while

8: output: ${\bf P}^{*},{\boldsymbol{\tau}}^{*}$

The optimization problem in (13) is a concave-convex fractional program since
the the numerator of the objective function is concave and the denominator is
affine. By introducing the parameter $\eta$, the fractional objective function
in (13) can be transformed to an equivalent non-fractional form, and the EE
optimization for a given $\eta$ can be reformulated as:

$\centering\begin{array}[]{l}F(\eta)=\mathop{\max}\limits_{{\bf
P},{\boldsymbol{\tau}}}{{\mathcal{A}}({\bf
P},{\boldsymbol{\tau}})}-\eta{{\mathcal{B}}({\bf P})}\\\
{\rm{s}}{\rm{.t}}{\rm{.}}\,\,\mathbf{C}_{1},\mathbf{C}_{2},\mathbf{C}_{3},\mathbf{C}_{4},\mathbf{C}_{5},\mathbf{C}_{6},\mathbf{C}_{7}.\end{array}\@add@centering$
(14)

The optimal solution ${\bf P}^{*}$ and ${\boldsymbol{\tau}}^{*}$ to (13) is
also optimal for (14) for a certain $\eta^{*}\geq 0$ that satisfies
$F(\eta^{*})=0$. For a given $\eta$, problem (14) is a convex problem in terms
of ${\bf P}$ and ${\boldsymbol{\tau}}$, therefore we employ an iterative
Dinkelbach-type algorithm to compute the optimal PA and SP by updating $\eta$
via (12) to find the roots of the equation $F(\eta)=0$. The proposed algorithm
is summarized in Algorithm 1. The EE problem is in fact a multi-objective
optimization problem for different values of $\eta$.

The Dinkelbach method has been shown to have a linear convergence rate [17].
In addition, it can be seen that since ${{\mathcal{A}}({\bf
P},{\boldsymbol{\tau}})}$ and ${{\mathcal{B}}({\bf P})}$ are concave and
convex, respectively, and the constraints form a convex set, Algorithm 1
requires to solve a convex problem in each iteration, which can be
accomplished with polynomial complexity.

Figure 2: Spectrum Partition as a function of sensing MI threshold ($R_{r}$)
and data rate threshold ($R_{c}$) when $R_{c}=R_{r}$.

Figure 3: Aggregate sensing MI and data rate as a function of sensing MI
threshold ($R_{r}$) and data rate threshold ($R_{c}$).

Figure 4: Sensing MI or data rate of each user as a function of the priority
of ISaC user and minimum QoS requirements.

##  IV Simulation results and Discussions

In this section, we perform comprehensive simulations to evaluate the
effectiveness of the proposed algorithms. A cell radius of 40 m is considered
with uniformly distributed users. Unless stated otherwise, the following
simulation parameters have been used, i.e., $W=100$ MHz [6], $T_{\rm
temp}=724$ K, $f_{c}=10$ GHz [7], $\sigma_{\rm RCS}=0.1$ $m^{2}$[7],
$\alpha_{r}=\alpha_{c}=2.5$ [7, 11], $P_{\max}=46$ dBm, $J=2$ [6],
$|\zeta_{1}|^{2}=0.01,|\zeta_{2}|^{2}=0.001$ [6], $\omega=33$ dBm [3],
$R_{r}=5$ Mbps, and $R_{c}=20$ Mbps. The users’ priorities are set as
$\Gamma_{1}=\Gamma_{2}=\Gamma_{3}=1/3$ [10]. The small-scale fading is modeled
according to Nakagami-$m$ fading, with $m=3$. The following benchmarks are
used for performance comparison. SP with equal power allocation (SP-EPA), PA
with equal spectrum partitioning (PA-ESP), and random allocation (RA) with
feasible solutions.

Fig. 4 depicts that the portion of the spectrum that is allocated to the
different users vary significantly between the proposed approach and the RA
scheme. For example, less spectrum is allocated to the ISaC user in RA scheme
which underscores the significance of allocating more spectrum to ISaC user
for optimal performance. The spectrum allocations for SP-EPA and PA-ESP are
not shown since the former follows the trend for the proposed scheme and the
latter has no spectrum optimization involved.

Figure 4 depicts the impact of varying the minimum communication and sensing
QoS requirements on the sum of MI and data rate considering the proposed
approach and benchmarks. It can be observed that the proposed joint approach
outperforms the baseline schemes by achieving an average of 10$\%$, 43$\%$,
and 67$\%$ performance gain as compared to the SP-EPA, PA-ESP, and RA schemes,
respectively. This highlights the significance of joint SP and PA. Moreover,
while the performance of the joint approach and the SP-EPA reduces with
increasing $R_{c}$ and $R_{r}$, the sum of MI and data rate for the PA-ESP and
RA benchmarks slightly improves. These trade-offs can be explained as follows.
Typically, the communication-only user contributes more to the objective
function (i.e., aggregate sensing MI and data rate) due to one-way PL effect.
Hence, more spectrum is allocated to communication-only user as compared to
sensing users (as shown in Fig. 4). Subsequently, the sensing performance
almost always achieves $R_{r}$ at most, whereas communication user almost
always exceeds $R_{c}$. When $R_{c}$ and $R_{r}$ increase, more spectrum is
needed to guarantee the higher sensing QoS requirement due to the two-way
path-loss. Thus, the spectrum allocated to communication-only user reduces
while that of sensing-only and ISaC users increases, as evident in Fig. 4.
This spectrum reduction leads to the reduction in the aggregate sensing MI and
data rate of all users for the proposed joint approach and the SP-EPA method
as shown in Fig. 3. Note that the performance of the PA-ESP and RA improves
with increasing $R_{c}$ and $R_{r}$ as no SP optimization is involved.

Figure 4 examines the effect of varying the priority of ISaC user $\Gamma_{2}$
on the performance of the semi-ISaC network for different QoS threshold
values. In this figure, the priority of the ISaC user $\Gamma_{2}$ is varied
by fixing $\Gamma_{1}=\Gamma_{3}$ such that
$\Gamma_{1}=\Gamma_{3}=(1-\Gamma_{2})/2$. As expected, the performance of the
communication-only user and the sensing-only user decreases, while the MI for
the ISaC user improves. However, the performance loss of sensing is evident
when sensing requirement $R_{r}$ is small and the loss is mostly insignificant
compared to communications-only user. Moreover, while setting $R_{r}>R_{c}$
reduces the data rate of the communication-only user, the MI of the sensing-
only user and ISaC user increases. This is due to the fact that increasing
$R_{r}$ from 5 Mbps to 30 Mbps while reducing $R_{c}$ from 20 Mbps to 5 Mbps
means allocating more resources to satisfy higher sensing QoS requirement. For
$\Gamma_{2}\geq 0.5$, the curves when $R_{r}<R_{c}$ and $R_{r}>R_{c}$ for the
ISaC user almost overlap since the user has very high priority.

Figure 7 depicts the impact of varying the values of the RCS and the maximum
transmit power on the sum of MI and data rate for the proposed approach.
First, it can be observed that increasing the value of the RCS leads to higher
sensing MI and data rate. This is because, the RCS value physically represents
the area into which targets reflects incoming radar signal. The greater the
RCS, the stronger the reflected radar echoes, leading to higher signal
strength for sensing purposes. Second, augmenting the transmit power budget
leads to higher aggregate sensing MI and data rate. This observation is the
same for different RCS values. The reason is, increasing the value of the
maximum transmit power results in a larger feasible region for the joint
optimization problem.

Figure 7 illustrates the EE and convergence behavior of Algorithm 1 under
different RCS and number of clutters for ISaC and sensing-only users. It can
be seen that the proposed algorithm converges within no more than 5 iterations
on average under different RCS and number of clutters, which verifies the fast
convergence of Algorithm 1.

Figure 5: Aggregate sensing MI and data rate as a function of transmit power
of Semi-ISaC BS and target RCS.

Figure 6: Convergence behavior of the proposed Algorithm 1.

Figure 7: EE as a function of sensing MI threshold ($R_{r}$) and data rate
threshold ($R_{c}$) when $R_{c}=R_{r}$.

Finally, Figure 7 compares the average EE of the proposed joint approach and
the considered benchmarks for different minimum communication and sensing QoS
requirements. It can be observed that the proposed scheme achieves superior EE
than the benchmarks, demonstrating the significance of joint SA and PA. The EE
reduces with increasing $R_{c}$ and $R_{r}$ since higher QoS threshold
restricts the feasible region.

##  V Conclusion

In this letter, a unified joint PA and SP optimization framework to maximize
the sum of MI and data rate, and EE in a semi-ISaC system has been proposed.
The sum of MI and data rate optimization problem was shown to be convex and
can be solved using standard convex optimization solvers. On the other hand,
the non-convex EE problem was transformed into into a concave-convex
fractional problem that can be solved via the Dinkelbach-method. Simulation
results have demonstrated the efficiency and the fast convergence of the
proposed joint schemes, the advantage of allocating more spectrum resources to
the ISaC user, and the impact of varying user priority on the performance of
the semi-ISaC network. Further research is required to consider the impact of
imperfect channel estimation and multi-antenna BSs.

## References

  * [1] F. Dong _et al._ , “Sensing as a service in 6G perceptive networks: A unified framework for ISAC resource allocation,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 5, pp. 3522–3536, May 2023.
  * [2] B. K. Chalise and B. Himed, “Performance tradeoff in a unified multi-static passive radar and communication system,” in _Proc. IEEE Radar Conf. (RadarConf)_ , Oklahoma City, OK, USA, Apr. 2018, pp. 653–658.
  * [3] J. Zou _et al._ , “Energy efficiency optimization for integrated sensing and communications systems,” in _Proc. IEEE Wireless Commun. Netw. Conf. (WCNC)_ , Austin, TX, USA, Apr. 2022, pp. 216–221.
  * [4] B. Zhao _et al._ , “Integrated sensing and communication aided dynamic resource allocation for random access in satellite terrestrial relay networks,” _IEEE Commun. Lett._ , vol. 27, no. 2, pp. 661–665, Feb. 2023\.
  * [5] F. Wang and H. Li, “Power allocation for coexisting multicarrier radar and communication systems in cluttered environments,” _IEEE Trans. Signal Process_ , vol. 69, pp. 1603–1613, 2021.
  * [6] M. Ashraf _et al._ , “Joint optimization of radar and communications performance in 6G cellular systems,” _IEEE Trans. Green Commun._ , vol. 7, no. 1, pp. 522–536, 2023.
  * [7] C. Zhang _et al._ , “Semi-integrated-sensing-and-communication (Semi-ISaC): From OMA to NOMA,” _IEEE Trans. Commun._ , vol. 71, no. 4, pp. 1878–1893, Apr. 2023.
  * [8] A. Bayesteh _et al._ , “Integrated sensing and communication (ISAC)—from concept to practice,” _Communications of Huawei Research_ , pp. 4–25, Sep. 2022.
  * [9] M. L. Rahman _et al._ , “Framework for a perceptive mobile network using joint communication and radar sensing,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 56, no. 3, pp. 1926–1941, Jun. 2020.
  * [10] Q. Qi _et al._ , “Integrating sensing, computing, and communication in 6G wireless networks: Design and optimization,” _IEEE Trans. Commun._ , vol. 70, no. 9, pp. 6212–6227, 2022.
  * [11] P. Kumari, S. A. Vorobyov, and R. W. Heath, “Adaptive virtual waveform design for millimeter-wave joint communication–radar,” _IEEE Trans. Signal Process._ , vol. 68, pp. 715–730, 2020.
  * [12] Z. He _et al._ , “Full-duplex communication for isac: Joint beamforming and power optimization,” _IEEE J. Sel. Areas Commun._ , vol. 41, no. 9, pp. 2920–2936, Sep. 2023.
  * [13] M. F. Keskin, V. Koivunen, and H. Wymeersch, “Limited feedforward waveform design for ofdm dual-functional radar-communications,” _IEEE Trans. Signal Process._ , vol. 69, pp. 2955–2970, 2021.
  * [14] Q. Zhang _et al._ , “Design and performance evaluation of joint sensing and communication integrated system for 5G mmWave enabled CAVs,” _IEEE J. Sel. Top. Signal Process._ , vol. 15, no. 6, pp. 1500–1514, 2021\.
  * [15] T. Tian _et al._ , “Mutual information-based power allocation and co-design for multicarrier radar and communication systems in coexistence,” _IEEE Access_ , vol. 7, pp. 159 300–159 312, 2019.
  * [16] T. M. Cover and J. A. Thomas, _Elements of Information Theory_. Hoboken, NJ, USA: Wiley, 2012.
  * [17] J. P. Crouzeix, J. A. Ferland, and S. Schaible, “An algorithm for generalized fractional programs,” _J. Optim. Theory Appl._ , vol. 47, no. 1, pp. 35–49, 1985.