# Distributed Optimization using Reduced Network Equivalents for Radial Power
Distribution Systems

Rabayet Sadnan , and Anamika Dubey Rabayet Sadnan and Anamika Dubey are with
the Department of EECS, WSU, Pullman, WA-99164, USA; e-mail:
<EMAIL_ADDRESS>
###### Abstract

The limitations of centralized optimization methods for power systems
operation have led to the distributed computing paradigm, particularly in
power distribution systems. The existing techniques reported in recent
literature for solving distributed optimization problems are not viable for
power distribution systems applications. The essential drawback remains a
large number of required communication rounds, i.e., macro-iterations among
the computing agents to solve one instance of the optimization problem; the
typical number of macro-iterations are in the order of $10^{2}\sim 10^{3}$. In
this paper, a new and scalable distributed optimization method based on
Equivalent Network Approximation (ENApp) is proposed to solve optimal power
flow (OPF) for a balanced radial distribution system. Specifically, the
distribution system’s radial topology is leveraged to reduce the decomposed
systems into upstream and downstream network equivalents. The proposed
innovations reduce the required number of macro-iterations/communication-
rounds for convergence by order of magnitude. The approach is validated using
IEEE 123-bus and IEEE 8500-node test systems and is shown to converge to the
same optimal solution as obtained using an equivalent centralized OPF (C-OPF)
model.

###### Index Terms:

Distributed optimization, Optimal Power Flow, Loss Minimization, Distributed
Energy Resources (DERs).

## Nomenclature

Sets

$\mathcal{G}=(\mathcal{N},\mathcal{E})$

Directed graph for the distribution system

$\mathcal{N}$

Set of all nodes

$\mathcal{E}$

Set of all distribution lines

$\mathcal{N}_{D}$

Set of all DER nodes in the system

$\mathcal{N}_{L}$

Set of all load nodes in the system

$A_{m}$

Directed graph representing the $m^{th}$ area.

$\mathcal{N}_{m}$

Set of all nodes in area $A_{m}$

$\mathcal{E}_{m}$

Set of all distribution lines in area $A_{m}$

$\mathcal{A}_{R}$

Set of all directed graphs defining all areas in a given distribution system.
Here, $\mathcal{A}_{R}=\\{A_{1},A_{2},...\\}$

$\mathcal{S}$

Set of operating constraints for the distribution system

$\mathcal{S}_{m}$

Set of operating constraints for the area $A_{m}$.

Functions

$F$

Cost function for the centralized problem

$\Phi_{m}$

Cost function for the distributed problem, $A_{m}$

$\psi_{1}$

Function projecting the network equivalence from upstream to downstream areas.

$\psi_{2}$

Function projecting the network equivalence from downstream to upstream areas.

Variables

$y_{1}$

Complicating variable representing voltage at the shared boundary bus.

$y_{2}$

Complicating variable representing power flow at the shared boundary bus.

$X$

Set of power flow and decision variables for the overall distribution system

$X_{m}$

Set of power flow, local, and complicating variables for area $A_{m}$

$x_{m}$

Set of local power flow and decision variables for area $A_{m}$, excluding
complicating variables.

$N_{A}=|\mathcal{A}_{R}|$

Total number of areas in the network

$x_{ij}$

Line reactance for line $(i,j)$

$r_{ij}$

Line resistance for ine $(i,j)$

$V_{i}$

Voltage magnitude for node $i$

$v_{i}=|V_{i}|^{2}$

Squared voltage magnitude for node $i$

$l_{ij}$

Squared magnitude of the current flow in line $(i,j)$

$P_{ij}$

Sending end active power for the line $(i,j)$

$Q_{ij}$

Sending end reactive power for the line $(i,j)$

$p_{L_{j}}$

Active load demand at the node $j$

$q_{L_{j}}$

Reactive load demand at node $j$

$p_{Dj}$

Active power generation by the DER at node $j$

$q_{Dj}$

Reactive power generation/demand by the DER at node $j$

$S_{DRj}$

kVA rating of the DER at node $j$

$I_{ij}^{rated}$

Thermal rating for the line $(i,j)$

##  I Introduction

The optimization of power systems operations, both at the transmission and
distribution levels, is crucial for an efficient and reliable grid.
Traditionally, the grid operations are centrally managed upon solving an
optimal power flow (OPF) problem using centralized optimization techniques;
the OPF problem is typically modeled as non-linear programming (NLP) problem
[1, 2]. With the increase in the number of controllable devices, especially in
power distribution systems, a centralized optimization paradigm suffers from
computational challenges and is susceptible to a single-point failure [3].
This has led to an interest in distributed optimization methods that (1)
reduce the computational requirements on a single decision-making unit by
distributing the problem into several smaller sub-problems that are
computationally simpler; (2) result in a decision-making paradigm with
multiple interacting agents that is robust to single-point failures; and (3)
relax the need for communication between the central controller and all
connected controllable/non-controllable assets.

The existing literature includes several methods to formulate the OPF as a
distributed optimization problem. In general, these methods adopt the
traditional distributed optimization techniques such as Augmented Lagrangian &
Method of Multipliers (ALMM), Alternating Direction Method of Multipliers
(ADMM), Auxiliary Problem Principle (APP), Predictor-Corrector Proximal
Multiplier method (PCPM), Analytical Target Cascading (ATC), etc., to model a
distributed optimal power flow problem (D-OPF) [3, 4]. For D-OPF applications,
the existing literature prefers methods based on ADMM and APP due to their
superior convergence properties [3]. Specific to radial distribution systems,
several methods based on the ADMM technique have been lately proposed [5, 6].
A D-OPF formulation includes macro-iterations to exchange the variables among
the decision-making agents and the micro-iterations to solve each agent’s
local subproblems. The macro-iterations dictate the communication systems’
requirements while micro and macro-iterations together decide the time-of-
convergence (ToC) for the algorithm. Unfortunately, both ADMM and APP
algorithms applied to OPF take a large number of macro-iterations ($\geq 100$
iterations), i.e., the message-passing between agents, to converge for
relatively small systems [7, 8]. While several recent methods attempted to
improve the ToC for local problems (micro-iterations) [9, 10, 6], they still
take hundreds if not thousands of macro-iterations to converge. Another
popular approach is the APP algorithm that linearizes the penalty term in the
augmented Lagrangian [11]. While the APP-based distributed OPF methods are
relatively faster and require sharing fewer variables, the number of macro-
iterations (i.e., communication rounds among distributed agents) needed to
reach convergence is of the same order as in ADMM [12].

In power distribution systems, having a large number of macro-iterations for a
single step of the optimization problem is not preferred as it will lead to
significant communication delays with the decision-making. Thus, a practical
implementation of such algorithms requires a very fast communication system to
reach a converged solution within a reasonable time. Deploying such a
communication system is expensive, and thus, they are typically not available
in power distribution systems. Besides requiring a large number of macro-
iterations, ADMM and APP based methods applied to D-OPF pose several other
limitations, including, but not limited to, parameter tuning and requiring at
least one grid-forming DER in each area. Lately, to address some of these
challenges, real-time feedback based online distributed algorithms have been
explored in the related literature for network optimization [13, 14, 15].
Generally, these algorithms do not wait to optimize for a time-step but
asymptotically arrive at an optimal decision over several steps of real-time
decision-making. For example, authors developed distributed optimal reactive
power control in [13], that can regulate the voltage and track the network-
level objective. However, they cannot achieve a global optimum solution [13].
Further, with 16 DERs in IEEE 123 bus system, this approach takes around 100s
of macro-iterations to converge, thus is slow in tracking the optimal solution
[13]. In [14], authors developed a unified approach for real-time optimization
of distribution grids. However, this approach requires a central coordinator
to update the dual variables. In [15], authors designed a distributed feedback
control to address some of the issues with the asynchronous and delayed
communications among smart agents; still, the algorithm takes hundreds of
iterations to converge/track the optimal solution for a mid-size feeder. This
raises further challenges to the performance of the algorithm for larger
feeders, especially during the fast varying phenomenon.

While distributed optimization methods are argued to be computationally
efficient for larger systems, the existing state-of-the-art methods are not
practically viable for power systems applications. Among several limitations,
the primary drawback remains the requirement for a large number of macro-
iterations requiring several hundred rounds of communication for information
exchange among distributed agents to solve one instance of the optimization
problem. Thus, implementing such algorithms requires expensive communication
infrastructure with high bandwidth and low communication delays. Also, the
existing methods need tuning of several parameters to achieve convergence.
Some techniques require intelligence at each node of the system, making it
challenging to scale for larger systems. The online distributed control
algorithms also suffer from the same challenges of slow convergence, resulting
in sub-optimal performance and slow tracking of the optimal solutions.

To address these challenges, in this paper, we propose a new distributed
optimization method to solve the distributed optimal power flow (D-OPF)
problem based on reduced Equivalent Network Approximation (ENApp) for radial
distribution systems. The novelty lies in the proposed decomposition approach
for the centralized optimization problem that actively leverages the radial
topology of the distribution systems. The proposed approach follows the primal
decomposition method combined with the block-coordinate descent algorithm for
distributed optimization. It does not require any central coordinator for
optimization and only involves communication among the neighbors, resulting in
a distributed algorithm as defined by the recent relevant literature [3].
Further, no parameter tuning is needed making the approach robust for any
operating conditions. The key contributions of this work are:

  * •

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;

  * •

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.

  * •

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.

##  II Modeling & Problem Formulation

In this paper, $(\cdot)^{T}$ represents matrix transpose; $(\cdot)^{*}$
represents the complex-conjugate; $|.|$ symbolizes the cardinality of a
discrete set or the absolute value of a number; $(\cdot)^{k}$ represents the
$k^{th}$ iteration; and $j=\sqrt{-1}$; $\overline{(.)}$ and $\underline{(.)}$
denotes the maximum and minimum limit of a given quantity.

###  II-A Network Modeling

Let us consider a radial power distribution network with $n$ nodes. Let
$\mathcal{N}$ be the set of all nodes in the system, i.e.,
$\mathcal{N}=\\{1,2,…,n\\}$. Let the directed graph
$\mathcal{G}=(\mathcal{N},\mathcal{E})$ is used to represent the distribution
system. Here $\mathcal{E}$ denotes the set of all distribution lines
connecting the ordered pair of buses $(i,j)$ i.e., from node $i$ to node $j$.
The distribution lines are modeled as series impedance, $r_{ij}+jx_{ij}$ for
$\forall(i,j)\in\mathcal{E}$. Let, $\mathcal{N}_{L}$ represents the set of
load nodes, and $\mathcal{N}_{D}$ represents the set of all DER nodes. For a
radial topology, the cardinality of sets $\mathcal{N}$ and $\mathcal{E}$ is
related as the following: $|\mathcal{N}|=|\mathcal{E}|+1$. Further,
$|\mathcal{N}|\geq|\mathcal{N}_{L}|\geq|\mathcal{N}_{D}|$. Let, for node $j$,
$i$ is the unique parent node, and $k$ be the set of all children nodes. Then,
in $k:j\rightarrow k$, $k$ represents the set of children nodes for the node
$j$. Next, we denote $v_{j}=|V_{j}|^{2}=V_{j}{V_{j}}^{*}$ as the squared
magnitude of voltage at node $j$ $(V_{j})$. Likewise, let $l_{ij}$ be the
squared magnitude of current flow in branch $(i,j)$.

The network is modeled using the branch flow equations [16] defined for each
line $(i,j)\in\mathcal{E}$ and $\forall j\in\mathcal{N}$.

$\displaystyle\small P_{ij}-r_{ij}l_{ij}-p_{L_{j}}+p_{Dj}=\sum_{k:j\rightarrow
k}P_{jk}$ (1a) $\displaystyle
Q_{ij}-x_{ij}l_{ij}-q_{L_{j}}+q_{Dj}=\sum_{k:j\rightarrow k}Q_{jk}$ (1b)
$\displaystyle
v_{j}=v_{i}-2(r_{ij}P_{ij}+x_{ij}Q_{ij})+(r_{ij}^{2}+x_{ij}^{2})l_{ij}$ (1c)
$\displaystyle v_{i}l_{ij}=P_{ij}^{2}+Q_{ij}^{2}$ (1d)

where, $P_{ij},Q_{ij}\in\mathbb{R}$ is the sending-end active and reactive
power flows for $(i,j)\in\mathcal{E}$. Complex power
$S_{L_{j}}=p_{L_{j}}+jq_{L_{j}}$ is the load connected and
$S_{D_{j}}=p_{Dj}+jq_{Dj}$ is the power output of DER connected at node $j$.
Note that, $\forall j\not\in\mathcal{N}_{D}$, $S_{D_{j}}=0$. Also, if
$j\not\in\mathcal{N}_{L}$, then $S_{L_{j}}=0$.

###  II-B DER Modeling

The DERs are modeled as balanced single-phase photovoltaic modules (PVs),
interfaced with the distribution system using smart inverters capable of two-
quadrant operation. The active power generation for the DER is specified by
$p_{Dj}$ and is assumed to be known and not-controllable. The reactive power
from DERs, $q_{Dj}$, is assumed to be controllable and modeled as the decision
variable for the resulting distribution-level optimal power flow problem. The
apparent power rating of the smart inverter limits the reactive power dispatch
($q_{Dj}$) from the DER.  Let the rating of the DER connected at node
$j\in\mathcal{N}_{D}$ be $S_{DRj}$. Then the bounds on $q_{Dj}$ $\forall$
$j\in N_{D}$ are expressed using (2).

$\displaystyle-\sqrt{S_{DRj}^{2}-p_{Dj}^{2}}\leq
q_{Dj}\leq\sqrt{S_{DRj}^{2}-p_{Dj}^{2}}$ (2a)

