\section{DOPF Development \& Integration}\label{overview}
In \cref{fig-flow-overview}, core features of DOPF development procedure are shown, along with the inputs required for and outputs generated from the successful demonstration of DOPF features. \\
\uline{\textit{Data Inputs and Outputs:}} For the distribution feeder of interest, its circuit model is the most important input for the DOPF development procedure, which allow to run power flows to validate that a viable solution of the DOPF exists. GridLAB-D and OpenDSS are common open source power flow simulators that can be utilized for this purpose. Other inputs include information on the equipment required to be controlled through DOPF. As penetration of DERs in distribution grids is increasing, multi-period optimization is of interest, and hence including time-series profiles can be considered essential for developing DOPF algorithms. From the core features of DOPF, expected outputs to be observed are the proposed changes in the controllable variables (e.g., taps positions, capacitor banks, controllable loads commands) and its impact on the grid states (e.g., node voltages and line currents).\\
\uline{\textit{DOPF Core Features Overview:}}
The core features of the DOPF development consist of 1) processing a grid model, 2) converting a highly non-convex distribution power flow problem into a tractable form, and 3) implementing an algorithm that can produce a reliable solution in terms of solution times and solution repeatability. The explanation of core features development and the rationale of the adopted steps in the procedure are presented next.
\subsection{DOPF Core Features Development}
\uline{\textit{Process Grid Model:}} The grid to be optimized by models needs to be converted to an appropriate vector and matrices format. This allows the DOPF developer to utilize packages/libraries which can improve computational performances (e.g., see Python's SciPy library implementation \cite{bressert2012scipy}). Based on the input circuit model format and the information available in the files, an appropriate topology reading algorithm may be written to develop such functionality as follows. First, an injection vector determining how much load/generation is injected/withdrawn at an appropriate node is developed. Second, incidence matrices mapping controllable equipment to the appropriate node in the injection vector is determined. Third, a verification model to map the injections to the grid states is developed to verify that the model obtained matches the real-world grid states. This step is basically translating the grid information into a form that is similar to power flow models, where admittance translates nodal power injections to voltage calculations. Consider a DOPF algorithm aiming to control load shedding $\mv{p}^{\textrm{fl}}$, capacitor banks $\mv{q}_{\textrm{cb}}$, tap-settings $ \mv{t}_{\textrm{rg}}$, PV active power shed $\mv{p}_{\text{pv}}$ and PV reactive power injections $\mv{q}_{\text{pv}}$, and power feed in from the bulk power system $p_0/q_0$. The injection vector to be controlled is then:
\begin{align}
 \mv{u}=\{\mv{p}^{\textrm{fl}}, \mv{q}_{\textrm{cb}}, \mv{t}_{\textrm{rg}}, \mv{p}_{\text{pv}}, \mv{q}_{\text{pv}}, p_0, q_0\}.
\end{align}
In this case, $\mv{x}$ is the set of state variables (like voltages and resultant current magnitude flows), $\mv{C_x}$ includes the impact of the controllable vector on state variables and $\mv{C_u}$ is the set of incidence matrices that maps the controllable equipment list onto the appropriate location in the network node list. An example of the incidence matrix for the regulator follows. An incidence matrix $\mv{C}_{\mv{r}}$ links regulator location to voltage changes such that $\mv{C}_{\mv{r}}(i,j)=1$ if tap regulator $j$ is upstream to node $i$, otherwise $\mv{C}_{\mv{r}}(i,j)=0$. The grid processing steps attempt to formulate the following equation to verify that the grid model is appropriate:
\begin{align}\label{eq:f}
 f(\mv{u}, \mv{x}, \mv{C_x}, \mv{C_u}) = \mv{0}
\end{align}
Different than the power flow problems, the DOPF requires non-equality constraints to be formulated too, such as voltage limitations, current flow limitations and controllable variables movement limitations. These constraints can be compacted in a similar manner to power flow model\footnote{where $\mv{D_x}$ is the inequality counterpart of $\mv{C_x}$, just as line flow utilizes different admittance quantities than nodal admittance.} as:
\begin{align}\label{eq:g}
 g(\mv{u}, \mv{x}, \mv{D_x}, \mv{C_u}) \leq \mv{0}
\end{align}

\uline{\textit{Obtain Tractable Power Flow Form:}} With the grid model processed and verified, the next step of the DOPF core feature development is to convert the non-convex distribution power flow model, along with the discrete controllable equipment functional space, into a tractable model form. The reason for performing this step is that distribution grids differ greatly in design and in the spirit of standardizing DOPF model development; a tractable power flow form must be made that can achieve a balance between the solution space and computational efforts. As a case of the simplest tractable form, a linear power flow model with continuous decision variables mapping can be picked. To ensure the tractable model form is able to be used in developing the DOPF algorithm, the model should be parameterized in the controllable equipment variables. This allows for validation of the errors on the approximation of the power flow as the change in the controllable equipment is proposed by the DOPF solution. From the most general grid model described above, a tractable approximation (e.g., linear model) can be expressed as:
\begin{align}\label{eq:lin_model}
    \tilde{\mv{x}} &= \mv{A}\mv{u} + \mv{C_u}\mv{u}.
\end{align}
where the above model is the linear mapping of control variables $\mv{u}$ to approximated grid states $\tilde{\mv{x}}$ using a linear sensitivity matrix $\mv{A}$ and the incidence matrix $\mv{C_u}$.
\vspace{-0.1cm}
\uline{\textit{Implement Reliable Solution Algorithm:}} This step utilizes the processed grid model and its approximation to develop a reliable solution format. By ``reliable,'' it is meant that the solution algorithm should provide acceptable results in a timely fashion. This can be achieved by implementing a solution algorithm using a standard off-the-shelf solver, which can then be benchmarked against the individual algorithms developed by the user. This step is more than just converting power flow formulation to an optimization problem. This step is introduced in the DOPF core features to allow for testing various different algorithms and techniques. These are necessary as a wide number of distribution grid operation techniques can be developed depending on the application and necessity of the operator, and such a solution algorithm development process can be used to obtain a reliable solution algorithm. Next, we present two examples where an approximate model \cref{eq:lin_model} is shown to be utilized for setting up 1) a benchmark method that utilizes a standard DOPF formulation to be solved by an off-the-shelf solver and 2) an alternate method that solves the standard DOPF formulation using a different algorithm.

\begin{mini!}|s|[2]
    {}{\mv{c}\big(\mv{u},\tilde{\mv{x}} \big) \label{eq:apprx_DSOprob_OF_1}}
    {\label{eqs:DSOprob_1}}
    {}\addConstraint{\mv{f} (\mv{u},\tilde{\mv{x}})  }{= 0   \label{eq:DSOprob_b_1}}{:\mv{\lambda}}
    \addConstraint{\mv{g} (\mv{u},\tilde{\mv{x}})  }{\leq 0   \label{eq:DSOprob_b_2}}{:\mv{\mu}}
\end{mini!}
where $\mv{c}(\cdot)$ is the objective function term to be minimized and \cref{eq:DSOprob_b_1} and \cref{eq:DSOprob_b_2} are the approximated counterparts of \cref{eq:f} and \cref{eq:g}, formulated in a similar manner as \cref{eq:lin_model}. Variables $\mv{\lambda}$ and $\mv{\mu}$ are the Lagrange multiplier of the respective constraints. To compliment the tractability of the optimization problem, similar to approximation of constraints, the objective function can also be made linear or quadratic.
The standard  DOPF formulation of \cref{eqs:DSOprob_1} could be solved using an off-the-shelf solver \cite{GurobiOptimization.2016}. However, below is an example of an alternate solution algorithm, which can be developed and benchmarked against \cref{eqs:DSOprob_1}. Consider the Lagrangian of \cref{eqs:DSOprob_1}:
\begin{align}
    &\mathcal{L}(\mv{u},\tilde{\mv{x}}, \mv{\lambda}, \mv{\mu})
    = \mv{c}\big(\mv{u},\tilde{\mv{x}}\big) - \mv{{\lambda}}(\mv{f} (\mv{u},\tilde{\mv{x}}))
     -{\mv{\mu}}(\mv{g} (\mv{u},\tilde{\mv{x}})),
\end{align}
which can be deployed to update primal and dual variables at each iteration $k$ using the following rules:
\begin{subequations}
\begin{align}\label{eq:algorithmupdate_1}
     \text{Primal Updates: \ \ } \mv{u}^{k+1} &= \argmin \mathcal{L}_{\mv{u}}(\cdot)  \\
     \mv{x}^{k+1} & = \argmin \mathcal{L}_{\mv{x}}(\cdot) \\
     \text{Dual Updates: \ \ } \mv{\lambda}^{k+1} &= \mv{\lambda}^{k} - \alpha\Delta \mathcal{L}_{\mv{\lambda}}(\cdot) \\
    \mv{\mu}^{k+1} &= \mv{\mu}^{k} - \beta \Delta \mathcal{L}_{\mv{\mu}}(\cdot),
\end{align}
\end{subequations}
until there are no variables updates. In \cref{eq:algorithmupdate_1}, $\Delta \mathcal{L}_{\cdot}(\cdot)$ is the first order derivative of $\mathcal{L}$ with respect to $(\cdot)$ and $(\alpha, \beta)$ are small positive scalars. The algorithm \cref{eq:algorithmupdate_1} is formally known as the gradient descent \cite{Boyd.2004}.

\begin{comment}

From the above-mentioned example of the linear grid model, a standard
\begin{mini!}|s|[2]
    {}{\mv{A}\mv{u}+\mv{u}^\top \mv{B}\mv{u} \label{eq:apprx_DSOprob_OF_2}}
    {\label{eqs:DSOprob_2}}
    {}\addConstraint{\mv{C}\mv{u} + \mv{D}\mv{x}  }{= 0   \label{eq:DSOprob_b_2}}{:\mv{\lambda}}
    \addConstraint{\mv{E}\mv{u} + \mv{F}\mv{x}  }{\leq 0   \label{eq:DSOprob_b_2}}{:\mv{\mu}}
\end{mini!}

