\section*{Nomenclature}

\vspace{-0.4cm}
\section{Introduction}
\begin{comment}

{\color{black}
\begin{itemize}[leftmargin=*]
    \item We proposed a new distributed optimization approach to solve the D-OPF problems in radial power distribution systems. Specifically, the radial topology of the distribution system is leveraged in developing a computationally tractable and practically viable D-OPF algorithm.
    The proposed D-OPF method converges to the solution of the same quality as C-OPF and takes a significantly fewer number of macro-iterations compared to the state-of-art distributed optimization methods;
    \item We provide a convergence analysis for the proposed D-OPF algorithm. Specifically, it is shown that the proposed D-OPF algorithm is guaranteed to converge at-most linearly at the boundary bus. Further, if the local subproblems converge and the centralized OPF admits a feasible solution, the proposed D-OPF method is guaranteed to converge.
    \item Extensive validation of the proposed approach is presented using IEEE 123-bus and IEEE 8500-node test system. It is demonstrated that the proposed D-OPF algorithm converges to the same solution as the C-OPF problem. Thus, the solution obtained from the D-OPF algorithm is of the same quality as obtained via solving an equivalent C-OPF algorithm. Further, compared to state-of-the-art D-OPF techniques, for example, ADMM, the proposed D-OPF algorithm takes an order of magnitude fewer communication rounds among distributed agents (macro-iterations) to converge. The applicability of the proposed approach is also demonstrated for an unbalanced distribution feeder.
\end{itemize}}

\end{comment}

Optimization of power distribution systems has gained significant attention with growing penetrations of controllable distributed energy resources (DERs). Traditionally, the distribution system is centrally managed using an advanced distribution management system (ADMS) that collects the measurement data and use it to coordinate DERs and other grid-edge technologies. With the increasing number of such controllable assets and growing system complexity, a centrally managed system suffers from several challenges including poor scalability, and susceptibility to a single-point failures \cite{momoh1999review1,castillo2013survey,molzahn2017survey}. The additional nonlinearity introduced by the three-phase and unbalanced power distribution systems may pose convergence issues for the centralized optimization problems and may result in failure to attain the optimal solution within a reasonable time \cite{molzahn2017survey}. Besides, a centrally managed system is impractical for emerging distribution system typically composed of several geographically separated areas controlled by different entities that may not want to share model and/or data. To address these challenges, distributed coordination and optimization paradigm has gained significant attention for grid-edge coordination.

The existing literature has shown the applicability of the distributed optimization paradigm in reducing the computational requirements for optimal power flow (OPF) problems and in mitigating other challenges including ownership boundaries and avoiding single-point failures \cite{molzahn2017survey}. To this end, general distributed optimization techniques, such as, Alternating Direction Method of Multipliers (ADMM), Auxiliary Problem Principle (APP), etc. have been adopted to solve distributed OPF problems for power distribution systems \cite{ kim2000comparison,zheng2015fully,boyd2011distributed}. Distributed OPF algorithms involve two iterative processes: (1) sub-problem optimization involving micro-iterations, where distributed agents individually solve optimization problems specific to their areas, (2) distributed coordination requiring macro-iterations, where distributed agents exchange the boundary variables typically with their immediate neighbors. The overall compute complexity of the algorithm depends upon the total number micro- and macro- iterations taken to achieve network-level convergence. The nature of the optimization problem dictates micro-iterations needed to solve each sub-problem; the added nonlinearity with three-phase OPF may pose convergence and compute challenges for this step. The distributed coordination approach dictates the number of macro-iterations needed for network-wide convergence; traditional algorithms take 100s of rounds of communication exchange or macro-iterations resulting in a slow converge. This calls for algorithmic advances in distributed OPF algorithms to realistically apply them to a real-world distribution system.

Although several recent methods attempted to reduce the computational requirements for local subproblems \cite{qu2017harnessing,mhanna2018adaptive,peng2016distributed}, challenges with convergence and speed persist due to added nonlinearities in three-phase unbalanced OPF models for power distribution system. Earlier {\em optimization with simulation} has been introduced for general complex optimization problems that can leverage the already developed complex models of the system, and avoid including the detailed system model in the optimization problem formulation \cite{april2004new, nakayama2002simulation,nguyen2014review}. Although this approach has not been directly applied to solve OPFs for power systems, the use of high-fidelity simulation models for design optimization has been reported \cite{zhou2019digital,pan2020digital}. This high-fidelity simulation model is referred as a digital twin (DT) of the physical power system and has been employed to aid with situational awareness and decision-making for power grid operations  \cite{brosinsky2018recent, he2019preliminary, darbali2021state}. The DT models provide an ability to perform high-fidelity simulations of the grid that can potentially help reduce the complexity of difficult OPF problems. However, the use of DT models has been limited to evaluating the physical system in advance for potential actions; their utility in directly assisting with difficult power system optimization problems has not been investigated.

Another challenge with traditional algorithms to solve distributed OPF algorithms emerge from a large number of macro-iterations ($\ge 100$ iterations) needed to converge for relatively small and/or balanced systems \cite{erseghe2014distributed, dall2013distributed,molzahn2017survey}.  Thus, these algorithms require an expensive communication infrastructure with high bandwidth and low communication delays. Although real-time feedback-based online distributed algorithms address some of the challenges \cite{bernstein2019real,magnusson2020distributed}, they still require thousands of time-steps to track the optimal solution, resulting in a sub-optimal system operation and constraint violations. Recently, authors proposed an algorithm based on the structure of the power flow problem for a radial distribution system that drastically reduces the number of macro-iterations  \cite{ sadnan2021distributed }. However, the sub-problems are modeled as nonlinear programming (NLP) problems that may pose convergence and compute challenges. Moreover, the performance of distributed OPF algorithms has not been evaluated in a {\em realistic cyber-physical environment} with appropriate models for  communication systems ~\cite {sadnan2021distributed, erseghe2014distributed,magnusson2020distributed}.

The objective of this paper is to develop a fast, computationally tractable, and communication-efficient distributed OPF (D-OPF) algorithm for {\em{unbalanced}} power distribution system. Our approach employs {\em optimization with simulation} concepts using DTs (at sub-problem level) and a distributed coordination algorithm to obtain a computationally efficient and near-optimal D-OPF approach. Specifically, we optimize a low-fidelity linear model of the sub-system and project the obtained solution to the high-fidelity simulation model using DTs resulting in a near-optimal but feasible power flow solution. Primarily, the use of DTs (at sub-problem level) excludes the necessity of solving complex nonlinear OPFs and reduces the overall solution time for the local subproblems. Then the distributed coordination algorithm is employed to exchange the projected power flow variables with the neighboring distributed areas. We carefully select shared variables based on the power flow in radial power distribution systems that help us reduce the required number of macro-iterations by an order of magnitude. Together these innovations result in highly scalable and near-optimal algorithms for solving difficult D-OPF problems.  The proposed approach is also extensively evaluated using a comprehensive cyber-power co-simulation environment.

\vspace{-0.3cm}
\section{Network Modeling \& Problem Formulation}
In this paper, $(\cdot)^T$ represents matrix transpose; $(\cdot)^*$ represents the complex-conjugate; $(\cdot)^{(n)}$ represents the $n^{th}$ iteration;  $\mathbb{Re}$ denotes the real part of the complex number; the superscript $p$ (without parenthesis) denotes the three-phases, i.e., $\{a,b,c\}$ of the system.
\vspace{-0.3cm}
\subsection{Power Flow Model}
\vspace{-0.1cm}
This section details the power flow models to formulate the OPF problems for an unbalanced power distribution system. The unbalanced power flow equations are based on our prior work \cite{jha2019bi}. Let us consider an unbalanced radial power distribution network of $n$ buses where, $\mathcal{N}$ denotes the set of buses in that system and $\mathcal{E}$ denotes the set of edges identifying distribution lines that connect the ordered pair of buses $\{ij\}$, $\forall$ $i,j \in \mathcal{N}$. Here, $\phi_j$ denotes the set of phases in the bus $j$. Let $v_j^p =|V_j^p|^2$ be the squared magnitude of voltage at bus $j\in \mathcal{N}$ for phase $p\in \phi_j$. Define $\phi_{ij} = \{pq : p\in \phi_i$ and $q\in \phi_j$  $\forall \{ij\}\in \mathcal{E}\} $. Let $l_{ij}^{pq} =(|I_{ij}^p||I_{ij}^q|)$ be the squared magnitude of the line current flowing in the phase ${pq}\in \phi_{ij}$ of line $(i,j)$. Also, $S^{pq}_{ij}=P^{pq}_{ij}+ Q^{pq}_{ij}$ and $z^{pq}_{ij}=r^{pq}_{ij}+j x^{pq}_{ij}$, where $pq \in \phi_{ij}$. Variable $p_{L,j}^p$ and $q_{L,j}^p$ denote the active and reactive load (respectively) connected at node $j$ of phase $p\in \phi_j$. Similarly, subscript $D,j$ denotes the DER generation at node $j$; e.g., $p^p_{D,j}$ denotes the active power generation at node $j$ for phase $p$. $\delta^{pq}_{ij}$ is the angle difference between the phase currents. First, we define the nonlinear power flow model and then the linearized model is detailed. For more details, please refer to \cite{jha2019bi}.

\subsubsection{Nonlinear Model}
With the approximated phase voltage and branch current angles, in \cite{jha2019bi} developed the nonlinear power flow model for an unbalanced radial power distribution system. The model is defined below in \eqref{nl_pf}. We model the loads in the systems as constant power loads.