###  II-C Centralized Non-linear and Relaxed OPF Problem

A centralized OPF (C-OPF) problem is defined to optimize the network for some
cost function characterized by a network-level problem objective, the power
flow constraints in (1a)-(1d), and the operating constraints on the power flow
variables. In this paper, we formulate the problem objective of minimizing
network power losses. The aim is to reduce the network losses by controlling
the reactive power output from DERs ($q_{Dj}$).

Let $X$ define the problem variables (including power flow and decision
variables), where $X=[P_{ij}$, $Q_{ij}$, $l_{ij}$, $v_{j}$, $q_{Dj}]^{T}$;
$\forall j\in\mathcal{N}$, and $\forall(i,j)\in\mathcal{E}$. Let $F(X)$ denote
the objective function representing the total power loss in the given
distribution system. Note that $F(X)$ is a function of both the power flow
variables and decision variables. Then, the centralized OPF problem is defined
as the following in (P1).

$\displaystyle\text{(P1)}\hskip 11.38092pt\min\hskip
5.69046pt{\color[rgb]{0,0,0}F(X)}=\sum_{(i,j)\in\mathcal{E}}l_{ij}r_{ij}\hskip
17.07182pt$ (3a) s.t.(1a) - (1d), and (2a) (3b)
$\displaystyle\underline{V}^{2}\leq v_{i}\leq\overline{V}^{2}\hskip
25.6073pt\forall i\in\mathcal{N}$ (3c) $\displaystyle
l_{ij}\leq\left(I^{rated}_{ij}\right)^{2}\hskip
14.22636pt\forall(i,j)\in\mathcal{E}$ (3d)

where, $\overline{V}=1.05$ and $\underline{V}=0.90$ are the limits on bus
voltages, and $(I^{rated}_{ij})^{2}$ is the thermal limit for the branch
$(i,j)$.

The NLP problem in (P1) can be reformulated as a convex optimization problem
upon relaxing the quadratic equality constraints (1d), (1d), as a convex
inequality [17]. This relaxation leads to a Second-Order Cone Programming
(SOCP) model, as detailed in (P2) using (4a)-(4c).

$\displaystyle\text{(P2)}\hskip 11.38092pt\min\hskip
5.69046pt{\color[rgb]{0,0,0}F(X)}=\sum_{(i,j)\in\mathcal{E}}l_{ij}r_{ij}$ (4a)
$\displaystyle\text{s.t.}\hskip
5.69046pt\text{\eqref{eqModel_nonlin1},\eqref{eqModel_nonlin2},
\eqref{eqModel_nonlin3}, \eqref{DG_lim}},\text{\eqref{nonlin_OPF1}, and
\eqref{nonlin_OPF2}}$ (4b) $\displaystyle v_{i}l_{ij}\geq
P_{ij}^{2}+Q_{ij}^{2}\hskip 17.07182pt\forall(i,j)\in\mathcal{E}$ (4c)

Remark 1: Note that for a balanced radial distribution system, the
aforementioned SOCP relaxation is exact under specific sets of conditions
[18]. Also, whenever the relaxed mode is exact, (4c) is active for the optimal
solutions of the relaxed problem (P2). That is, the obtained solution using
the relaxed model is also the optimal and power flow feasible for the original
NLP. However, there may be cases when the solutions obtained from the relaxed
SOCP model are not power flow feasible. These conditions are extensively
studies in recent literature [19, 20]. In those cases, the NLP in (P1) needs
to be solved to obtain the optimal and feasible power flow solution. For this
reason, in the proposed D-OPF method, we solve the NLP model in (P1). However,
we have introduced the SOCP model (P2) here to obtain the solution for the
centralized OPF problem for the large test feeder used in simulations. This is
because, for the large feeder, the NLP model did not converge. Thus, the
solutions from the SOCP model were used for the comparison and validation
purpose. Also, for the specific test case used in this paper, the SOCP
solutions were found to be power flow feasible for the original NLP model and
thus can be used for comparison of the solution quality.

Remark 2: It should be noted that the distribution system also includes legacy
voltage control devices such as voltage regulators and capacitor banks. The
legacy control devices, while not explicitly modeled here, can be easily
incorporated in the proposed D-OPF model. However, upon incorporating such
legacy devices, the optimization problem becomes Mixed-Integer Non-Linear
Programming (MINLP) problem. Previously, we proposed a bi-level approach to
manage the complexity of the MINLP problem resulting from the control of both
legacy devices and smart inverters [21]. This algorithm can be adopted to
solve the resulting MINLP problem in a computationally tractable manner.

Furthermore, note that in the above formulation we have only considered
reactive power from DGs ($q_{Dj}$) as controllable. However, the OPF
formulation can be easily extended to the cases when both active ($p_{Dj}$)
and reactive ($q_{Dj}$) power are controllable. In this case, the constraint
on DER power rating defined in (2) simply transforms into a SOCP constraint.

_Algorithm 1: BCD Method_

Initialize the variable blocks $(x_{1}^{0},x_{2}^{0},...,x_{s}^{0})$

And $k=1$;

while  _stopping criterion is not met_ do

for  _$i=1,2,...,s$_ do

$x_{i}^{k}:=\underset{x_{i}\in\mathcal{S}_{i}}{\text{argmin}}\hskip
5.69046pt{\color[rgb]{0,0,0}\Phi_{i}^{k}(x_{i})}$

end for

$X^{k}=(x_{1}^{k},x_{2}^{k},...,x_{s}^{k})$

$k\leftarrow k+1$

end while

Return Global Minimizer: $X^{*}\leftarrow X^{k}$

##  III Distributed Optimal Power Flow (D-OPF)

The loss minimization OPF problem defined in (P1) and (P2) is decomposable
into connected areas. In this paper, it is assumed that the distribution
system is composed of $N_{A}$ distributed areas. Let $\mathcal{A}_{R}$ be the
set of all decomposed areas in a given distribution system, where
$\mathcal{A}_{R}=\\{A_{1},A_{2},…,A_{N_{A}}\\}$. Also, let each area, $A_{m}$,
be defined as a directed graph
$A_{m}=\mathcal{G}(\mathcal{N}_{m},\mathcal{E}_{m})$, where
$m=\\{1,2,...,N_{A}\\}$ that is controllable using a Local Controller (LC).
The LCs solve a decomposed OPF problem to minimize the power losses for the
specific area by controlling the local decision variables. Next, the proposed
ENApp D-OPF framework is used to coordinate among the neighboring LCs and to
facilitate network-level optimization. The proposed distributed approach is
analogous to a primal-decomposition process that does an area-wise breakdown
of the centralized problem and requires exchanging only primal variables [22].
However, unlike the traditional primal-decomposition method, a new algorithm
is proposed that leverages a reduced equivalent network representation for the
power flow in radial distribution systems. The resulting sub-problems are
solved using an approach similar to the Block Coordinate Descent (BCD) method
described in [23]. The master problem requires solving only a set of non-
linear equations of the form $f=0$. Here, we use the Fixed Point Iteration
(FPI) method to solve the master problem.

###  III-A Block Coordinate Descent (BCD) Method

In the BCD algorithm, the original problem variables are partitioned into
several blocks. All but one block of the variables are fixed, and the problem
is minimized over the free blocks. In this way, the problem is minimized for
all the blocks. For example, the BCD method was used to solve the systems of
equations [24, 23]. The convergence properties of the BCD algorithm are
reported in [25, 26]. Generally, for a nonconvex problem, the BCD algorithm
admits a sub-linear convergence. While, for a strongly convex optimization
problem, the BCD algorithm is known to converge linearly.

A brief discussion of the BCD method to solve optimization problems is
described next. The objective is to solve the optimization problem detailed in
(5a) that aims at minimizing the cost function $F(X)$, subject to the
constraints defined by the constraint set $\mathcal{S}$. The BCD approach
solves the optimization problem defined in (5) in a distributed way. The
approach to solving (5) using the BCD method is detailed in Algorithm 1 [23].
Let variable $X$ be partitioned into $s$ blocks, where,
$X=(x_{1},x_{2},...,x_{s})$. We define the modified cost function,
$\Phi_{i}^{k}(x_{i})$, for the $i^{th}$ block at $k^{th}$ iteration, in (6a).
In the BCD approach, the decomposed cost function defined in (6) is minimized
over one block while keeping other blocks constant. This process is repeated
for other blocks until a prespecified stopping criterion is met. Thus, upon
convergence, the solution for the original optimization problem defined in (5)
is obtained in a distributed way.

$\small\min_{X}\hskip 5.69046ptF(X)\text{;}\hskip 8.5359pt\text{ s.t. }\hskip
8.5359pt\text{$X\in\mathcal{S}$}$ (5a)
$\small\Phi_{i}^{k}(x_{i})=F(x_{1}^{k},x_{2}^{k},...,x_{i-1}^{k},x_{i},x_{i+1}^{k-1},...,x_{s}^{k-1})\text{
s.t. }\hskip 2.84544pt\text{$x_{i}\in\mathcal{S}_{i}$}$ (6a)

Figure 1: An example of a two-area system.

###  III-B Proposed ENApp Algorithm for D-OPF

The developed D-OPF method leverages the radial topology of the power
distribution network. Note that for a radial distribution system, there is a
unique upstream area (UA) for any downstream area (DA). However, a UA (parent)
can be connected to multiple DAs (child). We model a dummy bus for each UA-DA
pairs; the dummy bus is shared by both the UA and DA regions. The dummy bus
acts as a fixed voltage source for the DA while as a load (PQ) bus for the
respective UA. This method is based on reduced network equivalents where the
DA and UA are replaced by a fixed load and a fixed voltage source at each
iteration, respectively. Each area solves its own OPF problem and then
exchanges the boundary variables only with its immediate neighbors.
Specifically, the DA shares complex power flow with the UA, while the UA
shares bus voltages with the DA. The macro-iterations defining the exchange of
variables among neighboring areas continue until both the UA and DA agree on
the same value for the shared bus variables within a prespecified tolerance.

Note that the proposed D-OPF method is a modified version of the BCD algorithm
for solving the sub-problems. However, in the proposed approach, the sub-
problems for the distributed optimization modules are executed in parallel as
opposed to the sequential execution adopted in the traditional BCD algorithm.
Further, unlike the primal decomposition described in [22], the
complicating/shared variables are split into two sets. One set of the shared
variables are minimized in one of the subproblems while keeping the other set
fixed. The process is then repeated for the other set of primal variables
while keeping the first set fixed this time, similar to the BCD method. The
iteration stops when the change in shared variables is within a prespecified
tolerance.

_Algorithm 2: ENApp D-OPF method for two Area System_

Initialize $\hat{y}^{0}_{1}\in S_{1}$ and $\hat{y}^{0}_{2}\in S_{2}$

Initialize $\epsilon=1,k=1$

while  _$|\epsilon| >0.001$_ do

Compute Parallely ($L_{1},L_{2}$) –

$L_{1}:$ $X_{1}^{k}$ := $\underset{X_{1}\in S_{1}}{\text{argmin}}\hskip
5.69046pt\Phi_{1}(X_{1},\hat{y}_{2}^{k-1})$

$L_{2}:\hskip 2.84544ptX_{2}^{k}:=\underset{{\color[rgb]{0,0,0}X_{2}\in
S_{2}}}{\text{argmin}}\hskip 5.69046pt\Phi_{2}(X_{2},\hat{y}_{1}^{k-1})$

Check residual

$\mathcal{R}$ = $\left[\begin{array}[]{c}y_{1}^{k}-\hat{y_{1}}^{k-1}\\\
y_{2}^{k}-\hat{y_{2}}^{k-1}\end{array}\right]$

FPI update:

$\hat{y}_{1}^{k}=\psi_{1}(\hat{y}_{2}^{k-1},x_{1}^{k})$

$\hat{y}_{2}^{k}=\psi_{2}(\hat{y}_{1}^{k-1},x_{2}^{k})$

Increase iteration count $k$:

$k\leftarrow k+1$

$\epsilon=max\hskip 3.0pt|\mathcal{R}|$

end while

Return Global Minimizer:
${\\{X^{\star}_{1},X^{\star}_{2}\\}=\\{X^{k}_{1},X^{k}_{2}\\}}$

####  III-B1 D-OPF for Two-Area System

First, we detail the proposed ENApp D-OPF algorithm using a two-area system.
The two-area model is then extended to obtain the generalized model for the
distribution system comprised of multiple areas. Suppose a network is composed
of two areas: Area 1 (UA) and Area 2 (DA), each with their own local
variables, $x_{1}$ and $x_{2}$ (see Fig. 1). Let
$\textbf{Y}=[y_{1},y_{2}]^{T}\in\mathcal{C}^{2}$ be the associated variable
for the shared bus. Note that Y is the complicating variable that couples the
optimization problems for the two areas. Here, $y_{1}$ and $y_{2}$ are the bus
voltage magnitude and the complex power flow through the bus shared between
Area 1 and Area 2, respectively.

Let the set of all local variables for Area 1 and Area 2 be $X_{1}$ and
$X_{2}$, respectively where, $X_{1}=\\{x_{1},y_{1}\\}{\hskip
5.69046pt}\text{\&{\hskip 5.69046pt}}X_{2}=\\{x_{2},y_{2}\\}$. Let $X$ be the
set of all problem variables where, $X=X_{1}\cup X_{2}$. The overall
centralized optimization problem is defined using (7), where, $S$ is the set
of constraints, and $F$ is a decomposable cost function.

$\small\min_{X\in S}{\hskip 5.69046ptF(X)}$ (7)

The problem defined in (7) can be written as (8).

$\small\min_{X_{1}\in S_{1},X_{2}\in S_{2}}{\hskip
5.69046pt}\Phi_{1}(X_{1},y_{2})+\Phi_{2}(X_{2},y_{1})$ (8)