where $\mv{A}$ and $\mv{B}$ are matrices that determine the cost of different decision variables in $\mv{u}$ and $\mv{C},
\mv{D}, \mv{E}, \mv{F}$ are the matrices that relate decision variables $\mv{u}$ with state variables $\mv{x}$ and depend upon the system parameters such as the incidence matrix and network impedance values.
\end{comment}
\vspace{-0.5cm}
\begin{figure}[h]
\centering
\includegraphics[width=0.4\textwidth]{DOPF_Federate_Runner2.png}
\caption{DOPF federate components}
 \label{fig-dopffedrun}
\end{figure}

\subsection{DOPF Federate Development}\label{sec-dopf-fed}

Co-simulation environments have been proposed for realistic assessment of the developed applications and services for the grid \cite{palmintier2017design}. In these co-simulation environments, each individual simulation environment is called a federate, which has a pre-defined set of inputs, outputs and simulation settings. This federate architecture is particularly inspired by the co-simulation framework called Hierarchical Engine for Large-Scale
Infrastructure Co-Simulation (HELICS) \cite{palmintier2017design}.
\Cref{fig-dopffedrun} shows components of the DOPF federate runner. In \cref{fig-dopffedrun}, basic structure of the DOPF federate is shown where it is set up to subscribe data from and publishes data to other federates.\\
 \uline{\textit{Simulation Setting:}} These settings include simulation duration, update intervals and specification on the types of core features to be selected for the federate.\\
 \uline{\textit{Publication/Subscription:}} These contain information on the data to be inputted (subscribed) from and outputted (published) to relevant feeders to achieve the desired DOPF federate functionality. For example, current tap positions of the installed regulators would be required to calculate an optimal setpoint from the DOPF.\\
 \uline{\textit{Component Definitions:}} These definitions decompose the publication and subscription requirements into frequency (e.g., static versus dynamic), data types (e.g., normal versus sparse matrices) and appropriate federate mappings (e.g., control variables to be sent to the feeder and current state to be obtained from the appropriate sensors). These definitions also allow to switch out the relevant federates so that different testing environments may be created (e.g., removing state estimator and replacing with actual quantities from feeder federate etc.) for the DOPF federate. \\
\uline{\textit{DOPF Federate Runner:}} The runner script implements the DOPF components such that it is able to 1) demonstrate the required DOPF functionality (e.g., ensuring calling the appropriate DOPF feature) and 2) is coordinated with other federates, i.e., stopping the DOPF federate at the appropriate times to input and output the data (e.g., generating control setpoints every 30 minutes and inputting voltage estimates every 25 minutes so that the DOPF functions are executed within the required time.)
\begin{comment}
\begin{mini!}|l|
	{}{ c(\mv{u}) }
	{\label{eqs:opti}}{}
	\addConstraint{f(\mv{x},\mv{u})}{= \mv{0}  \label{eq:const1}}{;\hspace{15pt}(\mv{\lambda})}
	\addConstraint{g(\mv{x},\mv{u})}{\leq \mv{0} \label{eq:const2}}{;\hspace{15pt}(\mv{\mu})}