\vspace{-0.4cm}
\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{nl_pf} \IEEEyessubnumber*
\nonumber  P_{ij}^{pp} - \sum_{q \in \phi_j}{l_{ij}^{pq} \left(r_{ij}^{pq} \cos(\delta_{ij}^{pq})- x_{ij}^{pq} \sin(\delta_{ij}^{pq})\right)}\\
   =\sum_{k:j \rightarrow k}P_{jk}^{pp} + p_{L,j}^p- p_{D,j}^p\\
    \nonumber  Q_{ij}^{pp} - \sum_{q \in \phi_j}{l_{ij}^{pq}\left(x_{ij}^{pq} \cos(\delta_{ij}^{pq})+ r_{ij}^{pq} \sin(\delta_{ij}^{pq})\right)}\\
 = \sum_{k:j \rightarrow k}Q_{jk}^{pp}+q_{L,j}^p  - q_{D,j}^p - q_{C,j}^p \\
   \nonumber v_j^p = v_i^p - \sum_{q \in \phi_j}{2 \mathbb{Re}\left[S_{ij}^{pq} (z_{ij}^{pq})^*\right]} + \sum_{q \in \phi_j}{z_{ij}^{pq} l_{ij}^{qq}}\\
   +\sum_{q1,q2 \in \phi_j, q1 \neq q2}{2\mathbb{Re}\left[ z_{ij}^{pq1} l_{ij}^{q1q2}\left(\angle(\delta_{ij}^{q1q2})\right)(z_{ij}^{pq2})^*\right]}\\
   (P_{ij}^{pp})^2 + (Q_{ij}^{pp})^2 = v_i^p  l_{ij}^{pp}\\
   (l_{ij}^{pq})^2 = l_{ij}^{pp}  l_{ij}^{qq}
\end{IEEEeqnarray}
\end{small}
\vspace{-0.45cm}

\subsubsection{Linear-approximated Model}
The computational complexity augmented by the nonlinear power flow models to the existing scalability issues associated with large-scale OPF problems can be reduced using the approximated linearized power flow models. The linearized power flow models can solve the power flow variables, such as nodal voltages, within reasonable accuracy and reduce the computation time by the order of magnitudes making the linear models scalable for large-scale power distribution systems. A linear-approximated power flow model, three-phase LinDistFlow model, is defined in \eqref{lin_pf}. The primary assumption is that the line loses are negligible compared to the power flow in the system.

\vspace{-0.4cm}
\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{lin_pf} \IEEEyessubnumber*
P_{ij}^{pp} =\sum_{k:j \rightarrow k}P_{jk}^{pp} + p_{L, j}^p- p_{D, j}^p\\
Q_{ij}^{pp} = \sum_{k:j \rightarrow k}Q_{jk}^{pp}+q_{L,j}^p  - q_{D,j}^p - q_{C,j}^p \\
v_j^p = v_i^p - \sum_{q \in \phi_j}{2 \mathbb{Re}\left[S_{ij}^{pq} (z_{ij}^{pq})^*\right]}
\end{IEEEeqnarray}
\end{small}
\vspace{-0.85cm}

\subsection{DER Modeling}
We define the DER models for different network objectives. In this paper, to achieve minimal cost, either reactive or active power generation of the DERs is set as the decision variables of the OPF problems. The active power generation, $p^p_{D,j}$ is assumed to measured and known, when the reactive power generation, $q^p_{D,j}$ is the decision variable. If $S^p_{DR,j}$ is defined as the nominal rating of the DER at node $j$ for phase $p\in \phi_j$, then the DER is modeled as \eqref{DER_Q}.

\vspace{-0.3cm}
\begin{small}
\begin{equation}\label{DER_Q}
-\sqrt{(S_{DR,j}^{p})^2 - (p_{D,j}^p)^2} \leq q_{D,j}^p \leq \sqrt{(S_{DR,j}^{p})^2 - (p_{D,j}^p)^2}
\end{equation}
\end{small}
\vspace{-0.4cm}

For network objectives, such as DER maximization, the active power generation, $p^p_{D,j}$ is set as decision variables, and we assume $q^p_{D,j} = 0$. The DER model is defined in \eqref{DER_P}. However, if both active and reactive power generation is considered as decision variables, then the DER is modeled as a second-order conic equation as defined in \eqref{DER_PQ}.

\vspace{-0.2cm}
\begin{small}
\begin{equation}\label{DER_P}
\small
0 \leq p_{D,j}^p \leq S_{DR,j}^{p}
\end{equation}
\vspace{-0.4cm}
\begin{equation}\label{DER_PQ}
\small
\left(p_{D,j}^p\right)^2+ \left(q_{D,j}^p\right)^2 \leq \left(S_{DR,j}^{p}\right)^2
\end{equation}
\end{small}
\vspace{-0.9cm}

\subsection{Centralized OPF Problem}
In this section, we define the centralized OPF problems for an unbalanced radial distribution system, that can be solved using a central controller; both nonlinear and linear power flow models have been used for such OPF problems with minimization and maximization objectives.

\subsubsection{Loss Minimization}
First, the active power loss minimization problem has been formulated where the line loss is reduced by generating optimal $q^p_{D,j}$, without violating operational limits. The optimization variable is $X = [P^{pp}_{ij}, Q^{pp}_{ij}, v^p_j, l^{pq}_{ij}, q^p_{D,j}]^T$, $ \forall j\in \mathcal{N}$ \& $\forall \{ij\} \in \mathcal{E}$. The OPF problem \textbf{(C1)} is defined in \eqref{loss_opf_nl} using nonlinear power flow models. Here, $\left(I^{rated}_{ij}\right)^2$ is the thermal limit for the line $\{ij\}\in E$, and $V_{max}$ \& $V_{min}$ denotes the nodal voltage limit.

\vspace{-0.4cm}
\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{loss_opf_nl} \IEEEyessubnumber*
\text{\textbf{(C1)}}  \hspace{0.4cm}  \text{Minimize} \sum_{p \in \phi_j, j:i \rightarrow j}{l^{pp}_{ij}r^{pp}_{ij}}\\
  \text{s.t.}  \hspace{0.3cm}  \text{equation \eqref{nl_pf} \& \eqref{DER_Q} } \label{C1b} \\
  l_{ij}^{p} \leq \left(I^{rated}_{ij}\right)^2  \text{\&} \hspace{0.3cm}
   {V_{min}}^2\leq v_j^p \leq {V_{max}}^2  \label{C1c}
\end{IEEEeqnarray}
\end{small}
\vspace{-0.4cm}

To take advantage of the computational efficacy and scalability, next we define the OPF problem with linear power flow models to minimize the line losses. To further ease the computations,
the nodal voltages are approximated as $1.00$ pu, and thus, the line loss has been approximated by a convex cost function $\sum \left(P^{pp}_{ij}\right)^2+\left(Q^{pp}_{ij}\right)^2$ for all edges $\{ij\}\in \mathcal{E}$. Here for the OPF problem with linear power flow model, the optimization variable is $X = [P^{pp}_{ij}, Q^{pp}_{ij}, v^p_j, q^p_{D,j}]^T$, $ \forall j\in \mathcal{N}$ \& $\forall \{ij\} \in \mathcal{E}$. The OPF problem with linear power flow model, \textbf{(C2)} is defined in \eqref{loss_opf_lin}.

\vspace{-0.4cm}
\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{loss_opf_lin} \IEEEyessubnumber*
\text{\textbf{(C2)}}  \hspace{0.4cm}  \text{Minimize} \sum_{p \in \phi_j, j:i \rightarrow j} {\left(P^{pp}_{ij}\right)^2}+{\left(Q^{pp}_{ij}\right)^2}\\
  \text{s.t.}  \hspace{0.3cm}  \text{equation \eqref{lin_pf} \& \eqref{DER_Q} }  \label{C2b} \\
   {V_{min}}^2\leq v_j^p \leq {V_{max}}^2  \label{C2c}
\end{IEEEeqnarray}
\end{small}
\vspace{-0.5cm}
\subsubsection{DER Maximization}
As a maximization problem, DER maximization objective has been formulated in this paper. Specifically, the network objective is to maximize active power generations, $p_{D,j}^p$, from all the DERs in all of the phases, without exceeding any operational limits. For the DER maximization OPF, the optimization variable is $X = [P^{pp}_{ij}, Q^{pp}_{ij}, v^p_j, l^{pq}_{ij}, p^p_{D,j}]^T$, $ \forall j\in \mathcal{N}$ \& $\forall \{ij\} \in \mathcal{E}$. The OPF problem \textbf{(C3)} is defined in \eqref{der_opf_nl} using the nonlinear power flow models. Similar to the previous OPF problems, $\left(I^{rated}_{ij}\right)^2$ is the thermal limit for the line $\{ij\}\in E$, and $V_{max}$ \& $V_{min}$ denotes the nodal voltage limit.

\vspace{-0.4cm}
\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{der_opf_nl} \IEEEyessubnumber*
\text{\textbf{(C3)}}  \hspace{0.4cm}  \text{Maximize} \sum_{\forall j \in \mathcal{N}_{}}\sum_{p \in \phi_j}{p^{p}_{D,j}}\\
  \text{s.t.}  \hspace{0.3cm}  \text{equation \eqref{nl_pf} \& \eqref{DER_P} }  \label{C3b} \\
  l_{ij}^{p} \leq \left(I^{rated}_{ij}\right)^2  \text{\&} \hspace{0.3cm}
   {V_{min}}^2\leq v_j^p \leq {V_{max}}^2  \label{C3c}
\end{IEEEeqnarray}
\end{small}
\vspace{-0.3cm}

Similar to the problem \textbf{(C2)}, the DER maximization objective problem can be formulated as a linear OPF model with linear approximated power flow model (eqn. \eqref{lin_pf}). Since the objective function for this problem is already a linear cost function, thus the whole OPF problem becomes linear optimization problem. The optimization problem is detailed in \eqref{der_opf_lin} and referred as problem \textbf{(C4)}. The optimization variable for \textbf{(C4)} is $X = [P^{pp}_{ij}, Q^{pp}_{ij}, v^p_j, p^p_{D,j}]^T$, $ \forall j\in \mathcal{N}$ \& $\forall \{ij\} \in \mathcal{E}$.