where, $S_{1}$ and $S_{2}$ are the set of constraints on local variables for
Area 1 and Area 2, respectively. Here, $\Phi_{1}$, $\Phi_{2}$ are the cost
functions for the respective local subproblems.

Thus, the original problem defined in (8) can be decomposed into the following
two sub-problems.

$\displaystyle\min_{X_{1}\in S_{1}}{\hskip 5.69046pt}\Phi_{1}(X_{1},y_{2})$
(9a) $\displaystyle\min_{X_{2}\in S_{2}}{\hskip
5.69046pt\Phi_{2}(X_{2},y_{1})}$ (9b)

The subproblems defined in (9) are solved as the following. The proposed
distributed algorithm first solves the subproblems individually to obtain
respective local and complicating variables. Here, $y_{2}$ and $y_{1}$ are
kept fixed to solve subproblem (9a) and sub-problem (9b), respectively for
$X_{1}$ and $X_{2}$. The complicating variables, $y_{1}$ and $y_{2}$, are then
exchanged among the two areas. The process is repeated until the master
problem converges. The master problem is described using interface equations
in (10) and solved using the FPI method. Note that the mathematical relation
between $y_{1}$ and $y_{2}$ is formulated based on nonlinear power flow
equations that are described as the interface equations in (10). The master
problem converges when the change in complicating variable in successive
iterations is within a pre-specified tolerance. The proposed algorithm for the
two-area system is detailed in Algorithm 2 and can be easily extended to
multiple areas as detailed later in this section.  An overview figure
representing the working of the proposed algorithm is shown in Fig. 2.

Figure 2: Descriptive picture of Algorithm 2: ENApp Method.

$\displaystyle y_{1}-\psi_{1}(y_{2},x_{1})=0$ (10a) $\displaystyle
y_{2}-\psi_{2}(y_{1},x_{2})=0$ (10b)

Remark 3 (Equivalent Network Approximation): The functions $\psi_{1}$ and
$\psi_{2}$ in (10) project/map the local variables ($x_{1},x_{2}$) of higher
dimensionality, to the shared variables ($y_{1},y_{2}$) of lower
dimensionality that represent network equivalents. The explicit form of these
functions for the D-OPF problem is simply power flow equations. More
explicitly, $\psi_{1}$ maps the boundary bus voltage with the power injections
in the local area and the shared bus, and $\psi_{2}$ maps the boundary bus
power flow with the power injections in the local area and voltage at the
shared bus. Thus, the functions $\psi_{1}$ and $\psi_{2}$ are simply mapping
the shared bus variables for the UA and DA areas using respective power flow
equations for the local areas. Exchanging boundary variables that correspond
to network equivalents help achieve a faster convergence in fewer macro-
iterations.

The network equivalents obtained in this work for UA and DAs are driven by the
common practice assumptions for network equivalents in radial power
distribution systems. Specifically, the aggregated load assumption to
represent DA is a common practice in power systems analysis. It states that
the downstream load of a bus can be approximated as a lumped load. This
enables us to represent the voltage at the shared bus as
$y_{1}=\psi_{1}(y_{2},x_{1})$ where $y_{2}$ is the lumped load representing
the network equivalent of the DA (Area 2) for the UA (Area 1). Similarly, as
the power flow solutions are unique for radial distribution systems [27],
$y_{1}$ represents the network equivalent for the UA (Area 1). Thus, the
equivalent load at the shared bus, $y_{2}$, can be expressed as
$y_{2}=\psi_{2}(y_{1},x_{2})$, where $y_{1}$ is the bus voltage representing
the network equivalent of the UA (Area 1) as seen from the DA (Area 2).  Note
that these network equivalents are not explicitly calculated. Rather, they are
updated based on local optimization solutions at each macro-iterations.

Thus, instead of all local variables, $x_{1}$ and $x_{2}$, we only need the
shared variables, $y_{1}$ and $y_{2}$, during the macro-iterations to capture
and share the information among distributed areas. Moreover, since the shared
variables $y_{1}$ and $y_{2}$ represent electrical network equivalents for the
respective UA and DA, it is able to provide more meaningful information to the
respective regions during macro-iterations. This leads to a fewer number of
variable exchanges to achieve convergence.

Figure 3: A simple radial distribution system.

####  III-B2 D-OPF for a System with Multiple Areas

For a system comprised of multiple areas, the OPF problem for the loss
minimization objective is defined  for each $A_{m}\in\mathcal{A}_{R}$ as shown
in (11a), where $m\in\\{1,2,...,N_{A}\\}$. Here, the problem variables for
area $A_{m}$ are defined as $X_{m}=[P_{ij}$, $Q_{ij}$, $l_{ij}$, $v_{j}$,
$q_{Dj}]^{T}$ $\forall j\in\mathcal{N}_{m}$ and
$\forall(i,j)\in\mathcal{E}_{m}$. Recall that the reactive power dispatch from
DERs, $q_{D_{j}}$, define the control/decision variables for the optimization
problem. The constraints for the OPF problem for area $A_{m}$ are defined in
(11b)-(11h). Here, (11b)-(11h) characterize the constraint set
$\mathcal{S}_{m}$, for the subproblem defined for $A_{m}$. Then, Algorithm 2
is implemented for each UA-DA pair. Note that we model a dummy bus between
each connected area, which enables us to identify unique UA-DA boundaries.
Thus Algorithm 2 can be extended for general radial distribution systems with
multiple laterals. The convergence is achieved when all the change in the
shared variables for all area pairs is within a pre-specified tolerance.

$\displaystyle\text{(P3)}\hskip 17.07182pt\min\hskip 5.69046pt\sum
l_{ij}r_{ij}$ (11a) $\displaystyle\text{s.t.}\hskip
11.38092ptP_{ij}-r_{ij}l_{ij}-p_{L_{j}}+p_{Dj}=\sum_{k:j\rightarrow k}P_{jk}$
(11b) $\displaystyle Q_{ij}-x_{ij}l_{ij}-q_{L_{j}}+q_{Dj}=\sum_{k:j\rightarrow
k}Q_{jk}$ (11c) $\displaystyle
v_{j}=v_{i}-2(r_{ij}P_{ij}+x_{ij}Q_{ij})+(r_{ij}^{2}+x_{ij}^{2})l_{ij}$ (11d)
$\displaystyle-\sqrt{S_{DRj}^{2}-p_{Dj}^{2}}\leq
q_{Dj}\leq\sqrt{S_{DRj}^{2}-p_{Dj}^{2}}$ (11e)
$\displaystyle\underline{V^{2}}\leq v_{i}\leq\overline{V^{2}}$ (11f)
$\displaystyle{\color[rgb]{0,0,0}l_{ij}\leq\left(I^{rated}_{ij}\right)^{2}}$
(11g) $\displaystyle v_{i}l_{ij}=P_{ij}^{2}+Q_{ij}^{2}\hskip 17.07182pt$ (11h)

###  III-C Convergence Analysis

The network parameters play a vital role in the fast convergence of the
proposed distributed algorithm. Here, for a simple system described in Fig. 3,
the convergence at the boundary for solving the master problem (nonlinear
power flow) is analyzed. In the system, load $S_{L}=P_{L}+jQ_{L}$ is connected
to the node $u$ through a branch with $z_{d}=r_{d}+jx_{d}$ series impedance. A
stiff voltage source at node $0$ with a squared voltage magnitude of $v_{0}$
is connected to the node $u$ through a series impedance of
$z_{u}=r_{u}+jx_{u}$. The power flow from the voltage source to the system is
$S_{0}=P_{0}+jQ_{0}$. Also, the power flow from node $u$ is denoted by
$S_{u}=P_{u}+jQ_{u}$. For the system (see Fig. 3), using the branch flow
equations, the squared voltage magnitude of the boundary bus ($u$ node) can be
written as (12a).

$\displaystyle
v_{u}=v_{0}-2(r_{u}P_{0}+x_{u}Q_{0})+z^{2}_{u}\frac{P_{0}^{2}+Q_{0}^{2}}{v_{0}}$
(12a)

Assumption 1: To show convergence at the boundary bus, we assume that the
$S_{0}$, $S_{u}$ and $S_{L}$ are in the order of $10^{0}\sim 10^{-1}$;
$|z_{u}|$ and $|z_{d}|$ are in the order of $10^{-2}$ or less. Note that these
are valid assumptions for a practical power distribution system.

Lemma 1: $P_{0}^{2}+Q_{0}^{2}$ can be approximated by $P_{u}^{2}+Q_{u}^{2}$.
Also, $P_{0}$ and $Q_{0}$ are in the same order as $P_{u}$ and $Q_{u}$,
respectively.

###### Proof.

By applying branch flow equation at the upstream bus, (13a) can be written in
terms of the line impedance and power injection at node $u$. The line loss is
relatively small compared to the real and reactive power flow in the branch.
By ignoring the line loss term, (13a) can be written as (14a).

$\displaystyle P_{0}^{2}+Q_{0}^{2}$ $\displaystyle=$ $\displaystyle
P_{u}^{2}+Q_{u}^{2}+\Bigg{\\{}\frac{P_{0}^{2}+Q_{0}^{2}}{v_{0}^{2}}[r_{u}^{2}+x_{u}^{2}]\Bigg{\\}}$
(13a) $\displaystyle+$
$\displaystyle\Bigg{\\{}\frac{2(P_{0}^{2}+Q_{0}^{2})}{v_{0}}[r_{u}P_{u}+x_{u}Q_{u}]\Bigg{\\}}$
(13b) $\displaystyle\frac{P_{0}^{2}+Q_{0}^{2}}{v_{0}}$ $\displaystyle=$
$\displaystyle\frac{P_{u}^{2}+Q_{u}^{2}}{v_{0}-2\\{r_{u}P_{u}+x_{u}Q_{u}\\}}$
(14a)

Since $v_{0}\approx 1$ and branch power losses are relatively smaller than the
branch power flow for a well-designed power system [18],
$v_{0}>>2\\{r_{u}P_{u}+x_{u}Q_{u}\\}$, and $P_{0}^{2}+Q_{0}^{2}\approx
P_{u}^{2}+Q_{u}^{2}$. ∎

Again, for the downstream area, assuming that the branch power losses are
relatively less compared to the branch power flow [18], we can express $P_{u}$
and $Q_{u}$ in terms of $P_{L}$ and $Q_{L}$ as (15).

$\displaystyle P_{u}=P_{L}+r_{d}\frac{P_{L}^{2}+Q_{L}^{2}}{v_{u}}\hskip
5.69046pt\text{and}\hskip
5.69046ptQ_{u}=Q_{L}+x_{d}\frac{P_{L}^{2}+Q_{L}^{2}}{v_{u}}$ (15)

Using Lemma 1, (12a), and (15), $v_{u}$ can be written as (16a). Here, at $k$
iteration, the voltage update step for the shared bus of UA uses the power
injection at the same bus obtained from the DA at $({k-1})^{th}$ iteration.
After incorporating the iteration indices in (16a), we obtain (LABEL:Con_3).

$\displaystyle
v_{u}=v_{0}-2(r_{u}P_{u}+x_{u}Q_{u})+z^{2}_{u}\frac{P_{u}^{2}+Q_{u}^{2}}{v_{0}}$
(16a) $\displaystyle v_{u}^{k}$ $\displaystyle=$ $\displaystyle
v_{0}-2(r_{u}P_{u}^{k-1}+x_{u}Q_{u}^{k-1})+z^{2}_{u}\frac{(P_{u}^{k-1})^{2}+(Q_{u}^{k-1})^{2}}{v_{0}}$
(17a) $\displaystyle=$ $\displaystyle
v_{0}-2\left\\{r_{u}P_{L}+x_{u}Q_{L}+(r_{u}r_{d}+x_{u}x_{d})\frac{S_{L}^{2}}{v_{u}^{k-1}}\right\\}$
$\displaystyle+\frac{z_{u}^{2}}{v_{0}}\left\\{S_{L}^{2}+\frac{z_{d}^{2}S_{L}^{4}}{(v_{u}^{k-1})^{2}}+\frac{2S_{L}^{2}}{v_{u}^{k-1}}\Big{(}r_{d}P_{L}+x_{d}Q_{L}\Big{)}\right\\}$
$\displaystyle\approx$ $\displaystyle
v_{0}-2(r_{u}P_{L}+x_{u}Q_{L})+\frac{z_{u}^{2}S_{L}^{2}}{v_{0}}+\frac{2S_{L}^{2}}{v_{u}^{k-1}}\Big{\\{}\frac{z_{u}^{2}}{v_{0}}(r_{d}P_{L}+...$
$\displaystyle x_{d}Q_{L})-(r_{u}r_{d}+x_{u}x_{d})\Big{\\}}\hskip
8.5359pt\text{ (ignoring $z_{u}^{2}z_{d}^{2}$ term) }$

Next, we define
$C_{1}=v_{0}-2(r_{u}P_{L}+x_{u}Q_{L})+\frac{z_{u}^{2}S_{L}^{2}}{v_{0}}$ and
$C_{2}=2S_{L}^{2}\big{\\{}\frac{z_{u}^{2}}{v_{0}}(r_{d}P_{L}+x_{d}Q_{L})-(r_{u}r_{d}+x_{u}x_{d})\big{\\}}$.
Then the voltage at the shared bus at $k^{th}$ iteration, $v_{u}^{k}$, can be
written as a function of the bus voltage at the previous iteration,
$v_{u}^{k-1}$, (18a).

$\displaystyle v_{u}^{k}=C_{1}+\frac{C_{2}}{v_{u}^{k-1}}$ (18a)