\end{mini!}
{\setlength{\parindent}{0cm}
The above problem can be cast as a convex or a non-convex optimization problem. A quadratic program variant of the above problem is given in \cref{sec:appendix}

To solve~\eqref{eqs:opti} as a central single optimization problem, the following complications exist.}
{\setlength{\parindent}{0cm}
\uline{\textit{Technical} Even for the case where feeder is not reconfigured through binary switches, the constraints \eqref{eq:const1}-\eqref{eq:const2} are usually highly nonlinear and non-convex. From the grid side, this non-convex contribution stems from power flow equations along with the legacy controls such as on-load tap changers, voltage regulators and feeder level switched capacitor banks. Similarly, from the technologies, state evolution (temperature of air conditioners) maybe non-convex.}
{\setlength{\parindent}{0cm}\uline{\textit{Diversity}Collecting and processing the information of the technology operations to be included in \eqref{eqs:opti} is a cumbersome task. For example, a BESS state-of-charge modeling is different than thermal modeling of flexible loads. Furthermore, the objectives to optimize technology operation also varies based on the preferences of its users.}
{\setlength{\parindent}{0cm}\uline{\textit{Privacy} The privacy issues may arise when sharing information between technology owners and grid operators exists. For example, the technology owners may not publicly make available their internal costs and constraints.}

\end{comment}
\subsection{DOPF Federate Integration}\label{sec-dopf-integration}
This section presents an example integration of the DOPF federate in a co-simulation environment. In \cref{fig-flow-overview2}, ``quantities'' are referred to complex nodal voltages and powers.\\
\uline{\textit{Feeder Federate:}} This represents an abstraction of a real-world grid in a co-simulation environment, which simulates the distribution grid power flow and generates appropriate quantities to be treated as the ``ground-truth.'' \\
\uline{\textit{Sensor Federate:}} For practicality purposes, information on the limited and practical sensor equipment of the distribution grid is simulated using a sensor federate. Inputs from the feeder federate on the nodes to be sensed are provided to this federate. Either less than total nodes, or noise on the ``ground truth'' quantities are assumed by the sensor federate to demonstrate the practical constraints of the grid.\\
\uline{\textit{State Estimate Federate:}} From the sensed quantities, the state estimator federate estimates the grid quantities.\\
\uline{\textit{Recorder Federate:}} To analyze performances of the federates, recorder federate is used to record the relevant quantities.\\
\uline{\textit{DOPF Federate:}} From the estimated grid quantities and grid topology information, the DOPF calculates control setpoints to be passed to the feeder federate. Note that as the model based DOPF core features are in discussion in this paper, grid information (e.g., topology and location of controllable equipment) is required. However, as data-driven (model-free) approaches are proposed in the literature, such requirements may be relaxed and hence may improve the overall data requirements for testing and developing DOPF federate. The DOPF federate runner described in \cref{sec-dopf-fed} coordinates the timing between different federates, i.e., obtaining relevant information (grid information and estimated quantities) and generating control setpoints in a timely fashion. Similar to the state estimate federate and feeder federate, variables determined by the DOPF federate are saved by the recorder federate. \Cref{sec:gadal-integration} gives an overview of the above-mentioned DOPF development and integration implementation \cite{Hanif_Distributed_Optimal_Power_2022}.
\begin{comment}
The proposed coordinating OPF in \cref{fig-methodology} sits in between the technology module and the grid module. The proposed algorithm obtains grid states and setpoints from external modules dedicated for performing grid feasibility and technology optimization, respectively. These states and setpoints are modified by the coordinating OPF and sent back to their respective modules. These external modules then provide the updated states and setpoints to the algorithm and the procedure repeats until the predefined criterion for the coordinating OPF is reached. In this way, the proposed method distributes information specific to the relevant modules i.e., setpoints for technology module and grid states/injections for grid module. This avoids the technical difficulties of solving (1b)-(1c) centrally. Also, the adopted information flow of the proposed method allows for diverse technologies to be included in the proposed method, while preserving their private information.
{\setlength{\parindent}{0cm}\paragraph{Grid Module}
Grid module only consists of solving a power flow problem along with the available legacy controllers such as voltage regulators, capacitor banks. Distribution grid software such as OpenDSS (an open source software) can perform these routines for a very large networks efficiently.}
{\setlength{\parindent}{0cm}\paragraph{Technology Module}
In this routine, technologies based on their individual preferences and constraints optimize their behavior. For example, a Battery Energy Storage System (BESS) can have an individual preference of participating in the wholesale market to maximize its revenue given its state of charge and power rating constraints.}

{\setlength{\parindent}{0cm}\paragraph{Coordinating OPF}
Consider at the beginning of the algorithm, at step~$k=0$, the algorithm collects setpoints $\mv{u}_{k}^{\textrm{tech}}$, grid states $\mv{x}_{k}$ and initialize LMs $\mv{\lambda}_k$/$\mv{\mu}_k$ and non-zero tuning parameters $\alpha/\beta$ to obtain modified set-points as follows:}
\begin{align} \label{eqs:coordinatingopf}
& { \footnotesize \mv{u}_{k+1}^{\textrm{tech}}:=   \frac{\partial  c(\mv{u}_{k}^{\textrm{tech}})}{\partial \mv{u}_{k}^{\textrm{tech}}} + (\frac{\partial g(\mv{x}_{k},\mv{u}_{k}^{\textrm{tech}})}{\partial {\mv{u}_{k}^{\textrm{tech}}}})^{\intercal}\mv{\mu}_k + (\frac{\partial f(\mv{x}_{k},\mv{u}_{k}^{\textrm{tech}})}{\partial \mv{u}_{k}^{\textrm{tech}}})^{\intercal}\mv{\lambda}_{k} = \mv{0}\nonumber}
\\
&{ \small g(\mv{x}_{k},\mv{u}_{k+1}^{\textrm{tech}}) :=  g(\mv{x}_{k},\mv{u}_{k}^{\textrm{tech}}) +   \frac{\partial g(\mv{x}_k,\mv{u}_{k+1}^{\textrm{tech}})}{\partial \mv{u}_{k}^{\textrm{tech}}} (\mv{u}_{k+1}^{\textrm{tech}} - \mv{u}_{k}^{\textrm{tech}})\nonumber} \\
&{ \small f(\mv{x}_{k},\mv{u}_{k+1}^{\textrm{tech}}) :=  f(\mv{x}_{k},\mv{u}_{k}^{\textrm{tech}}) +   \frac{\partial f(\mv{x}_k,\mv{u}_{k}^{\textrm{tech}})}{\partial \mv{u}_{k}^{\textrm{tech}}} (\mv{u}_{k+1}^{\textrm{tech}} - \mv{u}_{k}^{\textrm{tech}})\nonumber} \\
& { \small  \mv{\mu}_{k+1} = [\mv{\mu}_{k} + \alpha \cdot g(\mv{x}_{k},\mv{u}_{k+1}^{\textrm{tech}})]_{0}^{\infty}\nonumber} \\
&{ \small  \mv{\lambda}_{k+1} = \mv{\lambda}_{k} + \beta \cdot f(\mv{x}_{k},\mv{u}_{k+1}^{\textrm{tech}})}.
\end{align}
The new setpoints $\mv{u}_{k+1}^{\textrm{tech}}$ from \eqref{eqs:coordinatingopf} are in principle a gradient-descent step over penalizing setpoints which violate constraint space defined in \eqref{eqs:opti} \cite{Boyd.2014}, which come directly from taking the first-order derivative of the Lagrangian of \eqref{eqs:opti}\cite{Hanif.2018a}. These modified technology setpoints are then sent back to their individual modules in order for them to be verified against their internal requirements, i.e., to solve:
\begin{align}\label{eqs:feasibility-check}
    &    f(\mv{x}_{k},\mv{u}_{k+1}) = \mv{0} \nonumber\\
	&    g(\mv{x}_{k},\mv{u}_{k+1}) \leq \mv{0}.
\end{align}
The generated points above $(\mv{x}_{k+1},\mv{u}_{k+1})$, satisfy now both the internal technology constraints and grid feasibility, compactly represented as \eqref{eq:const1}-\eqref{eq:const2}. From the algorithm perspective, this essentially means that internal energy requirements are verified by the technology module, whereas power flows and legacy controls are executed by the grid module. With these alterations, a new states and setpoints~$(\mv{x}_{k+1},\mv{u}_{k+1}^{\textrm{tech}})$ are generated which are sent back to the coordinating OPF and the procedure is repeated until there is no modification of setpoints or the change in grid states. Note that throughout the algorithm, mappings are performed as needed. That is, the setpoints generated by the technology module are appropriately converted to active/reactive power injection vectors so that they correspond to the right order in the gradients computed in \eqref{eqs:coordinatingopf} and consequently in the grid module of \eqref{eqs:feasibility-check}.

The only computationally extensive part of the proposed algorithm is the gradient calculation, which requires admittance matrix inversion (see \cref{sensitivityderivations}). Ideally, gradients need to be updated for every new control state. One way to avoid this is to only update gradient if a large control variable change is observed within two consecutive iteration steps. Moreover, assuming grid topology does not change, the inversion of the admittance matrix (the most computational extensive part), can be pre-computed to speed up the gradient computation. It is going to be shown in \cref{Sec:Results} that these simplifying assumptions still yield the desired objective.

\end{comment}
\section{Appendix}\label{sec:gadal-integration}
The DOPF integration example implementation is shown in Fig. \ref{fig:dopf-tree}.
Each federated component must have its own directory with a minimum of the \texttt{component\_definitions.json} and the Python script to execute. The component definition file points the federate to the file to execute, maps the static inputs, and dynamic inputs/outputs. The DOPF implementation requires a \textit{scenarios} directory to initiate setup for the \texttt{LocalFeeder} models. The \texttt{LocalFeeder} directory generates the required DOPF information from the specified model and runs the co-simulation used by the DOPF federate. The \texttt{measuring\_federate} instantiates the federate variables to be published and \texttt{recorder} instantiates the federate variables to be subscribed to. The \texttt{wls\_federate} subscribes to the model voltages and powerflow for the given topology and publishes the estimated voltages. Figure \ref{fig:dopf-uml} outlines the DOPF federate and the primary utility functions used to solve. The \textit{OptimalPowerFlowFederate} class is primarily used to store the model information and establish the designated variable subscriptions and publications. The \texttt{opf\_grid\_utility\_script.py} contains the methods for federate matrices and vector ingestion.
\begin{figure}
    \centering
\includegraphics[width=0.5\textwidth]{opf_federate.png}
    \caption{UML Diagram for the DOPF Federate}
    \label{fig:dopf-uml}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.25\textwidth]{dopf_tree.png}
    \caption{DOPF Implementation File Structure}
    \label{fig:dopf-tree}
\end{figure}
\begin{comment}
\section{Numerical Modeling Details}\label{sec:appendix}
\sh{TODO: Validation of Tap Regulator Linearization Model}
{\setlength{\parindent}{0cm}\paragraph*{Notations} scalars are lower case $x$; vectors are bold and lower case~$(\mv{x})$; matrices are bold and upper case $(\mv{X})$. For complex~$(\mv{Y}) \in \mathbb{C} ^{\textrm{n} \times \textrm{n}}$ real part is $\Re(\mv{Y}) \in \mathbb{R} ^{\textrm{n} \times \textrm{n}}$, imaginary part is~$\Im(\mv{Y}) \in \mathbb{R} ^{\textrm{n} \times \textrm{n}}$ and complex conjugate is~$\conj{\mv{Y}}$.}
Operator $\diag(\mv{x})$ turns vector~$\mv{x}$ to a matrix with~$\mv{x}$ at its diagonal.

\section{Problem Statement}\label{sec:ps}
 The DSO solves the following multi-period constrained social welfare maximization problem
\begin{mini!}|s|[2]
    {}{\mv{f}\big(\mv{u},\mv{x} \big) \label{eq:apprx_DSOprob_OF_1}}
    {\label{eqs:DSOprob_1}}
    {}\addConstraint{\mv{g} (\mv{u},\mv{x}, \mv{W})  }{= 0   \label{eq:DSOprob_b_1}}{:\mv{\lambda}}
    \addConstraint{\mv{h} (\mv{u},\mv{x}, \mv{W})  }{\leq 0   \label{eq:DSOprob_b_1}}{:\mv{\mu}}
\end{mini!}

where $\mv{u}=\{\mv{p}^{\textrm{fl}}, \mv{q}_{\textrm{cb}}, \mv{t}_{\textrm{rg}}, p_0, q_0\}$ is the set of decision variables, $\mv{x}$ is the set of state variables (like voltages and currents) and $\mv{W}$ represents the network parameters such as the incidence matrix.

In the general case, the function $f(.)$, $g(.)$ and $h(.)$ are non-convex, however, for the linear OPF formulation, the above equations can be expressed as:

\begin{mini!}|s|[2]
    {}{\mv{A}\mv{u}+\mv{u}^\top \mv{B}\mv{u} \label{eq:apprx_DSOprob_OF_2}}
    {\label{eqs:DSOprob_2}}
    {}\addConstraint{\mv{C}\mv{u} + \mv{D}\mv{x}  }{= 0   \label{eq:DSOprob_b_2}}{:\mv{\lambda}}
    \addConstraint{\mv{E}\mv{u} + \mv{F}\mv{x}  }{\leq 0   \label{eq:DSOprob_b_2}}{:\mv{\mu}}
\end{mini!}

where $\mv{A}$ and $\mv{B}$ are matrices that determine the cost of different decision variables in $\mv{u}$ and $\mv{C},
\mv{D}, \mv{E}, \mv{F}$ are the matrices that relate decision variables $\mv{u}$ with state variables $\mv{x}$ and depend upon the system parameters such as the incidence matrix and network impedance values.
   \begin{mini!}|s|[2]
    {}{\mv{f}\big(\mv{p}_{\textrm{fl}}, \mv{q}_{\textrm{cb}}, \mv{t}_{\textrm{rg}}, p_0, q_0 \big) \label{eq:apprx_DSOprob_OF}}
    {\label{eqs:DSOprob}}
    {}\addConstraint{\mv{1}_{\textrm{L}}^\top (\mv{p}^{\textrm{cl}}  - \mv{I}_{\textrm{fl}}\mv{p}^{\textrm{fl}}) + p^{\text{l}} + \Delta p^l }{= p_0   \label{eq:DSOprob_b}}{:\lambda^{p}}
    \addConstraint{\mv{1}_{\textrm{L}}^\top (\mv{q}^{\textrm{cl}}  - \mv{I}_{\textrm{cb}}\mv{q}^{\textrm{cb}}) + q^{\text{l}} + \Delta q^l }{= q_0   \label{eq:DSOprob_b}}{:\lambda^{q}}
    \addConstraint{|\mv{v}^{-}| \leq |\mv{v}_\textrm{L}|}{\leq |\mv{v}^{+}| \label{eq:central_const3}}{:(\mv{\mu}_{\mv{v}}^{-},\mv{\mu}_{\mv{v}}^{+})}
	\addConstraint{ |{\mv{i}}_{ij}^{-}| \leq 	|\tilde{\mv{i}}_{ij}|}{\leq 	|{\mv{i}}_{ij}^{+}|\label{eq:central_const4}}{:(\mv{\mu}_{\mv{i}}_{ij}^{-},\mv{\mu}_{\mv{i}}_{ij}^{+})}
 	\addConstraint{ \mv{p}_{\textrm{fl}}^{-} \leq 	 \mv{p}_{\textrm{fl}}}{\leq 	\mv{p}_{\textrm{fl}}^{+}\label{eq:central_const5}}{:(\mv{\mu}_{\mv{p}_{\textrm{fl}}^{-}},\mv{\mu}_{\mv{p}_{\textrm{fl}}^{+}})}
 	\addConstraint{ \mv{q}_{\textrm{cb}}^{-} \leq 	 \mv{q}_{\textrm{cb}}}{\leq 	\mv{q}_{\textrm{cb}}^{+}\label{eq:central_const6}}{:(\mv{\mu}_{ \mv{q}_{\textrm{cb}}^{-}},\mv{\mu}_{ \mv{q}_{\textrm{cb}}^{+}})}
 		\addConstraint{ \mv{t}_{\textrm{rg}}^{-} \leq 	 \mv{t}_{\textrm{rg}}}{\leq 	\mv{t}_{\textrm{rg}}^{+}\label{eq:central_const7}}{:(\mv{\mu}_{\mv{t}_{\textrm{rg}}^{-}},\mv{\mu}_{\mv{t}_{\textrm{rg}}^{+}})}
\end{mini!}

The cost function consists of
\begin{align}
    \mv{f}&\big(\mv{p}_{\textrm{fl}}, \mv{q}_{\textrm{cb}}, \mv{t}_{\textrm{rg}}, p_0, q_0 \big) = \nonumber \\ &=\mv{c}^{a}_{\textrm{fl}}^{\intercal}(\mv{p}_{\textrm{fl}} - \mv{p}_{\textrm{fl}}^{r}) + \frac{1}{2} (\mv{p}_{\textrm{fl}} - \mv{p}_{\textrm{fl}}^{r})^{\intercal}\mv{C}^{b}_{\textrm{fl}}(\mv{p}_{\textrm{fl}} - \mv{p}_{\textrm{fl}}^{r}) \nonumber \\
    &+ \frac{1}{2} (\mv{q}_{\textrm{cb}}  - \mv{q}_{\textrm{cb}}^{r})^{\intercal}\mv{C}^{b}_{\textrm{cb}}(\mv{q}_{\textrm{cb}}  - \mv{q}_{\textrm{cb}}^{r}) \nonumber \\
    &+ \frac{1}{2} (\mv{t}_{\textrm{rg}}  - \mv{t}_{\textrm{rg}}^{r})^{\intercal}\mv{C}^{b}_{\textrm{rg}}(\mv{t}_{\textrm{rg}}  - \mv{t}_{\textrm{rg}}^{r}) \nonumber \\
    &+ c^{a}_{p0} p_0 + \frac{1}{2} p_0 c^{b}_{p0} p_0 + c^{a}_{q0} q_0 + \frac{1}{2} q_0 c^{b}_{q0} q_0
\end{align}
where $\mv{t}_{\textrm{rg}}$ and $\mv{t}_{\textrm{rg}}^{r}$ are the tap regulators and rated tap regulator positions ; $\mv{q}_{\textrm{cb}}$ and $\mv{q}_{\textrm{cb}}^{r}$ are the capacitor bank reactive power supply and the rated value of capacitor banks;  $\mv{p}_{\textrm{fl}}$ and $\mv{p}_{\textrm{fl}}^{r}$ are the flexible load reduction and their rated load values.
The global active and reactive power balance equation has an approximation of losses which is defined as:
\begin{align}
    \Delta p^l &= \mv{m}_{p}^{p^l}\mv{I}_{\textrm{fl}}\mv{p}^{\textrm{fl}} +  \mv{m}_{q}^{p^l}\mv{I}_{\textrm{cb}}\mv{q}^{\textrm{cb}} \nonumber \\
    \Delta q^l &= \mv{m}_{p}^{q^l}\mv{I}_{\textrm{fl}}\mv{p}^{\textrm{fl}} + \mv{m}_{q}^{p^l} \mv{I}_{\textrm{cb}}\mv{q}^{\textrm{cb}}
\end{align}
The voltage magnitudes are estimated using
\begin{align}
       |\tilde{\mv{v}}_\textrm{L}|  &= |\hat{\mv{v}}_\textrm{L}| + \mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{p}} \mv{I}_{\textrm{fl}}\mv{p}^{\textrm{fl}} + \mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{q}_{\textrm{L}}} \mv{m}_{q}^{p^l}\mv{I}_{\textrm{cb}}\mv{q}^{\textrm{cb}} \nonumber \\
       &+ \tau \diag(|\hat{\mv{v}}_\textrm{L}|)\mv{C}_{\mv{r}}(\mv{ \mv{t}_{\textrm{rg}}^{r}-\mv{t}_{\textrm{rg}}})