\vspace{-0.3cm}
\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{der_opf_lin} \IEEEyessubnumber*
\text{\textbf{(C4)}}  \hspace{0.4cm}  \text{Maximize} \sum_{\forall j \in \mathcal{N}_{}}\sum_{p \in \phi_j}{p^{p}_{D,j}}\\
  \text{s.t.}  \hspace{0.3cm}  \text{equation \eqref{lin_pf} \& \eqref{DER_P} }  \label{C4b} \\
   {V_{min}}^2\leq v_j^p \leq {V_{max}}^2  \label{C4c}
\end{IEEEeqnarray}
\end{small}
\vspace{-0.7cm}

\section{Distributed Algorithm}
The optimization problem for radial power distribution networks-- composed of several connected areas, are naturally decomposable into multiple sub-problems defined for the areas. In this section, first we detail the decomposition approach for the OPF problems of three-phase unbalanced systems, and then the local sub-problems are defined. Lastly, we briefly discuss the developed cyber-physical co-simulation platform, where both the physical and the communication layer are simulated{\color{black} to study the impact of the developed distributed OPF algorithm on a real life distribution network.}

\vspace{-0.3cm}
\subsection{Decomposition Method \& Network Approximations}
{\color{black}Let the distribution system be composed of $N$ connected areas/microgrids, and $\mathcal{A} = \{A_1,~A_2, …,~A_{N}\}$ be the set of all areas or microgrids. Any area $A_m \in \mathcal{A}$ is defined as the directed radial graph $A_m = \mathcal{G}(\mathcal{N}_m, \mathcal{E}_m)$; $\mathcal{N}_m$ and $\mathcal{E}_m$ defines the set of nodes and lines in the system. The set of central optimization variables is defined as $X = \bigcup_{m=1}^{N} X_m$, where $X_m $ represents the set of all local optimization variables for Area $A_m$. Similarly $\mathcal{S} = \bigcup_{m=1}^{N}\mathcal{S}_m$ represent the set of constraints for overall problem, and $\mathcal{S}_m$ represents the local constraints. If central objective $f(X)$ is a decomposable cost function, then the problem can be decomposed into several local sub-problems and written as \eqref{eqProposed_obj}. Here, $f_m$ is defined as the local objective function for Area $A_m$. The complicating variables, $y_{{m}^'}$ are shared by the neighboring areas of $A_m$, and are kept constant for solving the local optimization problem in this primal decomposition-based method.

\vspace{-0.3cm}
\begin{equation}\label{eqProposed_obj}
\small
\min_{X \in \mathcal{S}}{\hspace{0.1cm} f(X)} = {\hspace{0.1cm}} {\sum_{m=1}^N \min_{X_m \in \mathcal{S}_m} f_m(X_m,y_{{m}^'})}
\end{equation}
\vspace{-0.3cm}

As the overall system topology was radial to begin with, the decomposed areas of the networks are also connected radially and have unique upstream \& downstream sections. This specific structure of the system helps to identify the unique parent area and the child areas for any area $A_m$, which in turns associates the complicating/shared variables, $y_{{m}^'}$ for the local sub-problems -- exchanged among the neighboring computing agents to solve the overall problem. For the proposed decomposition approach, the complicating variables are the voltages and the power flows at the shared bus. Let for any Area $A_m$, bus $0$ is shared with its parent/Upstream area (UA) and node $k$ ($\{jk\} \in \mathcal{E}_m$) is shared with its child/Downstream Area (DA). Then the complicating variable, $y_{{m}^'} = [y^{ua},y^{da}]$, where $y^{ua} = [v_{0_{ua}}]$ , and $y^{da} = [P_{jk_{da}}, Q_{jk_{da}}]$. Here, each area solves related local sub-problems in parallel by assuming the constant complicating variable, i.e., a fixed voltage at the shared bus with the unique UA and constant loads at the shared buses with DAs. In terms of power system's physics, this approximates the whole upstream and downstream network segments of an area for radial power distribution systems.}

After solving local sub-problems, the respective complicating variables, i.e., the total power requirements in that area and the shared (with DA) bus voltages area are shared with UA and DAs, respectively. Then, the sub-problems are solved again with updated shared variables from the neighbors. When a consensus is attained at the shared boundaries, then the decision variables are dispatched within the area. The consensus at the boundary can be achieved using Fixed Point Iteration methods as described in \cite{sadnan2021online}.

\begin{figure*}[t]
    \centering
    \includegraphics[width=0.85\textwidth]{figs/Pictorial_desc_v2.eps}
        \vspace{-0.2cm}
    \caption{Use of Digital Twin Model in Distributed Optimization for Power System}
    \label{DTC}
\end{figure*}

\subsection{Local Distributed OPF Model}
In this section, the decomposed local sub-problems at $n^{th}$ macro-iteration, \textbf{(D1)-(D4)}, of the corresponding central OPFs, \textbf{(C1)-(C4)} are defined. For any Area $A_m$, let bus $0$ is shared with its UA and node $k$ ($\{jk\} \in \mathcal{E}_m$) is shared with its DAs. Also, the subscript of a node represents the area, that has solved the variable; e.g., $v_{0_m}$, $p_{k_m}$ are solved by area $A_m$, where $v_{0_{ua}}$, $P_{jk_{da}}$ are solved by UA and DA, respectively. In the proposed model, each local computing agents solves the corresponding distributed sub-problems, and over macro-iterations, convergence is achieved.

\subsubsection{Loss Minimization} The OPF problem \textbf{(C1)-(C2)} represent the central formulation of the loss minimization problem using nonlinear and approximated models, respectively, where the active power line losses are minimized by controlling the reactive power generations from the DERs. The corresponding decomposed sub-problems for Area $A_m$ , $\forall \{ij\}\in \mathcal{E}_m$ and $\forall j \in \mathcal{N}_m$  are defined in \eqref{D1}-\eqref{D2}, respectively. Constraint \eqref{eqvD1} and \eqref{eqvD2}  represents the approximated UA and DA with the help of complicating/shared variables, as described in the previous section.

\vspace{-0.3cm}
\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{D1} \IEEEyessubnumber*
\text{\textbf{(D1)}}\hspace{0.4cm} \min \hspace{0.2cm} f_m(X_m^{(n)}) =\sum_{p \in \phi_j, j:i \rightarrow j}{l^{pp}_{ij}r^{pp}_{ij}} \\
\text{s.t. equation \eqref{C1b} \& \eqref{C1c} } \label{conD1}\\
\hspace{-2pt}v_{0_m}^{p\hspace{1pt} (n) } = v_{0_{ua}}^{p\hspace{1pt}(n-1) }\label{eqvD1}; \hspace{0.1cm}
p_{k_m}^{p \hspace{1pt}(n)} = P_{jk_{da}}^{pp \hspace{1pt}(n-1)}; \hspace{0.1cm}
q_{k_m}^{p \hspace{1pt}(n)} = Q_{jk_{da}}^{pp \hspace{1pt}(n-1)}
\end{IEEEeqnarray}
\end{small}
\vspace{-0.3cm}

\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{D2} \IEEEyessubnumber*
\text{\textbf{(D2)}}\hspace{0.4cm} \min \hspace{0.2cm} f_m(X_m^{(n)}) =\sum_{p \in \phi_j, j:i \rightarrow j} {\left(P^{pp}_{ij}\right)^2}+{\left(Q^{pp}_{ij}\right)^2}\\
\text{s.t. equation \eqref{C2b} \& \eqref{C2c} } \label{conD2}\\
\hspace{-2pt}v_{0_m}^{p\hspace{1pt} (n) } = v_{0_{ua}}^{p\hspace{1pt}(n-1) }\label{eqvD2}; \hspace{0.1cm}
p_{k_m}^{p \hspace{1pt}(n)} = P_{jk_{da}}^{pp \hspace{1pt}(n-1)}; \hspace{0.1cm}
q_{k_m}^{p \hspace{1pt}(n)} = Q_{jk_{da}}^{pp \hspace{1pt}(n-1)}
\end{IEEEeqnarray}
\end{small}

\subsubsection{DER Maximization} The OPF problem \textbf{(C3)-(C4)} represent the central OPF formulation of the DER generation maximization problem using nonlinear and approximated models, respectively, where the real power generations for all the phases are maximized within the physical limits and without violating the operational constraints. Similar to loss minimization objective, the corresponding decomposed sub-problems for Area $A_m$ are defined in \eqref{D3}-\eqref{D4}, respectively. Similar to loss minimization problem, constraint \eqref{eqvD3} and \eqref{eqvD4}  represents the approximated UA and DA with the help of complicating/shared variables.

\vspace{-0.3cm}
\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{D3} \IEEEyessubnumber*
\text{\textbf{(D3)}}\hspace{0.4cm} \max \hspace{0.2cm} f_m(X_m^{(n)}) =\sum_{\forall j \in \mathcal{N}_{m}}\sum_{p \in \phi_j}{p^{p}_{D,j}}\\
\text{s.t. equation \eqref{C3b} \& \eqref{C3c} } \label{conD3}\\
\hspace{-2pt}v_{0_m}^{p\hspace{1pt} (n) } = v_{0_{ua}}^{p\hspace{1pt}(n-1) }\label{eqvD3}; \hspace{0.1cm}
p_{k_m}^{p \hspace{1pt}(n)} = P_{jk_{da}}^{pp \hspace{1pt}(n-1)}; \hspace{0.1cm}
q_{k_m}^{p \hspace{1pt}(n)} = Q_{jk_{da}}^{pp \hspace{1pt}(n-1)}
\end{IEEEeqnarray}
\end{small}
\vspace{-0.9cm}

\begin{small}
\begin{IEEEeqnarray}{C C}
\IEEEyesnumber\label{D4} \IEEEyessubnumber*
\text{\textbf{(D4)}}\hspace{0.4cm} \max \hspace{0.2cm} f_m(X_m^{(n)}) =\sum_{\forall j \in \mathcal{N}_{m}}\sum_{p \in \phi_j}{p^{p}_{D,j}}\\
\text{s.t. equation \eqref{C4b} \& \eqref{C4c} } \label{conD4}\\
\hspace{-2pt}v_{0_m}^{p\hspace{1pt} (n) } = v_{0_{ua}}^{p\hspace{1pt}(n-1) }\label{eqvD4}; \hspace{0.1cm}
p_{k_m}^{p \hspace{1pt}(n)} = P_{jk_{da}}^{pp \hspace{1pt}(n-1)}; \hspace{0.1cm}
q_{k_m}^{p \hspace{1pt}(n)} = Q_{jk_{da}}^{pp \hspace{1pt}(n-1)}
\end{IEEEeqnarray}
\end{small}
\vspace{-0.3cm}

\vspace{-0.6cm}
\subsection{Simulation Within Optimization: Use of Digital Twin}
To reduce the compute burdens to solve the NLP optimization problem, a method utilizing simulation within optimization has been proposed by incorporating the Digital Twin (DT) of the system. DTs are the detailed virtual model of the real-world system/network, linked to their physical counterpart by the continuous, real-time information flow. Specific to the power distribution systems, DTs are the software-based abstractions of complex physical networks of power distribution systems which are connected via a communication link to the actual physical system; real-time analysis can be conducted on these virtual models/DTs for decision support purposes. By using the DT models for solving optimization problems, we can examine the optimal solutions ahead of time, and without enacting the decision variables directly in the real systems, the feasibility of the systems can be ensured -- preventing the potential unwanted events for the network. This unfolds the opportunity to use approximated linear power flow models for solving OPFs -- defined in equation \eqref{lin_pf}, instead of using computationally complex nonlinear network models as described by equation \eqref{nl_pf}. {\color{black} Upon implementing the optimal decision variables in the DT models (e.g., optimal DER generations -- obtained from the approximated models), the corresponding nonlinear solutions of the state variables of the system, such as nodal voltages, branch power flows etc., of the system can be attained.

Mathematically, the optimization variable $X$ for the OPF problems is composed of state variables, $w$, and the optimization decision variables, $z$, i.e., $X = [w, z]^T$. For example, in the \textbf{(C1)} OPF problem, $w = [P^{pp}_{ij}, Q^{pp}_{ij}, v^p_j, l^{pq}_{ij}]$, and $z = q^p_{D,j}$ $ \forall j\in \mathcal{N}$ \& $\forall \{ij\} \in \mathcal{E}$. Similarly for approximated \textbf{(C4)} OPF problem, $w = [P^{pp}_{ij}, Q^{pp}_{ij}, v^p_j]$, and $z = p^p_{D,j}$ $ \forall j\in \mathcal{N}$ \& $\forall \{ij\} \in \mathcal{E}$. We use DT models for the simulation within optimization as such:
after solving the approximated OPF problem (equation \eqref{DT1}), the optimal decision variables, $z_{Lin}^*$ are implemented in the DT model -- defined by the function $\Phi_{DT}$, and then the nonlinear solutions are extracted. Thus, using \eqref{DT2}, the corresponding nonlinear state variables of the system, $w_{Nl}^*$ are obtained upon effectuating the linear optimizer $z_{Lin}^*$ in the system; the optimizer $X^*$ attained from this simulation within optimization approach is shown in equation \eqref{x_star}. Any nonlinear power-flow model of the system can be used to represent the functional form for the DT.}

\vspace{-0.3cm}
\begin{small}
\begin{equation}\label{DT1}
X_{Lin}^*=[w_{Lin}^*,z_{Lin}^*] := \underset{X_{Lin} \in \mathcal{S}}{\text{argmin}} {\hspace{0.1cm} f(X_{Lin})}
\end{equation}
\vspace{-0.3cm}
\begin{equation}\label{DT2}
w_{Nl}^* = \Phi_{DT} (z_{Lin}^*)
\end{equation}
\vspace{-0.5cm}
\begin{equation}\label{x_star}
X^* = [w^*_{Nl}, z_{Lin}^*]^T
\end{equation}
\end{small}
\vspace{-0.4cm}

Similarly, in the distributed decision-making process, the use of DTs for solving the OPF problems helps to build a routine that can leverage the computational speed and scalability of the approximated optimization problems with linear power flow models. Each distributed agent has the DT models ($\Phi_{DT}$) of the corresponding local networks. The agents first solve the local approximated sub-problems (similar to \eqref{DT1}), and then extract the nonlinear solutions by using \eqref{DT2} for the local system. During the macro-iteration, i.e., communication among these distributed agents, appropriate variables from nonlinear solution, $w_{Nl}^*$ are exchanged, instead of using the solutions from linear models, i.e., $w_{Lin}^*$.
Specifically, alternate to the direct exchange of linear solutions, the nonlinear solutions from DT simulations are shared among the distributed agents. This effectively helps to ensure the feasibility of the optimal solution of the overall problem at the end of the macro-iterations, even though the computationally simple and scalable approximated models are being solved. {\color{black} The block diagram of the proposed distributed algorithms with simulation within optimization using DT models is shown in Fig. \ref{DTC}}

\subsection{Co-simulation Platform}
A realistic evaluation of distributed algorithms require modeling the associated communication system. To this end, we developed a multi-agent cyber-physical co-simulation package using HELICS (Heirarchical Engine for Large-scale Infrastructure Co-Simulation), a generic co-simulation platform  \cite{palmintier2017design}. Our cyber-physical co-simulation package is composed of the following software: Python 3, GridLAB-D \cite{Chassin:2008fq}, NS-3 \cite{ns3}, which are used to simulate the computing agents, the power distribution networks, the communication networks, respectively. HELICS \cite{palmintier2017design} is then employed to synchronize the simulations and facilitates message passing among the three simulation software (see Fig.~\ref{cosim}). Each area has a controller agent shown at the top which has a single connection, via HELICS, to the communication network simulated in NS-3. In this way, the controllers are able to communicate with each other and with the device agents in the local area. The device agents  each have a connection, via HELICS, to the physical device simulated in GridLAB-D and to the communication network. The purpose of the device agent is twofold. It formats the raw data it receives from GridLAB-D so the control agent can understand it and it parses commands received from the control agent to the GridLAB-D. The control agents are capable of running their own GridLAB-D simulations using local area models as digital twins (DT). Each area employs the respective DT to obtain the actual power flow variables the system observes upon implementing the approximated OPF control variables.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.35\textwidth]{figs/Co-Sim_Diagram.eps}

    \caption{Co-Simulation Structure for Cyber-Physical System}
    \label{cosim}
\end{figure}

\vspace{-0.2cm}
\section{Results and Discussions}
\vspace{-0.2cm}
In this section, a detailed evaluation of the proposed D-OPF algorithm for a 3-phase unbalanced distribution system is presented. First, we have validated our distributed OPF approach for nonlinear optimization problems by comparing the solution with the central formulation. Then, we have simulated and verified the simulation within optimization approach that has been proposed in this paper with cyber-physical co-simulation platform and uses approximated OPF problems. Lastly, a time-series simulation has been performed to showcase the effect of network optimization using the proposed method, and different stressed communication network topologies have been tested to demonstrate the efficacy and robustness of the proposed distributed OPF algorithm.
\begin{figure}[t]
    \centering
    \includegraphics[width=0.43\textwidth]{figs/IEEE123_30DER.eps}
        \vspace{-0.2cm}
    \caption{Sectioned IEEE 123 Bus System}
    \label{123bus}
\end{figure}

\subsection{Simulated System}
To test and validate our proposed method, we have simulated two network objective problems with three-phase unbalanced IEEE-123 bus test system that has been assumed to be composed of 4 areas. Three different scenarios of DERs are simulated in this paper to test the proposed method -- (i) 10 three-phase DERs (at node 15, 23, 30, 37, 49, 50, 51, 67, 78, 107) with capacity of 60 kVA per phase, (ii) 30 three-phase DERs with capacity of 48 kVA per phase, and (iii) 30 three-phase DERs with capacity of 60 kVA per phase (Fig. \ref{123bus}). Please note, scenario (ii) is simulated for loss minimization problem to provide a more challenging problem for the algorithm in terms of required communication properties. Further, for the time series simulation, the loads are also varied randomly between 78\% to 98\% of the nominal values for 60 minutes. In addition to that, as a digital twin function, $\Phi_{DT}$, GridLab-D power flow simulator has been used in this study.

\begin{figure*}[t]
\centering
\begin{subfigure}{.32\textwidth}
  \centering  \includegraphics[width=1.0\linewidth]{figs/R_max.eps}
  \vspace{-0.6cm}
  \caption{Convergence at Boundary}
  \label{residual_cmp}
\end{subfigure}
\begin{subfigure}{.32\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{figs/V_ABC_30DER.eps}
  \vspace{-0.6cm}
  \caption{Node Voltages: Loss Minimization}
  \label{volt_lm}
\end{subfigure}
\begin{subfigure}{.32\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{figs/V_ABC_90DER.eps}
  \vspace{-0.6cm}
  \caption{Node Voltages: DER Maximization}
  \label{volt_pvc}
\end{subfigure}
\vspace{-0.2cm}
\caption{D-OPF Results for Loss Minimization and DER Maximization Objectives}
\label{D_nlp}
\end{figure*}

\subsection{Validation of Decomposition Method}
To validate the proposed D-OPF method for unbalanced power distribution systems, we have simulated two optimization problems,\textbf{ (D1)} \& \textbf{(D3)} for DER scenarios (i) \& (iii), respectively, and compared with the solutions of \textbf{(C1)} \& \textbf{(C3)} (see Table \ref{NL_val}). It is evident that the optimal solutions of D-OPF match with the corresponding Central OPF (C-OPF) solutions for both objectives. For example, the objective value for the DER maximization problem, i.e., total DER generation is 5.08 MW \& 5.18 MW, when solved using D-OPF \& C-OPF method, respectively. However, D-OPF takes significantly less time to solve, compared to the C-OPF method -- takes almost 10 times or more for some cases. From Fig. \ref{residual_cmp}, it can be seen that it only takes 4 macro-iterations to reach the global solution. In addition to that, the nodal voltages of D-OPF and C-OPF solutions are compared (Fig. \ref{volt_lm}-\ref{volt_pvc}). Later the nodal voltages of D-OPFs are compared with the GridLab-D/OpenDSS model -- upon enacting the decision variables from D-OPF, and the voltage difference is in the order of $10^{-4}$ or less (Fig. \ref{val_nlp}). These comparisons validate the feasibility and the optimality of the proposed D-OPF solutions for both minimization and maximization objectives.

\begin{figure}[h]
\centering
\begin{subfigure}{.24\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{figs/delV_validation_30DER.eps}

  \caption{Loss Minimization}
  \label{val_lm}
\end{subfigure}
\begin{subfigure}{.24\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{figs/delV_validation_90DER.eps}

  \caption{DER Maximization}
  \label{val_pvc}
\end{subfigure}
\caption{Validation for Decomposition Approach}
\label{val_nlp}
\end{figure}

\begin{small}
\begin{table}

    \centering

    \caption{Validation of D-OPF method}

    \setlength{\tabcolsep}{3pt}
    \begin{tabular}{c  c c  c c }
    \toprule
    \multirow{2}{*}{\textbf{OPF Problem}}  & \multicolumn{2}{c}{\textbf{Objective}} & \multicolumn{2}{c}{\textbf{Time (s)}} \\
    & D-OPF & C-OPF & D-OPF & C-OPF \\
    \midrule
    Loss Min.: Scenario (i)  & 43.976 kW & 43.975 kW & 20 & 520 \\

    DER Max.: Scenario (iii)    & 5.08 MW & 5.18 MW & 20 & 200 \\
    \bottomrule

    \end{tabular}
    \label{NL_val}
\end{table}
\end{small}

\subsection{Co-simulation \& Use of DT}
\vspace{-0.1cm}
Next, to showcase the robustness and the practicability of the proposed algorithm, the D-OPF method is implemented on a cyber-physical co-simulation platform. For this case, distributed problems \textbf{(D2) \& (D4)} are solved to leverage the simplicity of the linear network model for edge computing use-cases to solve the OPFs. In this paper, the DTs of the local areas are modeled using GridLab-D power flow simulator. For this paper, it is assumed that the distributed agents/local controllers communicate with neighbor every 2-seconds. And also, after the convergence is attained, it takes around $0.2$ seconds to dispatch the decision variables in the system.

\begin{small}
\begin{table}

    \centering
    \vspace{-0.15cm}
    \caption{Loss minimization Result Comparison}
    \vspace{-0.1cm}
    \label{loss_min10}
    \setlength{\tabcolsep}{3pt}
    \begin{tabular}{c c c c c}
    \toprule

    {\textbf{Method}}  & \multicolumn{2}{c}{\textbf{Loss (kW)}}  & \multicolumn{2}{c}{\textbf{Communication-round}}\\
    DER Scenario & Scenario (i) & Scenario (ii) &Scenario (i) &Scenario (ii)\\
    \midrule
    No OPF                     & 63.132 & 53.338 & 0 & 0 \\
    Linear C-OPF               & 44.875 & 26.508 & 0 & 0\\
    Linear D-OPF               & 44.959 & 26.546 & 5 & 5 \\
    Linear D-OPF+DT            & 44.959 & 27.238  & 5 & 5 \\
    Nonlinear C-OPF              & 43.975 & 27.160 & 0 & 0\\

    \bottomrule

    \end{tabular}
    \vspace{-0.6cm}
\end{table}
\end{small}

\subsubsection{Validation}
\begin{comment}
The implementation of the algorithm on the co-simulation platform using the linear approximate models are validated by comparing the resulting objective values.
Table \ref{loss_min10} and Table \ref{der_max}, show the comparisons between D-OPF and C-OPF results for three cases. Table~\ref{loss_min10} shows the loss minimization case for DER scenarios (i) and (iii), where real power generation is fix at $50$ kW and $20$ kW, respectively.
The approximated solution, i.e., the line loss is 44.959 kW, which is similar to the nonlinear solution of 43.975 kW. It takes around 5 communication-rounds among distributed agents to reach the convergence. For the DER scenario (iii) with 30 DERs, where the real power generation is fix at 20 kW, similar

remaining loss minimization studies we used a 30 DER case which provides a more challenging problem for the algorithm. Table~\ref{loss_min30} shows the C-OPF and D-OPF comparisons for the 30 DER case. Table~\ref{der_max} shows the results for DER maximization and uses the same physical system as in Table~\ref{loss_min30}, which has 30 DERs.
Each table has four columns. Column (A) is the objective value which is calculated from the results of the physical layer simulation in GridLAB-D. Column (B) is the time when the value in the previous column was achieved. Column (C) is the time the last dispatch from the MGCs was sent and will always be 0.2 seconds after the time the last MGC determined it was converged. The 0.2 second delay is meant to account for computation time inside the MGC. Column (D) is the number of communication rounds required before the final dispatch.
Between the time of final dispatch and the final solution objective from GridLAB-D there is between a 1.8 and 2.8 second delay despite having almost no communication delay (0.1ms). GridLAB-D is running on a one second timestep so it cannot implement the new dispatch until the next whole second and the results of the new dispatch will not show until at least one whole step later. Thus, a minimum 1.8 second delay is always added to the results.
In Table~\ref{der_max}, the times recorded in column (B) are earlier then the time of final dispatch in column (C). This is because, of the four MGCs, only MGC1 sends a dispatch that is different from the default values. In the D-OPF case, MGC1 dispatches at 12.2 seconds and only required 5 communication rounds. In the case with digital twins, MGC1 dispatches at 8.8 seconds and only required 3 communication rounds. Dispatches from other MGCs occurred later but did not impact the results of the simulation.
\end{comment}
Now we validate the co-simulation platform using the linear approximate models by comparing the resulting objective values. Table \ref{loss_min10} and Table \ref{der_max}, show the comparisons between D-OPF and C-OPF results for different cases. Table~\ref{loss_min10} shows the loss minimization case for DER scenarios (i) and (ii), where real power generations are fixed at $50$ and $20$ kW per phase, respectively. For both DER scenarios, line losses matches with the nonlinear solutions. For example, in scenario (i), the active power loss is $44.959$ kW for linear model with the DT model, and the nonlinear solution is $43.975$ kW. However, for scenario (ii), it can be seen that the use of DT model has be potential to get closer solutions to the nonlinear solution. In Table \ref{loss_min10}, the line loss for scenario (ii) is $26.546$ kW when DT models are not being used, however, if the DT models are used, the objective value is $27.23$ kW -- that is closer to the nonlinear solution $27.16$ kW. The use of DT with approximated model does not necessarily improve the results, but it guarantees the feasibility of the local solutions, before exchanging the boundary variables with the neighbors. Also, it takes around 5 communication-rounds among the agents to reach the consensus at the boundary. Similarly, Table \ref{der_max} validates the algorithm on the co-simulation platform for DER maximization problem that uses DER scenario (iii). When no OPF is performed, there were voltage limit violation in the system. However, if operational limits are considered, the optimal generation is $5.135$, $5.163$, and $5.18$ MW for linear D-OPF, linear D-OPF+DT, and nonlinear C-OPF, respectively. Thus the use of the DT model for DER maximization problem also brings the solution closer to the actual nonlinear solutions. Also, using DT model reduces the number of communication rounds from 7 to 5. Thus, Table \ref{loss_min10} and Table \ref{der_max} validate the co-simulation platform that uses computationally simple linear approximated OPF model with the DT model of the system.

\begin{small}
\begin{table}
    \vspace{-0.1cm}
    \centering
    \caption{DER Maximization Results Comparison}
    \vspace{-0.1cm}
    \label{der_max}
    \setlength{\tabcolsep}{3pt}
    \begin{tabular}{c  c  c}
    \toprule
    \textbf{Method}  & \textbf{Generation (MW)} & \textbf{Communication-round}\\
    DER Scenario & Scenario (iii) & Scenario (iii)\\
    \midrule
    No OPF                     & 5.400 & 0\\
    Linear C-OPF               & 5.135& 0\\
    Linear D-OPF               & 5.135 & 7\\
    Linear D-OPF+ DT           & 5.163 & 5\\
    Nonlinear C-OPF              & 5.180 & 0\\
    \bottomrule

    \end{tabular}
    \vspace{-0.7cm}
\end{table}
\end{small}

\subsubsection{Time-series Simulation}
In Fig.~\ref{dist_lp}, we show how the voltages in the system are impacted by the loss minimization (scenario (ii)) and DER maximization (scenario (iii)) objectives respectively. In the loss-minimization case, we expect that all of the bus voltages will be shifted towards the upper limit. In the DER maximization case we expect that bus voltages that exceed the voltage limit will be reduced to be within the voltage limits. These expectations are confirmed in Fig.~\ref{dist_lp}.

\begin{figure}[t]
\centering
\begin{subfigure}{.24\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{figs/loss_min_dist_30der_lowgen4x.eps}
  \vspace{-0.6cm}
  \caption{Loss Minimization}
  \label{dist_lm}
\end{subfigure}
\begin{subfigure}{.24\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{figs/der_max_dist.eps}

  \caption{DER Maximization}
  \label{dist_pvc}
\end{subfigure}
\caption{Voltage Distributions}
\label{dist_lp}
\end{figure}

\subsubsection{Communication Network Stress Test}
Next we test the robustness of the proposed algorithm, when subjected to poor communication network conditions; it is tested with two different communication network topologies (ideal and ring topology) for loss minimization problem \textbf{(D2)} with scenario (ii) and for DER maximization problem \textbf{(D4)} with scenario (iii). In the ideal topology, each area controller has a direct link with its neighbors and with devices in its own area. In the ring topology all devices and controllers are connected in a single large loop without regard to physical location. Each link between network nodes is a point-to-point link with a negligible delay and has different bandwidths (3kbps, 2kbps, or 1kbps) to show the worst-case performance of the algorithm (see Fig. \ref{comms_lp}). The lines marked with circles indicate results run with the ring topology and lines marked with a star indicate the ideal topology. The blue line shows the results if no communication is allowed between controllers and the controllers have perfect communication links with local devices.
There are several effects of poor communications that we have observed:
\begin{itemize}[leftmargin=*, noitemsep]
    \item  If a controller doesn't receive data from either UA or DA because of data delays, it will falsely assume convergence and dispatch the inverters prematurely. When the delayed data does arrive it will continue iterating again.

    \item Delayed communication between areas does not prevent convergence except in extreme cases.
    \item The quality of the resulting solution as measured after convergence is not impacted by the communication delays.

\end{itemize}

\begin{figure}[t]
\centering
\begin{subfigure}{.24\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{figs/lossmin_n1n2_prec.eps}
  \vspace{-0.6cm}
  \caption{Loss Minimization}
  \label{val_lm}
\end{subfigure}
\begin{subfigure}{.24\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{figs/pvcurt_n1n2.eps}

  \caption{DER Maximization}
  \label{val_pvc}
\end{subfigure}
\caption{Communication Stress Test}
\label{comms_lp}
\end{figure}

\subsubsection{Advantages of Using DT with LP}
The LP model can estimate the nodal voltages with moderate accuracy; however, for corner cases, such as when the system is operating at the voltage limits, the approximation error can lead to a higher optimality gap and even infeasible solutions. In Fig. \ref{fig:DT_corner_lm}, we have shown a snippet of voltage at phase C for a high loading condition; using a DT model yields more accurate nonlinear nodal voltages than LP solutions that approximate higher voltages. If the voltage lower bound constraint is active, an overestimation of the voltage can lead to infeasible space upon implementing the LP solutions to the actual system. However, by using the DT model, we can estimate the nonlinear voltages more accurately and can preemptively tackle any infeasibility. Again for the DER maximization case, when the voltage upper bound is active, the higher approximated voltage can lead us to sub-optimal DER generation solutions. Here, a low loading condition scenario has been simulated, and a snippet of phase A voltages are shown in Fig. \ref{fig:DT_corner_pv}. Upon implementing the LP solution in the OpenDSS, we obtain lower feeder voltages indicating that more DER can be integrated before the actual feeder limit is reached. Thus, LP model results in sub-optimal solutions. However, the LP+DT model does not overestimate the voltages and hence optimized for the maximum generation that can be integrated without causing any voltage violations. We observed that the optimal solutions for the two models differ by $300$ kW. This case further demonstrates the usability of the DT model in accurately solving the OPF problems without augmenting any complicated nonlinear formulation, but still resulting in more accurate optimal solutions.

\begin{figure}[t]
\centering
\begin{subfigure}{.24\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{figs/Loss_min_DT_Corner.eps}

  \caption{Loss Minimization}
  \label{fig:DT_corner_lm}
\end{subfigure}
\begin{subfigure}{.24\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{figs/PV_max_DT_Corner.eps}

  \caption{DER Maximization}
  \label{fig:DT_corner_pv}
\end{subfigure}
\caption{Voltage Approximation Error}
\label{fig:DT_corner}
\end{figure}

\vspace{-0.2cm}
\section{Conclusions}
\vspace{-0.05cm}
In this paper, a novel distributed optimal power flow method for radial three-phase, unbalanced power distribution system has been proposed. The proposed approach uses {\textit{simulation within optimization}} to improve the scalability and reduce solution time for each local subproblem. We employ detailed DT models of the local network in conjugation with a linear-approximated OPF to ensure a power-flow feasible and near-optimal solution for multiple OPF problems for power distribution system. Additionally, the proposed distributed coordination method achieves a converged network-level optimal solution using significantly fewer communication rounds. We also evaluated the performance of the proposed algorithm using a cyber-power co-simulation environment with various communication network parameters, to validate the robustness of the proposed distributed OPF algorithm under stressed communication. For various cases, it is shown that the proposed method successfully converges to the network-level optimal solution within reasonable time, even when the communication network is facing significant stress, such as, communication delays and a low bandwidth communication infrastructure.

\vspace{-0.3cm}
\bibliographystyle{ieeetr}


\section{Introduction}
Time evolution of quantum systems is one important class of problems where quantum computers are expected to have an exponential speed up relative to classical computers. One class of quantum field theories of particular interest are gauge theories, which are used to describe the physics of elementary particles, effective theories in nuclear physics, and condensed matter systems. Typically, gauge theories are studied using the Lagrangian formalism due to the fact that symmetries are manifest from this point of view. In order to study gauge theories using quantum computers however, one instead has to work with the Hamiltonian formalism. As with studying any physics problem computationally, it is of fundamental importance that the particular formulation chosen allows for efficient simulation using quantum computers.

One property of gauge theories that makes formulating the Hamiltonian non-trivial is that the Hilbert space of a given Hamiltonian spans multiple charge sectors. This fact can make quantum simulations of such a theory difficult for two reasons. Firstly, any noise from the environment during a quantum computation can cause a charge-violating transition and lead to incorrect results. Additionally, the naive basis of states for the Hilbert space is generally over-complete, requiring more qubits and gate operations than strictly necessary. Formulations have been developed that enforce gauge invariance including these unphysical states \cite{Halimeh_2020, Halimeh:2020ecg, Halimeh:2021lnv, Halimeh:2021vzf, Lamm_Scott_Yukari_2020, Tran_2021, PhysRevLett.109.125302, Banerjee_2012, https://doi.org/10.48550/arxiv.2012.08620, PhysRevX.3.041018, PhysRevA.90.042305, Stannigel_2014, Stryker_2019}, as well as removing them $\textit{a priori}$ \cite{Kaplan:2018vnj,Unmuth-Yockey:2018xak,Haase:2020kaj,Bender:2020ztu,PhysRevD.19.619,Bauer:2021gek}. Reviews of different approaches, including both analog and digital methods, can be found in Refs.~\cite{Wiese_2013, Zohar_2015, Dalmonte_2016, Aidelsburger_2021, Ba_uls_2020, Zohar_2021, Klco_2022, Bauer:2022hpo}.

In this work, we consider a formulation of pure U(1) gauge theory in 2+1 dimensions where unphysical states have been removed \textit{a priori} \cite{Bauer:2021gek}. We show that in this formulation, Suzuki-Trotter time evolution leads to an exponential volume scaling in the gate count, which is prohibitively expensive for realistic lattice sizes. We then show that this exponential volume scaling can be reduced to polynomial by performing a specific change of operator basis, and conclude by applying this change of basis to a specific lattice size.

\section{Gauge invariance and Gauss' Law}

The theory we consider is a pure U(1) gauge theory in 2+1 dimensions.  The Hamiltonian for the classical continuum theory is given by an integral over electric and magnetic fields
\begin{align}
    H = \int d^2x\, \left(\vec{E}(x)^2 + B(x)^2\right),
\end{align}
where the electric and magnetic fields are subject to their respective electric and magnetic Gauss' law, given by $\vec{\nabla} \cdot \vec{E}(x) = \rho(x)$ and $\vec{\nabla} \cdot B(x) = 0$, where $\rho(x)$ is the charge density. The quantum version of this theory is formulated on a finite lattice, and one has different options for what basis of states to use to represent the Hilbert space. Unless extra care is taken, the states in the Hilbert space of a given formulation do not automatically obey these Gauss law constraints. One example of such a formulation is the Kogut-Susskind formulation \cite{PhysRevD.11.395}.

As mentioned in the introduction, this property leads to complications with regards to quantum simulations. Firstly, because this naive basis of states is over-complete due to the gauge orbit degeneracy, more quantum resources are required for the calculation than strictly necessary. Furthermore, because the Hilbert space spans multiple charge sectors, quantum noise in a quantum computation can cause unphysical transitions between them. Figure~\ref{fig:hilbert_space} shows a schematic visualization of a Hilbert space. When performing a quantum simulation, one first prepares an initial state in a specific charge sector, which can be thought of as a specific arrangement of static charges. In a noise-free quantum simulation, the state would evolve in time while remaining in the same initial charge sector. However, the presence of noise could cause an unphysical transition from the initial charge sector to a different charge sector.
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{hilbert_space.png}
    \caption{A schematic visualization of a Hilbert space. Different values of $Q$ represent different gauge sectors, each with a different total charge. A gauge orbit is the set of physically equivalent gauge fields in a particular charge sector.}
    \label{fig:hilbert_space}
\end{figure}
One solution to these problems is to work with a formulation of the U(1) theory whose Hilbert space is restricted to only physical states; in other words, a formulation that is free of gauge redundancies. Such a formulation would not allow charge-violating transitions, and could also require less quantum resources due to the reduced Hilbert space size. Much work has been done studying gauge-redundancy free formulations of U(1) gauge theories \cite{Kaplan:2018vnj,Unmuth-Yockey:2018xak,Haase:2020kaj,Bender:2020ztu,PhysRevD.19.619,Bauer:2021gek}, which in general have become known as ``dual basis'' formulations.

The dual basis is formulated in terms of electric ``rotor'' operators and the usual magnetic field operators, both defined on the plaquettes of the lattice. The rotors are defined such that the transverse component of the electric field is given by $\vec{E}_T = \vec{\nabla} \times R$. By writing the longitudinal component of the electric field as $\vec{E}_L = \vec{\nabla} \cdot \rho$, one can see that the electric field $\vec{E} = \vec{E}_T + \vec{E}_L$ automatically satisfies electric Gauss' law. In these proceedings, we work in the charge-zero sector with $\rho(x)=0$. One important property of the rotor and magnetic field operators is that they are conjugate operators, satisfying $[B_p, R_{p'}] = i \delta_{pp'}$. This relation implies that the electric basis, where the rotor operators are diagonal, is related to the magnetic basis, where the magnetic field operators are diagonal, by the usual Fourier transformation. The Hamiltonian is given by a sum of the electric Hamiltonian $H_E$ and magnetic Hamiltonian $H_B$. For a periodic lattice with $N_x(N_y)$ sites in the $\hat{x}(\hat{y})$ directions, $H_E$ and $H_B$ are given by
\begin{equation}
    H_E = 2g^2 \sum_{p=1}^{N_x N_y} \left( \vec{\nabla} \times R_p \right)^2, \qquad H_B = \frac{1}{g^2} \left( N_x N_y - \sum_{p=1}^{N_x N_y} \cos B_p \right),
    \label{eq:HE_and_HB_pre_gauge_fix}
\end{equation}
where $g$ is the gauge coupling. If working in a sector with non-zero charge, the electric Hamiltonian will be modified appropriately. Note that we use lattice units throughout this work and therefore set the lattice spacing to one.

While the local Gauss' law constraints are automatically satisfied, there is one remaining global constraint that is not automatically satisfied. This can be seen by counting the degrees of freedom in this formulation. Working in the magnetic basis, an eigenstate of the Hamiltonian in Eq.~\eqref{eq:HE_and_HB_pre_gauge_fix} is completely specified by $N_x N_y$ quantum numbers, one for each $B_p$ operator. However, it was shown in Ref.~\cite{Kaplan:2018vnj} that a gauge-fixed U(1) gauge theory in two spatial dimensions requires only $N_x N_y-1$ quantum numbers, which implies there is one redundant degree of freedom in the current formulation. This redundancy corresponds to the constraint that the product of plaquettes around a close surface must be the identity, which is the lattice version of the integral form of magnetic Gauss' law $\sum_{p=1}^{N_x N_y} B_p = 0$. In order to constrain the Hilbert space to contain only physical states, we choose to remove the redundancy by writing $B_{N_x N_y} = - \sum_{p=1}^{N_x N_y-1} B_p$, and setting the associated rotor to zero, \textit{i.e.} $R_{N_x N_y}=0$. In this way, the system now contains $N_p \equiv N_x N_y - 1$ independent plaquettes. While the new electric Hamiltonian has the same form as before, the magnetic Hamiltonian becomes highly non-local. In particular, the gauge-redundancy free magnetic Hamiltonian is given by (up to an overall constant)
\begin{equation}
    H_B = -\frac{1}{g^2} \left[ \sum_{p=1}^{N_p} \cos B_p + \cos \left(  \sum_{p=1}^{N_p} B_p \right) \right],
\end{equation}
where the term that is the cosine of the sum of all plaquettes $\cos(\sum_{p=1}^{N_p} B_p)$ is the source of the non-locality. We will show in Sec.~\ref{ssec:gateCount} that this cosine is the source of the exponential volume scaling in the gate count.

\section{Quantum simulation}
\label{sec:quantumSimulation}
In this section we describe the strategy used to perform the quantum simulation. We start by outlining the time evolution strategy used, and then review the digitization scheme of the $R_p$ and $B_p$ operators given in Ref.~\cite{Bauer:2021gek}. We conclude this section by studying how the gate count scales with the volume $N_p$. Note that we only provide the details of the digitization required to perform a gate count study. Further details can be found in Ref.~\cite{Bauer:2021gek}.

\subsection{Time evolution strategy}
\label{ssec:timeEvolutionStrategy}
To study the gate count, it is sufficient to consider first order Suzuki-Trotter methods. In particular, the time evolution operator is approximated as
\begin{equation}
    U(t) = \left( e^{-i \delta t H_E} e^{-i \delta t H_B} \right)^{N_\text{steps}} + \mathcal{O}(\delta t),
\end{equation}
where $\delta t = t/N_\text{steps}$. We work in the magnetic basis, where $H_B$ is diagonal. Because rotors and magnetic field operators satisfy $[B_p, R_{p'}]=i \delta_{pp'}$, the electric and magnetic basis are related by a Fourier transform. Using this fact, the procedure we choose to perform a single Trotter step is to first implement the diagonal operator $e^{-i \delta t H_B}$, rotate to the electric basis using the quantum Fourier transform (QFT), implement diagonal operator $e^{-i \delta t H_E}$, and then rotate back to the magnetic basis using the inverse QFT. This procedure involves two classes of operators, namely the QFT, and diagonal operators. We perform a detailed gate count study of each step in Sec.~\ref{ssec:gateCount}.

\subsection{Digitization scheme}
\label{ssec:digitizationScheme}

The digitization scheme developed in Ref.~\cite{Bauer:2021gek} represents the rotors and magnetic field operators by diagonal matrices with evenly spaced eigenvalues. Each lattice site is represented by $n_q$ qubits, and so each operator is sampled $2^{n_q}$ times. Using the gate representation of such operators given in Ref.~\cite{PhysRevA.99.052335}, the $B$ and $R$ operators are written as\footnote{Note that Eq.~\eqref{eq:opRep} differs from Ref.~\cite{PhysRevA.99.052335} by a sign flip, a constant shift, and an overall multiplicative factor. This is due to a difference in the ordering of the eigenvalues, and the fact that we sample the operators asymmetrically about the zero eigenvalue.}
\begin{align}
    B = -\frac{b_\text{max}}{2^{n_q}} \left(\mathds{1} + \sum_{j=0}^{n_q-1} 2^j \sigma^z_j \right), \qquad R = -\frac{r_\text{max}}{2^{n_q}} \left(\mathds{1}+\sum_{j=0}^{n_q-1} 2^j \sigma^z_j\right),
    \label{eq:opRep}
\end{align}
where $\sigma^z_j$ is a Pauli-Z matrix acting on the $j^\text{th}$ qubit and $\mathds{1}$ is the $2^{n_q}\times 2^{n_q}$ identity matrix. The $B$ and $R$ operators are sampled in the ranges $[-b_\text{max}, b_\text{max}]$ and $[-r_\text{max}, r_\text{max}]$, respectively. Furthermore, the constant $b_\text{max}$ is chosen to minimize the digitization errors for a given value of the gauge coupling $g$, and the constant $r_\text{max}$ is a known function of $b_\text{max}$. Choosing a good value of $b_\text{max}$ can be done via a known analytic function and does not require performing a scan~\cite{Bauer:2021gek}. It was shown in Ref.~\cite{Bauer:2021gek} that sampling each operator a number of times corresponding to $n_q=3$ reproduces the low-lying spectrum to per-mille level accuracy. This fact will be important for the gate count study that follows.

\subsection{Gate count study}
\label{ssec:gateCount}

We are now in a position to study the gate count scaling for a single Trotter step. As previously mentioned, while $N_p$ will have to be taken large to simulate realistic lattice volumes, per-mill level accuracy can be achieved while keeping the number of qubits per lattice site $n_q$ small.

Each component of the Hamiltonian is implemented on a quantum computer separately, and we therefore consider the gate count of each component separately. Changing the basis requires performing a quantum Fourier transform. This can be done at each lattice site using a Fast Fourier Transform, which requires $O(n_q^2)$ gates per site \cite{nielsen_chuang_2010}. Therefore, the total gate count for changing the basis and back again for a single step is $O(n_q^2 N_p)$. Moving on to the electric Hamiltonian, because the lattice version of the curl is composed of differences, the electric Hamiltonian in Eq.~\eqref{eq:HE_and_HB_pre_gauge_fix} is composed of biliear terms of the form $R_p^2$ or $R_{p} R_{p'}$. Using the gate representation for $R$ in Eq.~\eqref{eq:opRep}, each binlinear term can be written schematically as $R^2 \sim \sum_{i,j=0}^{n_q-1} 2^{i+j} \sigma_i^z \sigma_j^z$, which is a sum of $n_q^2$ Pauli strings of length two. This implies that each bilinear term requires $\mathcal{O}(n_q^2)$ gates to exponentiate. The number of bilinear terms scales as $\mathcal{O}(N_p)$~\cite{Haase:2020kaj}, and so $e^{-i \delta t H_E}$ requires $\mathcal{O}(n_q^2 N_p)$ gates to implement.

The number of gates required to implement the complex exponential of each $\cos(B)$ term in the magnetic Hamiltonian can be found by writing $B$ in terms of Pauli-Z operators and then expanding the cosine using the Taylor series. Doing so we find
\begin{equation}
    \cos(B) = \sum_{m=0}^{\infty} \frac{(-1)^m}{(2m)!} (B)^{2m} = \sum_{i_0=0}^1 \sum_{i_1=0}^1 \dots \sum_{i_{n_q-1}=0}^1 a_{i_0 i_1\dots i_{n_q-1}} (\sigma_0^z)^{i_0} (\sigma_1^z)^{i_1} \dots (\sigma_{n_q-1}^z)^{i_{n_q-1}},
    \label{eq:cosBPauliZ}
\end{equation}
where $a_{i_0 i_1\dots i_{n_q-1}}$ is the coefficient a given Pauli string and is in general non-zero. Equation~\eqref{eq:cosBPauliZ} indicates that each $\cos(B)$ term is a sum $2^{n_q}$ Pauli strings of length $n_q$, and therefore requires $\mathcal{O}(2^{n_q})$ gates to exponentiate. While this is exponential in $n_q$, recall that $n_q$ can be kept small and still reproduce the low-lying spectrum of the theory. More importantly, because there are $N_p$ such terms in $H_B$, exponentiation of all single cosine terms requires $\mathcal{O}(2^{n_q} N_p)$ gates, which is linear in the volume. Furthermore, because each $\cos(B)$ term acts on a different set of $n_q$ qubits, these terms can implemented in parallel, removing the runtime volume dependence.

Following the same logic, using the Taylor series, the $\cos(\sum_{p=1}^{N_p} B_p)$ term is composed of a sum of $2^{n_q N_p}$ Pauli strings of length $n_q N_p$. The number of gates required to exponentiate this term is therefore $\mathcal{O}(2^{n_q N_p})$, which is exponential in the volume. For a modest sized $20 \times 20$ lattice, this term will require $\mathcal{O}(2^{400 n_q})$ gates per Trotter step.

\section{Operator basis change to break exponential volume scaling}
In the previous section we showed that the cosine of a single $B$ operator requires $\mathcal{O}(2^{n_q})$ gates to exponentiate, and that the cosine of the sum of every $B$ operator requires $\mathcal{O}(2^{n_q N_p})$ gates. This discussion generalizes such that the cosine of a sum of $M$ magnetic field operators requires $\mathcal{O}(2^{n_q M})$ gates to exponentiate. Therefore, an equivalent theory where any given cosine in the magnetic Hamiltonian contains a sum of no more than $M \sim \log_2 N_p$ magnetic field operators would require only a polynomial number of gates. In this section, we explain how this can be done using an operator change of basis put forth in Ref.~\cite{Grabowska:2022uos}. A detailed proof that such an operator basis change exists for any lattice volume, as well as an explicit efficient construction method, can be found in Ref.~\cite{Grabowska:2022uos}.

The strategy is to perform a change of operator basis
\begin{equation}
    B_p \to \mathcal{W}_{pp'} B_{p'}, \quad R_p \to \mathcal{W}_{pp'} R_{p'},
\end{equation}
where $\mathcal{W}$ is an orthogonal matrix. We will refer to this new operator basis as the weaved basis. Because $\mathcal{W}$ is orthogonal, the electric and magnetic basis are still related by a Fourier transform, and the time evolution strategy in the weaved basis is the same as outlined in Sec.~\ref{ssec:timeEvolutionStrategy}. The matrix $\mathcal{W}$ is chosen to be block diagonal with $N_s$ sub-blocks, written as $\mathcal{W} = \text{diag}(W_{d_{(1)}},  W_{d_{(2)}}, \dots,  W_{d_{(N_s)}})$. The matrices on the diagonal $W_d$ are so-called weaved matrices of dimension $d$. Because each $W_d$ is chosen to be orthogonal, the resulting matrix $\mathcal{W}$ is also orthogonal.

Note that any basis change that reduces the number of magnetic field operators in the $\cos(\sum_{p=1}^{N_p} B_p)$ term will increase the number of magnetic field operators appearing in the local $\cos(B_p)$ terms. These effects must be balanced in such a way as to ensure no term contains more than $\mathcal{O}(\log_2 N_p)$ magnetic field operators. The properties of $\mathcal{W}$ that result in this are \cite{Grabowska:2022uos}
\begin{enumerate}
    \item $\mathcal{W}$ is block diagonal with $N_s \sim \log_2 N_p$ sub-blocks
    \item Each sub-block $W_d$ has dimension $d \sim N_p / \log_2 N_p$
    \item First column of any $W_d$ has entries all equal to $1/\sqrt{d}$
    \item Each row of $W_d$ has no more than $\lceil \log_2 d \rceil + 1$ non-zero entries.
\end{enumerate}
The first three properties reduce the number of terms that appear in the $\cos(\sum_{p=1}^{N_p} B_p)$ term from $N_p$ to $\mathcal{O}(\log_2 N_p)$ and therefore this term now requires $\mathcal{O}(N_p^{n_q})$ gates to implement. The last property ensures that each of the previously local cosine terms contain no more than $\mathcal{O}(\log_2(N_p/\log_2 N_p))$ operators, and exponentiation of one of these terms requires $\mathcal{O}(N_p/\log_2 N_p)^{n_q}$ gates. The number of gates required to implement the exponential of $H_B$ in the weaved basis is therefore
\begin{equation}
    {\rm Gates}(e^{-i \delta t H_B}) = \mathcal{O}(N_p^{n_q}) + \mathcal{O}(N_p(N_p/\log_2 N_p)^{n_q}).
\end{equation}
Note that the number of gates is now polynomial in the volume $N_p$. While the change of operator basis will introduce more terms in the electric Hamiltonian, because $H_E$ is bilinear, the maximum number of terms is $\mathcal{O}(N_p^2)$. The worse case scaling to exponentiate $H_E$ in the weaved basis is therefore quadratic in the volume, given by $\mathcal{O}(n_q^2 N_p^2)$.

We now work through an example of how to choose $\mathcal{W}$ for $N_p=16$. We first provide the matrix and then walk through each of the four conditions required. The matrix is given by
\begin{equation}
    \mathcal{W} = \begin{pmatrix}
        W_4 & 0 & 0 & 0 \\
        0 & W_4 & 0 & 0 \\
        0 & 0 & W_4 & 0 \\
        0 & 0 & 0 & W_4
    \end{pmatrix},
    \qquad W_4 =  \begin{pmatrix}
        \frac{1}{2} & -\frac{1}{\sqrt{2}} & -\frac{1}{2} & 0 \\
        \frac{1}{2} & \frac{1}{\sqrt{2}} & -\frac{1}{2} & 0 \\
        \frac{1}{2} & 0 & \frac{1}{2} & -\frac{1}{\sqrt{2}} \\
        \frac{1}{2} & 0 & \frac{1}{2} & \frac{1}{\sqrt{2}}
    \end{pmatrix}.
    \label{eq:weavedNp16}
\end{equation}
We see that the number of sub-blocks is $N_s = 4$, the dimension of each sub-block is $d=4$, and the first column of each sub-block is equal to $1/2$.  Additionally, each row of a given sub-block has no more than three non-zero entries. The weaved matrix in Eq.~\eqref{eq:weavedNp16} therefore satisfies the required properties, and the number of gates required to exponentiate the magnetic Hamiltonian in the weaved basis for $N_p=16$ is $\mathcal{O}(16^{n_q})$.

This exponential reduction in gate count can also be understood visually. Figure~\ref{fig:visual} shows how the magnetic field operators are coupled in both the original and weaved basis, where the operators $\mathcal{O}_p$ represent the magnetic field operators, and all operators within the same box appear in a single cosine term in the Hamiltonian. The left image in Fig.~\ref{fig:visual} shows the connectivity in the original basis. The blue dashed boxes represent the $\cos(B)$ terms, and the red solid box represents the $\cos(\sum_{p=1}^{N_p} B_p)$ term. After performing the change of operator basis with the matrix $\mathcal{W}$ in Eq.\eqref{eq:weavedNp16}, the connectivity is given in the right image of Fig.~\ref{fig:visual}. The red solid box indicates that the number of magnetic field operators in the previously maximally coupled term has been reduced to four. The blue dashed and dotted boxes indicate that the previously local terms now each contain three magnetic field operators. The exponential reduction of the gate count in the weaved basis can be understood visually by the reduction in the maximum number of operators appearing in any given box from $N_p$ to $\log_2 N_p$. The number of gates required using \textit{e.g.} $n_q=2$ in the original basis is $\mathcal{O}(10^9)$ while in the weaved basis is only $\mathcal{O}(10^2)$.

\begin{figure}
\begin{minipage}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{connectivity_before.pdf}
\end{minipage}
\hfill\vline\hfill
\begin{minipage}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{connectivity_after.pdf}
\end{minipage}
\caption{Visual representation of the cosine terms in $H_B$ for $N_p=16$. Each operator $\mathcal{O}_p$ represents a magnetic field operator $B_p$, where the left(right) figure shows the original(weaved) basis. Operators inside the same box also appear in a single cosine of $H_B$. Boxes with different colors or line-styles correspond to different cosine terms. Left: The blue dashed rectangles with a single operator correspond to the local $\cos(B_p)$ terms in the original basis, and the red solid rectangle that contains every operator corresponds to the $\cos(\sum_{p=1}^{N_p} B_p)$ term in the original basis. Right: The blue dashed and dotted rectangles show, in the weaved basis, the increased connectivity of the previously local $\cos(B_p)$ terms. The red solid rectangle shows how the previously maximal connectivity of the $\cos(\sum_{p=1}^{N_p} B_p)$ term is reduced in the weaved basis.}
\label{fig:visual}
\end{figure}

\section{Conclusions}
Performing lattice gauge theory calculations on near-term quantum computers will require formulations that are robust to noise. It is also necessary, however, that the number of gates required to implement a given formulation does not scale exponentially with the system size. In these proceedings, we have shown that Suzuki-Trotter time evolution of a formulation of 2+1 U(1) lattice gauge theory that is fixed to a chosen charge sector, and therefore more robust to charge-violating transitions induced by quantum noise, requires a number of gates that scales exponentially with the volume. This exponential volume scaling will make even a single Trotter step for realistic lattice volumes prohibitively expensive. We demonstrated however, that this exponential volume scaling could be broken by performing an operator basis change to the weaved basis. This exponential reduction in the gate count will allow the possibility of simulating formulations of U(1) lattice gauge theories that are more robust to charge-violating transitions resulting from quantum noise.

Furthermore, while we applied the change of operator basis to a specific formulation of U(1) gauge theory in 2+1 dimensions, it can be applied to any bosonic theory. Moving forward, it will be interesting to see how it can be applied to formulations of non-abelian gauge theories.

\acknowledgments
The authors would like to thank Hank Lamm, Natalie Klco, Jesse Stryker and Tobias Osborne for useful discussions during the preparation of this work. CFK is supported by the U.S. Department of Energy (DOE) Computational Science Graduate Fellowship under award number DE-SC0020347. This work was supported by the DOE, Office of Science under contract DE-AC02-05CH11231, through Quantum Information Science Enabled Discovery
(QuantISED) for High Energy Physics (KA2401032), by the Office of Advanced Scientific Computing Research (ASCR) through the Accelerated Research for Quantum Computing Program, by the U.S. Department of Energy grant DE-FG02-97ER-41014 (Farrell), and the U.S. Department of Energy, Office of Science, Office of Nuclear Physics and \href{https://iqus.uw.edu}{\color{black}}{InQubator for Quantum Simulation (IQuS)} under Award Number DOE (NP) Award DE-SC0020970. This work is also supported, in part, through the \href{https://phys.washington.edu}{\color{black}}{Department of Physics} and \href{https://www.artsci.washington. edu}{\color{black}}{the College of Arts and Sciences}  at the University of Washington. This research benefited from the resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725.