The speed of the convergence of the sequence $\\{v_{u}^{k}\\}$, defined by
(18a), depends on the $\gamma=\frac{C_{2}}{C_{1}}$. Here, for any power
distribution system, $|\gamma|$ is bounded between $(0,1)$. The closer the
$|\gamma|$ is to the zero, the faster the sequence converges. The convergence
is guaranteed regardless of the value of $\gamma$.

The ratio between the difference of two consecutive terms in the sequence
$\\{v_{u}^{k}\\}$ can be expressed as (19a) and reduced to (20a) after some
algebraic manipulations.

$\displaystyle\frac{\left|{v_{u}^{k+2}-v_{u}^{k+1}}\right|}{\left|v_{u}^{k+1}-v_{u}^{k}\right|}=\frac{\left|C_{1}+\frac{C_{2}}{v_{u}^{k+1}}-v_{u}^{k+1}\right|}{\left|v_{u}^{k+1}-\frac{C_{2}}{v_{u}^{k+1}-C_{1}}\right|}=\frac{\left|v_{u}^{k+1}-C_{1}\right|}{\left|v_{u}^{k+1}\right|}$
(19a)

$\frac{\left|{v_{u}^{k+2}-v_{u}^{k+1}}\right|}{\left|v_{u}^{k+1}-v_{u}^{k}\right|}=\frac{\left|C_{1}+\frac{C_{2}}{v_{u}^{k}}-C_{1}\right|}{\left|C_{1}+\frac{C_{2}}{v_{u}^{k}}\right|}=\frac{\left|\frac{C_{2}}{C_{1}}\right|}{\left|v_{u}^{k}+\frac{C_{2}}{C_{1}}\right|}$
$\displaystyle\hskip
22.76228pt=\frac{\left|\gamma\right|}{\left|v_{u}^{k}+\gamma\right|}<1$ (20a)

From (20a), it is clear that the convergence is guaranteed for (18a) since
$v_{u}^{k}\approx 1$ and $|\gamma|$ is bounded between $(0,1)$. The
convergence is linear as there exists some $\delta\in(0,1)$ for which
$\frac{\left|\gamma\right|}{\left|v_{u}^{k}+\gamma\right|}\leq\delta$. Also,
for any practical power distribution system system, the line parameters such
as $r_{d}$, $x_{d}$ etc. are much smaller (in the order of $10^{-1}$ pu or
smaller) than the allowable voltage magnitude ($\sim 1$ pu). Thus, $|\gamma|$
is at-most 0.1 for power distribution network. The rate of convergence depends
on $|\gamma|$; the lower the $|\gamma|$, the higher the rate of convergence.
This property results in a faster convergence of the proposed algorithm at the
boundary bus requiring a significantly fewer number of macro-iterations or
communication-rounds among distributed agents.

Fig. 3 shows an equivalent representation of the two-area problem. The
upstream and downstream branch can be thought as an equivalent circuit model
of the UA and DA, respectively. Since in the ENApp D-OPF model, convergence is
searched at the boundary bus of the physically neighboring areas, this
analysis is valid for a multi-area system as well.

##  IV Results and Discussions

The proposed distributed OPF approach is thoroughly validated using IEEE
123-bus and IEEE-8500 node test systems. All experiments were simulated in
Matlab 2018b on a machine with 16GB memory and Core i7-8550U CPU @1.80GHz. The
NLP problems in both C-OPF and D-OPF models are solved using a commercial NLP
solver, Artelys Knitro [28].

###  IV-A Simulated System

The simulations are conducted using the following four test systems: (1) IEEE
123 bus system with 3 DERs (TS-1), (2) IEEE 123 bus with 9 DERs (TS-2), (3)
IEEE 8500 node test system with 113 DERs - about $10\%$
($|\mathcal{N}_{D}|/|\mathcal{N}_{L}|=0.1$) DER penetration (TS-3),and (4)
IEEE 8500 node test system with 544 DERs - about $50\%$
($|\mathcal{N}_{D}|/|\mathcal{N}_{L}|=0.5$) DER penetration (TS-4).  For IEEE
123-bus system, the kVA rating of the DERs is 60 kVA. For IEEE 8500 node test
system, the DER sizes are chosen randomly with a kVA rating ranging from 1.3
to 5.8 kVA. The simulation parameters associated with the respective test
systems are shown in Table I. Here, TS-3 and TS-4 are used to test the
scalability of the proposed ENApp D-OPF method for larger distribution
systems. We have simulated a single-phase equivalent circuit models for both
IEEE 123-bus and IEEE 8500-bus test systems. It is assumed that the test
systems are already decomposed and comprised of multiple areas. The
decomposition of both test systems and the connections between areas are shown
in Fig. 4 and 5. TS-1 and TS-2 are comprised of four areas (Fig. 4) and TS-3
and TS-4 are comprised of 17 areas (see Fig. 5).

Figure 4: TS-1 and TS-2 Test System Figure 5: TS-3 and TS-4 Test System

The proposed distributed OPF algorithm is validated for optimality upon
comparing against an equivalent centralized OPF model. We also compare the
compute time requirements for the proposed distributed OPF and the centralized
OPF models. Centralized OPF is formulated using the nonlinear branch flow
model as detailed in (P1) (see Section II.B).

TABLE I:  Summary for Simulated Systems

Test System | Area | Node | Total Variables | Constraints | DERs
---|---|---|---|---|---
TS-1 | 4 | 128 | 512 | 898 | 3
TS-2 | 4 | 128 | 518 | 910 | 9
TS-3 | 17 | 2522 | 10198 | 17876 | 113
TS-4 | 17 | 2522 | 10692 | 18738 | 544

###  IV-B IEEE 123-bus Test System

In this section, we present the results for the IEEE 123-bus system using TS-1
(with 3 DERs) and TS-2 (with 9 DERs). The OPF problem is solved for the loss
minimization problem. The results obtained from the proposed D-OPF method are
compared with the C-OPF problem defined in (P1). Here, we use a small test
system to demonstrate the convergence properties of the proposed D-OPF
algorithm.

For the 123-bus test system, the value of the objective function for the loss
minimization optimization problem and time required for C-OPF and ENApp D-OPF
to converge are compared in Table II. The table also shows the number of
macro-iterations it takes for the proposed algorithm to converge. For the
D-OPF approach, the compute time is the total time it takes for the D-OPF to
converge given each area is computing the local optimization solutions in
parallel. Note that the proposed approach converges to the same optimal
solution as that obtained using the C-OPF algorithm (see Table II). Thus the
proposed approach is able to achieve the optimal solution. Also, C-OPF takes
4.3 seconds to converge. D-OPF takes 4 macro-iterations and 3 seconds to
converge. For the small system, while the computational advantages are not
apparent, the major advantage is in the number of macro-iterations that ENApp
D-OPF takes to converge to the optimal solution; it takes 4 macro-iterations
to solve the non-linear OPF problem for TS-1 and TS-2. Notice that the
existing literature reports that the state-of-the-art distributed optimization
methods take hundreds of macro-iterations to converge for the IEEE 123-bus
test system [5, 6, 13].  We further provide additional comparison case studies
to better justify this statement (see Section IV.D).

TABLE II: C-OPF vs. ENApp D-OPF for IEEE 123-bus Test System Test System | Total Loss (kW) | Time (sec) | #Iterations
---|---|---|---
C-OPF | D-OPF | C-OPF | D-OPF | D-OPF
TS-1 | 21.7 | 21.9 | 4.3 | 3 | 4
TS-2 | 12.6 | 12.6 | 7.94 | 6.47 | 4

Next, the value of the objective function and the maximum residuals for
boundary variables for each macro-iteration of the D-OPF are shown for each
test system in Fig. 6. It can be observed that the D-OPF takes less than 4
iterations to converges to the optimal solutions for the IEE 123-bus test
system, see Fig. 6(a). The residual for the boundary variables also rapidly
decreases, see Fig. 6(b). The nodal voltages obtained using C-OPF and ENApp
D-OPF method are also compared and shown in Fig. 7. It can be observed that
the ENApp D-OPF converges to the same optimal solution as obtained by solving
the centralized OPF problem. Also, the difference in the nodal voltages
obtained using C-OPF and the proposed D-OPF is $<0.1\%$.

Figure 6: ENApp D-OPF Convergence Properties. At each macro-iteration (a)
Objective function; and (b) Maximum residual of boundary variables.

Figure 7: Comparison of Node Voltages for proposed ENApp D-OPF and C-OPF
Method, for (a) TS-1, (b) TS-2

###  IV-C IEEE 8500-node System

We demonstrate the scalability of the proposed approach using the single-phase
equivalent of IEEE 8500-node system with 2522 buses. We simulate two test
scenarios using the IEEE 8500-node test system: TS-3 with 113 DERs and TS-4
with 544 DERs. Note that this test case comprises of a large system for OPF
studies with a couple of thousand decision variables. Typically, the existing
distributed OPF algorithms have not been tested for a system of this size and
with these many DERs. Further, note that while the non-linear C-OPF (P1)
problem is easily solved for the small test system such as IEEE 123 bus
system, the commercial NLP solvers such as Artelys Knitro [28] fail to
converge for the selected large test system. In fact, one of the major
disadvantages of the large-scale centralized optimal power flow (C-OPF)
problems is the long compute time or failure to achieve convergence. Thus, for
validation, we compute optimal results using the centralized SOCP relaxation
of the NLP problem as detailed in (P2), denoted by C-OPF∗. Note that while, in
general, SOCP solutions may not be power flow feasible [19], they are found to
exact for the loss minimization problems under certain conditions [29]. In
this case, we explicitly validate that the SOCP solution is also power flow
feasible. Thus, for this case, the SOCP solution provides a suitable optimal
solution for comparison. We use cvx package to solve the relaxed SOCP problem
[30].

TABLE III: C-OPF vs. ENApp D-OPF for IEEE 8500-node Test System Test System | Total Loss (kW) | Time (sec) | #Iterations
---|---|---|---
C-OPF∗ | D-OPF | C-OPF | D-OPF | D-OPF
TS-3 | 249 | 249.1 | NA | 180 | 17
TS-4 | 152 | 152.8 | NA | 700 | 17

The results are shown in Table III. It can be observed that both C-OPF and
D-OPF converge to the same optimal solution; the error for the value of the
minimized objective function for the ENApp D-OPF method is less than 1%. The
compute time is not reported for the C-OPF problem as the NLP model did not
converge for the 8500-node test system. Therefore, we report the optimal
values obtained using the relaxed SOCP model for the C-OPF problem for
comparison. On the contrary, the proposed D-OPF algorithm is able to solve the
nonlinear programming model for the OPF detailed in (P3) and is able to
achieve convergence relatively fast for the large system. Furthermore, the
major advantage lies in the number of macro-iterations it takes to converge
for the 8500-node test system. Even for the large test system, the proposed
D-OPF achieves convergence using a relatively smaller number of macro-
iterations, here, 17 macro-iterations. Thus, unlike other approaches in the
literature, the number of macro-iterations increases sub-linearly with the
number of decomposed areas in the test system. To the knowledge of the
authors, no other distributed optimization methods reported in the literature
can solve the non-linear OPF problem and reach the optimum solution in fewer
macro-iterations for large-scale distribution systems (TS-3 or TS-4). The
existing methods take hundreds of rounds of communication at each step of the
optimization problem for large-scale distribution networks [5, 6]. Thus, this
test case demonstrates the scalability of the proposed distributed
optimization method for non-linear OPF problems.

The value of the objective function and the maximum residuals for boundary
variables with respect to the number of macro-iterations for the simulated
test systems are shown in Fig. 6. Note that the computation starts with a flat
start and still converges within a reasonable number of iterations. Further,
the maximum residual vector rapidly decreases with the iteration count. Next,
the nodal voltages obtained using C-OPF and ENApp D-OPF method are compared in
Fig. 8. From the figure, it is clear that the proposed ENApp D-OPF can solve
the original optimization problem regardless of the network size, and the
error in voltage magnitude is very low ($<0.1\%$).

Figure 8: Comparison of Node Voltages for proposed ENApp D-OPF and C-OPF
Method, for (a) TS-3, (b) TS-4

For completeness, we validate the feasibility of the SOCP solutions for the
simulation test cases TS-3 and TS-4. Note that the SOCP relaxation is exact if
the solution for the relaxed problem satisfies the equality conditions for the
relaxed equality constraint in (P2). That is,
$e_{ij}=v_{i}l_{ij}-(P_{ij}^{2}+Q_{ij}^{2})\approx 0,$
${\forall{ij}\in\mathcal{E}}$. It is observed that $e_{ij}$ is in the order of
$10^{-4}$ upon solving SOCP relaxed C-OPF problem for the IEEE 8500-node test
system (for both test systems, TS-3 & TS-4). Thus, the SOCP solution is exact
for the NLP OPF problem and can be used to benchmark the optimum solutions
obtained from the proposed D-OPF algorithm.

Figure 9: Convergence Properties for ADMM method [5] for TS-1 and TS-3 TABLE
IV: Comparison of ENApp D-OPF and ADMM Based Method

System | Total Loss (kW) | Time (sec) | #Iterations
---|---|---|---
ENApp | ADMM | ENApp | ADMM | ENApp | ADMM
TS-1 | 21.9 | 21.8 | 3 | 145 | 4 | 272
TS-3 | 249.1 | 249.1 | 180 | 6120 | 17 | 756

###  IV-D Comparison with ADMM