\end{align}
where $\tau$ is the regulator step-change ratio ($\tau = .00625$) and $\mv{C}_{\mv{r}} \in \mathbb{R}^{n\times r}$ is the matrix that links regulator location to voltage changes and is given by: $\mv{C}_{\mv{r}}(i,j)=1$ if tap regulator $j$ is upstream to node $i$, otherwise $\mv{C}_{\mv{r}}(i,j)=0$.
The current magnitude is estimated using
\begin{align}
    |\tilde{\mv{i}}_{ij}| = |\hat{\mv{i}}_{ij}| + \mv{M}^{|\mv{i}_{ij}|}_{\mv{p}}\mv{I}_{\textrm{fl}}\mv{p}^{\textrm{fl}} + \mv{M}^{|\mv{i}_{ij}|}_{\mv{q}_{\textrm{L}}}\mv{I}_{\textrm{cb}}\mv{q}^{\textrm{cb}}
\end{align}

\subsection{Distributed Optimization Problem Using Projected Gradient Descent Algorithm}
Consider the Lagrangian of \eqref{eqs:apprx_DSOprob_OF}:

\begin{align}
    &\mathcal{L}(\mv{u},\mv{x}, \mv{\lambda}, \mv{\mu})
    = \mv{f}\big(\mv{u},\mv{x}\big) - \mv{{\lambda}}(\mv{g} (\mv{u},\mv{x},\mv{W}))
     -{\mv{\mu}}(\mv{h} (\mv{u},\mv{x},\mv{W})),
\end{align}

\begin{subequations}
\begin{align}\label{eq:algorithmupdate_1}
     \text{Primal Updates: \ \ } \mv{u}^{k+1} &= \argmin \mathcal{L}_{\mv{u}}(\cdot)  \\
     \mv{x}^{k+1} & = \mv{k}(\mv{u}^{k+1},\mv{\lambda}^k,\mv{\mu}^k) \\
     \text{Dual Updates: \ \ } \mv{\lambda}^{k+1} &= \mv{\lambda}^{k} - \alpha\Delta \mathcal{L}_{\mv{\lambda}}(\cdot) \\
    \mv{\mu}^{k+1} &= \mv{\mu}^{k} - \beta \Delta \mathcal{L}_{\mv{\mu}}(\cdot),
\end{align}
\end{subequations}

\begin{align}
    &\mathcal{L}(\mv{p}_{\textrm{fl}}, \mv{q}_{\textrm{cb}}, \mv{t}_{\textrm{rg}}, p_0, q_0, \lambda^{(\cdot)}, \mv{\mu}^{(\cdot)}) \nonumber \\
    &= \mv{c}\big(\mv{p}_{\textrm{fl}}, \mv{q}_{\textrm{cb}}, \mv{t}_{\textrm{rg}}, p_0, q_0 \big) \nonumber \\
    & - {\lambda}_{p}(\mv{1}_{\textrm{L}}^\top (\mv{p}^{\textrm{cl}}  - \mv{I}_{\textrm{fl}}\mv{p}^{\textrm{fl}}) + p^{\text{l}} + \Delta p^l - p_0) \nonumber \\
    & -{\lambda}_{q}(\mv{1}_{\textrm{L}}^\top (\mv{q}^{\textrm{cl}}  - \mv{I}_{\textrm{cb}}\mv{q}^{\textrm{cb}}) + q^{\text{l}} + \Delta q^l - q_0) \nonumber \\
    &  -(\mv{\mu}_{\mv{v}}^{+})^{\intercal}(-|\mv{v}_\textrm{L}| + |\mv{v}^{+}|) -  (\mv{\mu}_{\mv{v}}^{-})^{\intercal}(|\mv{v}_\textrm{L}| - |\mv{v}^{-}|) \nonumber \\
    & -(\mv{\mu}_{\mv{i}}_{ij}^{+})^{\intercal}(-|\tilde{\mv{i}}_{ij}| + 	|{\mv{i}}_{ij}^{+}|)  - (\mv{\mu}_{\mv{i}}_{ij}^{-})^{\intercal}(|\tilde{\mv{i}}_{ij}| - 	|{\mv{i}}_{ij}^{-}|) \nonumber \\
    &  - (\mv{\mu}_{\mv{p}_{\textrm{fl}}^{+}})^{\intercal}(-\mv{p}_{\textrm{fl}} + \mv{p}_{\textrm{fl}}^{+}) -  (\mv{\mu}_{\mv{p}_{\textrm{fl}}^{-}})^{\intercal}(\mv{p}_{\textrm{fl}} -  \mv{p}_{\textrm{fl}}^{-}) \nonumber \\
    & - (\mv{\mu}_{\mv{q}_{\textrm{cb}}^{+}})^{\intercal}(-\mv{q}_{\textrm{cb}} + \mv{q}_{\textrm{cb}}^{+}) -  (\mv{\mu}_{\mv{q}_{\textrm{cb}}^{-}})^{\intercal}(\mv{q}_{\textrm{cb}} -  \mv{q}_{\textrm{cb}}^{-}) \nonumber \\
    &- (\mv{\mu}_{\mv{t}_{\textrm{rg}}^{+}})^{\intercal}(-\mv{t}_{\textrm{rg}} + \mv{t}_{\textrm{rg}}^{+}) -  (\mv{\mu}_{\mv{t}_{\textrm{rg}}^{-}})^{\intercal}(\mv{t}_{\textrm{rg}} -  \mv{t}_{\textrm{rg}}^{-}),
\end{align}
where $\lambda^{(\cdot)}$ and $\mv{\mu}^{(\cdot)}$ respectively contains all equality and inequality Lagrange multipliers listed in \eqref{eqs:opti_central}. Now using the above Lagrangian, and denoting $\Delta \mathcal{L}_{x}(\cdot)$ as the partial derivative of Lagrangian $\mathcal{L}$ with respect to variable $x$ , the distributed dual subgradient to solve \eqref{eqs:opti_central} can be written as:
\begin{subequations}
\begin{align}\label{eq:algorithmupdate}
     \text{Primal Updates: \ \ } \mv{p}_{\textrm{fl}}^{k+1} &= \argmin \mathcal{L}_{\mv{p}_{\textrm{fl}}}(\cdot)  \\
     \mv{q}_{\textrm{cb}}^{k+1} &= \argmin \mathcal{L}_{\mv{q}_{\textrm{cb}}}(\cdot)  \\
     \mv{t}_{\textrm{rg}}^{k+1} &= \argmin \mathcal{L}_{\mv{t}_{\textrm{rg}}}(\cdot)  \\
     \text{Dual Updates: \ \ } \lambda^{k+1} &= \lambda^{k} - \alpha\Delta \mathcal{L}_{\lambda}(\cdot) \\
    \mv{\mu}^{k+1} &= \mv{\mu}^{k} - \beta \Delta \mathcal{L}_{\mv{\mu}}(\cdot),
\end{align}
\end{subequations}
where $\alpha$ and $\beta$ are the update rates. The above distributed algorithm also incorporates the following conditions which are necessary to be for a feasible solution of \eqref{eqs:opti_central}. First, the first order optimality conditions states
\begin{align}\label{eq:firstorderoptimality}
    \Delta \mathcal{L}_{\mv{p}_{\textrm{fl}}}(\cdot) =  \Delta \mathcal{L}_{\lambda}(\cdot) = \Delta \mathcal{L}_{\mv{\mu}}(\cdot) = 0,
