# Capturing Logical Structure of Visually Structured Documents with Multimodal
Transition Parser

Yuta Koreeda  Christopher D. Manning  Stanford University, Stanford, CA, USA

###### Abstract

While many NLP pipelines assume raw, clean texts, many texts we encounter in
the wild, including a vast majority of legal documents, are not so clean, with
many of them being visually structured documents (VSDs) such as PDFs.
Conventional preprocessing tools for VSDs mainly focused on word segmentation
and coarse layout analysis, whereas fine-grained logical structure analysis
(such as identifying paragraph boundaries and their hierarchies) of VSDs is
underexplored. To that end, we proposed to formulate the task as prediction of
_transition labels_ between text fragments that maps the fragments to a tree,
and developed a feature-based machine learning system that fuses visual,
textual and semantic cues. Our system is easily customizable to different
types of VSDs and it significantly outperformed baselines in identifying
different structures in VSDs. For example, our system obtained a paragraph
boundary detection F1 score of 0.953 which is significantly better than a
popular PDF-to-text tool with an F1 score of 0.739.

##  1 Introduction

Despite recent motivation to utilize NLP for wider range of real world
applications, most NLP papers, tasks and pipelines assume raw, clean texts.
However, many texts we encounter in the wild, including a vast majority of
legal documents (e.g., contracts and legal codes), are not so clean, with many
of them being visually structured documents (VSDs) such as PDFs. For example,
of 7.3 million text documents found in Panama Papers (which arguably
approximates the distribution of data one would see in the wild),
approximately 30% were PDFs111Calculated from Obermaier et al. (2016) by
regarding their emails, PDFs and text documents as the denominator.. Good
preprocessing of VSDs is crucial in order to apply recent advances in NLP to
real world applications.

Figure 1: Overview of the logical structure analysis for VSDs and its
formulation.

Thus far, the most micro and macro extremes of VSD preprocessing have been
extensively studied, such as word segmentation and layout analysis (detecting
figures, body texts, etc.; Soto and Yoo, 2019; Stahl et al., 2018),
respectively. While these two lines of studies allow extracting a sequence of
words in the body of a document, neither of them accounts for local, logical
structures such as paragraph boundaries and their hierarchies.

These structures convey important information in any domain, but they are
particulary important in the legal domain. For example, Figure 1(1) shows raw
text extracted from a non-disclosure agreement (NDA) in PDF format. An
information extraction (IE) system must be aware of the hierarchical structure
to successfully identify target information (e.g., extracting “definition of
confidential information” requires understanding of hierarchy as in Figure
1(2)). Furthermore, we must utilize the logical structures to remove debris
that has slipped through layout analysis (“Page 1 of 5” in this case) and
other structural artifacts (such as semicolons and section numbers) for a
generic NLP pipeline to work properly.

Yet, such logical structure analysis is difficult. Even the best PDF-to-text
tool with a word-related error rate as low as 1.0% suffers from 17.0% newline
detection error Bast and Korzen (2017) that is arguably the easiest form of
logical structure analysis.

The goal of this study is to develop a fine-grained logical structure analysis
system for VSDs. We propose a transition parser-like formulation of logical
structure analysis, where we predict a _transition label_ between each
consecutive pair of text fragments (e.g., two fragments are in a same
paragraph, or in different paragraphs of different hierarchies). Based on such
formulation, we developed a feature-based machine learning system that fuses
multimodal cues: visual (such as indentation and line spacing), textual (such
as section numbering and punctuation), and semantic (such as language model
coherence) cues. Finally, we show that our system is easily customizable to
different types of VSDs and that it significantly outperforms baselines in
identifying different structures in VSDs. For example, our system obtained a
paragraph boundary detection F1 score of 0.953 that is significantly better
than PDFMiner222https://euske.github.io/pdfminer/, a popular PDF-to-text tool,
with an F1 score of 0.739. We open-sourced our system and
dataset333https://github.com/stanfordnlp/pdf-struct.

##  2 Problem Setting and Our Formulation

In this study, we concentrate on logical structure analysis of VSDs. The input
is a sequence of text blocks (Figure 1(3)) that can be obtained by utilizing
existing coarse layout analysis and word-level preprocessing tools. We aim to
extract paragraphs and identify their relationships. This is equivalent to
creating a tree with each block as a node (Figure 1(4)).

We propose to formulate this tree generation problem as identification of a
_transition label_ between each consecutive pair of blocks (Figure 1(5)) that
defines their relationship in the tree. We define the transition
$\mathit{trans}_{i}$ between $i$-th block (hereafter $b_{i}$) and $b_{i+1}$ as
one of the following:

continuous

$b_{i}$ and $b_{i+1}$ are continuous in a single paragraph (Figure 1(6)).

consecutive

$b_{i+1}$ is the start of a new paragraph at the same level as $b_{i}$ (Figure
1(7)).

down

$b_{i+1}$ is the start of a new paragraph that is a child (a lower level) of
the paragraph that $b_{i}$ belongs to (Figure 1(6)).

up

$b_{i+1}$ is the start of a new paragraph that is in a higher level than the
paragraph that $b_{i}$ belongs to (Figure 1(8)).

omitted

$i$-th block is debris and omitted (Figure 1(9)). $\mathit{trans}_{i-1}$ is
carried over to the relationship betwen $b_{i-1}$ and $b_{i+1}$.

While down is well-defined (because we assume a tree), up can be ambiguous as
to how many levels we should raise. To that end, we also introduce a pointer
to each up block, which points at $b_{j}$ whose level $b_{i}$ belongs to
($\mathit{ptr}_{i}=b_{j}$, where $j<i$; Figure 1(8)).

##  3 Dataset

In this study, we target four types of VSDs in different file formats and
languages:

Contract pdfen

English NDAs in PDF format.

Law pdfen

English executive orders from local authorities.

Contract txten

English NDAs in visually structured plain text format.

Contract pdfja

Japanese NDAs in PDF format.

Examples of each type of VSDs are shown in footnote 5.

(a) Contract pdfen

(b) Law pdfen

(c) Contract txten

(d) Contract pdfja

Figure 2: Examples of VSDs in our dataset555(a)
http://www.astho.org/Programs/Infectious-Disease/Healthcare-Associated-
Infections/Electronic-Health-Records/Toolkit/Data-Use-Agreement-New-York-
City/, (b)
https://www2.illinois.gov/IISNews/21288-Gov._Pritzker_Stay_at_Home_Order.pdf,
(c)
https://www.sec.gov/Archives/edgar/data/86115/0000930661/0000930661-99-001321-index.htm,
and (d) http://www.septima.co.jp/contracts/27_himitsuhoji.pdf

For PDFs, we downloaded PDFs from Google.com search result. Since our focus is
not on coarse layout analysis or word-level preprocessing, we selected single
column documents and extracted blocks with an existing software. Specifically,
we utilized PDFMiner and extracted each LTTextLine, which roughly corresponds
to each line of text, as a block. We merged multiple LTTextLines where
LTTextLines are vertically overlapping.

For plain texts, we searched documents filed at
EDGAR666https://www.sec.gov/edgar.shtml. We simply used each non-blank line of
a plain text as a block.

We annotated all documents by hand. We describe more details of the data
collection and annotation in Section A.1.

The data statistics are given in Table 1. While the number of documents is
somewhat limited, we note that each document comes with many text blocks and
evaluations were stable. Furthermore, it was enough to reliably show the
difference between our system and baselines in our experiments.

##  4 Proposed System

###  4.1 Transition Parser

In this work, we propose to employ handcrafted features and a machine
learning-based classifier as the transition parser. This strategy is more
suited to our task than utilizing deep learning because  (1) we can
incorporate visual, textual and semantic cues, and  (2) it only requires a
small number of training data which is critical in the legal domain where most
data is proprietary  .

| Contract pdfen | Law pdfen | Contract txten | Contract pdfja
---|---|---|---|---
Format | PDF | PDF | Text | PDF
Language | English | English | English | Japanese
#Documents | 40 | 40 | 22 | 40
#Text blocks | 137.9 |  | 165.9 |  | 142.0 |  | 73.7 |
Max. depth | 3.4 |  | 3.9 |  | 3.1 |  | 3.0 |
#continuous | 95.4 | (68%) | 110.6 | (67%) | 109.9 | (77%) | 33.9 | (44%)
#consecutive | 20.3 | (17%) | 30.8 | (20%) | 15.3 | (12%) | 14.9 | (20%)
#up | 8.5 | ( 6%) | 7.1 | ( 4%) | 4.8 | ( 3%) | 11.0 | (15%)
#down | 9.4 | ( 6%) | 9.9 | ( 6%) | 4.6 | ( 3%) | 12.1 | (17%)
#omitted | 4.4 | ( 3%) | 7.6 | ( 3%) | 7.4 | ( 4%) | 1.8 | ( 2%)

  * A number in the second set of rows indicates an average count over documents. A percentage represents an average ratio of each label.

Table 1: Dataset information