As stated before, the state-of-the-art distributed optimization approaches
such as, ADMM based methods require a large number of communication rounds
(macro-iterations) among distributed agents, which is in the order of
$~{}10^{2}$ for mid-size networks. Next, we compare the proposed ENApp D-OPF
method with another ADMM based technique adopted from [5] for TS-1 and TS-3.
The convergence of the ADMM approach is demonstrated using the maximum
residual at the boundary bus for TS-1 and TS-3 (see Fig. 9). The two
distributed methods are also compared for solution quality and compute time in
Table IV. Although the ADMM-based method converges to the same solution as the
proposed ENApp, it takes around 272 macro-iterations for the mid-sized IEEE
123-bus system (TS-1) and 756 macro-iterations for the large IEEE 8500-node
feeder (TS-3). Whereas, the proposed method requires 4 and 17 macro-iterations
for IEEE 123-bus (TS-1) and IEEE 8500-node (TS-3), respectively to converge to
a solution of the same quality. The required compute time to converge is also
very high for the ADMM-based method.

###  IV-E Extension to a Three-Phase Unbalanced Radial Distribution Feeder

To showcase the efficacy of the proposed ENApp D-OPF algorithm for the
unbalanced distribution network, we simulated a 3-phase unbalanced IEEE-13 bus
system (TS-5) using our proposed method (See Fig. 10 in Appendix A.). The
local nonlinear unbalanced OPF model (P4) is described in Appendix A. The
network is assumed to be composed of 2 areas with 3 DERs at each area. The
result of the 3-phase unbalanced distribution system is shown in Table V. From
the table, we can see that all the 6 DERs have the same optimal Q injections
for both proposed ENApp D-OPF and an equivalent C-OPF algorithm. Thus, the
proposed approach is able to obtain an optimal solution of the same quality as
the C-OPF model for an unbalanced distribution feeder. Also, both methods
achieve the same minimum value for the active power losses of are 44.5 kW. The
ENApp D-OPF converges within 2 iterations for this unbalanced system. Thus,
the proposed ENApp D-OPF is also as efficient as for single-phase systems.

TABLE V: C-OPF vs. ENApp D-OPF for Unbalanced System (TS-5) | Bus |  C-OPF† |  D-OPF†
---|---|---|---
Q Injection | 3 (A) | 66.3 kVAr | 66.3 kVAr
3 (B) | 66.3 kVAr | 66.3 kVAr
3 (C) | 66.3 kVAr | 66.3 kVAr
7 (A) | 132.7 kVAr | 132.7 kVAr
7 (B) | 132.7 kVAr | 132.7 kVAr
7 (C) | 132.7 kVAr | 132.7 kVAr
Total Loss | 44.5 kW | 44.5 kW
Macro-Iteration for D-OPF† | 2

###  IV-F Discussions

The advantages of the proposed ENApp D-OPF approach are briefly summarized
here. First, the fewer rounds of communication required among neighbors in the
proposed D-OPF model make it practically feasible and relatively less
expensive to deploy than popular D-OPF methods (that take several hundred
macro-iterations). Besides, the fewer rounds of communication among
neighboring areas introduce lesser and manageable communication delays. The
sub-linear increase in the iteration count with the increase in the number of
the decomposed area also bodes well for scaling the approach for any size of
the distribution system. Further, the proposed method does not require any
parameter tuning to achieve convergence. If the original problem has the
convergence guarantees, then the distributed version is also guaranteed to
converge. This statement is also true for nonlinear optimization problems.
Moreover, lower observability requirements improve overall system security.
Finally, the proposed method is easily scalable to larger distribution systems
without losing any of the advantages related to convergence.

Next, we briefly detail the limitations of the proposed algorithm. As detailed
before, the proposed approach actively leverages the radial topology of the
distribution feeder. Thus, the developed algorithm, in its current form, is
not applicable to general meshed feeders. A particular case of meshed
topology, where each area is radially connected while sub-areas are meshed,
can be solved using the proposed approach. The expansion of the proposed
methods for weekly meshed systems is a part of future work.

Finally, it is to be noted that advanced communication infrastructure is
required to implement the proposed D-OPF method for the distribution systems.
Although the current distribution feeders do not support such advanced
communication capabilities, moving forward with the integration of ADMS and
push from the utilities to better manage their distribution systems and
distributed energy resources (DERs), it is only a matter of time when such
information will be available and can be used in everyday optimization of the
distribution systems [31, 32].

##  V Conclusions

This paper presents a novel distributed optimization algorithm to solve the
optimal power flow (OPF) problem for a radial distribution system. The
proposed approach is based on the fundamental properties of power flow in
radial power distribution systems that allow us to represent the decomposed
radial distribution systems as equivalent voltage sources (for the upstream
node) and aggregated loads (for the downstream nodes). We term this as the
equivalent network approximation based distributed OPF (ENApp D-OPF)
algorithm. Upon leveraging this property, the number of macro-
iterations/communication-rounds required for the convergence of the
distributed algorithm is significantly reduced. The proposed approach is
validated using simulations on the single-phase equivalent of IEEE 123-bus and
IEEE 8500-node test systems and shown to converge to the same optimal
solutions as obtained by the centralized OPF (C-OPF) method. Further, for the
IEEE 8500-node test system, while the nonlinear model for the C-OPF cannot be
solved using commercial solvers, the proposed D-OPF algorithm solves the
nonlinear programming problem and converges to the optimum solution within a
reasonable time.  The proposed method’s applicability for a radial unbalanced
three-phase system is also demonstrated using a small unbalanced distribution
feeder. The extension of the proposed method for weakly meshed topologies is
currently under investigation and is a part of our future research.

##  Appendix A Unbalanced 3-phase Distributed OPF Model

This section details the D-OPF formulation for an unbalanced power
distribution system. The power flow equations are based on our prior work on
unbalanced D-OPF [21]. A part of radial power distribution network (Area
$A_{m}$) of $n$ buses is considered where, $\mathcal{N}_{m}=\\{1,2,…,n\\}$
denotes the set of buses in that system. Assume $\mathcal{E}_{m}$ denotes set
of edges identifying distribution lines that connect the ordered pair of buses
$(i,j)$, $\forall$ $i,j\in\mathcal{N}_{m}$. The three-phases of the system are
denoted by $\\{a,b,c\\}$ and $\phi_{i}$ denotes the set of phases in the bus
$i$. Let $v_{j}^{p}=|V_{j}^{p}|^{2}$ be the squared magnitude of voltage at
bus $j$ for phase $p\in\phi_{j}$. Define $\phi_{ij}=\\{pq:p\in\phi_{i}$ and
$q\in\phi_{j}$ $\forall(i,j)\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}+jx^{pq}_{ij}$, where $pq\in\phi_{ij}$.
$S_{Lj}^{p}=p_{Lj}^{p}+jq_{Lj}^{p}$ denote the load connected at node $j$ of
phase $p\in\phi_{i}$. $\delta^{pq}_{ij}$ is the angle difference between the
phase currents. $\mathbb{Re}$ denotes the real part of the complex number. For
detail on power flow model please refer to [21]. The OPF problem is detailed
next. The representative IEEE 13-bus feeder decomposed into two areas is given
in Fig. 10.

$\small\text{(P4)}\hskip
17.07182pt\text{Minimize}\sum_{p\in\phi_{j},j:i\rightarrow
j}{l^{p}_{ij}r^{p}_{ij}}$ s.t. $\displaystyle
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)}$
$\displaystyle=\sum_{k:j\rightarrow k}P_{jk}^{pp}+p_{Lj}^{p}-p_{Dj}^{p}\hskip
91.04872pt$ $\displaystyle
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)}$
$\displaystyle=\sum_{k:j\rightarrow
k}Q_{jk}^{pp}+q_{Lj}^{p}-q_{Dj}^{p}-q_{Cj}^{p}\hskip 91.04872pt$
$\displaystyle
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}}$
$\displaystyle+\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]}$
$\displaystyle(P_{ij}^{pp})^{2}+(Q_{ij}^{pp})^{2}=v_{i}^{p}l_{ij}^{pp}$
$\displaystyle(l_{ij}^{pq})^{2}=l_{ij}^{pp}l_{ij}^{qq}$ $\displaystyle
q_{Dj}^{p}\leq\sqrt{(s_{DR}^{p})^{2}-(p_{Dj}^{p})^{2}}$ $\displaystyle
q_{Dj}^{p}\geq-\sqrt{(s_{DR}^{p})^{2}-(p_{Dj}^{p})^{2}}$
$\displaystyle\underline{V}^{2}\leq v_{j}^{p}\leq\overline{V}^{2}$

Figure 10: TS-5: IEEE 13-bus test system with two areas.

## References

  * [1] J. A. Momoh, R. Adapa, and M. El-Hawary, “A review of selected optimal power flow literature to 1993. i. nonlinear and quadratic programming approaches,” IEEE transactions on power systems, vol. 14, no. 1, pp. 96–104, 1999.
  * [2] A. Castillo and R. P. O’Neill, “Survey of approaches to solving the ACOPF,” US Federal Energy Regulatory Commission, Tech. Rep, 2013.
  * [3] D. K. Molzahn, F. Dörfler, H. Sandberg, S. H. Low, S. Chakrabarti, R. Baldick, and J. Lavaei, “A survey of distributed optimization and control algorithms for electric power systems,” IEEE Transactions on Smart Grid, vol. 8, no. 6, pp. 2941–2962, 2017.
  * [4] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al., “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends® in Machine learning, vol. 3, no. 1, pp. 1–122, 2011\.
  * [5] W. Zheng, W. Wu, B. Zhang, H. Sun, and Y. Liu, “A fully distributed reactive power optimization and control method for active distribution networks,” IEEE Transactions on Smart Grid, vol. 7, no. 2, pp. 1021–1033, 2015.
  * [6] Q. Peng and S. H. Low, “Distributed optimal power flow algorithm for radial networks, i: Balanced single phase case,” IEEE Transactions on Smart Grid, vol. 9, no. 1, pp. 111–121, 2016.
  * [7] T. Erseghe, “Distributed optimal power flow using admm,” IEEE transactions on power systems, vol. 29, no. 5, pp. 2370–2380, 2014.
  * [8] E. Dall’Anese, H. Zhu, and G. B. Giannakis, “Distributed optimal power flow for smart microgrids,” IEEE Transactions on Smart Grid, vol. 4, no. 3, pp. 1464–1475, 2013.
  * [9] G. Qu and N. Li, “Harnessing smoothness to accelerate distributed optimization,” IEEE Transactions on Control of Network Systems, vol. 5, no. 3, pp. 1245–1260, 2017.
  * [10] S. Mhanna, G. Verbič, and A. C. Chapman, “Adaptive admm for distributed ac optimal power flow,” IEEE Transactions on Power Systems, vol. 34, no. 3, pp. 2025–2035, 2018.
  * [11] B. H. Kim and R. Baldick, “A comparison of distributed optimal power flow algorithms,” IEEE Transactions on Power Systems, vol. 15, no. 2, pp. 599–604, 2000.
  * [12] W. Lu, M. Liu, S. Lin, and L. Li, “Fully decentralized optimal power flow of multi-area interconnected power systems based on distributed interior point method,” IEEE Transactions on Power Systems, vol. 33, no. 1, pp. 901–910, 2017.
  * [13] G. Cavraro and R. Carli, “Local and distributed voltage control algorithms in distribution networks,” IEEE Transactions on Power Systems, vol. 33, no. 2, pp. 1420–1430, 2017.
  * [14] A. Bernstein and E. Dall’Anese, “Real-time feedback-based optimization of distribution grids: A unified approach,” IEEE Transactions on Control of Network Systems, vol. 6, no. 3, pp. 1197–1209, 2019.
  * [15] S. Magnússon, G. Qu, and N. Li, “Distributed optimal voltage control with asynchronous and delayed communication,” IEEE Transactions on Smart Grid, 2020.
  * [16] M. E. Baran and F. F. Wu, “Optimal capacitor placement on radial distribution systems,” IEEE Transactions on power Delivery, vol. 4, no. 1, pp. 725–734, 1989.
  * [17] D. K. Molzahn, I. A. Hiskens, et al., “A survey of relaxations and approximations of the power flow equations,” Foundations and Trends® in Electric Energy Systems, vol. 4, no. 1-2, pp. 1–221, 2019.
  * [18] M. Farivar and S. H. Low, “Branch flow model: Relaxations and convexification—part i,” IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 2554–2564, 2013.
  * [19] Q. Li and V. Vittal, “Non-iterative enhanced sdp relaxations for optimal scheduling of distributed energy storage in distribution systems,” IEEE Transactions on Power Systems, vol. 32, no. 3, pp. 1721–1732, 2016.
  * [20] R. R. Jha and A. Dubey, “Exact distribution optimal power flow (d-opf) model using convex iteration technique,” in 2019 IEEE Power Energy Society General Meeting (PESGM), pp. 1–5, 2019.
  * [21] R. R. Jha, A. Dubey, C.-C. Liu, and K. P. Schneider, “Bi-level volt-var optimization to coordinate smart inverters with voltage control devices,” IEEE Transactions on Power Systems, vol. 34, no. 3, pp. 1801–1813, 2019\.
  * [22] C. Langbort, L. Xiao, R. D’Andrea, and S. Boyd, “A decomposition approach to distributed analysis of networked systems,” in 43rd IEEE Conference on Decision and Control, vol. 4, pp. 3980–3985, IEEE, 2004.
  * [23] Y. Xu and W. Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” SIAM Journal on imaging sciences, vol. 6, no. 3, pp. 1758–1789, 2013.
  * [24] J. M. Ortega and W. C. Rheinboldt, Iterative solution of nonlinear equations in several variables, vol. 30.  Siam, 1970.
  * [25] A. Beck and L. Tetruashvili, “On the convergence of block coordinate descent type methods,” SIAM journal on Optimization, vol. 23, no. 4, pp. 2037–2060, 2013.
  * [26] Y. Xu and W. Yin, “A globally convergent algorithm for nonconvex optimization based on block coordinate update,” Journal of Scientific Computing, vol. 72, no. 2, pp. 700–734, 2017.
  * [27] H. D. Chiang and M. E. Baran, “On the existence and uniqueness of load flow solution for radial distribution power networks,” IEEE Transactions on Circuits and Systems, vol. 37, no. 3, pp. 410–416, 1990.
  * [28] J. Nocedal, “Knitro: An integrated package for nonlinear optimization,” in Large-Scale Nonlinear Optimization, pp. 35–60, Springer, 2006.
  * [29] M. Farivar and S. H. Low, “Branch flow model: Relaxations and convexification—part ii,” IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 2565–2572, 2013.
  * [30] M. Grant and S. Boyd, “Cvx: Matlab software for disciplined convex programming, version 2.1,” 2014.
  * [31] Y. Ngo, H. Arant, C. Bilby, and A. Grice, “Investing for the future: How small utilities are finding success with advanced distribution management systems,” IEEE Power and Energy Magazine, vol. 18, no. 1, pp. 34–42, 2020\.
  * [32] E. Boardman, “Advanced applications in an advanced distribution management system: Essentials for implementation and integration,” IEEE Power and Energy Magazine, vol. 18, no. 1, pp. 43–54, 2020.