\end{align}
non-negative inequality LMs $\mv{\mu} \geq \mv{0}$ with complimentary slackness condition on inequality constraints \eqref{eq:central_const3}--\eqref{eq:central_const5}, and all feasible primal constraints  \eqref{eq:central_const1}--\eqref{eq:central_const5}. These conditions are also generalized and called KKT conditions. We do not formally define KKT conditions here for the sake of brevity.
For primal updates, we obtain the update rule using the partial derivative as:
\begin{align}
    &\Delta \mathcal{L}_{\mv{p}_{\textrm{fl}}}(\cdot) =   \mv{c}^{a}_{\textrm{fl}} + \mv{C}^{b}_{\textrm{fl}}(\mv{p}_{\textrm{fl}} - \mv{p}_{\textrm{fl}}^{r}) - \lambda_p(\mv{m}^{p^{l}}_{p}\mv{I}_{\textrm{fl}})^{\intercal} - \lambda_p(-\mv{1}_{\textrm{L}}^{\intercal}\mv{I}_{\textrm{fl}})^{\intercal}    \nonumber \\
   & - \lambda_q(\mv{m}_{p}^{q^l}\mv{I}_{\textrm{fl}})^{\intercal}  + \mv{I}_{\textrm{fl}}^{\intercal}\mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{p}_{\textrm{L}}}(\mv{\mu}_{\mv{v}}^{+}-\mv{\mu}_{\mv{v}}^{-}) +
	\mv{I}_{\textrm{fl}}^{\intercal}\mv{M}^{|\mv{i}_{ij}|}_{\mv{p}_{\textrm{L}}}(\mv{\mu}_{\mv{i}}_{ij}^{+}-\mv{\mu}_{\mv{i}}_{ij}^{-}) \nonumber \\
& + (\mv{\mu}_{\mv{p}_{\textrm{fl}}^{+}}-\mv{\mu}_{\mv{p}_{\textrm{fl}}^{-}}) = \mv{0}, \\
&\mv{p}_{\textrm{fl}}(k+1) = \Bigg[\mv{p}_{\textrm{fl}}^{r} + \Big(-\mv{c}^{a}_{\textrm{fl}} + \lambda_p(k)(\mv{m}^{p^{l}}_{p}\mv{I}_{\textrm{fl}})^{\intercal} +  \lambda_p(k)(-\mv{1}_{\textrm{L}}^{\intercal}\mv{I}_{\textrm{fl}})^{\intercal}    \nonumber \\
   &+ \lambda_q(\mv{m}_{p}^{q^l}\mv{I}_{\textrm{fl}})^{\intercal} - \mv{I}_{\textrm{fl}}^{\intercal}\mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{p}_{\textrm{L}}}(\mv{\mu}_{\mv{v}}^{+}(k)-\mv{\mu}_{\mv{v}}^{-}(k))  \nonumber \\
   &-\mv{I}_{\textrm{fl}}^{\intercal}\mv{M}^{|\mv{i}_{ij}|}_{\mv{p}_{\textrm{L}}}(\mv{\mu}_{\mv{i}}_{ij}^{+}(k)-\mv{\mu}_{\mv{i}}_{ij}^{-}(k))\Big)/\mv{C}^{b}\Bigg]^{\mv{p}_{\textrm{fl}}^{+}}_{\mv{p}_{\textrm{fl}}^{-}}
\end{align}
where notice that projection on to $[{\mv{p}_{\textrm{fl}}^{+}},{\mv{p}_{\textrm{fl}}^{-}}]$ removes the need to include and update later on the associated dual variables $[\mv{\mu}_{\mv{p}_{\textrm{fl}}^{-}},\mv{\mu}_{\mv{p}_{\textrm{fl}}^{+}}]$. We proceed with capacitor bank primal updates as:
\begin{align}
    &\Delta \mathcal{L}_{\mv{q}_{\textrm{cb}}}(\cdot) = \nonumber \\
    &  \mv{C}^{b}_{\textrm{cb}}(\mv{q}_{\textrm{cb}} - \mv{q}_{\textrm{cb}}^{r}) - \lambda_p(\mv{m}^{p^{l}}_{q}\mv{I}_{\textrm{cb}})^{\intercal} - \lambda_q(-\mv{1}_{\textrm{L}}^{\intercal}\mv{I}_{\textrm{cb}})^{\intercal}    \nonumber \\
   & - \lambda_q(\mv{m}_{q}^{q^l}\mv{I}_{\textrm{cb}})^{\intercal}  + \mv{I}_{\textrm{cb}}^{\intercal}\mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{q}_{\textrm{L}}}(\mv{\mu}_{\mv{v}}^{+}-\mv{\mu}_{\mv{v}}^{-}) \nonumber \\
& +\mv{I}_{\textrm{cb}}^{\intercal}\mv{M}^{|\mv{i}_{ij}|}_{\mv{q}_{\textrm{L}}}(\mv{\mu}_{\mv{i}}_{ij}^{+}-\mv{\mu}_{\mv{i}}_{ij}^{-}) + (\mv{\mu}_{\mv{q}_{\textrm{cb}}^{+}}-\mv{\mu}_{\mv{q}_{\textrm{cb}}^{-}}) = \mv{0}, \\
&\mv{q}_{\textrm{cb}}(k+1) = \Bigg[\mv{q}_{\textrm{cb}}^{r} +\Big( \lambda_p(\mv{m}^{p^{l}}_{q}\mv{I}_{\textrm{cb}})^{\intercal} +  \lambda_q(-\mv{1}_{\textrm{L}}^{\intercal}\mv{I}_{\textrm{cb}})^{\intercal}     \nonumber \\
   &+ \lambda_q(\mv{m}_{q}^{q^l}\mv{I}_{\textrm{cb}})^{\intercal} - \mv{I}_{\textrm{cb}}^{\intercal}\mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{q}_{\textrm{L}}}(\mv{\mu}_{\mv{v}}^{+}(k)-\mv{\mu}_{\mv{v}}^{-}(k))  \nonumber \\
   &-\mv{I}_{\textrm{cb}}^{\intercal}\mv{M}^{|\mv{i}_{ij}|}_{\mv{q}_{\textrm{L}}}(\mv{\mu}_{\mv{i}}_{ij}^{+}(k)-\mv{\mu}_{\mv{i}}_{ij}^{-}(k))\Big)/\mv{C}^{b}_{\textrm{cb}}\Bigg]^{\mv{q}_{\textrm{cb}}^{+}}_{\mv{q}_{\textrm{cb}}^{-}}
\end{align}
and for tap regulators as:
\begin{align}
    &\Delta \mathcal{L}_{\mv{t}_{\textrm{rg}}}(\cdot) =
      \mv{C}^{b}_{\textrm{rg}}(\mv{t}_{\textrm{rg}} - \mv{t}_{\textrm{rg}}^{r})   + \tau\mv{C}_{\mv{r}}^{\intercal} diag(|\hat{\mv{v}}_\textrm{L}|)(\mv{\mu}_{\mv{v}}^{+}-\mv{\mu}_{\mv{v}}^{-}) \nonumber \\
& - (\mv{\mu}_{\mv{t}_{\textrm{rg}}^{+}}-\mv{\mu}_{\mv{t}_{\textrm{rg}}^{-}}) = \mv{0}, \\
&\mv{t}_{\textrm{rg}}(k+1) = \Bigg[\mv{t}_{\textrm{rg}}^{r} -\Big(\tau\mv{C}_{\mv{r}}^{\intercal} diag(|\hat{\mv{v}}_\textrm{L}|)(\mv{\mu}_{\mv{v}}^{+}-\mv{\mu}_{\mv{v}}^{-}) \Big)/\mv{C}^{b}_{\textrm{rg}}\Bigg]^{\mv{t}_{\textrm{rg}}^{+}}_{\mv{t}_{\textrm{rg}}^{-}}
\end{align}
The net active and reactive power exports to the distribution grid can also be obtained using
\begin{align}
    &\Delta \mathcal{L}_{p_{0}}(\cdot) = \lambda_p -c_{p0}^a - c_{p0}^b p_0 = 0 \nonumber \\
    & p_0(k+1) =  (\lambda_p - c_{p0}^a)/c_{p0}^b \\
    &\Delta \mathcal{L}_{q_{0}}(\cdot) = \lambda_q -c_{q0}^a - c_{q0}^b q_0 = 0 \\
    & q_0(k+1) =  (\lambda_q - c_{q0}^a)/c_{q0}^b
\end{align}

The dual updates are performed as follows.
\begin{subequations}
\begin{align}
& \lambda_p(k+1) = \lambda_p(k)  \nonumber \\
&-\alpha_{\lambda_p}(\mv{1}_{\textrm{L}}^\top (\mv{p}^{\textrm{cl}}  - \mv{I}_{\textrm{fl}}\mv{p}^{\textrm{fl}}) + p^{\text{l}} + \Delta p^l - p_0) \\
& \lambda_q(k+1) = \lambda_q(k)  \nonumber \\
&- \alpha_{\lambda_q}(\mv{1}_{\textrm{L}}^\top (\mv{q}^{\textrm{cl}}  - \mv{I}_{\textrm{cb}}\mv{q}^{\textrm{cb}}) + q^{\text{l}} + \Delta q^l - q_0) \\
& \mv{\mu}_{\mv{v}}^{-}(k+1) = \Big[\mv{\mu}_{\mv{v}}^{-}(k) - \beta_{\mv{\mu}_{\mv{v}}^{-}} (|\mv{v}_\textrm{L}|(k) - |\mv{v}^{-}|)\Big]^{\infty}_{0} \\
& \mv{\mu}_{\mv{v}}^{+}(k+1) = \Big[\mv{\mu}_{\mv{v}}^{+}(k) - \beta_{\mv{\mu}_{\mv{v}}^{+}}(|\mv{v}^{+}|-|\mv{v}_\textrm{L}|(k))\Big]^{\infty}_{0}  \\
& \mv{\mu}_{\mv{\mu}_{\mv{i}_{ij}^{-}}}(k+1) = \Big[\mv{\mu}_{\mv{\mu}_{\mv{i}_{ij}^{-}}}(k) - \beta_{\mv{\mu}_{\mv{i}}^{-}}(|\tilde{\mv{i}}_{ij}|(k) - |{\mv{i}}_{ij}^{-}|)\Big]^{\infty}_{0}  \\
& \mv{\mu}_{\mv{\mu}_{\mv{i}_{ij}^{+}}}(k+1) = \Big[\mv{\mu}_{\mv{\mu}_{\mv{i}_{ij}^{+}}}(k) - \beta_{\mv{\mu}_{\mv{i}}^{+}}(|{\mv{i}}_{ij}^{+}| - |\tilde{\mv{i}}_{ij}|(k) )\Big]^{\infty}_{0}
\end{align}
\end{subequations}
where the projection $[0,\infty]$ enforces the non-negative inequality Lagrange multiplier constraint. Note that with the increase in the primal variables and constraints of the original problem, for examples addition of reactive power procurement in the objective function and including taps limitations in the constraints, the primal and dual updates are to be augmented accordingly.
\subsection{Algorithm Interpretation}
The proposed OPF has a physically intuitive flow. The LMs associated with the respective constraints can be imagined as penalizing variables. These variables become non-zero, once the respective constraints become binding (non-zero). The value of the LMs increases until the constraints are satisfied. In this case, the LMs will stop changing as well as the control variables, making the algorithm reach an equilibrium. For implementing soft constraints, the respective LMs can still serve as a signal to penalize control variables to follow in the direction of adhering to the imposed limits. However, the algorithm termination criterion shall be adjusted to reflect that. An example of such a scenario can be imagine when a technology's available flexibility is not enough to bring all grid voltages within a specified limit.