|  |  |  | Document type
---|---|---|---|---
|  |  |  | Contract pdfen |  |
|  | Description |  Blocks | /Law pdfen | Contract txten | Contract pdfja
Visual features |
|  V1 | Indentation (up, down or same) |  1-2, 2-3 | ✓ | ✓ | ✓
|  V2 | Indentation after erasing numbering |  1-2, 2-3 |  | ✓ |
|  V3 | Centered |  2, 3 | ✓ | ✓ | ✓
|  V4 | Line break before right margin* |  1, 2 | ✓ | ✓ | ✓
|  V5 | Page change |  1-2, 2-3 | ✓ |  | ✓
|  V6 | Within top 15% of a page |  2 | ✓ |  | ✓
|  V7 | Within bottom 15% of a page |  2 | ✓ |  | ✓
|  V8 | Larger line spacing* |  1-2, 2-3 | ✓ | ✓ | ✓
|  V9 | Justified with spaces in middle |  2, 3 | ✓ | ✓ | ✓
|  V10 | Similar text in a similar position* |  2 | ✓ |  | ✓
|  V11 | Emphasis by spaces between characters |  1, 2 |  |  | ✓
|  V12 | Emphasis by parentheses |  1, 2 |  |  | ✓
Textual features |
|  T1 | Numbering transition* |  2 | ✓ | ✓ | ✓
|  T2 | Punctuated |  1, 2 | ✓ | ✓ | ✓
|  T3 | List start (/[-;:,]$/) |  1, 2 | ✓ | ✓ | ✓
|  T4 | List elements (/(;|,|and|or)$/) |  2 | ✓ | ✓ |
|  T5 | Page number (strict) |  1, 2, 3 | ✓ | ✓ | ✓
|  T6 | Page number (tolerant) |  1, 2, 3 | ✓ | ✓ | ✓
|  T7 | Starts with “whereas” |  3 | ✓ | ✓ |
|  T8 | Starts with “now, therefore” |  3 | ✓ | ✓ |
|  T9 | Dictionary-like (includes “:” & not V4) |  2, 3 | ✓ | ✓ |
|  T10 | All capital |  2, 3 | ✓ | ✓ |
|  T11 | Contiguous blank field (underbars) |  1-2, 2-3 | ✓ | ✓ | ✓
|  T12 | Horizontal line (“*-=#%_+” only) |  1, 2, 3 |  | ✓ |
Semantic features |
|  S1 | Language model coherence* |  1-2-3 | ✓ | ✓ | ✓

  * The “Blocks” columns list blocks used to extract features for $\mathit{trans}_{2}$ (e.g. “1-2, 2-3” means $\left[b_{i-1},b_{i}\right]$ and $\left[b_{i},b_{i+1}\right]$ are used to extract two sets of features). Features with a similar intended functionality are assigned the same feature name and implementations may vary for different document types. *: Explained in detail in Section 4.1.

Table 2: List of features for each feature extractor

For each block, our parser extracts features from a context of four blocks and
performs multi-class classification over the five transition labels. Since
omitted changes targets of transition, we also omit omitted blocks in feature
extraction. For $\mathit{trans}_{i}\neq\texttt{omitted}$, we extract features
from $\left[b_{i-1},b_{i},b_{j},b_{j+1}\right]$ where $b_{j}$ is the first
block after $b_{i}$ with $\mathit{trans}_{j}\neq\texttt{omitted}$. For
$\mathit{trans}_{i}=\texttt{omitted}$, we extract features from
$\left[b_{i-1},b_{i},b_{i+1},b_{i+2}\right]$. At test time, since we need to
know the presence of omitted before feature extraction, we run a first pass of
predictions to identify blocks with omitted, then use that information to
dynamically extract features to identify other labels.

Our system can be customized to different types of documents by modifying the
features. We have designed a feature set for each document type by visually
inspecting the training dataset (Table 2). For Contract txten, we regarded
space characters as horizontal spacing and blank lines as vertical spacing,
which allowed us to define features that are analogous to those for PDFs.

While readers can reference our open-sourced code for the concrete
implementation, we will discuss some of the features that have important
implementation details. For a target block $b_{i}$:

Numbering transition (T1)

A categorical feature that itself is a heuristic transition parser. It
identifies a numbering in each block and keeps a memory of the largest
numberings by their types (i.e., its alphanumeric type and styling, such as
IV. and (a)). It outputs  (1) continuousif no numbering is found,  (2)
consecutiveif the numbering in $b_{i+1}$ is contiguous to the numbering in
$b_{i}$,  (3) upif not consecutive and there is a corresponding number in the
memory, and  (4) downif it is none of above and it is the first number in its
numbering type  . For example, B0 in Figure 1 is down as 1. is the first
numbering type that it sees and “1” will be added to the memory. B1 and B2 are
continuous as no numbering is found and B3 is consecutive as a number “2” is
found in the same type as 1.. B4 is down as it contains a new numbering type.

Language model coherence (S1)

To determine if $b_{i}$ should be classified as omitted, it utilizes language
model to classify whether it is more natural to have $b_{i}$ or $b_{i+1}$
after $b_{i-1}$. Specifically, we use GPT-2 Radford et al. (2019) to calculate
language model loss $\ell(i,i-1)$ for $b_{i}$ given $b_{i-1}$ as a context
(i.e., fed into the model but not used in the loss calculation). We then
calculate $\ell(i,i-1)-\ell(i+1,i-1)$ as the feature. If it is more coherent
to have $b_{i}$ after $b_{i-1}$, $\ell(i,i-1)$ will be smaller than
$\ell(i+1,i-1)$ and the feature value will be negative. We also utilize
$\ell(i+1,i)-\ell(i+1,i-1)$.

Similar text in similar position (V10)

Headers and footers tend to appear in similar positions across different pages
with similar texts. For example, a contract may have the contract’s title on
every pages at the same position. This feature is $1$ if there exists a block
$b_{j}$ such that blocks’ overlapping area is larger than 50% of their
bounding box (treating as if they are on the same page), and their edit
distance is small777$d(b_{i},b_{j})/\max(len(b_{i}),len(b_{j}))<0.1$, where
$d$ gives the Levenshtein distance and $len$ gives the length of text..

Line break before right margin (V4)

A Boolean feature that is $0$ if the block spans to the right margin and $1$
otherwise (i.e., breaks before the right margin). To distinguish the body and
the margin of the document, we apply 1D clustering888We utilized a naïve 1D
clustering, where it greedily adds elements from a sorted list to a cluster
while the maximum difference of the elements is within a user-defined
threshold. on the right positions of the blocks and extract the rightmost
cluster with minimum members of six per page (to ignore headers/footers) as
the right margin (Figure 3). This margin information is used in other features
(V3, V6 and V7).

Larger line spacing (V8)

A Boolean feature that is $0$ if line spacing is _normal_ and $1$ otherwise.
To determine the _normal_ line spacing, we apply 1D clustering on line
spacings and pick a cluster with the largest number of members.

Figure 3: A sketch of how we determine right margin. We apply 1D clustering on
the right positions of the blocks and choose the rightmost cluster with at
least a user-defined number of members. If we choose to have a minimum of two
members, the right margin would be the cluster with three members.

###  4.2 Pointer Identification

We also implement the pointer identification with handcrafted features and a
machine learning-based classifier. Since a down transition creates a new level
that a block can point back to, we extract all pairs of
$\left[b_{j},b_{i}\right]$ ($b_{j}\in\mathcal{C}_{i}$) with
$\mathit{trans}_{i}=\texttt{up}$, $\mathit{trans}_{j}=\texttt{down}$ and
$j<i$. We then extract features from $\left[b_{j},b_{i}\right]$ and train a
binary classifier to predict
$p\left(\mathit{ptr}_{i}=b_{j}\middle|b_{j},b_{i}\right)$. In training, we use
ground truth down labels to extract candidates $\mathcal{C}_{i}$. At test
time, we aggregate $\mathcal{C}_{i}$ from predicted transition labels and
predict the pointer by
$\mathit{ptr}_{i}=\operatorname{argmax}_{b_{j}\in\mathcal{C}_{i}}p\left(\mathit{ptr}_{i}=b_{j}\middle|b_{j},b_{i}\right)$.

While our pointer points at a block with down ($b_{j}$), it is sometimes
important to extract features from the first block in the paragraph that
$b_{j}$ belongs to, which we will hereafter refer as $b_{head(j)}$. Using
$b_{head(j)}$, we extract the following features from
$\left[b_{j},b_{i}\right]$:

Consecutive numbering

Boolean features on whether a numbering in $b_{i}$ is contiguous to a
numbering in $b_{j}$ and $b_{head(j)}$, respectively.

Indentation

Categorical features on whether indentation gets larger, smaller or stays the
same from $b_{j}$ to $b_{i}$ and from $b_{head(j)}$ to $b_{i+1}$,
respectively.

Left aligned

Binary features on whether $b_{j}$, $b_{i+1}$ and $b_{head(j)}$ are left
aligned, respectively.

Transition counts

We count numbers of blocks $\\{b_{k}\\}_{j<k<i}$ with down and with up,
respectively. We use these two numbers along with their difference as
features. This is based on an intuition that a closer block with down tends to
be more important.

Pointer features are also customizable, but we used the same features999More
precisely, the pointer features are implemented slightly different for
different document types, such as numbering being modified to Japanese for
Contract pdfja, but they are intended to have similar functionalities. for all
the document types.

⬇

class PDFFeatureExtractor(BaseFeatureExtractor):

def __init__(self, text_boxes):

bboxes = np.array(

[tb.bbox for tb in text_boxes])

page_top = bboxes[:, 3].max()

page_bottom = bboxes[:, 1].min()

self.header_thresh = \ page_top - 0.15 * (page_top - page_bottom)

…

\par@single_input_feature([1])

def header_region(self, tb):

return tb.bbox[3] > self.header_thresh

\par@pairwise_feature([(0, 1), (1, 2)])

def page_change(self, tb1, tb2):

if tb1 is None or tb2 is None:

return True

return tb1.page != tb2.page

\par@pointer_feature()

def pointer_left_aligned(

self, head_tb, tb1, tb2, tb3):

return self.left_aligned(tb1)

\par…

Figure 4: The Python implmenetation of a feature extractor Figure 5:
Evaluation from IE perspective. For each of ground truth and predicted trees,
we extract a relationship matrix (right) that describes all the pairwise
relationships and calculate F1 scores/accuracy by comparing the matrices.

While we call our system a “transition parser”, we do not employ a stack and
instead employ the graph-based parser-like formulation for the pointer
identification. We selected this strategy because of the recent success of
graph-based parsers Dozat and Manning (2017); Zhang et al. (2019).

##  5 Implementation and Customization

In this section, we briefly describe the implementation of our system that
allows easy customization to different types of VSDs. Our system employs
modular and customizable design and is implemented in Python. A user may
implement a new feature extractor simply by writing a new feature extractor
class where each feature is implemented as its class function (Figure 4). For
example, @single_input_feature([1]) denotes that the subsequent function
should be applied to the second block of each context (thus corresponding to
feature V6). Likewise, the features for pointer identification can be
implemented by marking a function with @pointer_feature(), which takes a
candidate block $b_{j}$ (tb1), a target block $b_{i}$ (tb2), the block next to
the target block $b_{i+1}$ (tb3) and $b_{head(j)}$ (head_tb) as an input.

A feature extractor object is instantiated for each document where all feature
functions are automatically aggregated to produce the feature vector. A new
feature extractor can inherit from an existing feature extractor (e.g.,
feature extractors for Contract pdfen and Contract pdfja both inherit from a
base PDF feature extractor), which makes it easy to reuse implementations.

|  |  | Contract pdfen |  | Law pdfen |  | Contract txten |  | Contract pdfja
---|---|---|---|---|---|---|---|---|---
Relationship |  |  | Visual | Number | Ours |  | Visual | Number | Ours |  | Visual | Number | Ours |  | Visual | Number | Ours
Same paragraph |  Micro | P | 0.982 | 0.484 | 0.944 |  | 0.891 | 0.219 | 0.858 |  | 0.993 | 0.540 | 0.983 |  | 0.446 | 0.402 | 0.973
R | 0.683 | 0.947 | 0.951 |  | 0.681 | 0.969 | 0.957 |  | 0.708 | 0.917 | 0.978 |  | 0.552 | 0.985 | 0.966
F | 0.806 | 0.641 | 0.947 |  | 0.772 | 0.357 | 0.905 |  | 0.826 | 0.680 | 0.980 |  | 0.494 | 0.571 | 0.969
Macro | P | 0.980 | 0.644 | 0.955 |  | 0.906 | 0.328 | 0.936 |  | 0.990 | 0.595 | 0.969 |  | 0.481 | 0.478 | 0.971
R | 0.670 | 0.966 | 0.951 |  | 0.634 | 0.974 | 0.951 |  | 0.746 | 0.934 | 0.976 |  | 0.527 | 0.985 | 0.956
F | 0.782 | 0.736 | 0.948 |  | 0.731 | 0.452 | 0.933 |  | 0.847 | 0.687 | 0.971 |  | 0.450 | 0.617 | 0.955
Siblings |  Micro | P | 0.332 | 0.677 | 0.841 |  | 0.430 | 0.647 | 0.849 |  | 0.397 | 0.780 | 0.784 |  | 0.106 | 0.151 | 0.699
R | 0.323 | 0.765 | 0.736 |  | 0.283 | 0.504 | 0.712 |  | 0.481 | 0.763 | 0.723 |  | 0.506 | 0.571 | 0.691
F | 0.328 | 0.718 | 0.785 |  | 0.341 | 0.567 | 0.774 |  | 0.435 | 0.772 | 0.752 |  | 0.176 | 0.238 | 0.695
Macro | P | 0.443 | 0.678 | 0.791 |  | 0.598 | 0.493 | 0.793 |  | 0.482 | 0.677 | 0.814 |  | 0.347 | 0.237 | 0.719
R | 0.427 | 0.691 | 0.751 |  | 0.417 | 0.379 | 0.696 |  | 0.557 | 0.603 | 0.758 |  | 0.506 | 0.536 | 0.663
F | 0.337 | 0.650 | 0.748 |  | 0.410 | 0.385 | 0.724 |  | 0.435 | 0.605 | 0.754 |  | 0.292 | 0.283 | 0.671
Descendants |  Micro | P | 0.381 | 0.184 | 0.502 |  | 0.627 | 0.132 | 0.456 |  | 0.239 | 0.190 | 0.541 |  | 0.536 | 0.125 | 0.577
R | 0.123 | 0.879 | 0.807 |  | 0.303 | 0.881 | 0.858 |  | 0.048 | 0.888 | 0.771 |  | 0.340 | 0.580 | 0.826
F | 0.186 | 0.304 | 0.619 |  | 0.409 | 0.229 | 0.596 |  | 0.080 | 0.313 | 0.635 |  | 0.416 | 0.205 | 0.679
Macro | P | 0.295 | 0.242 | 0.655 |  | 0.438 | 0.173 | 0.581 |  | 0.193 | 0.269 | 0.639 |  | 0.462 | 0.122 | 0.737
R | 0.194 | 0.848 | 0.798 |  | 0.314 | 0.764 | 0.837 |  | 0.072 | 0.859 | 0.735 |  | 0.358 | 0.519 | 0.834
F | 0.203 | 0.340 | 0.641 |  | 0.327 | 0.230 | 0.617 |  | 0.096 | 0.367 | 0.625 |  | 0.372 | 0.195 | 0.739
Accuracy | Micro | 0.772 | 0.778 | 0.914 |  | 0.827 | 0.685 | 0.908 |  | 0.587 | 0.674 | 0.828 |  | 0.618 | 0.623 | 0.940
Macro | 0.686 | 0.679 | 0.889 |  | 0.732 | 0.427 | 0.840 |  | 0.571 | 0.580 | 0.841 |  | 0.623 | 0.492 | 0.899
Average F1 | Micro | 0.440 | 0.555 | 0.784 |  | 0.507 | 0.384 | 0.758 |  | 0.447 | 0.588 | 0.789 |  | 0.362 | 0.338 | 0.781
Macro | 0.441 | 0.576 | 0.779 |  | 0.489 | 0.356 | 0.758 |  | 0.459 | 0.553 | 0.783 |  | 0.372 | 0.365 | 0.788

  * “Micro”: Micro-average, “Macro”: Macro-average, “P”: Precision, “R”: Recall, “F”: F1 score

Table 3: Results for evaluation on IE perspective

##  6 Experiments

###  6.1 Evaluation Metrics

While we do report transition prediction accuracy, it is not a true task
metric since it is rooted on our formulation of the task. Looking back at our
initial motivation in Section 1, we introduce two sets of evaluation metrics.

The first set of metrics is rooted on IE perspective. For IE, it is important
to identify ancestor-descendant and sibling relationships because it allows,
for example, identifying a subject (in an ancestral block) and its objects (a
decendant block and its siblings). Thus, we evaluate F1 scores for identifying
pairs of blocks in  (1) same paragraph,  (2) sibling, and  (3) ancestor-
descendant  relationships, respectively (Figure 5). Note that we do not
include cousin blocks in the sibling relationship, because it is not clear
whether cousin blocks have any meaningful information in the context of IE.

We use the second set of metrics to evaluate a system’s efficacy as a
preprocessing tool for more general NLP pipelines. We evaluate paragraph
boundary identification metrics since paragraph boundaries can be used to
determine appropriate chunks of text to be fed into the NLP pipelines. We also
report accuracy for removing debris with omitted.

We used five-folds cross validation for the evaluation.

###  6.2 Baselines

We compared our system against the following baselines:

Numbering baseline

Hatsutori et al. (2017) This baseline detects numberings using a set of
regular expressions and identifies dropping in hierarchy when the type of
numberings has changed. Adopting Hatsutori et al. (2017) to our problem
formulation, our implementation is the same as the feature “numbering
transition (T1).”

Visual baseline

This baseline relies purely on visual cues; i.e., indentation and line
spacing. For each pair of consecutive blocks, this baseline outputs  (1)
continuouswhen indentation does not change and line spacing is normal (as in
feature V8),  (2) consecutivewhen indentation does not change and line spacing
is larger than normal,  (3) downwhen indentation gets larger, and  (4) upwhen
indentation gets smaller  . On up, it points back at the closest block with
the same indentation.

PDFMiner

We use this popular open-source project to detect paragraph boundaries as in
Bast and Korzen (2017). PDFMiner relies purely on geometric heuristics to
detect paragraph breaks.

|  |  | Contract pdfen |  | Law pdfen |  | Contract txten |  | Contract pdfja
---|---|---|---|---|---|---|---|---|---
Criteria |  |  | PDFMiner | Visual | Number | Ours |  | PDFMiner | Visual | Number | Ours |  | Visual | Number | Ours |  | PDFMiner | Visual | Number | Ours
Paragraph boundary |  Micro | P | 0.672 | 0.563 | 0.914 | 0.958 |  | 0.546 | 0.536 | 0.911 | 0.948 |  | 0.465 | 0.783 | 0.955 |  | 0.531 | 0.603 | 0.961 | 0.970
R | 0.822 | 0.968 | 0.700 | 0.948 |  | 0.858 | 0.916 | 0.637 | 0.948 |  | 0.989 | 0.637 | 0.945 |  | 0.850 | 0.663 | 0.627 | 0.991
F | 0.739 | 0.712 | 0.793 | 0.953 |  | 0.667 | 0.676 | 0.750 | 0.948 |  | 0.633 | 0.702 | 0.950 |  | 0.653 | 0.632 | 0.759 | 0.980
Macro | P | 0.698 | 0.598 | 0.921 | 0.958 |  | 0.632 | 0.565 | 0.866 | 0.946 |  | 0.527 | 0.840 | 0.953 |  | 0.585 | 0.645 | 0.964 | 0.970
R | 0.798 | 0.964 | 0.703 | 0.945 |  | 0.874 | 0.930 | 0.522 | 0.943 |  | 0.984 | 0.633 | 0.944 |  | 0.867 | 0.653 | 0.624 | 0.988
F | 0.722 | 0.729 | 0.772 | 0.947 |  | 0.703 | 0.692 | 0.620 | 0.940 |  | 0.673 | 0.693 | 0.947 |  | 0.661 | 0.627 | 0.745 | 0.976
Block elimination |  Micro | P | — | — | — | 0.969 |  | — | — | — | 0.979 |  | — | — | 0.865 |  | — | — | — | 1.000
R | — | — | — | 0.897 |  | — | — | — | 0.755 |  | — | — | 0.914 |  | — | — | — | 0.849
F | — | — | — | 0.932 |  | — | — | — | 0.852 |  | — | — | 0.889 |  | — | — | — | 0.919
Macro | P | — | — | — | 0.948 |  | — | — | — | 0.929 |  | — | — | 0.815 |  | — | — | — | 0.929
R | — | — | — | 0.906 |  | — | — | — | 0.858 |  | — | — | 0.816 |  | — | — | — | 0.866
F | — | — | — | 0.913 |  | — | — | — | 0.874 |  | — | — | 0.800 |  | — | — | — | 0.888

  * “Micro”: Micro-average, “Macro”: Macro-average, “P”: Precision, “R”: Recall, “F”: F1 score

Table 4: Results for evaluation on preprocessing perspective

| Contract pdfen | Law pdfen | Contract txten | Contract pdfja
---|---|---|---|---
# | Forward | Backward | Forward | Backward | Forward | Backward | Forward | Backward
| All | (0.914) | All | (0.914) | All | (0.908) | All | (0.908) | All | (0.828) | All | (0.828) | All | (0.940) | All | (0.940)
1 | T1, 2 | (0.763) | T1, 2 | (0.855) | T1, 2 | (0.685) | T1, 2 | (0.854) | V8, 2-3 | (0.333) | T2, 2 | (0.820) | T1, 2 | (0.596) | T1, 2 | (0.934)
2 | V10, 2 | (0.796) | T10, 3 | (0.794) | V8, 2-3 | (0.883) | V10, 2 | (0.859) | T10, 2 | (0.465) | V9, 2 | (0.811) | V12, 2 | (0.686) | V9, 2 | (0.909)
3 | T10, 3 | (0.818) | T10, 2 | (0.794) | V10, 2 | (0.893) | T2, 2 | (0.836) | T9, 2 | (0.716) | T3, 2 | (0.805) | V1, 2-3 | (0.821) | V8, 2-3 | (0.882)
4 | T7, 3 | (0.853) | V1, 2-3 | (0.796) | V8, 1-2 | (0.885) | V8, 2-3 | (0.800) | T3, 2 | (0.727) | V5 | (0.785) | T2, 2 | (0.813) | S1‡, 2 | (0.865)
5 | T10, 2 | (0.813) | S1‡, 2 | (0.808) | V5, 2-3 | (0.858) | V1, 2-3 | (0.747) | T6, 2 | (0.721) | T10, 2 | (0.781) | V8, 2-3 | (0.887) | T2, 2 | (0.856)
6 | V1, 2-3 | (0.844) | T8, 2-3 | (0.801) | T7, 3 | (0.881) | V4, 2-3 | (0.716) | T4, 2 | (0.723) | V8, 2-3 | (0.752) | V9, 3 | (0.906) | V4, 2 | (0.824)
7 | V4, 2 | (0.868) | T2, 2 | (0.774) | V1, 2-3 | (0.898) | V9, 2 | (0.676) | T2, 2 | (0.721) | S1†, 2 | (0.749) | T2, 1 | (0.913) | V1, 2-3 | (0.787)
8 | T2, 2 | (0.886) | V4, 2-3 | (0.692) | V3, 2 | (0.904) | S1†, 2 | (0.717) | V4, 2 | (0.722) | V9, 3 | (0.751) | V4, 2-3 | (0.926) | V9, 2 | (0.799)

  * Numbers in parentheses show micro-average transition label prediction accuracy. The first line shows the results with all features. †: $\ell(i,i-1)-\ell(i+1,i-1)$ variant. ‡: $\ell(i+1,i)-\ell(i+1,i-1)$ variant.

Table 5: Eight most important features chosen by greedy forward selection and
backward elimination.

###  6.3 Implementation Details

We used Random Forest Breiman (2001) as the transition and pointer
classifiers, which is suited for categorical features that occupy the majority
of our features. We did not tune hyperparameters of the Random Forest
classifier and used default values of scikit-learn Pedregosa et al. (2011).

For language model coherence feature S1, we used GPT-2
medium101010https://huggingface.co/gpt2 for English documents and japanese-
gpt2-medium111111https://huggingface.co/rinna/japanese-gpt2-medium) for
Japanese documents.

###  6.4 Results

Structure and preprocessing evaluations are shown on Table 3 and Table 4,
respectively. Our system obtained micro-average structure prediction accuracy
of 0.914 for Contract pdfen, 0.908 for Law pdfen, 0.828 for Contract txten and
0.940 for Contract pdfja, significantly outperforming the best baselines with
0.778, 0.827, 0.674 and 0.623, respectively. Our system performed the best
with respect to F1 scores for all but one structure relationships.

The difference was even more significant for paragraph boundary detection. For
Contract pdfen, our system obtained a micro-average paragraph boundary
detection F1 score of 0.953 that is significantly better than PDFMiner with an
F1 score of 0.739. PDFMiner performed on par with our visual baseline and
generally performed worse than our numbering baseline. This shows the
importance of incorporating textual information to preprocess VSDs.

Micro-average transition label prediction accuracies were 0.951 (Contract
pdfen), 0.938 (Law pdfen), 0.955 (Contract txten) and 0.923 (Contract pdfja).

We investigated the importance of each feature with greedy forward selection
and greedy backward elimination of the features (Table 5). We can observe that
our system makes a balanced use of the visual and textual cues. “Indentation
(V1)”, “Larger line spacing (V8)” and “numbering hierarchy (T1)”, which
partially represent the baselines, were ranked high in many cases. At the same
time, other features such as “all capital (T10)” and “punctuated (T2)” were
also contributing significantly to the accuracy, which made our system much
superior to the baselines.

The feature importance revealed that the semantic cue (S1) was no more
important than other cues. We suspect that the feature (which compares whether
adjacent or non-adjacent block is more likely given a context) had fallen back
to mere language model with the context being ignored in some cases, possibly
due to GPT-2 not being fine-tuned on the legal domain.

We also conducted a qualitative error analysis. For Contract pdfen, we found
that our system was performing poorly on documents where they had bold or
underlined section titles, followed by paragraphs without any indentation
(predicted continuous instead of down). We believe incorporating typographic
features would improve our system as implied by the success of the “all
capital (T10)” feature.

For Contract txten, we found that blocks that are all capitals or are all
underbars were misclassified as omitted. All capital words and underbars are
frequently used to denote headers and footers, but they were used as section
titles and input fields in these examples. Unlike for Contract pdfen, we
attribute this problem to lack of training data, as those should have been
classified correctly with other features (such as T4 and T8) if the system had
seen similar patterns in the training data.

Interestingly, we observed that the system tends to do better in documents
that are hierarchically more complex. This may be because hierarchically
complex documents tend to incorporate more cues to support humans comprehend
the documents.

##  7 Related Works

As discussed in Section 1, previous works mainly focused on word segmentation
and layout analysis, whereas fine-grained logical structure analysis of VSDs
is less addressed. Nevertheless, there exist some studies that focus on
similar goals.

Abreu et al. (2019) and Ferrés et al. (2018) have tried to deal with logical
structure analysis by identifying specific structures in VSDs such as
subheadings. However, these studies are too coarse-grained and cannot handle
paragraph-level logical structure, thus they are unable to satisfy the need we
have discussed in Section 1. FinSBD-3 shared task Au et al. (2021) is more
fine-grained than those works and incorporates extraction of list items.
However, its main focus is not on analysis of logical structures; it has only
four static levels for list hierarchies and does not consider hierarchies in
non-list paragraphs.

Hatsutori et al. (2017) proposed a rule-based system that purely relies on
numberings. We compared our system against it in Section 6 and showed that our
system, which also incorporates textual and semantic cues, is superior to
their method.

Sporleder and Lapata (2004) proposed a paragraph boundary detection method for
plain texts that purely relies on textual and semantic cues. While their
method is not intended for VSDs, some of their ideas could be incorporated to
our work as additional features. We leave use of more advanced semantic cues
for a future work.

While the goal is different, our textual features have some similarity to
those used in sentence boundary detection Gillick (2009). Since our goal is to
predict structures as well as boundaries, we employ richer textual and visual
features that they do not utilize.

LayoutLM Xu et al. (2020, 2021) incorporates multimodal self-supervised
learning to utilize deep learning for form understanding. While it may
alleviate the need for a large training dataset, it is not trivial to adopt
the same method for logical structure analysis as text blocks would not fit
onto the LayoutLM’s context. Furthermore, it is easier to diagnose and to
improve our system as it utilizes a combination of hand-crafted features,
while deep learning systems tend to be completely black box.

##  8 Conclusions

We proposed a transition parser-like formulation of the logical structure
analysis of VSDs and developed a feature-based machine learning system that
fuses visual, textual and semantic cues. Our system significantly outperformed
baselines and an existing open-source software on different types of VSDs. The
experiment revealed that incorporating both the visual and textual cues is
crucial in successfully conducting logical structure analysis of VSDs. As a
future work, we will incorporate typographic and more advanced semantic cues.

## Acknowledgements

We used computational resource of AI Bridging Cloud Infrastructure (ABCI)
provided by the National Institute of Advanced Industrial Science and
Technology (AIST) for the experiments.

## References

  * Abreu et al. (2019) Carla Abreu, Henrique Cardoso, and Eugénio Oliveira. 2019.  FinDSE@FinTOC-2019 Shared Task.  In _Proceedings of the Second Financial Narrative Processing Workshop_.
  * Au et al. (2021) Willy Au, Abderrahim Ait-Azzi, and Juyeon Kang. 2021.  FinSBD-2021: The 3rd Shared Task on Structure Boundary Detection in Unstructured Text in the Financial Domain.  In _Companion Proceedings of the Web Conference 2021_.
  * Bast and Korzen (2017) Hannah Bast and Claudius Korzen. 2017.  A Benchmark and Evaluation for Text Extraction from PDF.  In _2017 ACM/IEEE Joint Conference on Digital Libraries_.
  * Breiman (2001) Leo Breiman. 2001.  Random Forests.  _Machine Learning_ , 45(1):5–32.
  * Dozat and Manning (2017) Timothy Dozat and Christopher D. Manning. 2017.  Deep Biaffine Attention for Neural Dependency Parsing.  In _5th International Conference on Learning Representations_.
  * Ferrés et al. (2018) Daniel Ferrés, Horacio Saggion, Francesco Ronzano, and Àlex Bravo. 2018.  PDFdigest: an Adaptable Layout-Aware PDF-to-XML Textual Content Extractor for Scientific Articles.  In _Proceedings of the Eleventh International Conference on Language Resources and Evaluation_.
  * Gillick (2009) Dan Gillick. 2009.  Sentence Boundary Detection and the Problem with the U.S. In _Proceedings of Human Language Technologies: The 2009 Annual Conference of the North American Chapter of the Association for Computational Linguistics_.
  * Hatsutori et al. (2017) Yoichi Hatsutori, Katsumasa Yoshikawa, and Haruki Imai. 2017.  Estimating Legal Document Structure by Considering Style Information and Table of Contents.  In _New Frontiers in Artificial Intelligence_ , pages 270–283. Springer International Publishing.
  * Obermaier et al. (2016) Frederik Obermaier, Bastian Obermayer, Vanessa Wormer, and Wolfgang Jaschensky. 2016\.  About the Panama Papers.  _Süddeutsche Zeitung_.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. 2011.  Scikit-learn: Machine Learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830.
  * Radford et al. (2019) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019.  Language Models are Unsupervised Multitask Learners.  _OpenAI blog_ , 1(8):9.
  * Soto and Yoo (2019) Carlos Soto and Shinjae Yoo. 2019.  Visual Detection with Context for Document Layout Analysis.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing_.
  * Sporleder and Lapata (2004) Caroline Sporleder and Mirella Lapata. 2004.  Automatic Paragraph Identification: A Study across Languages and Domains.  In _Proceedings of the 2004 Conference on Empirical Methods in Natural Language Processing_.
  * Stahl et al. (2018) Christopher Stahl, Steven Young, Drahomira Herrmannova, Robert Patton, and Jack Wells. 2018.  DeepPDF: A Deep Learning Approach to Extracting Text from PDFs.  In _Proceedings of the Eleventh International Conference on Language Resources and Evaluation_.
  * The U.S. Securities and Exchange Commission (2018) The U.S. Securities and Exchange Commission. 2018.  _EDGAR® Public Dissemination Service Technical Specification_.
  * Xu et al. (2021) Yang Xu, Yiheng Xu, Tengchao Lv, Lei Cui, Furu Wei, Guoxin Wang, Yijuan Lu, Dinei Florencio, Cha Zhang, Wanxiang Che, Min Zhang, and Lidong Zhou. 2021.  LayoutLMv2: Multi-modal pre-training for visually-rich document understanding.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing_.
  * Xu et al. (2020) Yiheng Xu, Minghao Li, Lei Cui, Shaohan Huang, Furu Wei, and Ming Zhou. 2020.  Layoutlm: Pre-training of text and layout for document image understanding.  In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_.
  * Zhang et al. (2019) Sheng Zhang, Xutai Ma, Kevin Duh, and Benjamin Van Durme. 2019.  AMR Parsing as Sequence-to-Graph Transduction.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_.

##  Appendix A Appendix

###  A.1 Details of Data Collection and Annotation

In this section, we provide supplemental information regarding the data
collection and the annotation discussed in Section 3.

For PDFs, we queried Google search engines and downloaded the PDF files that
the search engines returned. We used the following queries and the domains:

Contract pdfen

“ “non-disclosure” agreement filetype:pdf” on seven domains from countries
where English is widely spoken (US “.com”, UK “.co.uk”, Australia “.com.au”,
New Zealand “.co.nz”, Singapore “.com.sg”, Canada “.ca”, South Africa
“.co.za”).

Law pdfen

“site:*.gov “order” filetype:pdf” on “google.com”.

Contract pdfja

““秘密保持契約書” filetype:pdf” on “google.co.jp”.

For the collection of Contract txten, we first download all the documents
filed at EDGAR from 1996 to 2020 in a form of daily
archives121212https://www.sec.gov/Archives/edgar/Oldloads/. We uncompressed
each archive and deserialized files using regular expressions by referencing
to the EDGAR specificationsThe U.S. Securities and Exchange Commission (2018),
which gave us 12,851,835 filings each of which contains multiple documents. We
then extracted NDA candidates from the documents by a rule-based filtering.
Using meta-data obtained during the deserialization, we extracted documents
whose file type starts with “EX” (denotes that it is an exhibit), its file
extension is one of “.pdf”, “.PDF”, “.txt”, “.TXT”, “.html”, “.HTML”, “.htm”
or “HTM”, and its content is matched by a regular expression
“(?<![a-zA-Z,̇"()]␣*)([Nn]on[-␣][Dd]isclosure)|(NON[-␣]DISCLOSURE)”.

We then randomly selected documents that fulfill following criteria:

  * •

it is an NDA or an executive order,

  * •

it has embedded texts (for PDFs),

  * •

it is a single column document, and

  * •

a similar document is not yet in the dataset.

The last criterion mainly targets contracts from same organizations and
executive orders from same authorities. It ensures that we get a wide variety
of documents in our dataset.

The datasets were annotated by one of the authors. We did not employ majority
vote to improve annotation consistency, because labels can be easily
determined by a brief inspection of the document.


[acronym]long-short

# Tensor Power Flow Formulations for Multidimensional Analyses in Distribution
Systems

Edgar Mauricio Salazar Duque1, Juan S. Giraldo2, Pedro P. Vergara3, Phuong H.
Nguyen1, Han (J.G.) Slootweg1  1 Eindhoven University of Technology,
Netherlands, 3 Delft University of Technology, Netherlands  2 Energy
Transition Studies Group, Netherlands Organisation for Applied Scientific
Research (TNO), Netherlands

###### Abstract

In this paper, we present two multidimensional power flow formulations based
on a fixed-point iteration (FPI) algorithm to efficiently solve hundreds of
thousands of power flows (PFs) in distribution systems. The presented
algorithms are the base for a new TensorPowerFlow (TPF) tool and shine for
their simplicity, benefiting from multicore central processing unit (CPU) and
graphics processing unit (GPU) parallelization. We also focus on the
mathematical convergence properties of the algorithm, showing that its unique
solution is at the practical operational point, which is the solution of high-
voltage and low-current. The proof is validated using numerical simulations
showing the robustness of the FPI algorithm compared to the classical Newton-
Raphson (NR) approach. In the case study, a benchmark with different PF
solution methods is performed, showing that for applications requiring a
yearly simulation at 1-minute resolution the computation time is decreased by
a factor of 164, compared to the NR in its sparse formulation.

###### Index Terms:

Power flow, fixed-point iteration, tensor power flow, mixed computer
resources.

##  I Introduction

The power flow study is crucial for different technical analyses of electrical
distribution systems. For example, power flows are widely used in time series
simulation (TSS), where long-term analysis depends on high granularity in
time, e.g., integration of distributed energy resources, Volt/Var control, and
hosting capacity [1]. Another widespread use of multiple power flows is the
probabilistic power flow (PPF), in which exogenous uncertainties are modeled
using scenarios, evaluated using power flows, and its impact evaluated using
stochastic analysis [2]. Many other applications require multiple power flow
executions, such as metaheuristic-based optimization, contingency analysis,
and machine learning in power systems [3].

Notably, the common factor between the aforementioned applications is the
execution of multiple (thousands or even millions) power flows to provide
insightful results. These applications are multidimensional regardless of the
technique used to reduce the number of timestamps, scenarios, or training size
[1, 2, 3]. The multidimensionality motivates efficient power flow formulations
and new techniques, providing fast and accurate results. Multidimensional
problems have been tackled in the past, such as in [2] for the PPF problem
using the embedded holomorphic power flow, or in [4] where tensors are used to
formulate the three-phase power flow. However, they are not scalable or even
computationally more exhaustive than traditional alternatives.

Advances in computer hardware, such as the increase in the number of cores in
CPUs and the evolution of GPU designs, transitioning from simple graphics
processors to highly parallel multiprocessors of many cores [5], opened a new
paradigm of programming and rethinking PF algorithms. New lines of research
look at reducing computational time by combining CPU and GPU resources to
improve the speed of convergence of Newton-Raphson algorithms [6, 7, 8, 9],
which is seen as a preferred method to improve due to its quadratic
convergence. However, the formulation of these approaches is designed
specifically to solve one power flow in the least amount of time.
Unfortunately, these algorithms are heavily focused on GPU-based
architectures, making their implementation a tedious process and requiring
specific GPU programming knowledge.

In this paper, we advocate for using a fixed point iteration algorithm to
solve multidimensional PF formulations. This fixed point algorithm has shown
robust performance [10]; it has a simple formulation for the case of
distribution system analyses, it is suitable and scalable for multidimensional
applications in the form of a tensor, and can also benefit from multicore CPU
and GPU parallelization. In this work, we also focus on the mathematical
convergence properties of the algorithm, showing its unique point of solution
at the practical operational point (if the solution exists), which is the
high-voltage, low-current solution. The contributions of the paper are as
follows:

  * •

Present a practical tool for multidimensional power flow analysis in
distribution systems, named TensorPowerFlow (TPF) 111The public link to the
repository will be provided after the acceptance of the paper in its dense and
sparse versions, which are based on a fixed-point iteration (FPI) algorithm.
The TPF opens new possibilities for using mixed computing resources (CPUs or
GPUs) to increase performance.

  * •

Give a geometric and physical interpretation of the existence of the power
flow solution using an FPI iterative algorithm, showing the robustness and
numerical stability of the FPI algorithm.

##  II Single Dimension Fixed Point Power Flow

We take a network with one substation, $b=|\Omega_{\mathrm{d}}|$ demand nodes,
and $\phi=|\Omega_{\phi}|$ phases. Nodal voltages and currents are related by
the admittance matrix
$\boldsymbol{Y}\in\mathbb{C}^{(b+1)\,\phi\times(b+1)\,\phi}$ as follows:

$\displaystyle\small\begin{bmatrix}\boldsymbol{i}_{s}\\\
-\boldsymbol{i}_{d}\end{bmatrix}=\begin{bmatrix}\boldsymbol{Y}_{ss}&\boldsymbol{Y}_{sd}\\\
\boldsymbol{Y}_{ds}&\boldsymbol{Y}_{dd}\end{bmatrix}\begin{bmatrix}\boldsymbol{v}_{s}\\\
\boldsymbol{v}_{d}\end{bmatrix}$ (1)

where vectors ${{\boldsymbol{i}}_{{s}}}\in\mathbb{C}^{\phi\times 1}$ and
${{\boldsymbol{i}}_{{d}}}\in\mathbb{C}^{b\phi\times 1}$ represent complex
injections of nodal current at the substation and at the demand nodes
$d\in\Omega_{\mathrm{d}}$, while
${{\boldsymbol{v}}_{{s}}}\in\mathbb{C}^{\phi\times 1}$ and
${\boldsymbol{v}}:={{\boldsymbol{v}}_{{d}}\in\mathbb{C}^{b\phi\times 1}}$ are
complex components of the respective nodal voltages. It must be noted that (1)
is general and can represent single-phase, polyphase, radial, meshed networks
and distributed generation of constant power.

In this paper, we call a single-dimensional power flow to the solution of (1)
for a single vector of nominal complex power
$\boldsymbol{s}\in\mathbb{C}^{b\phi\times 1}$ at the demand nodes, which is a
snapshot of consumption in the grid. The single-dimensional power flow
formulation, namely the FPI algorithm, is the basis of the tensorized version
and it is based on the notation of the equivalent version of the successive
approximation method (SAM) presented in [10] as:

$\displaystyle\boldsymbol{v}_{(n+1)}={\boldsymbol{F}\boldsymbol{v}_{(n)}^{*}}^{\circ(-1)}+\boldsymbol{w}$
(2)

with

$\displaystyle\boldsymbol{A}=\textbf{{diag}}\left(\alpha_{P}\\!\odot\\!\boldsymbol{s}^{*}\right);\hskip
15.0pt\boldsymbol{B}=\textbf{{diag}}\left({\alpha_{Z}\odot\boldsymbol{s}^{*}}\right)+\boldsymbol{Y}_{dd};$
$\displaystyle\boldsymbol{c}=\boldsymbol{Y}_{ds}\boldsymbol{v}_{s}+{\alpha_{I}\odot\boldsymbol{s}^{*}};$
$\displaystyle\boldsymbol{F}=-\boldsymbol{B}^{-1}\boldsymbol{A};\hskip
20.0pt\boldsymbol{w}=-{\boldsymbol{B}}^{-1}\boldsymbol{c}$ (3)

where $\odot$ is the Hadamard product,
${\boldsymbol{v}_{(n)}^{*}}^{\circ(-1)}$ the element-wise reciprocal of the
conjugated vector containing the nodal voltages at iteration $n$, and
$\alpha_{Z},\,\alpha_{I}$, and $\alpha_{P}$ represent the coefficients of the
ZIP load model per demand node and phase. It is worth mentioning that
$\boldsymbol{A}\in\mathbb{C}^{b\phi\times b\phi}$,
$\boldsymbol{B}\in\mathbb{C}^{b\phi\times b\phi}$, and
$\boldsymbol{c}\in\mathbb{C}^{b\phi\times 1}$ are constant matrices and vector
within the iterative process for a particular operating point. Thus,
$\boldsymbol{F}\in\mathbb{C}^{b\phi\times b\phi}$ and
$\boldsymbol{w}\in\mathbb{C}^{b\phi\times 1}$ can be calculated once used at
each iteration.

The algorithm in (2) can be rearranged as an iterative mapping using
$\boldsymbol{v}_{(n)}^{*\circ{-1}}=\boldsymbol{v}_{(n)}\oslash\|\boldsymbol{v}_{(n)}\|^{2}_{\circ}$,
where $\oslash$ is the Hadamard division and $\|\cdot\|_{\circ}$ is the
element-wise euclidean norm, as

$\displaystyle\boldsymbol{v}_{(n+1)}$ $\displaystyle=$
$\displaystyle-\boldsymbol{B}^{-1}\;\textbf{{diag}}(\alpha_{p}\odot\boldsymbol{s}^{*}\oslash\|\boldsymbol{v}_{(n)}\|^{2}_{\circ})\boldsymbol{v}_{(n)}+\boldsymbol{w}$
(4) $\displaystyle=$ $\displaystyle T(\boldsymbol{v}_{(n)})$

Here, the Banach fixed-point theorem is used to prove that (4) is a
contraction mapping in order to have a solution. This is the case where
$T(\cdot)$ satisfies $\|\boldsymbol{v}_{(n+1)}-\boldsymbol{v}_{(n)}\|_{1}\leq
k\|T(\boldsymbol{v}_{(n+1)})-T(\boldsymbol{v}_{(n)})\|_{1}$, where the norm-1
distance is used in our case ($\|\cdot\|_{1}$), that is, for vectors
$\|\boldsymbol{x}\|_{1}=\sum_{i}|x_{i}|$, where $|\cdot|$ is the absolute
value, and for matrices $\|A\|_{1}=\max_{j}(\sum_{i}|a_{ij}|)$, which is the
column norm. When $k<1$, the contraction mapping $T(\cdot)$ has a unique point
of convergence. This can be shown as

$\displaystyle\IEEEeqnarraymulticol{3}{l}{\|\boldsymbol{v}_{(n+1)}-\boldsymbol{v}_{(n)}\|_{1}}$
$\displaystyle=$
$\displaystyle\|\boldsymbol{B}^{-1}\;\textbf{{diag}}(\alpha_{p}\odot\boldsymbol{s}^{*}\oslash\|\boldsymbol{v}_{(n+1)}\|^{2}_{\circ})\boldsymbol{v}_{(n+1)}$
(5)
$\displaystyle-\>\boldsymbol{B}^{-1}\;\textbf{{diag}}(\alpha_{p}\odot\boldsymbol{s}^{*}\oslash\|\boldsymbol{v}_{(n)}\|^{2}_{\circ})\boldsymbol{v}_{(n)}\|_{1}$

where the conjugate load power in each bus-phases is expressed by its
equivalent load impedance solution $\boldsymbol{z}_{l}$. i.e.,
$\boldsymbol{s}^{*}\oslash\|\boldsymbol{v}_{(n)}\|^{2}_{\odot}=\mathbbm{1}\oslash\boldsymbol{z}_{l}$,
where $\mathbbm{1}$ is a vector of ones of dimension $b\phi$. Assuming that
$\boldsymbol{s}$ is a feasible load consumption in the grid (details of this
feasibility are given in Section III), the iterative algorithm reaches closer
voltage values for each iteration, i.e.,
$\boldsymbol{v}_{(n+1)}\mapsto\boldsymbol{v}_{(n)}$, then (5) can be reduced
using Hölder’s inequality as

$\displaystyle\|\boldsymbol{v}_{(n+1)}-\boldsymbol{v}_{(n)}\|_{1}$
$\displaystyle=$
$\displaystyle\|\boldsymbol{B}^{-1}\;\textbf{{diag}}(\alpha_{p}\oslash\boldsymbol{z}_{l})(\boldsymbol{v}_{(n+1)}-\boldsymbol{v}_{(n)})\|_{1}$
(6) $\displaystyle\leq$
$\displaystyle\|\boldsymbol{B}^{-1}\;\textbf{{diag}}(\alpha_{p}\oslash\boldsymbol{z}_{l})\|_{1}\|(\boldsymbol{v}_{(n+1)}-\boldsymbol{v}_{(n)})\|_{1}$
$\displaystyle\leq$ $\displaystyle
k\|(\boldsymbol{v}_{(n+1)}-\boldsymbol{v}_{(n)})\|_{1}$

From (6), for purely constant power $\alpha_{p}=1$, and noticing that
$\boldsymbol{B}^{-1}$ is the grid impedance matrix
$\boldsymbol{Z}_{\mathrm{B}}$. Then, the contraction scalar is

$\displaystyle k$ $\displaystyle=$
$\displaystyle\|\boldsymbol{Z}_{\mathrm{B}}\textbf{{diag}}(\mathbbm{1}\oslash\boldsymbol{z}_{l})\|_{1}=|\hat{z}_{jj}|/|\hat{z}_{l,j}|$
(7)

where the hat notation represents the solution values of the biggest ratio
between impedances. The only possible solution of (7) in order to be a
contraction mapping, i.e., $k<1$, is when $|z_{jj}|<|z_{l,j}|$. The diagonal
entries of $\boldsymbol{Z}_{\mathrm{B}}$, i.e., $|z_{jj}|$, are the Thevenin
impedance equivalent of bus $j$ [11], meaning that the only solution for (7)
is for the operational point of the network with high impedance (high voltage,
low current), which is the feasible operational state of the network. The
physical interpretation of (7) can be shown by analysis of an equivalent two-
bus system, which is discussed in detail in the next section.

Figure 1: Geometry of the solutions of the power flow problem for a two-bus
system with parameters $z_{s}=1.0+\boldsymbol{j}0.5$ and $\|v_{o}\|=1.0$. (a)
and (b) are the contour graphs for the circular formulation (9) and (10), for
$p_{l}$ and $q_{l}$, respectively; black circles highlight a feasible
combination of active and reactive power. (c) The red circles highlight an
infeasible value of the load power (no crossing between the circles), while
the feasible values have two points of solution that form a line that passes
through the origin. (d) All the solutions for the critical load power, which
has only one point of contact between the circles, form a circle with radius
$\|z_{s}\|$.

##  III Geometric interpretation of the existence of a power flow solution

Consider a grid formed by two nodes: node 0 as a source node, acting as the
slack bus with known voltage $v_{0}$ and reference angle $\theta_{0}=0$, and
node 1 with load $s_{l}$, operating as a load bus. Then, the power of the load
is given by $s_{l}=v_{l}i_{l}^{*}$, defining
$v_{l}=(z_{l}/(z_{s}+z_{l}))v_{0}$, $i_{l}=(v_{0}/(z_{s}+z_{l}))$, where the
source load is $z_{s}=r_{s}+\boldsymbol{j}x_{s}$, and the equivalent load
impedance is $z_{l}=r_{l}+\boldsymbol{j}x_{l}$. The power on the load defined
in terms of the source voltage and impedances is

$\displaystyle s_{l}$ $\displaystyle=$ $\displaystyle
p_{l}+\boldsymbol{j}q_{l}=\big{(}\frac{r_{l}}{a^{2}+b^{2}}+\boldsymbol{j}\frac{x_{l}}{a^{2}+b^{2}}\big{)}\|v_{0}\|^{2},$
(8)

where $a=(r_{s}+r_{l})$ and $b=(x_{s}+x_{l})$. Rearranging the expressions for
$p_{l}$ and $q_{l}$ into circle equations as (9) and (10) respectively, we
have that

$\displaystyle(r_{l}-c_{1,p})^{2}+(x_{l}-c_{2,p})^{2}$ $\displaystyle=$
$\displaystyle r_{p}^{2}$ (9)
$\displaystyle(r_{l}-c_{1,q})^{2}+(x_{l}-c_{2,q})^{2}$ $\displaystyle=$
$\displaystyle r_{q}^{2}$ (10) $\displaystyle c_{1,p}$
$\displaystyle=\frac{\|v_{0}\|^{2}}{2p_{l}}-r_{s}$ $\displaystyle c_{1,q}$
$\displaystyle=-r_{s}$ $\displaystyle c_{2,p}$ $\displaystyle=-x_{s}$
$\displaystyle c_{2,q}$ $\displaystyle=\frac{\|v_{0}\|^{2}}{2q_{l}}-x_{s}$
$\displaystyle r_{p}$
$\displaystyle=\frac{\|v_{0}\|}{2p_{l}}\sqrt{\|v_{o}\|^{2}-4r_{s}p_{l}}$
$\displaystyle r_{q}$
$\displaystyle=\frac{\|v_{0}\|}{2q_{l}}\sqrt{\|v_{o}\|^{2}-4x_{s}q_{l}}$

The contour plots of (9) and (10) are shown in Fig. 1(a,b), highlighting as an
example the circles for the power values $p_{l}=0.18$ [p.u] and $q_{l}=0.11$
[p.u]. The intersections of the circles are the two possible impedance
solutions for $z_{l}$, i.e., $z_{l,1}$ and $z_{l,2}$, where the values for
$p_{l}$ and $q_{l}$ can exist simultaneously. For example, Fig. 2(c) shows two
cases: the first is where the power flow has a solution and the circles
intersect (black line) at two points, and the second is where the red circle
does not have an intersection, meaning an infeasible combination of power
values. The intersection points of the circles form a line that passes through
the origin in all cases. This can be shown using the intercept in the line
equation $y=mr_{l}+\beta$, shown as a blue line in Fig. 1(c). The intercept is
defined in terms of the parameters of the two circles as

$\displaystyle\beta$ $\displaystyle=$
$\displaystyle\frac{\overbrace{(c_{1,p}^{2}+c_{2,p}^{2}+r_{q}^{2})}^{B_{0}}-\overbrace{(c_{1,q}^{2}+c_{2,q}^{2}+r_{p}^{2})}^{B_{1}}}{2(c_{2,p}-c_{2,q})}$
(11)

where the numerator of (11) is zero because $B_{0}=B_{1}$.

The worst-case loading scenario occurs when the two circles have only one
point of contact (Fig. 1(d)). The tangent line of the point of contact between
the circles is perpendicular to their radii and also passes through the origin
according to (11). Therefore, the altitude ($h$) of the scalene triangle
formed by the centers of the circles and the origin, that is, points
$\langle(c_{1,p},c_{2,p}),(c_{1,q},c_{2,q}),(0,0)\rangle$, can be calculated
using the triangle inequalities as

$\displaystyle h^{2}$ $\displaystyle=$
$\displaystyle(c_{1,p})^{2}+(c_{2,p})^{2}-r_{p}^{2}=r_{s}^{2}+x_{s}^{2}=\|z_{s}\|^{2}.$
(12)

This means that all unique points of contact between the two circles lie in
another circle defined by $\|z_{s}\|$, which is the magenta circle shown in
Fig. 1(c,d). This also means that the critical points where the power flow is
feasible and with only one impedance solution $z_{l}$, $z_{l,1}=z_{l,2}$, are
when $\|z_{s}\|=\|z_{l}\|$, which is the point of maximum power transfer. From
this, it is clear and important to notice that there are two possible
solutions for the two-bus system that are not critical; one lives inside and
the other outside the circle defined by $\|z_{s}\|$, meaning that
$\|z_{l,1}\|<\|z_{s}\|$ and $\|z_{l,2}\|>\|z_{s}\|$.

A region of convergence in the complex power plane $(p_{l},q_{l})$ can be
defined for the power in the load $s_{l}$ using the critical point of contact
of the two circles in the impedance plane $(r_{l},x_{l})$. When the circles
have one point of contact, then the distance between the centers of the
circles equals the sum of its radii, that is,
$r_{p}+r_{q}=\sqrt{(c_{1,p}-c_{1,q})^{2}+(c_{2,p}-c_{2,q})^{2}}$. Rearranging
this expression we have the following quadratic equation, named
$f_{\alpha}(p_{l},q_{l})|_{(r_{s},x_{s})}$, for $p_{l}$ and $q_{l}$

$\displaystyle Gq_{l}^{2}+Hp_{l}q_{l}+Ip_{l}^{2}+Jq_{l}+Kp_{l}+L$
$\displaystyle=$ $\displaystyle 0$ (13) $\displaystyle G$
$\displaystyle=r_{s}^{2}$ $\displaystyle H$ $\displaystyle=-2r_{s}x_{s}$
$\displaystyle I$ $\displaystyle=\|v_{0}\|^{2}x_{s}$ $\displaystyle J$
$\displaystyle=x_{s}^{2}$ $\displaystyle K$ $\displaystyle=\|v_{0}\|^{2}r_{s}$
$\displaystyle L$ $\displaystyle=-\|v_{0}\|^{2}/4$

Figure 2: Example of feasibility regions depicted by (13) and (14). (a) Each
point of the maximum power transfer circle $\|z_{s}\|$ in the impedance plane
$(r_{s},x_{s})$ parametrize the parabola in (13) that defines the feasible
load power values. (b) The vertexes of the parametrized parabolas are rotating
around the circle defined by (14) in the load power plane $(p_{l},q_{l})$.
Only points of the first quadrant are shown, highlighting three example points
($z_{s,1}$, $z_{s,2}$, $z_{s,3}$) with their respective parabolas. Feasibility
region for the case of $z_{s,2}$ is emphasised with red and yellow. The green
region is the union of all possible feasible regions for all parabolas.

The conic section defined by (13) is a parabola due to its discriminant
$\Delta=H^{2}-4GI=0$. Parabola is parametrized by the source impedance
$z_{s}$, i.e., $r_{s}$ and $x_{s}$. An example of this parabola where $z_{s}$
has only a resistive value is shown in blue in Fig. 2(a). The parabola curve
represents all $(p_{l},q_{l})$ values where the two-bus system has the maximum
power transfer, which means that it has only one equivalent load impedance
solution. The region defined on the left side of the parabola is all the power
values that have a valid solution, meaning that $z_{l,1}$ and $z_{l,2}$ exist,
and the right regions are values for which the system is infeasible. This
region changes as a result of the parameterization values of the parabola. The
blue parabola is when $z_{l}$ has only a reactance value. If we compute all
sets of parabolas
$\mathcal{S}=\\{f_{\alpha}(p_{l},q_{l})|_{(r_{s},x_{s})}\>|\>r_{s}^{2}+x_{s}^{2}=\|z_{s}\|^{2}\\}$
and have the union of all feasible power flow regions defined by $\mathcal{S}$
we have the circle delimiting such a region by
$\|s_{l}\|\leq\frac{\|v_{0}\|^{2}}{4\|z_{s}\|}$. This circle is the collection
of all the vertices of the rotating parabolas. This means that the condition
to have a feasible solution for the power flow is

$\displaystyle\|v_{0}\|^{2}$ $\displaystyle\geq$ $\displaystyle
4\|s_{l}\|\|z_{s}\|.$ (14)

This condition that is in the norm-form, is used to prove the existence of the
power flow solution in [12, 13, 14], and here it shows a geometrical
derivation and interpretation. It should be noted that this region (14) is a
conservative estimate, since the actual regions depend on specific
parameterization of $r_{s}$ and $x_{s}$ [14], as shown in Fig. 2(b) for the
$z_{s,2}$ example.

Therefore, if the power consumption in the grid follows at least the condition
(14)[13], and the power flow converges using the iterative fixed point
algorithm in (2), and following that the condition (7) must be less than one,
this implies that the unique point of convergence is the high-impedance
solution.

Figure 3: Convergence analysis of NR and FPI algorithms. (a) and (c)
corresponds to the regions of convergence for different initial values of
voltages $v_{0}$ for NR and FPI algorithms, respectively. The green region
corresponds to the high voltage (high impedance) and the red region to the low
voltage (low impedance) solutions. (b) and (d) correspond to the number of
iterations required for convergence, for the NR and FPI algorithms,
respectively.

The two-bus system is solved using the NR and FPI algorithms, with different
starting points $\boldsymbol{v}_{0}$ to numerically confirm the unique
convergence point of (2). The first row of subplots in Fig. 3 corresponds to
the NR and the second FPI solutions. The green and red colors in Fig. 3(a,c)
represent the attraction regions for high-voltage (high-impedance) and low-
voltage solutions, respectively. The number of iterations for each algorithm
is shown as colored contour plots in Fig. 3(b,d) for the NR and FPI
algorithms, respectively. The results confirm the robustness of the FPI
algorithm, which converges to the operational point irrelevant from the
initial voltage estimate, unlike NR, on which for a simple two-bus system, the
risk of falling into a nonoperational point (but mathematically valid) still
exists and is dependent on the $\boldsymbol{v}_{0}$ starting point.
Additionally, the FPI algorithm has a consistent number of iterations in the
entire complex domain for $v_{l}$. The interested reader is referred to [10]
for further numerical comparisons.

The NR algorithm in polar coordinates needs fewer iterations than the FPI
approach [10, 15]. Moreover, NR benefits from system sparsity, making it ideal
for sequential computers. Still, its per-iteration computation cost is
significantly higher due to the need for Jacobian inverse calculations. When
dealing with hundreds to millions of power flows, the cumulative cost of NR
iterations adds up, leading to longer processing times. Thanks to advancements
in computer hardware, especially faster matrix multiplications, the FPI
algorithm is an excellent choice for extensive power flow simulations on
modern computers for distribution networks. Its reliance on successive matrix
multiplications, along with its simple formulation, makes it easy to program.
Moreover, it can naturally extend to a tensor setting.

##  IV Multidimensional Fixed Point Power Flow

###  IV-A Tensor Power Flow - Dense Formulation

Consider a study that requires the analysis of an extensive number of cases of
load consumption. A tensor of the power can be built as
$\boldsymbol{S}\in\mathbb{C}^{\ldots\times p\times r\times t\times b\phi}$,
where $p$, $r$, and $t$ could mean a number of experiments, scenarios, and
time steps, with the possibility of extending the tensor to more dimensions.
The FPI algorithm in (2) in its tensor form is described as

$\displaystyle\boldsymbol{V}_{(n+1)}$ $\displaystyle=$
$\displaystyle\boldsymbol{\mathcal{F}}{\boldsymbol{V}^{*}_{(n)}}^{\circ(-1)}+\boldsymbol{\mathcal{W}}$
(15)

where dimensions of the matrices are
$\boldsymbol{\mathcal{F}}\in\mathbb{C}^{\ldots\times p\times r\times t\times
b\phi\times b\phi}$,
$\boldsymbol{V}_{(n+1)},\,{\boldsymbol{V}^{*}_{(n)}}^{\circ(-1)},\boldsymbol{\mathcal{W}}\in\mathbb{C}^{\ldots\times
p\times r\times b\phi\times t}$. An example of the structure of the tensors is
shown in Fig. 4(a). Recall that the submatrices $\boldsymbol{F}$ in (II),
which form the tensor $\boldsymbol{\mathcal{F}}$, are composed by a matrix
multiplication that involves $\boldsymbol{Z}_{B}$, meaning that
$\boldsymbol{\mathcal{F}}$ is dense. It should be noted that, in the specific
case of purely constant power and fixed topology, building a tensor
$\boldsymbol{\mathcal{F}}$ only requires the repetition of the constant matrix
$\boldsymbol{Z}_{B}$. Furthermore, with fixed voltage in the distribution
transformer, the tensor $\boldsymbol{\mathcal{W}}$ is also constant and is the
duplication of the vector $\boldsymbol{w}$ for all the cases under study. To
reduce the number of operations, the repetitions of the matrix/vector, and for
generalization for any number of dimensions, we define the dimensional tensor
elements, $\tau$, as ${\tau=\ldots\times p\times r\times t}$, in order to
reshape the tensors in (15) for efficient computation.

The reshaped form of the power and voltage tensors in their two-dimensional
matrix form is shown in Fig. 4(b). Their dimensions are
$\boldsymbol{\dot{S}^{*}},\,\boldsymbol{\dot{V}}_{(n+1)},\,{\boldsymbol{\dot{V}}^{*}_{(n)}}\in\mathbb{C}^{b\phi\times\tau}$,
where the dot notation stands for the reshaped version of the tensors. The
power matrix $\boldsymbol{\dot{S}^{*}}$ is the concatenation of the power
vectors $\boldsymbol{s}^{*}$ for all the cases under study along the secondary
axis. It should be mentioned that reshaping the tensor does not invalidate the
convergence of the FPI algorithm, as the update of the voltage values is the
same as (2). The advantage of the tensor form is that matrix operations can be
accelerated via parallelization on the CPU using standard low-level basic
linear algebra subprograms (BLAS), e.g. OpenBLAS, LAPACK, IntelMKL, or
exploiting the use of the multicore architectures from GPUs, which are
specifically optimized for the parallel processing matrix multiplications.
Algorithm 1 shows the implementation of the re-shaped version of (15).

Algorithm 1 : Tensor power flow - Dense

1:$b\phi$ = number of bus-phases, $\tau$ = number of power flows.

2:Input parameters: $\boldsymbol{\dot{S}}\in\mathbb{C}^{b\phi\times\tau}$,
$\boldsymbol{Z}_{B}\in\mathbb{C}^{b\phi\times b\phi}$,
$\boldsymbol{W}\in\mathbb{C}^{b\phi\times 1}$, tolerance, iterations.

3:Output parameters:
$\boldsymbol{\dot{V}}_{n}\in\mathbb{C}^{b\phi\times\tau}$, $n$.

4:$\boldsymbol{\dot{V}}_{0}=\mathbb{1}+\boldsymbol{j}\mathbb{0},\quad\boldsymbol{\dot{V}}_{0}\in\mathbb{C}^{b\phi\times\tau}$;
$n=0$; $\mathrm{tol}=\infty$

5:while $\mathrm{tol}\geq\mathrm{tolerance}$ $\boldsymbol{\mathrm{and}}$ $n<$
iterations do

6: for $i=1$ to ${\tau}$ do

7:
$\boldsymbol{\dot{V}}_{(n+1)}[i]=\boldsymbol{Z}_{B}(\boldsymbol{\dot{S}}[i]\odot\boldsymbol{\dot{V}}_{(n)}^{\circ(-1)}[i])^{*}+\boldsymbol{W}$

8: $\triangleright$ Iterate over the dimension $\tau$ of
$\boldsymbol{\dot{S}}$, and $\boldsymbol{\dot{V}}_{(n)}$.

9: end for

10:
$\mathrm{tol}=\mathrm{max}(\lVert\boldsymbol{s}_{(n+1)}-\boldsymbol{s}_{(n)}\rVert^{2})$;
$n=n+1$

11:end while

12:return ($\boldsymbol{\dot{V}}_{(n)}$, $n$)

Although the dense formulation is simple to implement, the tensor
$\boldsymbol{\mathcal{F}}$ could take a considerable amount of memory because
it is dense. Therefore, for cases of networks with large $b\phi$ the sparse
formulation is proposed.

$b\phi$
---
$r$
---
$r$
---
$r$
---
$t$
---
$p$
---
$b\phi$
---
$\boldsymbol{\mathcal{F}}$
---
$\boldsymbol{\mathcal{W}}$
---
${\boldsymbol{V}^{*}_{(n)}}^{\circ(-1)}$
---
$\boldsymbol{V}_{(n+1)}$
---
$r$
---
$\boldsymbol{Z}_{B}$
---
$\boldsymbol{\dot{S}}^{*}$
---
${\boldsymbol{\dot{V}}^{*}_{(n)}}^{\circ(-1)}$
---
$\boldsymbol{W}$
---
$\boldsymbol{\dot{V}}_{(n+1)}$
---
$\tau$
---
$\tau$
---
$\tau$
---
$b\phi$
---
$b\phi$
---
$r$
---
$r$
---
$\boldsymbol{\mathcal{M}}$
---
$\boldsymbol{\dot{\mathcal{M}}}$
---
$\boldsymbol{\mathcal{H}}$
---
$\boldsymbol{\mathcal{H}}$
---
$b\phi$
---
$r$
---
$t$
---
$p$
---
$r$
---
${\boldsymbol{V}^{*}_{(n)}}^{\circ(-1)}$
---
${\boldsymbol{\dot{V}}^{*}_{(n)}}^{\circ(-1)}$
---
$\boldsymbol{V}_{(n+1)}$
---
$\boldsymbol{\dot{V}}_{(n+1)}$
---
$\tau\cdot b\phi$
---
$\tau\cdot b\phi$
---
$b\phi$
---
Figure 4: Example of the tensor power flow formulations for a simulation with
a four-dimensional power tensor $\boldsymbol{S}\in\mathbb{C}^{p\times r\times
t\times b\phi}$, for $p=r=2$ and $t=b\phi=3$. (a) Visualization of the tensor
dense formulation of (15). (b) In the case of constant power, the reshaped
tensors use resources efficiently as the operations highlighted in red are
performed concurrently. (c) Visualization of the tensor sparse formulation in
(16), where tensors $\boldsymbol{\mathcal{M}}$ and $\boldsymbol{\mathcal{H}}$
are sparse. (d) Reshaped sparse formulation (17), which is solved iteratively
using a direct sparse solver.

###  IV-B Tensor Power Flow - Sparse Formulation

The tensor algorithm in (15) can be reformulated to exploit the sparsity of
$\boldsymbol{Y}_{dd}$ and $\boldsymbol{Y}_{ds}$. The sparse formulation is
defined as

$\displaystyle\boldsymbol{\mathcal{M}}\boldsymbol{V}_{(n+1)}$ $\displaystyle=$
$\displaystyle{\boldsymbol{V}^{*}_{(n)}}^{\circ(-1)}+\boldsymbol{\mathcal{H}}.$
(16)

where ${\boldsymbol{\mathcal{M}}\in\mathbb{C}^{\ldots\times p\times r\times
t\times b\phi\times b\phi}}$ is the tensor containing the resulting tensor
operation of
$\boldsymbol{\mathcal{M}}=-\boldsymbol{\mathcal{A}}^{\circ(-1)}\boldsymbol{\mathcal{B}}$,
where
$\boldsymbol{\mathcal{A}},\,\boldsymbol{\mathcal{B}}\in\mathbb{C}^{\ldots\times
p\times r\times t\times b\phi\times b\phi}$ and
${\boldsymbol{\mathcal{H}}\in\mathbb{C}^{\ldots\times p\times r\times
b\phi\times t}}$ a tensor containing the resulting vectors
$\boldsymbol{\mathcal{H}}=\boldsymbol{\mathcal{A}}^{\circ(-1)}\boldsymbol{\mathcal{C}}$,
where $\boldsymbol{\mathcal{C}}\in\mathbb{C}^{\ldots\times p\times r\times
b\phi\times t}$. An example of (16) is shown in Fig. 4(c). The sparse
formulation does not include the inverse of $\boldsymbol{Y}_{dd}$ in any form,
and similarly to the dense formulation, the expression in (16) can be re-
shaped in a two-dimensional matrix form to construct a sparse linear system of
the type $Ax=b$ as

$\displaystyle\underbrace{\vphantom{\boldsymbol{\dot{V}}_{(n+1)}}\boldsymbol{\dot{\mathcal{M}}}}_{A}\underbrace{\boldsymbol{\dot{V}}_{(n+1)}}_{x}$
$\displaystyle=$
$\displaystyle\underbrace{{\boldsymbol{\dot{V}}^{*}_{(n)}}^{\circ(-1)}+\boldsymbol{\dot{\mathcal{H}}}}_{b},$
(17)

where the reshaped sparse matrix $\boldsymbol{\dot{\mathcal{M}}}$ is the
diagonal concatenation of the submatrices of $\boldsymbol{\mathcal{M}}$.
Vectors
$\boldsymbol{\dot{\mathcal{V}}}_{(n+1)},\,\boldsymbol{\dot{\mathcal{V}}}_{(n)}$,
and $\boldsymbol{\dot{\mathcal{H}}}$ are vertically arranged. A visual example
of the form of this reshape is shown in Fig. 4(d). The system in (17) can be
solved iteratively by a sparse direct solver.

The sparse direct solvers have basically three steps: (i) analysis of the
matrix $\boldsymbol{\dot{\mathcal{M}}}$ to reduce fill-in (via Cholesky, LU,
or QR decomposition) and symbolic factorization, (ii) numerical factorization,
and (iii) solving the system [16]. It is critical to note that the main
advantage of (17) is that the sparse matrix $\boldsymbol{\dot{\mathcal{M}}}$
is constant and does not change during iterations. Therefore, steps (i) and
(ii) are performed only once, in the first iteration, which significantly
reduces the computing time for subsequent steps. The implementation of (17) is
shown in Algorithm 2. It is worth recalling that in the NR algorithm, the
sparse Jacobian matrix needs to be updated, requiring the first two steps in
every iteration to calculate its inverse, increasing the computational time.
Notice that rearranging the formulation from (15) to (16) does not invalidate
the convergence of the algorithm, as the sequence of the voltage values, that
is, $\\{\boldsymbol{V}_{(n)}\\}_{n=0}^{\infty}$, is still the same.

Algorithm 2 : Tensor power flow - Sparse

1:Input parameters:

2:$\boldsymbol{\dot{S}}\in\mathbb{C}^{b\phi\times\tau}$,
$\boldsymbol{Y_{dd}}\in\mathbb{C}^{b\phi\times
b\phi},\boldsymbol{Y_{ds}}\in\mathbb{C}^{b\phi\times 1}$, tolerance,
iterations. $\triangleright$ $\boldsymbol{Y_{dd}}$ and $\boldsymbol{Y_{ds}}$
are sparse.

3:Output parameters:
$\boldsymbol{\dot{V}}_{(n)}\in\mathbb{C}^{b\phi\times\tau}$, $n$.

4:$\boldsymbol{\dot{\mathcal{M}}}=-\mathrm{diag}(\boldsymbol{\dot{S}}_{[0]}^{*(-1)})\boldsymbol{Y_{dd}}$
$\triangleright$ [0] first entry over dimension $\tau$.

5:$\boldsymbol{\dot{\mathcal{H}}}=-\mathrm{diag}(\boldsymbol{\dot{S}}_{[0]}^{*(-1)})\boldsymbol{Y_{ds}}$

6:for $i=1$ to ${\tau}$ do

7:
$\boldsymbol{m}=-\mathrm{diag}(\boldsymbol{\dot{S}}_{[i]}^{*(-1)})\boldsymbol{Y_{dd}}$

8:
$\boldsymbol{h}=-\mathrm{diag}(\boldsymbol{\dot{S}}_{[i]}^{*(-1)})\boldsymbol{Y_{ds}}$

9:
$\boldsymbol{\dot{\mathcal{M}}}=\left(\begin{array}[]{c|c}\boldsymbol{\dot{\mathcal{M}}}&\mathbb{0}\\\
\hline\cr\mathbb{0}&\boldsymbol{m}\end{array}\right)$ $\triangleright$ Concat.
sparse matrix diagonally

10:
$\boldsymbol{\dot{\mathcal{H}}}=\left(\begin{array}[]{c}\boldsymbol{\dot{\mathcal{H}}}\\\
\hline\cr\boldsymbol{h}\end{array}\right)$ $\triangleright$ Concatenate sparse
vectors

11:end for

12:$\boldsymbol{\dot{V}}_{0}=\mathbb{1}+\boldsymbol{j}\mathbb{0},\quad\boldsymbol{\dot{V}}_{0}\in\mathbb{C}^{(b\phi\cdot\tau)\times
1}$; $n=0$; $\mathrm{tol}=\infty$

13:while $\mathrm{tol}\geq\mathrm{tolerance}$ $\boldsymbol{\mathrm{and}}$ $n<$
iterations do

14:
$\boldsymbol{\dot{V}}_{(n+1)}=\mathrm{solve}(\boldsymbol{\dot{\mathcal{M}}},{\boldsymbol{\dot{V}}^{*}}^{\circ(-1)}_{(n)}+\boldsymbol{\dot{\mathcal{N}}})$
$\triangleright$ Sparse system (17)

15:
$\mathrm{tol}=\mathrm{max}(\lVert\boldsymbol{s}_{(n+1)}-\boldsymbol{s}_{(n)}\rVert^{2})$;
$n=n+1$

16:end while

17:$\boldsymbol{V}_{n}=\mathrm{reshape}(\boldsymbol{\dot{V}}_{n},(b\phi\times\tau))$

18:return ($\boldsymbol{V}_{(n)}$, $n$)

##  V Simulation results

In this section, we discuss the comparison of Algorithms 1 and 2, labeled
Tensor (Dense) and Tensor (Sparse), respectively, against current methods,
such as SAM [10], NR with its sparse formulation in polar coordinates (NR
(Sparse)) [17], and the backward and forward sweep method (BFS) [18].
Additionally, the tensor-dense formulation is programmed to use a GPU to
quantify computational speed improvements; this implementation is named Tensor
(GPU). The implementation of Algorithms 1 and 2 is publicly available to the
community as a Python package named TensorPowerFlow222The public link to the
repository will be provided after the acceptance of the paper. The experiments
were conducted on a conventional laptop with an Intel(R) Core(TM) i7-7700HQ
CPU @ 2.80GHz, with 8 Logical Processor(s), and a NVIDIA Quadro M1200 with 4GB
(GDDR5) of memory.

###  V-A Computational Performance

The goal of the experiments is to test the computation performance of the
different methods for two aspects: (i) grid size and (ii) dimensionality. For
the first aspect, different networks with bus-phases $b\phi$, which range from
9 to 5k, are generated using a $k$-ary tree generative model, with random
$k$-child between 1 and 5, in order to simulate the radial structure of the
distribution system [19]. For the second aspect, the dimensional tensor
elements $\tau$, which are the number of PFs to compute, range from 10 to
525k, and the power values are generated using a multivariate elliptical
copula [20] to simulate realistic scenarios of consumption profiles. The same
power scenarios and networks are used for all methods. Additionally, for
methods that use sparse solvers, the Intel oneMKL - PARDISO is used in all of
them for a fair comparison.

Figure 5: Comparison of the performance of the algorithms. Computational wall-
times ($t_{c}$) increase the bus-phases size, $b\phi$, for (a) 1 power flow,
and (b) 500 power flows. (c) Test for increased dimensional tensor elements
$\tau$, for a grid size of 500 bus-phases ($b\phi$). (d) Computational
complexity fit with asymptotic complexity model $t_{c}=c\cdot n^{k}$ (solid
lines are $c$, dotted $k$). Tensor dense, sparse, and GPU are the proposed
algorithms in this paper.

The performance of the methods for the increase of $b\phi$ is shown for one PF
in Fig. 5(a), and for 500 PF Fig. 5(b). The best methods are SAM and Tensor
(Dense) for one PF with an execution time of 0.1 [ms] for the smallest
network. However, both methods reduce its performance for an increased grid
size, where around 900 bus-phases, the preference is for the Tensor (Sparse)
formulation. In the case of 500 PF (Fig.5(b)), the SAM lags behind because it
requires more matrix multiplications in its original formulation; the Tensor
(Dense) continues as the top performer for smaller grids, seconded by the
Tensor (Sparse), where after a grid size of 2k is better than the Tensor
(Dense). An interesting behavior is seen for the Tensor (GPU) implementation,
which is significantly lower for a small amount of power flows (Fig. 5(a)),
which is explained due to the overhead time required to transfer the data from
the host computer to the GPU. However, when the number of matrix operations is
large, either for increase $\tau$ or $b\phi$, the Tensor (GPU) is the best
performer, as shown in Fig. 5(b)(c). The computational complexity of the
algorithms with a model of the form $t_{c}=c\,n^{k}$, is shown Fig. 5(d),
confirming that SAM and BFS scale poorly for larger grids.

TABLE I: Computing Wall-times for solving one year of power flows at different
time resolutions (table values in minutes).

Algorithm |  | Grid size
---
($b\phi$)
| @1 hr
---
(8760 pf)
| @30 min
---
(17520 pf)
| @15 min
---
(35040 pf)
| @1 min
---
(525600 pf)
BFS | 100 | 1.98 | 3.96 | 7.91 | 118.71
NR (Sparse) | 100 | 1.03 | 2.06 | 4.12 | 61.85
Tensor (Sparse) | 100 | 0.17 | 0.34 | 0.68 | 10.14
SAM | 100 | 0.02 | 0.04 | 0.09 | 1.29
Tensor (GPU) | 100 | 0.02 | 0.04 | 0.09 | 1.30
Tensor (Dense) | 100 | 0.01 | 0.01 | 0.02 | 0.37
BFS | 5000 | 186.49 | 372.99 | 745.97 | -
NR (Sparse) | 5000 | 8.77 | 17.54 | 35.08 | 526.23
Tensor (Sparse) | 5000 | 2.43 | 4.86 | 9.71 | 145.71
SAM | 5000 | 108.09 | 216.18 | 432.36 | -
Tensor (GPU) | 5000 | 0.92 | 1.84 | 3.69 | 55.34
Tensor (Dense) | 5000 | 2.60 | 5.19 | 10.38 | 155.71

###  V-B Yearly Time Series Simulation

A comparative test is performed for a yearly TSS study for grid sizes ranging
from 100 to 5000 bus-phases ($b\phi$), for time resolution of 1, 15, 30, and
60 minutes. The results are shown in Table I, where the maximum waiting time
for the results is set to 9 hours. In general, all the FPI methods are
preferred and best suited for large amounts of PFs instead of conventional NR
and BFS techniques. Among the iterative methods, the proposed Tensor (Dense)
algorithm stood out as the fastest for the smaller grid size case. It
completed the 1-minute resolution case (approximately 525k PFs calculations)
in only 22 seconds. To put this in perspective, when compared with the NR
(Sparse) method –which took more than one hour to accomplish the same task–
the Tensor (Dense) revealed a remarkable speedup factor of 164. In the case of
larger grids (5k bus-phases), the proposed methods are still more efficient
than NR (Sparse), and the Tensor (Sparse) is the second best for all time
resolutions. However, the GPU implementation of Tensor (Dense) is the fastest,
showing the speed up from GPU, which is capable of massively parallelizing the
matrix multiplications automatically in a multi-thread setting. It is worth
noting that a GPU-based sparse solver could be utilized to enhance the Tensor
(Sparse), though this area requires further exploration.

##  VI Conclusions

This paper presented practical PF formulations for analyses of distribution
systems in a multidimensional scope. The formulations are based on a FPI
algorithm, and two algorithms are presented in their tensor forms for dense
and sparse versions. The convergence proof of the algorithms and the existence
of a solution are presented with its geometrical and physical interpretation.
The mathematical analysis shows that the FPI algorithm converges to the high
impedance value in the case of the existence of a solution, and it can be
extended to the multidimensional formulation. The performance of the proposed
algorithms compared to conventional methods such as NR and BFS was evaluated
in two aspects: the size of the grid (number of bus-phases, $b\phi$) and
dimensionality (number of PFs, $\tau$). The algorithms were tested in a
practical TSS application for different grid sizes and time resolutions.
Tensor (Dense) emerged as the fastest algorithm, providing a significant
speedup over traditional methods like NR (Sparse), especially for smaller grid
sizes. The study also introduced a GPU implementation, Tensor (GPU), which
showed reduced performance for a small number of power flows due to data
transfer overhead but excelled when the number of matrix operations was
substantial, as seen when the grid size and dimensionality increased.

## References

  * [1] M. U. Qureshi, S. Grijalva, M. J. Reno, J. Deboever, X. Zhang, and R. J. Broderick, “A fast scalable quasi-static time series analysis method for PV impact studies using linear sensitivity model,” _IEEE Trans. Sust. Energy_ , vol. 10, no. 1, pp. 301–310, 2019.
  * [2] C. Liu, K. Sun, B. Wang, and W. Ju, “Probabilistic power flow analysis using multidimensional holomorphic embedding and generalized cumulants,” _IEEE Trans. Power Syst._ , vol. 33, no. 6, pp. 7132–7142, 2018.
  * [3] E. M. S. Duque, J. S. Giraldo, P. P. Vergara, P. Nguyen, A. van der Molen, and H. Slootweg, “Community energy storage operation via reinforcement learning with eligibility traces,” _Electr. Power Syst. Res._ , vol. 212, p. 108515, 2022.
  * [4] A. Garcés, J. J. Mora, and M.-A. Useche, “Putting tensors back in power systems analysis,” in _2019 Int. Conf. on Smart Energy Systems and Technologies (SEST)_ , 2019, pp. 1–5.
  * [5] J. Nickolls and W. J. Dally, “The GPU Computing Era,” _IEEE Micro_ , vol. 30, no. 2, pp. 56–69, Mar. 2010.
  * [6] G. Zhou, R. Bo, L. Chien, X. Zhang, F. Shi, C. Xu, and Y. Feng, “GPU-Based Batch LU-Factorization Solver for Concurrent Analysis of Massive Power Flows,” _IEEE Trans. Power Syst._ , vol. 32, no. 6, pp. 4975–4977, Nov. 2017.
  * [7] X. Li, F. Li, H. Yuan, H. Cui, and Q. Hu, “GPU-Based Fast Decoupled Power Flow With Preconditioned Iterative Solver and Inexact Newton Method,” _IEEE Trans. Power Syst._ , vol. 32, no. 4, pp. 2695–2703, Jul. 2017.
  * [8] G. Zhou, R. Bo, L. Chien, X. Zhang, S. Yang, and D. Su, “GPU-Accelerated Algorithm for Online Probabilistic Power Flow,” _IEEE Trans. Power Syst._ , vol. 33, no. 1, pp. 1132–1135, Jan. 2018.
  * [9] G. Zhou, Y. Feng, R. Bo, L. Chien, X. Zhang, Y. Lang, Y. Jia, and Z. Chen, “GPU-Accelerated Batch-ACPF Solution for N-1 Static Security Analysis,” _IEEE Trans. Smart Grid_ , vol. 8, no. 3, pp. 1406–1416, May 2017.
  * [10] J. S. Giraldo, O. D. Montoya, P. P. Vergara, and F. Milano, “A fixed-point current injection power flow for electric distribution systems using Laurent series,” _Electr. Power Syst. Res._ , vol. 211, p. 108326, 2022\.
  * [11] J. J. Grainger and W. D. Stevenson, _Power system analysis_ , ser. McGraw-Hill series in electrical and computer engineering Power and energy. New York, NY St. Louis San Francisco: McGraw-Hill, Inc, 1994.
  * [12] U. Sur and G. Sarkar, “Existence of Explicit and Unique Necessary Conditions for Power Flow Insolvability in Power Distribution Systems,” _IEEE Syst. Jour._ , vol. 13, no. 1, pp. 702–709, Mar. 2019.
  * [13] S. Bolognani and S. Zampieri, “On the Existence and Linear Approximation of the Power Flow Solution in Power Distribution Networks,” _IEEE Trans. Power Syst._ , vol. 31, no. 1, pp. 163–172, Jan. 2016.
  * [14] S. Yu, H. D. Nguyen, and K. S. Turitsyn, “Simple certificate of solvability of power flow equations for distribution systems,” in _2015 IEEE Power & Energy Society General Meeting_. Denver, CO, USA: IEEE, Jul. 2015, pp. 1–5.
  * [15] A. Ahmadi, M. C. Smith, E. R. Collins, V. Dargahi, and S. Jin, “Fast Newton-Raphson Power Flow Analysis Based on Sparse Techniques and Parallel Processing,” _IEEE Trans. Power Syst._ , vol. 37, no. 3, pp. 1695–1705, May 2022.
  * [16] T. A. Davis, _Direct methods for sparse linear systems_ , ser. Fundamentals of algorithms. Philadelphia: SIAM, 2006\.
  * [17] B. Sereeter, C. Vuik, and C. Witteveen, “On a comparison of Newton–Raphson solvers for power flow problems,” _J. Comput. Appl. Math._ , vol. 360, pp. 157–169, Nov. 2019.
  * [18] E. Bompard, E. Carpaneto, G. Chicco, and R. Napoli, “Convergence of the backward/forward sweep method for the load-flow analysis of radial distribution systems,” _Int. J. Electr. Power Energy Syst._ , vol. 22, no. 7, pp. 521–530, Oct. 2000.
  * [19] J. A. Storer, _An Introduction to Data Structures and Algorithms_. Boston, MA: Birkhäuser Boston, 2002.
  * [20] E. M. S. Duque, P. P. Vergara, P. H. Nguyen, A. Van Der Molen, and J. G. Slootweg, “Conditional Multivariate Elliptical Copulas to Model Residential Load Profiles From Smart Meter Data,” _IEEE Trans. Smart Grid_ , vol. 12, no. 5, pp. 4280–4294, Sep. 2021.