#  p-forms from syzygies

Alan Muniz  IMECC - UNICAMP
Departamento de Matemática
Rua Sérgio Buarque de Holanda, 651
13083-970 Campinas-SP, Brazil<EMAIL_ADDRESS>
(Date: December, 2022)

###### Abstract.

We describe how to use syzygies to produce homogeneous polynomial differential
$p$-forms with prescribed vanishing locus.

###### Key words and phrases:

Differential forms, Distributions, Singular scheme, Vanishing Locus, Syzygy

###### 2020 Mathematics Subject Classification:

14D20, 14J60, 14F06, 13D02

##  1\. Introduction

In Mathematics, being able to compute explicit examples is of great
importance. In the particular study of distributions on projective spaces,
such examples are governed by homogeneous polynomial differential forms which
induce exact sequences of the form

$0\longrightarrow F\longrightarrow{\rm
T}{\mathbb{P}^{n}}\stackrel{{\scriptstyle\omega}}{{\longrightarrow}}N\longrightarrow
0$

where $F$ is a rank $n-p$ reflexive sheaf and $N$ is a rank $p$ torsion-free
sheaf. Properties of the vanishing locus of the $p$-form $\omega$ reflect on
properties of the sheaves $F$ and $N$; thus providing the $p$-form leads to
the understanding of the sheaves. In this direction we prove the following
result. Fix $k$ a field.

######  Theorem A.

Let $Z\subset\mathbb{P}^{n}$ a closed subscheme with (saturated) homogeneous
ideal $I_{Z}\subset R=k[x_{0},\dots,x_{n}]$ and let $\mathcal{A}^{p}(Z)$ be
the $R$-module of polynomial $p$-forms vanishing on $Z$. Then we have an exact
sequence of graded $R$-modules:

$0\longrightarrow I_{Z}\otimes\langle x_{i}dx_{j}-x_{j}dx_{i}\rangle_{0\leq
i<j\leq
n}\longrightarrow\mathcal{A}^{p}(Z)\longrightarrow\operatorname{Tor}_{p}^{R}(I_{Z},k)(p+1)\longrightarrow
0.$

Moreover, if $(I_{Z})_{d}=0$ then
$\mathcal{A}^{p}(Z)_{d}\cong\operatorname{Tor}_{p}^{R}(I_{Z},k)_{d+p+1}$.

Recall that $\operatorname{Tor}_{p}^{R}(I_{Z},k)$ can be identified with the
space of $p$-th syzygies of a minimal set of generators for $I_{Z}$. Given a
(minimal) free resolution

$0\longrightarrow
F_{n}\stackrel{{\scriptstyle\phi_{n}}}{{\longrightarrow}}F_{n-1}\stackrel{{\scriptstyle\phi_{n-1}}}{{\longrightarrow}}\dots\longrightarrow
F_{0}\stackrel{{\scriptstyle\phi_{0}}}{{\longrightarrow}}I_{Z}\longrightarrow
0,$

we construct a map
$\xi_{p}\colon\operatorname{Tor}_{p}^{S}(I_{Z},k)\rightarrow\mathcal{A}^{p}(Z)$
given by differentiating and combining the entries of
$\phi_{1},\dots,\phi_{p}$, i.e., $\xi_{p}$ depends on syzygies up to order
$p$, see Proposition 3.2.

The use of syzygies to describe distributions goes back at least to the work
of Campillo and Olivares [CO01], see also [CO21] and references therein. In
[CJM22, §4] the case of $1$-forms is essentially described, serving as a
prelude to the present work. Note that for $p=1$ Theorem A gives a slightly
more complete version of [CJM22, Proposition 4.5].

After recalling the basic notions pertaining to this work in §2, we prove in
§3 Proposition 3.1, from which Theorem A follows. Finally we provide some
examples in §4. Our computations come from an implementation of these routines
in Macaulay2 [GS].

### Acknowledgements

This work is supported by INCTmat/MCT/Brazil, CNPq grant number 160934/2022-2.
It originates from discussions with Maurício Corrêa, Marcos Jardim and Renato
Vidal to whom I am grateful.

##  2\. Preliminaries and notation

We begin by recalling some basic facts and establish the notation to be used
throughout the paper. Let $k$ be a field, that we may assume algebraically
closed of characteristic zero. Fix $V$ a $k$-vector space, $n:=\dim V-1$, and
let ${\mathbb{P}^{n}}=\mathbb{P}(V)$ the projective space of lines in $V$
through the origin. Let ${\rm T}{\mathbb{P}^{n}}$ and
$\Omega^{1}_{\mathbb{P}^{n}}$ denote, respectively, the tangent and cotangent
bundles/sheaves ${\mathbb{P}^{n}}$. We have the the Euler sequence

(1)
$0\longrightarrow{\mathcal{O}}_{\mathbb{P}^{n}}\stackrel{{\scriptstyle\operatorname{rad}}}{{\longrightarrow}}{\mathcal{O}}_{\mathbb{P}^{n}}(1)\otimes
V\longrightarrow{\rm T}{\mathbb{P}^{n}}\longrightarrow 0.$

Then we may identify $V=H^{0}({\rm T}{\mathbb{P}^{n}}(-1))$ the space of
constant vector fields (on the affine space $\mathbb{A}^{n+1}_{k}$). Fixing
homogeneous coordinates $(x_{0}:\dots:x_{n})$, we have that $V$ is spanned by
the derivations $\frac{\partial}{\partial
x_{0}},\dots,\frac{\partial}{\partial x_{n}}$ so that the map
$\operatorname{rad}$ in (1) is written as the inclusion of the radial vector
field:

$\operatorname{rad}=x_{0}\frac{\partial}{\partial
x_{0}}+\cdots+x_{n}\frac{\partial}{\partial x_{n}}.$

Dualizing (1) we get

(2)
$0\longrightarrow\Omega_{\mathbb{P}^{n}}^{1}\longrightarrow{\mathcal{O}}_{\mathbb{P}^{n}}(-1)\otimes
V^{\ast}\stackrel{{\scriptstyle\iota_{\operatorname{rad}}}}{{\longrightarrow}}{\mathcal{O}}_{\mathbb{P}^{n}}\longrightarrow
0,$

and considering $\\{dx_{0},\dots,dx_{n}\\}$ the basis of $V^{\ast}$ dual to
$\\{\frac{\partial}{\partial x_{0}},\dots,\frac{\partial}{\partial x_{n}}\\}$
we see that $\iota_{\operatorname{rad}}$ is the contraction – or interior
product – of (local) differential $1$-forms with the radial vector field:
$\iota_{\operatorname{rad}}(\omega)=\omega(\operatorname{rad})$. Furthermore,
we take exterior powers of (2) to arrive at

(3)
$0\longrightarrow\Omega_{\mathbb{P}^{n}}^{p}\longrightarrow{\mathcal{O}}_{\mathbb{P}^{n}}(-k)\otimes{\textstyle\bigwedge}^{p}V^{\ast}\stackrel{{\scriptstyle\iota_{\operatorname{rad}}}}{{\longrightarrow}}\Omega_{\mathbb{P}^{n}}^{p-1}\longrightarrow
0,$

where
$\iota_{\operatorname{rad}}\omega(v_{1},\dots,v_{p-1})=\omega(\operatorname{rad},v_{1},\dots,v_{p-1})$
is the contraction of $p$-forms with $\operatorname{rad}$.

Note that from (3) we have that global sections of
$\Omega_{\mathbb{P}^{n}}^{p}(d+p+1)$ are in bijection with homogeneous
differential $p$-forms

$\omega=\sum_{0\leq i_{1}<\dots<i_{p}\leq n}A_{i_{1}\dots
i_{p}}dx_{i_{1}}\wedge\dots\wedge dx_{i_{p}}$

satisfying $\iota_{\operatorname{rad}}\omega=0$ and $\deg A_{i_{1}\dots
i_{p}}=d+1$.

To fix notation, let $R$ denote the polynomial ring $R=k[x_{0},\dots,x_{n}]$
and let $\Omega^{p}_{R}=R\otimes{\textstyle\bigwedge}^{p}V^{\ast}$ be the free
$R$-module of polynomial differential $p$-forms; $\Omega^{0}_{R}=R$ and
$\Omega^{l}_{R}=0$ for $l<0$. Then the radial vector field
$\operatorname{rad}$ defines a $R$-linear map
$\iota_{\operatorname{rad}}\colon\Omega^{p}_{R}\rightarrow\Omega^{p-1}_{R}$ so
that its kernel is $\bigoplus_{r}H^{0}(\Omega_{\mathbb{P}^{n}}^{p}(r))$.

###  2.1. Distributions

Given $1\leq k\leq n-1$, a codimension $p$ distribution $\mathscr{D}$ on
${\mathbb{P}^{n}}$ is defined by a short exact sequence

(4) $0\longrightarrow{\rm
T}_{\mathscr{D}}\stackrel{{\scriptstyle\phi}}{{\longrightarrow}}{\rm
T}{\mathbb{P}^{n}}\stackrel{{\scriptstyle\psi}}{{\longrightarrow}}{\rm
N}_{\mathscr{D}}\longrightarrow 0$

such that ${\rm N}_{\mathscr{D}}$ is a rank $p$ torsion-free sheaf; whence
${\rm T}_{\mathscr{D}}$ is a rank $n-p$ reflexive sheaf. The distribution
$\mathscr{D}$ is integrable, i.e. a foliation, if $\phi({\rm
T}_{\mathscr{D}})$ is closed under the Lie bracket of vector fields.

Taking exterior powers of (4) yields a differential $p$-form $\omega\in
H^{0}(\Omega_{\mathbb{P}^{n}}^{p}(d+p+1))$, where $d:=c_{1}({\rm
T}_{\mathscr{D}}(-1))$ is called the degree of $\mathscr{D}$. It follows that
the map $\psi$ is given by contraction: $\psi(v)=\iota_{v}\omega$. The
coefficients of $\omega$ generate the singular scheme
$\operatorname{Sing}(\mathscr{D})\subset{\mathbb{P}^{n}}$, supported on the
set of points where ${\rm N}_{\mathscr{D}}$ is not free. As ${\rm
N}_{\mathscr{D}}$ is torsion-free,
$\operatorname{{codim}}\operatorname{Sing}(\mathscr{D})\geq 2$.

Therefore to study degree $d$ codimension $k$ distributions on
${\mathbb{P}^{n}}$ we may focus on homogeneous $p$-forms representing sections
of $H^{0}(\Omega_{\mathbb{P}^{n}}^{p}(d+p+1))$. But first notice that not
every such $p$-form induces a distribution.

######  Example 2.1.

Consider the $2$-form $\omega\in H^{0}(\Omega_{\mathbb{P}^{n}}^{2}(3))$ given
by

$\omega=x_{0}(dx_{1}\wedge dx_{2}+dx_{3}\wedge
dx_{4})-dx_{0}\wedge(x_{1}dx_{2}-x_{2}dx_{1}+x_{3}dx_{4}-x_{4}dx_{3}).$

One can readily check that the contraction map
$\iota_{\bullet}\omega\colon{\rm
T}{\mathbb{P}^{n}}\rightarrow\Omega_{\mathbb{P}^{n}}^{p-1}(3)$ is injective;
hence $\omega$ does not define a distribution. For instance, on the affine
chart $U_{0}=\\{x_{0}=1\\}$ we have natural local coordinates
$(x_{1},x_{2},x_{3},x_{4})$ and, for any local vector field
$v=\sum_{j=1}^{4}a_{j}\frac{\partial}{\partial x_{j}}$, we have

$\iota_{v}\omega|_{U_{0}}=\iota_{v}(dx_{1}\wedge dx_{2}+dx_{3}\wedge
dx_{4})=a_{1}dx_{2}-a_{2}dx_{1}+a_{4}dx_{3}-a_{3}dx_{4}$

which vanishes if and only if so does $v$.