In order to implement the proposed algorithm, there is an implicit requirement that the gradients exist for each iteration in \eqref{eqs:coordinatingopf}. Assuming the generated control variables are feasible in the technology optimization space, this guarantee then comes down to power flow solution existence. There exists large number of literature dedicated to analyze this feasibility~\cite{Bernstein.2017a}, which is based on
formulating power flows with a fixed-point form. Also, fixed-point form allows to obtain gradients of grid states, yielding a global behavior of power flow equations. This is different to the usual first-order gradients as they are only accurate around the linearization point. This work also utilizes fixed-point form gradients~\cref{sensitivityderivations}. Also, the grid feasibility conditions to ensure that the new control variable ensures grid feasibility can also be included in~\eqref{eqs:coordinatingopf}. However, this may increase the computation per iteration step, as these conditions are nonlinear. Since practical grid control points yield a feasible power flow solution \cite{Bruno2011}, we did not explore this research path in this paper.

\subsection{Sensitivity Calculations}\label{sensitivityderivations}
This section describes the calculation of sensitivities which are utilized in the above explained distributed optimization algorithm.

\vspace{-0.05cm}{\setlength{\parindent}{0cm}\paragraph{Grid Model}
Assume steady state at bus~$i$ described by; complex voltage~$v_{i} = |v|_{i}e^{j\theta_{i}} \in \mathbb{C}$, with magnitude~$|v|_i \in \mathbb{R}$ and angle~$\theta_i \in \mathbb{R}$; and complex injections~$s_i = p_i + j q_i$, where $p_i \in \mathbb{R} $ and~$q_i \in \mathbb{R}$ are the active and reactive powers. For the whole grid, with buses numbered as:~$\{0,1,\ldots, \textrm{n}\} \in \textrm{N}$, where~$0$ is for the root-bus and the rest are~$\textrm{n}$ constant active reactive injection (PQ) buses, and complex voltages as: $\mv{v} := (v_{0},(\mv{v}_{\textrm{L}})^\intercal)^{\intercal} \in \mathbb{C}^{\textrm{N}}$, with $\mv{v}_{\textrm{L}} := (v_{1}, \ldots, v_{\textrm{n}})^{\intercal} \in \mathbb{C}^{\textrm{n}}$ and~$v_{0} \in \mathbb{C}$ are the voltages for~$\textrm{n}$~buses and the root-bus, respectively;  complex power vector~$\mv{s} := (s_{0},(\mv{s}_{\textrm{L}})^{\intercal})^\intercal \in \mathbb{C}^{\textrm{N}}$, where~$s_{0} \in \mathbb{C}$ and $\mv{s}_\textrm{L} := \mv{p}_\textrm{L} + j\mv{q}_\textrm{L} \in \mathbb{C}^{\textrm{n}}$. The derivation to follow can also be generalized to include ZIP loads~\cite{Bazrafshan.2017} along with delta connection and missing phases~\cite{Bernstein.2017a}.
Using above, power flow equations can be written as:}
{\normal\begin{align}\label{eq:pf_PQ_vonly}\vspace{-0.02cm}
{\mv{s}_\textrm{L}}  &= \diag(\mv{v}_\textrm{L}) (\conj{\mv{Y}}_{\text{L0}} \conj{v}_{0} + 	\conj{\mv{Y}}_{\text{LL}} \conj{\mv{v}}_\textrm{L}),
\end{align}}
where the goal is to find~$\mv{v}_\textrm{L}$ given the specified loading~$\mv{s}_\textrm{L}$ and the desired voltage at the root-bus~$v_{0}$. Above uses a split version of the nodal admittance matrix~${\mv{Y}} \in \mathbb{C}^{(\textrm{n}+1) \times (\textrm{n}+1)}$ as:
{\normal\begin{bmatrix}
\mv{Y}_{\textrm{00}} & \ \ \ 	\mv{Y}_{\textrm{0L}} \\
\mv{Y}_{\textrm{L0}} & \ \ \ 	\mv{Y}_{\textrm{LL}}
\end{bmatrix}}where~$\mv{Y}_{\text{00}} \in \mathbb{C}^{1}$,  $\mv{Y}_{\text{0L}} \in \mathbb{C}^{1 \times \textrm{n}}$,  $\mv{Y}_{\text{L0}}  \in \mathbb{C}^{\textrm{n} \times 1}$ and  $\mv{Y}_{\text{LL}} \in \mathbb{C}^{\textrm{n} \times \textrm{n}}$.
Now solving for~$\mv{v}_\textrm{L}$~\cite{Wang.2017}:
\begin{align}\label{eq:pf_fixedpoint}
\mv{v}_\textrm{L}   = \mv{{w}} + {\mv{Y}_{\text{LL}}}^{-1}{\diag(\mv{\conj{v}}_\textrm{L})}^{-1}\conj{\mv{s}}_\textrm{L},
\end{align}
where~$\mv{{w}} := -{\mv{Y}_{\text{LL}}}^{-1}\mv{Y}_{\text{L0}} v_0$ is the no-load voltage. Equation~\cref{eq:pf_fixedpoint} is of the fixed-point form, i.e., can be solved for $\hat{\mv{v}}_\textrm{L}$ at each iteration~$k$, given an initial estimate $\mv{v}^{0}_\textrm{L}$ and known $\hat{\mv{s}}_\textrm{L}$. Relevant for constraints in \eqref{eqs:opti}, we derive expressions for line currents $\mv{i}_{ij}$ $\in \mathbb{C}^{\textrm{m}}$ and system losses $s_l$ $\in \mathbb{C}$.
\begin{align}
	\mv{i}_{ij} &= \mv{Y}^{\textrm{ij}}_{0} v_0 + \mv{Y}^{\textrm{ij}}_\textrm{L} \hat{\mv{v}}_\textrm{L}, \label{eq:linecurrentwholegrid}	  \\
	s^{\textrm{l}} &= \conj{v}_0(\conj{\mv{Y}}_{\text{00}}\conj{v}_0 +	\conj{\mv{Y}}_{\text{0L}}\conj{\hat{\mv{v}}}_\textrm{L}) + \hat{\conj{\mv{v}}}_\textrm{L}^{\intercal}(\conj{\mv{Y}}_{\text{L0}}\conj{v}_0 +	\conj{\mv{Y}}_{\text{LL}}\hat{\conj{\mv{v}}}_\textrm{L})
\end{align}
where {\normal{$\mv{Y}^{\textrm{ij}} := (\mv{Y}^{\textrm{ij}}_{0},\mv{Y}^{\textrm{ij}}_\textrm{L}) \in \mathbb{C}^{\textrm{m} \times (\textrm{n}+1)}$}} are the line admittance matrix. For \eqref{eqs:coordinatingopf}, we deploy following approximations:

{\normal
\begin{subequations}\label{eqs:approx}
	\begin{align}
	|\tilde{\mv{v}}_\textrm{L}|  &= \mv{\hat{a}} + \mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{p}} \mv{p}_\textrm{L} + \mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{q}_{\textrm{L}}} \mv{q}_\textrm{L}, \label{eq:linvm} \\
	|\tilde{\mv{i}}_{ij}|  &= \mv{\hat{b}} + \mv{M}^{|\mv{i}_{ij}|}_{\mv{p}}\mv{p}_\textrm{L} + \mv{M}^{|\mv{i}_{ij}|}_{\mv{q}_{\textrm{L}}}\mv{q}_\textrm{L}, \label{eq:linpf_f} \\
	\Re(\tilde{s}^\textrm{l}) := \tilde{p}^{\textrm{l}}  &= \hat{c} + \mv{m}^{p^{\textrm{l}}}_{\mv{p}_\textrm{L}}\mv{p}_\textrm{L} + \mv{m}^{p^{\textrm{l}}}_{\mv{q}_\textrm{L}}\mv{q}_\textrm{L}, \label{eq:loss_p}. \\
	\tilde{q}^{\textrm{l}}  &= \hat{d} + \mv{m}^{q^{\textrm{l}}}_{\mv{p}}\mv{p}_\textrm{L} + \mv{m}^{q^{\textrm{l}}}_{\mv{q}_\textrm{L}}\mv{q}_\textrm{L} \label{eq:loss_q}.
	\end{align}
\end{subequations}}
Consider the first iteration from the satisfied ~\eqref{eq:pf_fixedpoint}, the complex voltage sensitivity with respect to PQ grid injections is:
{\normal\begin{equation}\label{eq:sens_v_cplx}
(\mv{M}^{\mv{v}_\textrm{L}}_{\mv{p}_\textrm{L}}, \mv{M}^{\mv{v}_\textrm{L}}_{\mv{q}_\textrm{L}}) := \Big(\frac{\partial \mv{v}_\textrm{L}}{\partial \mv{p}_\textrm{L}}, j\frac{\partial \mv{v}_\textrm{L}}{\partial \mv{q}_{\textrm{L}}}\Big) = \Big( {\mv{Y}_{\text{LL}}}^{-1}{\diag(\mv{\conj{v}}_\textrm{L})}^{-1}, -j {\mv{Y}_{\text{LL}}}^{-1}{\diag(\mv{\conj{v}}_\textrm{L})}^{-1} \Big). \nonumber
\end{equation}}
From this, the linear complex voltage~$\tilde{\mv{v}}_{\textrm{L}}$ follows:
{\normal \begin{equation} \label{eq:cmplxvoltage}
\tilde{\mv{v}}_{\textrm{L}} = \mv{w} + \mv{M}^{\mv{v}_\textrm{L}}_{\mv{p}_\textrm{L}} \mv{p}_\textrm{L} + \mv{M}^{\mv{v}_\textrm{L}}_{\mv{q}_\textrm{L}} \mv{q}_\textrm{L}
\end{equation}}
\noindent Now, linear voltage magnitudes in \eqref{eq:linvm} follows with {\normal{\normal$\hat{\mv{a}} = \mv{|w|}$,
$\mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{p}_\textrm{L}} := \Re(|\mv{W}|(\mv{W}^{-1}\mv{M}^{\mv{v}_\textrm{L}}_{\mv{p}_\textrm{L}}))$ and $\mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{q}_\textrm{L}} := \Re(|\mv{W}|(\mv{W}^{-1}\mv{M}^{\mv{v}_\textrm{L}}_{\mv{q}_\textrm{L}}))$, with $\mv{W} := \diag(\mv{w})$}}.
The complex linear current is: \vspace{-0.05cm}
\begin{align}
\tilde{\mv{i}}_{ij} = \mv{Y}^{\textrm{ij}}_{0} v_0 + \mv{Y}^{\textrm{ij}}_\textrm{L} \mv{w} + \mv{M}^{\mv{i}_{ij}}_{\mv{p}_{\textrm{L}}}\mv{p}_\textrm{L} + \mv{M}^{\mv{i}_{ij}}_{\mv{q}_{\textrm{L}}}\mv{q}_\textrm{L} \label{eq:losswholegrid}
\end{align}
where {\normal{$(\mv{M}^{\mv{i}_{ij}}_{\mv{p}_\textrm{L}},\mv{M}^{\mv{i}_{ij}}_{\mv{q}_\textrm{L}}) := (\mv{Y}^{\textrm{ij}}_\textrm{L}\mv{M}^{\mv{v}_\textrm{L}}_{\mv{p}_\textrm{L}}, \mv{Y}^{\textrm{ij}}_\textrm{L}\mv{M}^{\mv{v}_\textrm{L}}_{\mv{q}_\textrm{L}} \mv{q}_\textrm{L})$}}, allowing us to obtain the required parameters of~\cref{eq:linpf_f} as
{\normal
\begin{align*}
    \mv{M}^{|\mv{i}_{ij}|}_{\mv{p}_\textrm{L}} &:= \big(\diag(|\hat{\mv{i}}_{ij}| )\big)\Re\big(\diag(\conj{\hat{\mv{i}}}_{ij})\mv{M}^{\mv{i}_{ij}}_{\mv{p}_\textrm{L}}\big),
\end{align*}
\begin{align*}
    \mv{M}^{|\mv{i}_{ij}|}_{\mv{q}_\textrm{L}} &:= \big(\diag(|\hat{\mv{i}}_{ij}| )\big)\Im\big(\diag(\conj{\hat{\mv{i}}}_{ij})\mv{M}^{\mv{i}_{ij}}_{\mv{q}_\textrm{L}}\big),
\end{align*}
\begin{align}
\hat{\mv{b}} &:= |\mv{Y}^{\textrm{ij}}_{0} v_0 + \mv{Y}^{\textrm{ij}}_\textrm{L} \mv{w}|
\end{align}}
Finally, we have the approximate active power loss \cref{eq:loss_p}:\vspace{-0.05cm}
\begin{equation}\label{sensitivitiy_lossterm}
\frac{\partial \mv{s}^{\text{l}}}{\partial \mv{v}_\textrm{L}}  := v_{0}\mv{\conj{Y}}_{\text{0L}} +  (\mv{Y}_{\text{L0}}v_{0} + \mv{Y}_{\text{LL}}\mv{v}_\textrm{L})^{\intercal}  +  {\mv{v}_\textrm{L}}^{\intercal} \mv{\conj{Y}}_{\text{LL}}.
\end{equation}
Now following the chain rule and sensitivities {\normal{$\mv{M}^{\mv{v}_\textrm{L}}_{\mv{s}_\textrm{L}} := (\mv{M}^{\mv{v}_\textrm{L}}_{\mv{p}_\textrm{L}}, \mv{M}^{\mv{v}_\textrm{L}}_{\mv{q}_\textrm{L}})$}}~\eqref{eq:cmplxvoltage}, we get
{\normal \begin{equation}
(\mv{M}^{p^{\textrm{l}}}_{\mv{p}_\textrm{L}},\mv{M}^{p^{\textrm{l}}}_{\mv{q}_\textrm{L}}) := (\frac{\partial \mv{s}^{\text{l}}}{\partial \mv{p}_\textrm{L}},\frac{\partial \mv{s}^{\text{l}}}{\partial \mv{q}_\textrm{L}}) = \Re(\frac{\partial \mv{s}^{\text{l}}}{\partial \mv{v}_\textrm{L}}\mv{M}^{\mv{v}_\textrm{L}}_{\mv{s}_\textrm{L}}) + \Im(\frac{\partial \mv{s}^{\text{l}}}{\partial \mv{v}_\textrm{L}}\mv{M}^{\mv{v}_\textrm{L}}_{\mv{q}_\textrm{L}}) \nonumber
\end{equation}}
In order to complete the loss approximation~\eqref{eq:loss_p}, the required no-load condition \eqref{sensitivitiy_lossterm} can now be written as:
\begin{equation}\label{sensitivitiy_lossterm_noload}
\frac{\partial \mv{s}^{\text{l}}}{\partial \mv{w}}  := v_{0}\mv{\conj{Y}}_{\text{0L}} +  (\mv{Y}_{\text{L0}}u_{\text{0}} + \mv{Y}_{\text{LL}}\mv{w})^{\intercal}  +  {\mv{w}}^{\intercal} \mv{\conj{Y}}_{\text{LL}}.
\end{equation}
with the required no load coefficient as:
\begin{equation}
    \hat{c} = \conj{v}_0(\conj{\mv{Y}}_{\text{00}}\conj{v}_0 +	\conj{\mv{Y}}_{\text{0L}}\conj{\mv{w}}) + \hat{\conj{\mv{w}}}^{\intercal}(\conj{\mv{Y}}_{\text{L0}}\conj{v}_0 +	\conj{\mv{Y}}_{\text{LL}}\hat{\conj{\mv{w}}})
\end{equation}
Refer to~\cite{Hanif.2018c} for more details on the above approximations.

\subsection{Tap Regulator Linearization Model}
Tap regulators are discrete devices that can be used to manage the voltages in a distribution network. Let $\mv{u} \in \mathbb{Z}^r$ be the vector of tap regulator set-points, with $r$ being the number of tap regulators. Then because of the tap regulator actions the voltage equations in~\eqref{eq:linvm} is modified as:
\begin{align}
    |\tilde{\mv{v}}_\textrm{L}|  &= \mv{\hat{a}} + \mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{p}} \mv{p}_\textrm{L} + \mv{M}^{|\mv{v}_\textrm{L}|}_{\mv{q}_{\textrm{L}}} \mv{q}_\textrm{L}+\mv{M}_{\mv{r}}\mv{\Delta u}
\end{align}
where $\mv{\Delta u} \in \mathbb{Z}^r$ is the change in the regulator tap positions from nominal values and $\mv{M}_{\mv{r}} \in \mathbb{R}^{n\times r}$ is the sensitivity matrix of voltage magnitude with respect to the changes in tap positions, given by:
   $ \mv{M}_{\mv{r}}= \frac{\partial |{\mv{v}}_\textrm{L}|}{\partial \mv{u}}$
\\
To derive the expression for $\mv{M}_{\mv{r}}$ consider a system with tap regulators between each node, then the voltage magnitudes are related as:
\begin{align}
    |v_1|&=|v_0|+\tau u_{1}|v_0|\\
    |v_2|&=|v_1|+\tau u_{2}|v_1|\\
    \vdots\\
    |v_n|&=|v_{n-1}|+\tau u_{n}|v_{n-1}|
\end{align}
Calculating sensitivity of voltage magnitudes to tap changes in $u_1$ and assuming that the nominal tap positions are at zero, one gets:
\begin{align}
    \frac{\partial |v_1|}{\partial u_1}&=\tau |v_0|\\
    \frac{\partial |v_2|}{\partial u_1}&=\tau |v_0|\\
    \vdots\\
    \frac{\partial |v_n|}{\partial u_1}&=\tau |v_0|
\end{align}
Similarly sensitivities w.r.t $u_2$ can be calculated as:
\begin{align}
    \frac{\partial |v_1|}{\partial u_2}&=0\\
    \frac{\partial |v_2|}{\partial u_2}&=\tau |v_1|\\
    \vdots\\
    \frac{\partial |v_n|}{\partial u_2}&=\tau |v_1|
\end{align}
Such sensitivities can similarly be calculated w.r.t all $\mathbf{u}$ to obtain the sensitivity matrix $\mv{M}_{\mv{r}}$ as:

\begin{align}
    \begin{bmatrix} \frac{\partial |v_1|}{\partial a_1} \ \frac{\partial |v_1|}{\partial a_2} \ \hdots \ \frac{\partial |v_1|}{\partial a_n}\\
    \frac{\partial |v_2|}{\partial a_1} \ \frac{\partial |v_2|}{\partial a_2} \ \hdots \ \frac{\partial |v_2|}{\partial a_n}\\
    \vdots\\
    \frac{\partial |v_n|}{\partial a_1} \ \frac{\partial |v_n|}{\partial a_2} \ \hdots \ \frac{\partial |v_n|}{\partial a_n}
    \end{bmatrix}=\tau \begin{bmatrix}
    |v_0| & 0 & \hdots 0\\
    0 & |v_1| & \hdots 0\\
    \vdots\\
    0 & 0 & \hdots |v_{n-1}|\end{bmatrix}
    \begin{bmatrix}
    1 & 0 & \hdots 0\\
    1 & 1 & \hdots 0\\
    1 & 1 & \hdots 1
    \end{bmatrix}
\end{align}