In order to ensure that a $p$-form $\omega\in
H^{0}(\Omega_{\mathbb{P}^{n}}^{p}(d+p+1))$ defines a distribution, we must
impose some conditions. First we require that the coefficients of $\omega$
have no common factor, so that the singular scheme has codimension at least
two. The other condition to impose is that $\omega$ be locally decomposable
off the singular set (LDS for short), see [dM00]. This means that each point
in ${\mathbb{P}^{n}}\setminus\operatorname{Sing}(\mathscr{D})$ has a Zariski
open neighborhood $U$ such that $\omega|_{U}$ decomposes as a product of
$1$-forms. It follows from the analysis of associated systems, see [God69,
Chapter I], that if the ${\mathcal{O}}_{\mathbb{P}^{n}}(U)$-module ${\rm
N}_{\mathscr{D}}(U)$ is free of rank $p$, then
$\omega|_{U}=\eta_{1}\wedge\dots\wedge\eta_{p}$, for some
$\eta_{j}\in\Omega_{\mathbb{P}^{n}}^{1}(U)$. Conversely, if $\omega|_{U}$ is
decomposable then, up to reducing $U$, ${\rm N}_{\mathscr{D}}(U)$ is free of
rank $p$. In order to have an integrable distribution, we must require further
that $d\eta_{i}=\sum_{j}\theta_{ij}\wedge\eta_{j}$ for some
$\theta_{ij}\in\Omega_{\mathbb{P}^{n}}^{1}(U)$. Fortunately, there is a more
hands on characterization for LDS and integrable forms.

######  Remark 2.2.

To simplify our notation we hereafter conventionalise
${\textstyle\bigwedge}^{0}V^{\ast}={\textstyle\bigwedge}^{0}V=k$ and
$\iota_{v}\omega:=\omega$ for $v\in{\textstyle\bigwedge}^{0}V$.

######  Lemma 2.3.

A homogeneous $p$-form $\omega$ on ${\mathbb{P}^{n}}$ is locally decomposable
off the singular set (LDS) if

$(\iota_{v}\omega)\wedge\omega=0,\quad\text{ for every
}v\in{\textstyle\bigwedge}^{p-1}V;$

here $\iota_{v_{1}\wedge\dots\wedge
v_{p-1}}\omega:=\iota_{v_{p-1}}\cdots\iota_{v_{1}}\omega$. Moreover, a LDS
form $\omega$ is integrable if

$(\iota_{v}\omega)\wedge d\omega=0,\quad\text{ for every
}v\in{\textstyle\bigwedge}^{p-1}V.$

The proof is an iterated application of de Rham-Saito Division Lemma [Sai76]
after localizing to the principal open subset
$D_{+}(f)\subset{\mathbb{P}^{n}}$, for $f$ a coefficient of $\omega$, and we
leave it to the reader.

Given an LDS $p$-form $\omega$ defining a distribution $\mathscr{D}$, we want
to compute the tangent and normal sheaves of $\mathscr{D}$; to do so we may
use a suitable complex. This was observed in [CACJ18, p.13] for codimension
$p=1$. The general case is similar. Taking exterior powers of the Euler
sequence (2) we get a natural inclusion
$\Omega_{\mathbb{P}^{n}}^{p-1}(d+p+1)\hookrightarrow{\mathcal{O}}_{\mathbb{P}^{n}}(d+2)\otimes{\textstyle\bigwedge}^{p-1}V^{\ast}$.
On the other hand, as in (1), ${\rm T}{\mathbb{P}^{n}}$ is the cokernel of
$\operatorname{rad}\colon{\mathcal{O}}_{\mathbb{P}^{n}}\rightarrow{\mathcal{O}}_{\mathbb{P}^{n}}(1)\otimes
V$, induced by the radial vector field. Hence we consider the composition

$C_{\omega}\colon{\mathcal{O}}_{\mathbb{P}^{n}}(1)\otimes
V\twoheadrightarrow{\rm
T}{\mathbb{P}^{n}}\stackrel{{\scriptstyle\omega}}{{\longrightarrow}}\Omega_{\mathbb{P}^{n}}^{p-1}(d+p+1)\hookrightarrow{\mathcal{O}}_{\mathbb{P}^{n}}(d+2)\otimes{\textstyle\bigwedge}^{p-1}V^{\ast}.$

Note that ${\rm N}_{\mathscr{D}}$ is isomorphic to the image of $C_{\omega}$
and we also get:

(5)
${\mathcal{O}}_{\mathbb{P}^{n}}\stackrel{{\scriptstyle\operatorname{rad}}}{{\longrightarrow}}{\mathcal{O}}_{\mathbb{P}^{n}}(1)\otimes
V\stackrel{{\scriptstyle
C_{\omega}}}{{\longrightarrow}}{\mathcal{O}}_{\mathbb{P}^{n}}(d+2)\otimes{\textstyle\bigwedge}^{p-1}V^{\ast}.$

This complex is interesting because the associated complex of free $R$-modules
is computationally convenient to describe ${\rm T}_{\mathscr{D}}$. The
following is straightforward.

######  Lemma 2.4.

Let $\mathscr{D}$ be a codimension $p$ distribution on ${\mathbb{P}^{n}}$ of
degree $d$ given by a homogeneous $p$-form $\omega$. Then ${\rm
T}_{\mathscr{D}}$ is the middle cohomology of the complex (5) and ${\rm
N}_{\mathscr{D}}$ is the image of $C_{\omega}$.

##  3\. Forms with prescribed vanishing locus

Now we turn to the main objective of this work that is to describe the module
of homogeneous $p$-forms, not necessarily LDS, vanishing along some given
subscheme. To provide a description of distributions one may further apply
Lemma 2.3.

Let $Z\subset{\mathbb{P}^{n}}$ be a closed subscheme with ideal sheaf
$\mathscr{I}_{Z}$ and consider

$\mathcal{A}^{p}(Z):=\bigoplus_{d\geq
0}H^{0}(\Omega_{\mathbb{P}^{n}}^{p}(d+p+1)\otimes\mathscr{I}_{Z}).$

the $R$-module of twisted differential $p$-forms that vanish on $Z$. Let
$I_{Z}$ denote the saturated homogeneous ideal of $Z$, i.e.,
$I_{Z}:=\bigoplus_{j}H^{0}(\mathscr{I}_{Z}(j))\subset R$.

######  Proposition 3.1.

Let $Z\subset\mathbb{P}^{n}$ be a closed subscheme then

(6) $0\longrightarrow
I_{Z}\otimes\iota_{\operatorname{rad}}\Omega^{p+1}_{R}\longrightarrow\mathcal{A}^{p}(Z)\longrightarrow\operatorname{Tor}_{p}^{S}(I_{Z},k)(p+1)\longrightarrow
0.$

###### Proof.

Consider the $(p+1)$-st exterior power of Euler sequence tensorized with the
sheaf $\mathscr{I}_{Z}(d+p+1)$:

$0\longrightarrow\Omega^{p+1}_{\mathbb{P}^{n}}(d+p+1)\otimes\mathscr{I}_{Z}\longrightarrow\mathscr{I}_{Z}(d)\otimes{\textstyle\bigwedge}^{p+1}V^{\ast}\stackrel{{\scriptstyle\iota_{\operatorname{rad}}}}{{\longrightarrow}}\Omega^{p}_{\mathbb{P}^{n}}(d+p+1)\otimes\mathscr{I}_{Z}\longrightarrow
0,$

which is exact since $\Omega^{p}_{\mathbb{P}^{n}}$ is locally free. From the
long sequence of cohomology we get

$H^{0}(\mathscr{I}_{Z}(d))\otimes{\textstyle\bigwedge}^{p+1}V^{\ast}\stackrel{{\scriptstyle\iota_{\operatorname{rad}}}}{{\longrightarrow}}H^{0}(\Omega^{p}_{\mathbb{P}^{n}}(d+p+1)\otimes\mathscr{I}_{Z})\stackrel{{\scriptstyle\phi}}{{\longrightarrow}}H^{1}(\Omega^{p+1}_{\mathbb{P}^{n}}(d+p+1)\otimes\mathscr{I}_{Z})$

is exact. From [Eis05, Theorem 5.8] the image of $\phi$ is precisely
$\operatorname{Tor}^{S}_{p}(I_{Z},k)_{d+p+1}$.

Taking the direct sum over $d\geq 0$ we get the desired sequence of
$R$-modules. ∎

In most cases of interest one wants to describe degree $d$ distributions
singular along a $Z$ such that $H^{0}(\mathscr{I}_{Z}(d))=0$, so that

$H^{0}(\Omega^{p}_{\mathbb{P}^{n}}(d+p+1)\otimes\mathscr{I}_{Z})\simeq\operatorname{Tor}^{S}_{p}(I_{Z},k)_{d+p+1}.$

It is an interesting open question to decide whether
$H^{0}(\mathscr{I}_{Z}(d))=0$ holds for $Z=\operatorname{Sing}(\mathscr{D})$.
This is true, for instance, if $k=1$ and $\dim Z=0$, see [CJM22, Lemma 4.2].

###  3.1. $p$-forms and syzygies

Note that the sequence of graded $k$-vector spaces underlying (6) must split
and we derive such a splitting from the syzygies of $I_{Z}$. Consider the
minimal graded free resolution

(7) $0\longrightarrow
F_{n}\stackrel{{\scriptstyle\phi_{n}}}{{\longrightarrow}}F_{n-1}\stackrel{{\scriptstyle\phi_{n-1}}}{{\longrightarrow}}\dots\longrightarrow
F_{0}\stackrel{{\scriptstyle\phi_{0}}}{{\longrightarrow}}I_{Z}\longrightarrow
0,$

where $F_{i}=\bigoplus_{j}S(-j)^{\beta_{i,j}}$. Recall that, since the
resolution is minimal,

$\operatorname{Tor}^{S}_{p}(I_{Z},k)\simeq F_{p}\otimes
k=\bigoplus_{j}k^{\beta_{p,j}},$

with $k^{\beta_{p,j}}$ in degree $j$. Also, fixed the minimal generators given
by $\phi_{0}$, the module of $p$-th syzygies of $I_{Z}$ is the image of
$\phi_{p}$. Note that if we tensorize (7) with the free module
$\Omega^{l}_{R}$ we get an exact sequence

(8) $0\longrightarrow
F_{n}\otimes\Omega^{l}_{R}\stackrel{{\scriptstyle\phi_{n}}}{{\longrightarrow}}F_{n-1}\otimes\Omega^{l}_{R}\stackrel{{\scriptstyle\phi_{n-1}}}{{\longrightarrow}}\dots\longrightarrow
F_{0}\otimes\Omega^{l}_{R}\stackrel{{\scriptstyle\phi_{0}}}{{\longrightarrow}}I_{Z}\otimes\Omega^{l}_{R}\longrightarrow
0$

where $\phi_{j}=\phi_{j}\otimes 1$ by abuse of notation.

We then define a $k$-linear map
$\delta\colon\Omega^{p}_{R}\rightarrow\Omega^{p+1}_{R}$ by setting

$\delta\omega=\frac{d\omega}{\deg\omega},$

on a homogeneous $\omega$. Here $\deg\omega$ is the total degree of $\omega$
considering $\deg dx_{i}=\deg x_{i}=1$. The important property $\delta$ has is
that

$\iota_{\operatorname{rad}}\delta\omega=\omega,\quad\text{for}\quad\omega\in\ker\iota_{\operatorname{rad}}.$

Below, for a matrix of $p$-forms $G=(g_{ij})$ we denote $\delta G=(\delta
g_{ij})$ and similarly for $\iota_{\operatorname{rad}}$.

######  Proposition 3.2.

The morphism
$\xi_{p}\colon\operatorname{Tor}_{p}^{S}(I_{Z},k)\rightarrow\mathcal{A}^{p}(Z)$
defined by

$\xi_{p}(t)=\phi_{0}\cdot\delta(\phi_{1}\cdot\delta(\cdots\delta\phi_{p}))\cdot
t.$

is injective and its image does not intersect the image of
$I_{Z}\otimes\iota_{\operatorname{rad}}\Omega^{p+1}_{R}$.

###### Proof.

First note that, due to $R$-linearity,
$\iota_{\operatorname{rad}}(\phi_{j}\cdot\delta\phi_{j+1})=\phi_{j}\cdot\phi_{j+1}=0$.
It is then straightforward to show that
$\iota_{\operatorname{rad}}\xi_{p}(t)=0$.

We will assume by contradiction that $\xi_{p}(t)\in
I_{Z}\otimes\iota_{\operatorname{rad}}\Omega^{p}_{R}$ and conclude that $t=0$,
hence proving both claims at once. From the assumption there exists $\eta_{0}$
a column vector of $(p+1)$-forms such that

$\xi_{p}(t)=\phi_{0}\cdot\iota_{\operatorname{rad}}\eta_{0}.$

Hence

$\phi_{0}\cdot\left(\delta(\phi_{1}\cdot\delta(\cdots\delta\phi_{p}))\cdot
t-\iota_{\operatorname{rad}}\eta_{0}\right)=0$

and, due to the exactness of (7) twisted by $\Omega^{p}_{R}$, there exists
$\eta_{1}$ a column vector of $p$-forms such that

$\delta(\phi_{1}\cdot\delta(\cdots\delta\phi_{p}))\cdot
t-\iota_{\operatorname{rad}}\eta_{0}=\phi_{1}\cdot\eta_{1}.$

Applying $\iota_{\operatorname{rad}}$ we get, due to $R$-linearity,

$\phi_{1}\cdot\delta(\cdots\delta\phi_{p}))\cdot
t=\phi_{1}\cdot\iota_{\operatorname{rad}}\eta_{1}\Longrightarrow\phi_{1}\cdot\left(\delta(\cdots\delta\phi_{p}))\cdot
t-\iota_{\operatorname{rad}}\eta_{1}\right)=0.$

Thus there exists $\eta_{2}$, a column vector of $(k-1)$-forms, such that
$\delta(\cdots\delta\phi_{p}))\cdot
t-\iota_{\operatorname{rad}}\eta_{1}=\phi_{2}\cdot\eta_{2}$. Iterating this
process we arrive at

$\delta\phi_{p}\cdot
t-\iota_{\operatorname{rad}}\eta_{p-1}=\phi_{p}\cdot\eta_{p}$

where $\eta_{p}$ is a column vector of $1$-forms. Hence there exists a matrix
of polynomials $A$ such that

$t-\iota_{\operatorname{rad}}\eta_{p}=\phi_{p+1}\cdot A.$

Since the resolution (7) is minimal, each entry of $\phi_{p+1}$ is a
homogeneous polynomial; the same is true for
$\iota_{\operatorname{rad}}\eta_{p}$. On the other hand, the entries of
$t\in\operatorname{Tor}_{p}^{S}(I_{Z},k)$ are constants. Thus, comparing
degrees we see that $t=0$.

∎

Note that from (6) we expect that $f\xi_{p}(t)\in
I_{Z}\otimes\iota_{\operatorname{rad}}\Omega^{p}_{R}$, for any homogeneous
polynomial $f\in S$. Indeed, we can write

$f\xi_{p}(t)=\phi_{0}\cdot\iota_{\operatorname{rad}}\left(\delta
f\wedge\delta(\phi_{1}\cdot\delta(\cdots\delta\phi_{p}))\cdot t\right).$

Also we have that $\xi_{p}(t)$ is not homogeneous unless
$t\in\operatorname{Tor}^{S}_{p}(I_{Z},k)_{m}$ for some $m$; in this case the
total degree of $\xi_{p}(t)$ is $m$. Moreover, if
$h^{0}(\mathscr{I}_{Z}(d))=0$ we can pass to the first linear strand of (7):

$0\longrightarrow
F_{n}^{0}\stackrel{{\scriptstyle\phi_{n}^{0}}}{{\longrightarrow}}F_{n-1}^{0}\stackrel{{\scriptstyle\phi_{n-1}^{0}}}{{\longrightarrow}}\dots\longrightarrow
F_{0}^{0}\stackrel{{\scriptstyle\phi_{0}^{0}}}{{\longrightarrow}}I_{Z}$

where $F_{j}^{0}=S(-j-d-1)^{\beta_{j,j+d+1}}$ and $\phi_{j}^{0}$ are the
corresponding linear blocks. Thus
$H^{0}(\Omega^{p}_{\mathbb{P}^{n}}(d+p+1)\otimes\mathscr{I}_{Z})\simeq\operatorname{Tor}^{S}_{p}(I_{Z},k)_{d+p+1}$
may be computed from

$\xi_{p}^{0}(t)=\phi_{0}^{0}\cdot d\phi_{1}^{0}\wedge\dots\wedge
d\phi_{p}^{0}\cdot t,$

which involves only the degree $d+1$ generators of $I_{Z}$.

##  4\. Examples

In this section we compute some examples. We will focus on degree $d$
distributions singular along $Z$ such that $(I_{Z})_{d}=0$. This is expected
to always hold for $Z$ the full singular scheme, see the introduction to
[GJM22].

######  Example 4.1 ($n=2,p=1,d=1$).

Let us start with a simple example. Let $Z\subset\mathbb{P}^{2}$ be a reduced
subscheme of length $3$ not contained in a line. Then we may suppose

$I_{Z}=(x_{0},x_{1})\cap(x_{0},x_{2})\cap(x_{1},x_{2})=(x_{0}x_{1},x_{0}x_{2},x_{1}x_{2}).$

The resolution is given by

$0\longrightarrow R^{2}\xrightarrow{\begin{pmatrix}x_{2}&0\\\ -x_{1}&x_{1}\\\
0&-x_{0}\end{pmatrix}}R^{3}\xrightarrow{\begin{pmatrix}x_{0}x_{1}&x_{0}x_{2}&x_{1}x_{2}\end{pmatrix}}I_{Z}\longrightarrow
0$

Thus $H^{0}(\Omega_{\mathbb{P}^{2}}^{1}(3)\otimes\mathscr{I}_{Z})\cong k^{2}$
spanned by

$\omega_{1}=x_{0}x_{1}dx_{2}-x_{0}x_{2}dx_{1}\quad\text{and}\quad\omega_{2}=x_{0}x_{2}dx_{1}-x_{1}x_{2}dx_{0}.$

Note that both $\omega_{1}$ and $\omega_{2}$ vanish along a line and a point,
but a general linear combination of them vanish precisely at $Z$.

######  Example 4.2 ($n=3,p=1,d=1$).

Also for $d=1$ consider $Z\subset\mathbb{P}^{3}$ a twisted cubic:

$0\longrightarrow R^{2}\xrightarrow{\begin{pmatrix}x_{0}&x_{1}\\\
x_{1}&x_{2}\\\
x_{2}&x_{3}\end{pmatrix}}R^{3}\xrightarrow{\begin{pmatrix}x_{1}x_{3}-x_{2}^{2}&x_{1}x_{2}-x_{0}x_{3}&x_{0}x_{2}-x_{1}^{2}\end{pmatrix}}I_{Z}\longrightarrow
0$

Then we also get
$H^{0}(\Omega_{\mathbb{P}^{2}}^{1}(3)\otimes\mathscr{I}_{Z})\cong k^{2}$. A
general element vanish only on $Z$.

######  Example 4.3 ($n=3,p=1,d=1$).

Next we describe a pathological example, still for codimension one and degree
one on $\mathbb{P}^{3}$. Consider $Z$ given by

$I_{Z}=(x_{0}^{2},x_{1}^{2},x_{0}x_{2},x_{1}x_{2},x_{2}^{2}-x_{0}x_{1}).$

It has length $5$ and is supported on a single point. It is a simple example
of a point that is not a local complete intersection. Thus there is no local
$1$-form vanishing only on $Z$. Due to Theorem A, any $1$-form vanishing on
$Z$ can be written as $\omega=Adx_{0}+Bdx_{1}+Cdx_{2}$ where

$\displaystyle A$
$\displaystyle=t_{0}x_{1}x_{2}+t_{2}x_{1}^{2}+t_{3}x_{0}x_{2}+t_{4}(x_{0}x_{1}-x_{2}^{2})$
$\displaystyle B$
$\displaystyle=-t_{0}x_{0}x_{2}+t_{1}x_{1}x_{2}-t_{2}(x_{0}x_{1}-x_{2}^{2})-t_{4}x_{0}^{2}$
$\displaystyle C$
$\displaystyle=-t_{1}x_{1}^{2}-t_{2}x_{1}x_{2}-t_{3}x_{0}^{2}+t_{4}x_{0}x_{2}$

and $t_{0},\dots,t_{4}\in\mathbb{C}$. Note that $\omega$ does not depend on
$x_{3}$ so it is a linear pullback of a $1$-form $\eta$ on $\mathbb{P}^{2}$
and the singular locus is thus a cone over the singular locus of $\eta$.
Therefore any $\omega$ vanishing on $Z$ must vanish along $3$ lines concurring
at $Z_{\rm red}$.

######  Example 4.4 ($n=3,p=2,d=2$).

Any $(n-1)$-form $\omega\in H^{0}(\Omega^{n-1}_{\mathbb{P}^{n}}(d+n)$ can be
written as
$\omega=\iota_{\operatorname{rad}}\iota_{v}\,dx_{0}\wedge\dots\wedge dx_{n}$
for some vector field $v\in H^{0}({\rm T}{\mathbb{P}^{n}}(d-1)$; in particular
it is LDS. The distribution, actually the foliation, described by $\omega$ is
often better described by $v$, and to get this vector field from $\omega$ we
just note that $\iota_{v}\,dx_{0}\wedge\dots\wedge
dx_{n}=\frac{1}{d+n}d\omega$. If
$v=\sum_{j=0}^{n}a_{j}\frac{\partial}{\partial x_{j}}$ then the singular
scheme is defined by the maximal minors of the matrix

$\begin{pmatrix}x_{0}&\cdots&x_{n}\\\ a_{0}&\cdots&a_{n}\end{pmatrix},$

which coincides with the ideal generated by the coefficients of $\omega$, up
to saturation.

Now we specialize to $\mathbb{P}^{3}$. In [CJM19] foliations by curves on
$\mathbb{P}^{3}$ are studied with a special focus on those having locally-free
conormal sheaf ${\rm N}_{\mathscr{D}}^{\vee}$. Nonetheless estimates on Chern
classes predict a foliation of degree $2$ with conormal sheaf satisfying
$c_{1}({\rm N}_{\mathscr{D}}^{\vee})=-5,c_{2}({\rm N}_{\mathscr{D}}^{\vee})=9$
and $c_{3}({\rm N}_{\mathscr{D}}^{\vee})=3$; here $c_{3}>0$ implies non-
locally-free. Then [CJM19, Theorem 4.1] translates it to predicting a
foliation singular along $Z=C\cup P$ where $C$ is a curve of degree $2$ and
genus $-2$ and $P$ is zero dimensional of length $3$. Then consider, for
instance,

$\displaystyle C$
$\displaystyle=V(x_{0}^{2},x_{0}x_{1},x_{1}^{2},x_{0}(x_{2}^{2}-x_{3}^{2})-x_{1}x_{3}x_{2}),$
$\displaystyle P$ $\displaystyle=V(x_{2}-x_{0},x_{2}+x_{1},x_{3})\cup
V(x_{1}-x_{0},x_{2},x_{3}+x_{0})$ $\displaystyle\qquad\cup
V(x_{1}-2x_{0},x_{2}+x_{0},x_{3}-x_{0}).$

Computing the syzygies we can construct, inside a space of dimension $4$, the
vector field $v=\sum_{j=0}^{3}a_{j}\frac{\partial}{\partial x_{j}}$ where

$\displaystyle a_{0}$
$\displaystyle=-4\,x_{0}^{2}-50\,x_{0}x_{1}+20\,x_{1}^{2},$ $\displaystyle
a_{1}$ $\displaystyle=-40\,x_{0}^{2}+16\,x_{0}x_{1}-10\,x_{1}^{2},$
$\displaystyle a_{2}$
$\displaystyle=40\,x_{0}^{2}-45\,x_{0}x_{1}+35\,x_{1}^{2}-4\,x_{0}x_{2}+50\,x_{1}x_{2}+30\,x_{0}x_{3},$
$\displaystyle a_{3}$
$\displaystyle=50\,x_{0}^{2}+40\,x_{0}x_{1}-40\,x_{1}^{2}+30\,x_{0}x_{2}-4\,x_{0}x_{3}+20\,x_{1}x_{3}.$

To check that $v$ is singular precisely along $Z$ one may follow the Macaulay2
routine below.

i1 : R = QQ[x_0..x_3]; i2 : C = ideal(x_0^2, x_0*x_1, x_1^2, x_0*(x_2^2-x_3^2)
- x_1*x_3*x_2); o2 : Ideal of R i3 : P = intersect(ideal(x_2-x_0,
x_2+x_1,x_3),ideal(x_1-x_0, x_2,x_3+x_0), ideal(x_1-2*x_0, x_2+x_0,x_3-x_0));
o3 : Ideal of R i4 : Z = intersect(C,P); o4 : Ideal of R i5 : a0 =
-4*x_0^2-50*x_0*x_1+20*x_1^2; i6 : a1 = -40*x_0^2+16*x_0*x_1-10*x_1^2; i7 : a2
= 40*x_0^2-45*x_0*x_1+35*x_1^2-4*x_0*x_2+50*x_1*x_2+30*x_0*x_3; i8 : a3 =
50*x_0^2+40*x_0*x_1-40*x_1^2+30*x_0*x_2-4*x_0*x_3+20*x_1*x_3; i9 : singD =
saturate minors(2, matrix{{x_0,x_1,x_2,x_3},{a0,a1,a2,a3}}); o9 : Ideal of R
i10 : Z == singD o10 = true

## References

  * [CACJ18] Omegar Calvo-Andrade, Maurício Corrêa, and Marcos Jardim, _Codimension One Holomorphic Distributions on the Projective Three-space_ , International Mathematics Research Notices (2018), 9011–9074, rny251.
  * [CJM19] Maurício Corrêa, Marcos Jardim, and Simone Marchesi, _Classification of foliations by curves of low degree on the three-dimensional projective space_ , 2019.
  * [CJM22] Maurício Corrêa, Marcos Jardim, and Alan Muniz, _Moduli of distributions via singular schemes_ , Math. Z. 301 (2022), no. 3, 2709–2731. MR 4437336
  * [CO01] Antonio Campillo and Jorge Olivares, _Polarity with respect to a foliation and Cayley-Bacharach theorems_ , J. Reine Angew. Math. 534 (2001), 95–118. MR 1831632
  * [CO21] by same author, _Cayley-Bacharach and singularities of foliations_ , Bull. Braz. Math. Soc. (N.S.) 52 (2021), no. 3, 477–498. MR 4290357
  * [dM00] Airton S. de Medeiros, _Singular foliations and differential $p$-forms_, Ann. Fac. Sci. Toulouse Math. (6) 9 (2000), no. 3, 451–466. MR 1842027
  * [Eis05] David Eisenbud, _The geometry of syzygies_ , Graduate Texts in Mathematics, vol. 229, Springer-Verlag, New York, 2005, A second course in commutative algebra and algebraic geometry. MR 2103875
  * [GJM22] Hugo Galeano, Marcos Jardim, and Alan Muniz, _Codimension one distributions of degree 2 on the three-dimensional projective space_ , J. Pure Appl. Algebra 226 (2022), no. 2, Paper No. 106840, 32. MR 4287792
  * [God69] Claude Godbillon, _Géométrie différentielle et mécanique analytique_ , Hermann, Paris, 1969. MR 0242081
  * [GS] Daniel R. Grayson and Michael E. Stillman, _Macaulay2, a software system for research in algebraic geometry_ , Available at http://www.math.uiuc.edu/Macaulay2/.
  * [Sai76] Kyoji Saito, _On a generalization of de-Rham lemma_ , Ann. Inst. Fourier (Grenoble) 26 (1976), no. 2, vii, 165–170. MR 413155