Based on this expression, the sensitivity matrix can be expressed as:
\begin{align}
    \mv{M}_{\mv{r}}=\tau diag(\mv{w})\mv{C}_{\mv{r}}
\end{align}
where $\tau$ is the regulator step-change ratio ($\tau = .00625$) and $\mv{C}_{\mv{r}} \in \mathbb{R}^{n\times r}$ is the matrix that links regulator location to voltage changes and is given by: $\mv{C}_{\mv{r}}(i,j)=1$ if tap regulator $j$ is upstream to node $i$, otherwise $\mv{C}_{\mv{r}}(i,j)=0$.
\end{comment}
\section{Case Studies and Analysis}\label{Sec:Results}
This section presents results from the implemented DOPF core features. For the simulation, we have used the modified IEEE 123-bus test system. First we validate the linear power model, and then we compare the own implemented gradient-descent based distributed DOPF solution with its central variant which uses the off-the-shelf solver \cite{GurobiOptimization.2016} to showcase the scalability and computational efficiency.\\
\uline{\textit{Linear Model Analysis:}}
First, we showcase the validation of the linearized power flow model by changing the (i) flexible load, (ii) PV active power generations, (iii) regulator taps and (iv) kVars of the capacitor banks by introducing a multiplier -- denoted as \textit{``Load Multiplier''} -- and changing the value of the variable \textit{``Load Multiplier''} $\in [0,~1]$. Then we compare the nonlinear nodal voltage, current in the lines and line losses with approximated linear solutions, where OpenDSS is used to obtain the nonlinear solutions \cite{openDSS.2018}. In Fig. \ref{fig:Volt_flow_linearizations_123}, the linear approximated voltage and current for different nodes and lines are shown, where it can be seen relatively accurate linearized model is obtained around which can be used to generate DOPF solutions.
In addition to that, the approximated line losses have been compared in Fig. \ref{fig:Loss_linearization_123}, a common feature in DOPF minimization objectives. Similar to the node voltage and line current approximations, both the total kW and total kVar line losses can be approximated with the linearized model.\\
\uline{\textit{Solution Benchmark:}}
We benchmark the distributed DOPF solution (implemented using gradient descent algorithm) with its centralized variants implemented using off-the-shelf solver. We utilize (i) flexible loads, (ii) PV active power sheds, (iii) PV reactive power sheds, (iv) regulator tap settings and (v) capacitor banks switches for the optimal operation of the system. The problem objective is to reduce the cost of operation by minimally changing the current setpoints when subjected to a sudden change in the system and bring the system to a state in which the pre-specified operational limit has been maintained. The benchmark is shown in Fig. \ref{fig:Benchmark_opf_123} and Fig. \ref{fig:CvD_Stat}. We compare the major optimal control variables -- the shed of flexible loads and PV power outputs in Fig. \ref{fig:Benchmark_opf_123}; distributed control variables reaches the central solution in 79 iterations. Other control variables, such as cap bank switches and tap positions, remain the same for the system changes and thus produce trivial solutions, which are therefore not presented in this paper.
\begin{figure}[ht]
\centering
\includegraphics[width=0.5\textwidth]{C_v_D_PQF_123.png}
\caption{Benchmarking Distributed Solution across Iterations with the Central Solution Obtained from an Off-the-Shelf Solver for 123-Bus.}
 \label{fig:Benchmark_opf_123}
\end{figure}
To demonstrate that the DOPF solution is reliable and scalable, we create 100 different instances of control variable changes and compare both the distributed with central solutions. In the left figure of Fig. \ref{fig:CvD_Stat}, we see that the distributed solution is almost 100 times faster than the central solution -- $\sim$ 0.4 seconds  compared to $\sim$40 seconds -- and thus scalable for larger power distribution systems. However, the solution quality is not compromised as the objective value matches the central solutions (see the right figure of Fig. \ref{fig:CvD_Stat}).
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{C_v_D_123_Stat.png}
\caption{Statistics of 100 Randomly Initialized DOPF problem.}
 \label{fig:CvD_Stat}
\end{figure}

\section{Case Studies and Analysis}\label{Sec:Results}
This section presents results from the implemented DOPF core features.
\subsection{Linear Model Analysis}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/flow_plot_node12.png}
\caption{Line flow linearization results of three phases of one of the 13 Node circuit lines.}
 \label{fig-flow-linearizations}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/Flow_val_123.eps}
\caption{Line flow linearization results of 123 Bus System.}
 \label{fig-flow-linearizations_123I}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/v_plot_node0.png}
\caption{Node voltage linearization results for three phases of one of the 13 Node circuit node.}
 \label{fig-flow-linearizations}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/Volt_val_123.eps}
\caption{Node voltage linearization results for 123 Bus System.}
 \label{fig-flow-linearizations_123V}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/losses_plot.png}
\caption{Total active and reactive power loss linearization results.}
 \label{fig-flow-linearizations}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/Loss_val_123.eps}
\caption{Total active and reactive power loss linearization results for 123 bus.}
 \label{fig-flow-linearizations_123L}
\end{figure}

\subsection{Solution Benchmark}
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/solution_plot.png}
\caption{Bench-marking distributed solution across iterations (solid lines) with the central LP/QP solution (dotted line) obtained from an off-the-shelf solver.}
 \label{fig-flow-linearizations}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/C_v_D_123.eps}
\caption{Bench-marking distributed solution across iterations (solid lines) with the central LP/QP solution (dotted line) obtained from an off-the-shelf solver for 123 Bus.}
 \label{ooo}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/C_v_D_123_Q.eps}
\caption{Bench-marking distributed solution across iterations (solid lines) with the central LP/QP solution (dotted line) obtained from an off-the-shelf solver for 123 Bus.}
 \label{ooo}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/statistics_plot.png}
\caption{Statistics of 100 randomly initialized flexibility constraints and primal variable initial values.}
 \label{fig-flow-linearizations}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{figs/C_v_D_123_Stat.eps}
\caption{Statistics of 100 randomly initialized flexibility constraints and primal variable initial values for 123 bus.}
 \label{fig-flow-linearizations}
\end{figure}

\section{Conclusion}\label{Sec:Conclusions}
This paper presented a DOPF development and integration method. The intention of the paper was to help bring DOPF from a research tool to a more widely accepted solution in the industry. To that end, the paper presented a systematic way of developing and integrating a DOPF solution. For the development of DOPF, the three steps of a reliable solution development were: 1) processing the grid, 2) obtaining a tractable solution, and 3) implementing multiple solution methodologies for benchmarking purposes. Tests were conducted on IEEE 123-bus system where multiple solution algorithms implementation allowed reliably testing the developed DOPF scheme.

\section{Introduction}
\label{Sec:Introduction}

Electric distribution grids are experiencing significant changes in their infrastructure, such as installation of distributed energy resources (DERs). This infrastructure update must also accompany a renewed operation strategy for the distribution grids. Distribution grid optimal power flow (DOPF) could provide a pathway to implement such innovation in the distribution grids. In bulk power systems, optimal power flow (OPF) has received tremendous attention, as it has proved to dispatch flexible resources in order to balance supply and demand, while at the same time maintaining grid and device security constraints and limits~\cite{carpentier1962contribution}. Indeed, OPF is included in the most bulk power system Energy Management System (EMS) tools. However, DOPF is inherently different than OPF, as distribution grids have higher R/X ratio, have different phases in different portions of the grid and have a different topology. Hence, the same bulk system implementations are not simply transferable to DOPF and consequently to Distribution Management Systems (DMSs).

The theoretical research on the DOPF formulation and solution algorithms is an active research area. Traditionally \textit{DistFlow} based algorithms that are derived from the branch power flow models are used to represent the DOPF problem~\cite{Baran.1989,nazir2020optimal}. To deal with the non-convexity of the \textit{DistFlow} equations, several methods have been proposed in the literature, including using linear approximations~\cite{Baran.1989}, convex relaxation methods~\cite{gan2014convex} and convex restriction methods~\cite{nazir2021grid}. In addition, in distribution systems, the scheduling of discrete mechanical devices such as transformers, switches and capacitor banks needs to be optimized in conjunction with continuous flexible loads~\cite{nazir2018receding}.
With the emergence of these theoretical advancements, open-source example implementations of DOPF have been provided in the~\textit{PowerModels} package~\cite{coffrin2018powermodels} and the {\it{PowerModelsDistribution}}, specific for distribution systems~\cite{FOBES2020106664}. An OpenDSS driven OPF example is given in~\cite{rigoni2020open}.

Even though DOPF algorithms and their examples are receiving theoretical attention, their inclusion in Advanced Distribution Management System (ADMS) platforms \cite{agalgaonkar2016adms} has still not taken place. Some examples of the state-of-the-art ADMS applications are Volt/Var control set-point generation and Conservative Voltage Reduction (CVR) mechanisms, which can be considered the derivatives of the possible DOPF problem formulations. Implementing a generalized DOPF solution to eradicate all possible issues in a distribution grid is a challenging task. This is because distribution grids vary greatly in terms of size, locations and their unique operational challenges (e.g., remote versus island grids). Hence, the one-size-fits all solution of bulk systems, which usually aim to minimize the total cost of energy, may not be easily transferable to distribution grids. As a result, the distribution grids realm requires additional research on efficient methods to design, develop, integrate and test DOPF algorithms and solutions. The present paper aims to address this critical gap in research and development.

This paper presents a step-by-step procedure for developing the DOPF solution and testing it in an integrated simulation environment. The paper demonstrates the core feature development procedure for the DOPF. This procedure is important because it helps the developer to implement a reliable DOPF solution. The paper proposes implementation is achieved by processing the grid models in
 a form suitable for adoption by computationally efficient computer libraries/packages, developing approximation to reduce complexity and comparing alternate solutions. Co-simulation platforms are proposed to provide realistic testing conditions, since they allow for detailed physical simulations that mimic real-world communication infrastructure \cite{huang2018simulation}. Using the proposed core features, the paper demonstrates the development of a DOPF federate that can be integrated in a co-simulation environment. In the end, the paper presents an example of a fully integrated DOPF application in a co-simulation environment.

\input{DOPF_Development.tex}
\input{DOPF_Integration_Example.tex}
\input{case_studies}
\input{conclusion_outlook}
\input{appendix.tex}
\tiny{
\bibliographystyle{myIEEEtran}