DER

    distributed energy resource
DG

    distributed generator
ADN

    active distribution network
FRT

    fault ride through

# Probabilistic Stability Assessment for Active Distribution Grids

Sebastian Liemannゼ*, Lia Strengeリ, Paul Schultzパ, Holm Hinnersホ,
Johannis Porstヨ, Marcel Sarstedtマ, Frank Hellmannフ ゼInstitute of Energy
Systems, Energy Efficiency and Energy Economics (ie3),
Technical University of Dortmund, Dortmund, Germany *sebastian.liemann@tu-
dortmund.de
リControl Systems Group, Technische Universität Berlin, Berlin, Germany
フパRD4, Potsdam Institute for Climate Impact Research (PIK), Potsdam, Germany
ホInstitute of Automation, University of Bremen, Bremen, Germany
ヨInstitute of Electrical Energy Systems, Friedrich-Alexander University
Erlangen-Nürnberg, Germany
マInstitute of Electric Power Systems, Leibniz University Hannover, Hannover,
Germany

###### Abstract

This paper demonstrates the concept of probabilistic stability assessment on
large-signal stability in the use case of short circuits in an active
distribution grid. Here, the concept of survivability is applied, which
extends classical stability assessments by evaluating the stability and
operational limits during transients for a wide range of operating points and
failures. For this purpose, a free, open-source, and computationally efficient
environment (Julia) for dynamic simulation of power grids is used to
demonstrate its capabilities. The model implementation is validated against
established commercial software and deviations are minimal with respect to
power flow and dynamic simulations. The results of a large-scale survivability
analysis reveal i) a broad field of application for probabilistic stability
analysis and ii) that new non-intuitive stability correlations can be
obtained. Hence, the proposed method shows strong potential to efficiently
conduct power system stability analysis in active distribution grids.

###### Index Terms:

active distribution grid, Julia, probabilistic stability, short-circuits,
survivability

††publicationid: pubid: 978-1-6654-3597-0/21/$31.00 ©2021 IEEE

##  I Introduction

The operation of electrical grids is getting much more complex, due to the
increasing coordination effort between distributed energy resources and sector
coupling interfaces to the gas, heat or public communication system. In
addition, grid operators of different voltage levels have to cooperate more
closely in order to keep the electrical system stable. In this context, active
distribution networks will play a crucial part in the coordination of future
grids, to increase the number of DERs that can be integrated into the system.
For example, ADNs can be used to provide ancillary services for higher voltage
levels [1]. Since the control system in an ADN often allows a high degree of
operational freedom, many different grid operating points are reachable by the
operator [1]. For a stable operation of the grid, the stability of these
operating points have to be assessed in advance by dynamic simulations in
order to initiate necessary countermeasures. In case of ADNs, large signal
stability becomes especially relevant, as their controls and nonlinearites
have to be considered for a correct stability assessment. However, due to the
large operational envelope, classical stability assessments, e.g. short-
circuit calculations or tripping of large generators, only cover a small part
of these operating points and their corresponding stability regions. To
counteract this shortcoming, probabilistic stability assessment methods in
combination with a high-performance simulation tool can help to evaluate the
stability for a large set of operating points and parameters [2].

In particular, stability assessments based on a linearisation of the dynamics
around a stable operating point are usually insufficient in the face of large
(i.e. finite) perturbations due to the common nonlinearities and
multistability of power systems. Hence, the extent of the corresponding basin
of attraction becomes the determining quantity to assess the effect of large
state deviations. Direct methods based on energy functions [3, 4] have a long
tradition in transient stability analysis and are used to approximate the
basin geometry. Geometric approaches, however, are subject to the curse of
dimensionality, requiring high computational effort and hence are limited to
small systems [5, 6]. Alternatively, the volume of the basin [7, 8] can be
considered to quantify the stability of an operating point with respect to
finite perturbations. Probabilistic stability measures such as basin stability
[9] identify the volume measure of the basin with the probability that a power
system remains synchronised, given a specific distribution of finite
perturbations. This extends the classic small and large signal stability
assessment by not only evaluating single scenarios but rather a wide range of
the state space around operating points. The related concept of survivability
[10] extends this notion to include operational limits on transient
deviations. The advantage of these approaches is that the stability measures
can be efficiently evaluated for large systems using Monte Carlo sampling [11,
12].

Probabilistic stability measures have already been successfully applied to
power systems like swing equation networks [13, 14, 15, 16, 17] or DC
microgrids [18, 19]. However, an application to a more realistic grid model of
e.g. ADNs is still missing. Therefore, this paper demonstrates the practical
applicability of survivability by means of transient simulations of an ADN. As
dynamic simulations with many variations of control parameters are very time-
consuming, a simulation tool with high performance is needed. Here, the
programming language Julia is explored with its package PowerDynamics.jl [20]
to perform these simulations. However, the focus in this paper is on
demonstrating the concept of survivability and not on its computational speed
of its implementation. Since PowerDynamics.jl is still under development and
many simulations have to be carried out, differences in model structures or
calculation methods could potentially lead to large discrepancies in the
simulation results, compared to commercial tools. Therefore, a simulation tool
comparison between PowerDynamics.jl, MATLAB Simulink and DIgSILENT
PowerFactory is done first to show the accuracy of the PowerDynamics.jl
package and validate the model.
The paper is structured as follows. In Section II, a description about the
structure and control of the implemented ADN is given. In Section III, the
utilised tools are briefly described, where also their power flows and dynamic
simulations results with fault ride through (FRT) are compared. The
probabilistic stability assessment is carried out in Section IV and the paper
closes with a conclusion and outlook in Section V.

#### Notation

In this paper, capital letters like $U$ denote SI units, whereas small letters
$u$ denote per units and underlined letters $\underline{u}$ denote complex
values.

##  II Model and active distribution grid control

MV-01MV-02MV-03MV-04MV-05MV-08MV-10MV-11MV-07MV-06MV-09OLTCDGDGDGDGDGDGDGDGDGDGDG$P_{\mathrm{ref}}-P_{\mathrm{meas}}=\Delta
P$$Q_{\mathrm{ref}}-Q_{\mathrm{meas}}=\Delta Q$HV-00

Figure 1: 12-bus Cigré MV test network (own visualization, based on [21])

The model of the active distribution grid is based on the 12-bus version of
the CIGRE medium-voltage (MV) benchmark grid for renewable energy integration
found in [21] and shown schematically in Fig. 1. Compared to the 14-bus
version it omits the second branch consisting of buses 12 to 14 and their
second transformer. The parameters for the transformer, lines and loads are
taken from the benchmark grid [21]. The transformer tap is set to 10 and no
phase shift is considered. In the following simulations, the benchmark grid is
adapted such that the switches for the connection of MV-06 to MV-07 and MV-04
to MV-11 remain open and no power supply by DG units is considered during the
power flow comparison (cf. Section III-A). The operational limits of the
distributed generators are set to $0$ and $1\text{\,}\mathrm{MW}$ for active
power generation and $\pm$1\text{\,}\mathrm{Mvar}$$ for reactive power
generation. Depending on the simulation, the loads are either modelled as
static loads without voltage dependence or as static impedance. More detailed
descriptions are available in the cited sources.

###  II-A Interconnection power flow tracking

In this paper, the grid is referred to as an ADN due to its capability of
following reference functions for the active and reactive power flow at the
grid interconnection between the medium- and high-voltage grid. The control
architecture which makes this controlled behaviour possible was developed and
described in [1]. This interconnection power flow controller is based on a
distributed architecture: a distributed generator (DG) is connected to each
bus which receives global power flow exchange reference values
$P_{\mathrm{ref}}$ and $Q_{\mathrm{ref}}$ as well as measurements
$P_{\mathrm{meas}}$ and $Q_{\mathrm{meas}}$ taken at the medium-voltage side
of the transformer (bus MV-01). The global control errors $\Delta
P(t)=P_{\mathrm{ref}}(t)-P_{\mathrm{meas}}(t)$ and $\Delta
Q(t)=Q_{\mathrm{ref}}(t)-Q_{\mathrm{meas}}(t)$ are then subjected to local
integrator blocks at each medium-voltage bus which in turn yield the local
power injection reference signals $P_{\mathrm{r,l}}$ and $Q_{\mathrm{r,l}}$
determining the individual generator set points at that bus (cf. Fig. 2). This
generator set point is then saturated using the generator limits
$P_{\mathrm{max}}/Q_{\mathrm{max}}$ and $P_{\mathrm{min}}/Q_{\mathrm{min}}$.
This implies that there are eleven I-controllers for active and reactive
power, one at each medium-voltage bus, which all respond to the same control
error measured at the interconnection. Fig. 1 shows this signal flow using
dotted lines and Fig. 2 shows it for a single generator, the power flow
controller being on the left of the figure. Note that the other power flow
controllers receive the same $\Delta P$ and $\Delta Q$ signals, but their
integrators might yield different setpoints for each generator.

Within the DG unit control, a $dq$-transformation is used together with a
simplified model of a phase-locked loop to generate reference values for a
controlled current source. An overview of the DG control layout is illustrated
in Fig. 2. The integrators are all frozen if safe voltage limits
($\pm$$0.1\text{\,}\mathrm{pu}$ is taken here) are violated at any bus and the
grid then awaits stepping by the on-load tap change transformer before
integration of the control error continues [1].

Generator$\Delta P$$\Delta Q$$P_{\mathrm{max}}$,
$Q_{\mathrm{max}}$$\underline{U}_{\mathrm{PCC}}$$I_{\mathrm{d}}(P_{\mathrm{r},g},\underline{U})$$I_{\mathrm{q}}(Q_{\mathrm{r},g},\underline{U})$$\theta\left(\underline{U}_{\mathrm{PCC}}\right)$$I_{0}(I_{\mathrm{d,r}},I_{\mathrm{q,r}},\theta)$$I_{\mathrm{ref},a}(I_{0})$$I_{\mathrm{ref},c}(I_{0})$$P_{\mathrm{r},g}$$Q_{\mathrm{r},g}$ControlledCurrentSourcesCurrent
ControlPLL$\frac{e^{-sT_{\mathrm{t}}}}{Ts+1}$$P_{\mathrm{r,l}}$$Q_{\mathrm{r,l}}$$-$FRT$\frac{1}{s}$$-$$P_{\mathrm{ref}}$$Q_{\mathrm{ref}}$$P_{\mathrm{meas}}$$Q_{\mathrm{meas}}$limits$P_{\mathrm{min}}$,
$Q_{\mathrm{min}}$$I_{\mathrm{d,r}}$$I_{\mathrm{q,r}}$$I_{\mathrm{d}}$$I_{\mathrm{q}}$$\theta$$I_{\mathrm{ref},b}(I_{0})$

Figure 2: Distributed Generator Internal Layout

###  II-B Fault ride-through control

FRT controllers are designed to support voltages in case of short-circuit
events and prevent cascading loss of generation units. They accomplish this by
injecting a reactive current proportional to the bus voltage error in order to
support the grid voltage [22]. In the distributed generator model described
above, such a control mechanism is developed and integrated as an additional
block which modifies the currents $I_{\mathrm{d}}$ and $I_{\mathrm{q}}$
depending on the voltage at the point of common coupling
$\underline{U}_{\mathrm{PCC}}$. A schematic of the FRT control blocks is given
in Fig. 3. The first control block observes the magnitude of the voltage
$|\underline{U}_{\mathrm{PCC}}|$ at the point of common coupling. The output
$e=1$ represents an error flag. The dead band $U_{\mathrm{dead}}$ for the FRT
control is set to $\pm 10~{}\%$ of $\underline{U}_{\mathrm{PCC}}$. In such an
event the block also calculates the additional reactive current injection

$\displaystyle
I_{q+}=-k_{\mathrm{FRT}}(|\underline{U}_{\mathrm{PCC}}|-U_{\mathrm{ref}}\pm
U_{\mathrm{dead}}),$ (1)

with the sign depending on whether the deviation is above or below the
reference voltage (positive for negative deviation, as in line faults). The
reference voltage $U_{\mathrm{ref}}$ is set to $1\text{\,}\mathrm{pu}$. The
block FRT Limit then curtails $I_{\mathrm{d}}$ and $I_{\mathrm{q}}$ to the
maximum generator current (i.e. $I_{\mathrm{max}}=$1\text{\,}\mathrm{pu}$$)
with preference given to $I_{\mathrm{q}}$ during voltage deviations and
$I_{\mathrm{d}}$ otherwise.

$|\underline{U}_{\mathrm{PCC}}|$$e$$I_{\mathrm{q}+}$FRT
Limit$I_{\mathrm{d}}$$I_{\mathrm{q}}$$I_{\mathrm{d,r}}$$I_{\mathrm{q,r}}$$e$$I_{\mathrm{q}+}$$\Delta
U$1

Figure 3: Details of the Fault Ride-Through control in Fig. 2.

##  III Model validation

For a comprehensive numerical analysis of the model developed in Sec. II, the
power system and ADN control implementation in PowerDynamics.jl is validated
against established environments like DIgSILENT PowerFactory and MATLAB
Simulink with respect to the power flow solution and the short-circuit
behaviour inside the dynamic simulation. Consistent results despite
unavoidable tool-dependent differences in the model implementation indicate
robustness of the numerical procedure.

First, PowerDynamics.jl [20] is a package for dynamic modeling and simulation
of electrical power systems written in the high-performance programming
language Julia. The power flow is calculated by determining the respective
steady state of the dynamical equations using the Newton-Raphson algorithm in
combination with forward-mode automatic differentiation to determine the
Jacobian [23]. Numerical RMS simulations are performed with the Rodas4 solver,
a 4th-order stiff-aware Rosenbrock method for the RMS simulations provided by
DifferentialEquations.jl [24].

Second, DIgSILENT PowerFactory is a commercial power system analysis software
that offers various models and functions for power system analysis depending
on the implemented study case. For the load flow solution, the Newton-Raphson
algorithm based on power equations are chosen, while the RMS simulations are
solved with a trapezoidal method [25].

Finally, MATLAB Simulink is a programming environment for modeling, simulating
and analyzing multi-domain dynamical systems. To model and simulate electrical
power systems, the Simscape Power Systems toolbox is used. Power flow
calculations are based on the Newton-Raphson algorithm in power_loadflow.m.
Dynamic investigations within the RMS domain are performed using the
trapezoidal Heun method integrated in the MATLAB suite.

###  III-A Power flow comparison

The initial values of the state variables within dynamic simulations are based
on the steady state of the system. In electric power systems, the steady state
is defined by the complex bus voltage vector
$\underline{\mathbf{\mathit{U}}}_{PF}$ which is the solution of the power flow
equation determined by a power flow calculation (cf. [26]). Differences in the
power flow solution between different environments may arise even from subtle
differences in the modelling of the components used. A comparison of the
simulation tools regarding the ADN power flow results reveals that the
absolute deviations of the bus voltages are less than
${10}^{-4}\text{\,}\mathrm{pu}$ for the magnitudes and
${10}^{-3}\text{\,}\mathrm{\SIUnitSymbolDegree}$ for the phase angles. Only
the phase angle results in Simulink differ slightly from the other simulation
tools while an implementation in pure MATLAB based on [26] does not show this
deviation. In this comparison, the DG units do not inject any power and only
the loads are active. For the next sections, this power flow is used as the
starting point for the subsequent dynamic simulations. Here, the measured
power flow over the transformer with
$P_{\mathrm{meas}}=$24.373\text{\,}\mathrm{MW}$$ and
$Q_{\mathrm{meas}}=$6.115\text{\,}\mathrm{Mvar}$$ are used as the base
reference points $P_{\mathrm{ref}}^{0}$ and $Q_{\mathrm{ref}}^{0}$.

00.511.522.5${\text{\it{}t}}\text{/s}\rightarrow$3.50.20.40.60.8${\text{\it{}u}}\text{/p.u.}\uparrow$1.2SimulinkPowerFactoryPowerDynamics.jl

2.982.993${\text{\it{}t}}\text{/s}\rightarrow$3.020.50.7${\text{\it{}u}}\text{/p.u.}\uparrow$1.13.143.153.16${\text{\it{}t}}\text{/s}\rightarrow$3.180.9751${\text{\it{}u}}\text{/p.u.}\uparrow$1.05

Figure 4: Comparison of ADN Control and FRT behaviour in Simulink,
PowerFactory and PowerDynamics.jl for all bus voltages

###  III-B Dynamic simulation and short-circuit comparison

The agreement in the power flow solution extends to the behaviour of the ADN
control during short-circuit events. The dynamic simulation results of each
simulation tool are plotted in Fig. 4, where all bus voltages (except HV-00)
of the test system are illustrated. The short-circuit takes place at bus 8 at
$t=$3\text{\,}\mathrm{s}$$ with a fault-impedance of
$R_{\mathrm{on}}=$3.5\text{\,}\mathrm{\SIUnitSymbolOhm}$$ and is cleared at
$t=$3.15\text{\,}\mathrm{s}$$. Before the short-circuit, the ADN is subjected
to a reference step at $t=$1\text{\,}\mathrm{s}$$ from $P_{\mathrm{ref}}^{0}$
and $Q_{\mathrm{ref}}^{0}$ to $P_{\mathrm{ref}}=$24\text{\,}\mathrm{MW}$$ and
$Q_{\mathrm{ref}}=$5\text{\,}\mathrm{Mvar}$$ at the low-voltage side of the
transformer to demonstrate the capabilities of the control architecture. As it
can be seen by Fig. 4 the deviations during the dynamic simulations are very
small between all simulation tools. Only in MATLAB Simulink, small differences
can be detected during the occurrence and disappearance of the short circuit,
which can be attributed to numerical inaccuracies of the solver at these
discontinuous events. In general, the results are consistent and demonstrate
the applicability of the implementation in Julia for dynamic ADN simulations.
Therefore, the following probabilistic stability analysis is carried out in
PowerDynamics.jl due to its computational efficiency.

##  IV Probabilistic stability analysis

In this section, a probabilistic stability assessment, in particular the
survivability, of the dynamic ADN model is performed with respect to short-
circuit clearance.111The code can be found at github.com/strangeli/adn-
survivability . It can be shown that this recent concept from dynamical
systems theory yields complementary information about the stability of an
ADN’s operating point. In contrast to classical stability concepts such as the
critical clearing time of generators, the concept of survivability extends
analyses by covering a wide range of stability regions due to the large
examination space. Especially for the survivability, the monitoring of
operational limits of assets during transients are also included to further
improve the assessment of the overall system stability. Thus, simulations
which are numerically stable but would not occur in reality, e.g. due to
protection devices, can be better evaluated in the stability assessment. To
evaluate the short circuit’s stability impact at one of the buses, the VDE-
AR-N 4110 FRT limiting curve [22] shown in Fig. 5 is relied on. In particular,
survival is defined as a successful fault-ride through according to the
limiting curve (cp. Fig 5 and [22]).

The fault model considered here is an ensemble of short-circuits at medium-
voltage buses with a random fault duration and fault resistance. The fault
duration is taken from a normal distribution $\Delta
t_{\mathrm{fault}}=\mathcal{N}($150\text{\,}\mathrm{ms}$,$10\text{\,}\mathrm{ms}$)$
with a mean of $150\text{\,}\mathrm{ms}$ and a standard deviation of
$10\text{\,}\mathrm{ms}$. The fault resistance is randomly drawn from a
uniform distribution
$R_{on}=\mathcal{U}($3\text{\,}\mathrm{\SIUnitSymbolOhm}$,$10\text{\,}\mathrm{\SIUnitSymbolOhm}$)$.
Such a fault leads to a transient voltage drop as it is visualized in Fig. 4.

0.153${\text{\it{}t}}\text{/s}\rightarrow$600.15${\text{\it{}u}}\text{/p.u.}\uparrow$0.850.90connected(”survive”)disconnected(”not
survive”)

Figure 5: Low voltage fault-ride through limiting curve (own visualization,
based on [22])

In the next subsections, firstly the survivability of the whole system is
estimated for short-circuits at single nodes to demonstrate its general
concept. Secondly, the focus is set so the ADN control where a large envelope
of possible operating points for $P_{\mathrm{ref}}$ and $Q_{\mathrm{ref}}$ are
assessed according to their short-circuit survivability. In particular, the
second study shows how the survivability in combination with an efficient
simulation tool helps to investigate broad stability regions using large-
signal stability simulations. Here, the focus is only on the identification of
safe operating regions and not how they are interpreted or used by the grid
operator.

###  IV-A Single-node survivability

The concept of single-node survivability [10] is an example for a
probabilistic stability assessment of deterministic systems. In this paper, it
is defined as the conditional probability of survival at all buses given a
local fault at a single bus according to the setup above. Following the
approach in [10], the probability of bus $i$ to survive a random fault is the
single-node survivability value $\mu(i|R_{on},\Delta t_{\mathrm{fault}})$
which is estimated by using Monte-Carlo sampling. For this, each bus of Fig. 1
is subjected to a sample of $N=1000$ short-circuits from the fault model and
the cases $T_{i}$ that survive according to our definition are counted. Since
the survival of each simulation run is a binary outcome, a Bernoulli
experiment is performed with empirical mean $\mu(i|R_{on},\Delta
t_{\mathrm{fault}})=\tfrac{T_{i}}{N}$ and confidence interval bounded by
$1/(2\sqrt{N})\approx 1.6\%$ [27] independent of the state space dimension.

Figure 6: Single-node survivability

Fig. 6 shows the estimated single-node survivability for short circuits at the
respective bus on the x-axis for voltage-independent (PQ) and static impedance
(Z) loads. It is impacted both by the load model and the location of the
fault. The results show that bus 3 has the lowest survivability in case of
constant loads, which is not obvious, as it is more centrally connected
compared to other nodes which is usually beneficial in case of short-circuits.
However, a reason might be the high line impedance between buses 2 and 3. As
the short-circuit at bus 3 leads to low voltages at the buses behind, the
loads are drawing high currents to keep their power constant. All these
currents have to flow over the line 2-3 which results into an even bigger
voltage drop and thus decreasing the probability of a successful low-voltage
ride through. At the same time, the DG units are not capable of lifting the
fault voltage further within their FRT limits, which could have reduced the
flow between buses 2 and 3. This also explains the high survivability for
buses at the end of the feeder (e.g. bus 6 and 11), as only at these buses the
voltage is low and the voltage drop over the line. For constant impedance
loads, this effect cannot be observed, as the load current decreases. Overall,
the modeling of loads has a strong impact on the survivability in view of the
different results. This suggests that a survivability analysis can discover
non-intuitive effects, but can also confirm known relationships in terms of
stability.

###  IV-B $P_{\mathrm{ref}}$ and $Q_{\mathrm{ref}}$ envelope survivability

0.850.5pt5pt

Figure 7: Short-circuit survivability at bus MV-03 for static loads (top) and
for voltage-dependent loads (bottom)

From Fig. 6, it is shown that bus 3 has the lowest single-node survivability
with $\mu(3|R_{on},\Delta t_{\mathrm{fault}})\approx 0.881\pm 0.016$ (PQ),
defining the worst-case outcome. Therefore, it is chosen for a more detailed
analysis. In particular, the flexibility options available to the grid
operator with respect to changes in the reference set points
$P_{\mathrm{ref}}$ and $Q_{\mathrm{ref}}$ are of interest. To do this, the
single-node survivability analysis at bus 3 is combined with a step change in
the reference signal as in Fig. 4. After each step, a short circuit at bus 3
is simulated for those samples where the relative global control error $\Delta
P/P_{meas}$ is below 5%. This numerical experiment is carried out with
$10^{6}$ samples for $P_{\mathrm{ref}}$ and $Q_{\mathrm{ref}}$, demonstrating
the computational efficiency provided by Julia which allows such large
simulations to be performed at manageable cost. The output of a single
simulation is a boolean value for its survivability. Fig. 7 shows the voltage
survivability for a random short circuit after a reference step at $t=0.25$
$\text{\,}\mathrm{s}$ starting from the former power flow
$\\{P_{\mathrm{ref}}^{0},Q_{\mathrm{ref}}^{0}\\}$ to a random value taken from
a uniform distribution of
$P_{\mathrm{ref}}\in[P_{\mathrm{ref}}^{0}-$15\text{\,}\mathrm{MW}$;P_{\mathrm{ref}}^{0}+$15\text{\,}\mathrm{MW}$]$
and
$Q_{\mathrm{ref}}\in[Q_{\mathrm{ref}}^{0}-$10\text{\,}\mathrm{Mvar}$;Q_{\mathrm{ref}}^{0}+$10\text{\,}\mathrm{Mvar}$]$.
The survivability is calculated by a clustering of the sampled region into
overlapping squared cells of size
$$0.5\text{\,}\mathrm{MW}$\times$0.5\text{\,}\mathrm{Mvar}$$. For each cell,
the survivability is calculated by the number of surviving samples – according
to the definition and fault model in Sec. IV-A – divided by the total number
of samples in the cell (which is ensured to be greater or equal 100).

Hence, the results in the Fig. 7 show the feasibility region of the reference
steps for a fixed OLTC tap position and the variation of the single-node
survivability at bus 3 within this region for voltage-independent (Fig. 7,
top) and voltage-dependent loads (Fig. 7, bottom). The black area of
$\mu(3|R_{on},\Delta t_{\mathrm{fault}})\approx 0$ can be reached with a
reference step, but in a fault-case the FRT limits are easily violated. This
is plausible since the DG units have to provide a lot of inductive reactive
power and therefore the voltage in the grid decreases. With increasing active
power infeed by the DG units, the voltage rises, increasing the survivability
in turn. Comparing the figures, voltage-dependent loads have a larger
feasibility region than voltage-independent loads and they have in general
higher survivability within the feasibility region. This can be explained by
the reduced power requirement of voltage-dependent loads during the fault and
therefore they help maintain the voltage. New “forbidden” zones arise for
voltage-independent loads (upper left) and the boundaries between low and high
survivability are very steep near the edges. These boundaries identified here
may be worth studying in more detail in future work.

##  V Conclusions

In this paper, the safe operating point region with respect to fault
tolerance, a task expected to be of high importance to future ADN operators,
is investigated. Concretely, it could be demonstrated that probabilistic
methods can be an effective tool to handle increasing complexity in grid
operation. Such methods are highly adaptive to specific research questions and
can be used iteratively for exploring stability regions and their relation to
system characteristics. Here, the study based on short-circuit survivability
reveals numerous open questions, especially regarding the origin of the
stability thresholds. Future research should target alternative failure models
and parametrisations, e.g. parameter studies for limits on the short-circuit.
The proposed method and environment also allow for an extension to other
probabilistic stability measures [28, 29, 30] and for a comparison of various
network-local perturbations [13] to gain a further deeper understanding of the
stability properties of ADNs. Remarkably enough, it could also be shown that
high performance open source implementations like PowerDynamics.jl are
sufficiently capable and accurate to correctly capture the systems behaviour,
making comprehensive transient stability studies accessible. However, detailed
statements about the computational efficiency need to be adressed in future
research.

## Acknowledgment

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) – KU 837/39-1 / RA 516/13-1, project numbers 360248793,
360290054, 360460668, 359921210. The authors acknowledge the support of BMBF,
Condynet2 FK. 03EK3055A. All authors gratefully acknowledge the European
Regional Development Fund (ERDF), the German Federal Ministry of Education and
Research and the Land Brandenburg for supporting this project by providing
resources on the high performance computer system at the Potsdam Institute for
Climate Impact Research.

## References

  * [1] D. G. Mayorga, L. Robitzky, S. Liemann, U. Hager, J. Myrzik, and C. Rehtanz, “Distribution network control scheme for power flow regulation at the interconnection point between transmission and distribution system,” _2016 IEEE Innovative Smart Grid Technologies - Asia (ISGT-Asia)_ , pp. 23–28, Nov. 2016.
  * [2] P. Brzeski, M. Lazarek, T. Kapitaniak, J. Kurths, and P. Perlikowski, “Basin stability approach for quantifying responses of multistable systems with parameters mismatch,” _Meccanica_ , vol. 51, no. 11, pp. 2713–2726, Nov. 2016.
  * [3] H.-D. Chiang, _Direct Methods for Stability Analysis of Electric Power Systems_. Hoboken, NJ, USA,: John Wiley & Sons, Inc., 2010.
  * [4] T. L. Vu and K. Turitsyn, “Lyapunov functions family approach to transient stability assessment,” _IEEE Transactions on Power Systems_ , vol. 31, no. 2, pp. 1269–1277, Mar. 2016.
  * [5] L. Lovász and S. Vempala, “Simulated annealing in convex bodies and an o*(n4) volume algorithm,” _Journal of Computer and System Sciences_ , vol. 72, no. 2, pp. 392–417, 2006.
  * [6] A. Gajduk, M. Todorovski, and L. Kocarev, “Stability of power grids: An overview,” _The European Physical Journal Special Topics_ , vol. 223, no. 12, pp. 2387–2409, 2014.
  * [7] P. Brzeski, J. Wojewoda, T. Kapitaniak, J. Kurths, and P. Perlikowski, “Sample-based approach can outperform the classical dynamical analysis - experimental confirmation of the basin stability method,” _Scientific Reports_ , vol. 7, no. 1, p. 6121, Dec. 2017.
  * [8] D. Dudkowski, K. Czołczyński, and T. Kapitaniak, “Multistability and basin stability in coupled pendulum clocks,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 29, no. October, p. 103140, 2019.
  * [9] P. J. Menck, J. Heitzig, N. Marwan, and J. Kurths, “How basin stability complements the linear-stability paradigm,” _Nature Physics_ , vol. 9, no. 2, pp. 89–92, Feb. 2013.
  * [10] F. Hellmann, P. Schultz, C. Grabow, J. Heitzig, and J. Kurths, “Survivability of Deterministic Dynamical Systems,” _Scientific Reports_ , vol. 6, no. 1, Sep. 2016. [Online]. Available: http://www.nature.com/articles/srep29654
  * [11] J. von Neumann, “Various techniques used in connection with random digits,” _National Bureau of Standards, Applied Math Series_ , no. 12, pp. 36–38, 1951\.
  * [12] M. Evans and T. Swartz, _Approximating integrals via Monte Carlo and deterministic methods_. OUP Oxford, 2000, vol. 20.
  * [13] J. Nitzbon, P. Schultz, J. Heitzig, J. Kurths, and F. Hellmann, “Deciphering the imprint of topology on nonlinear dynamical network stability,” _New Journal of Physics_ , vol. 19, no. 3, p. 033029, Mar. 2017.
  * [14] M. F. Wolff, P. G. Lind, and P. Maass, “Power grid stability under perturbation of single nodes: Effects of heterogeneity and internal nodes,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 28, no. 10, p. 103120, 2018.
  * [15] Y. Feld and A. K. Hartmann, “Large-deviations of the basin stability of power grids,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 29, no. 11, p. 113103, Nov. 2019.
  * [16] H. Kim, M. J. Lee, S. H. Lee, and S.-W. Son, “On structural and dynamical factors determining the integrated basin instability of power-grid nodes,” _Chaos_ , vol. 29, no. 10, Jun. 2019.
  * [17] Z. Liu, X. He, Z. Ding, and Z. Zhang, “A basin stability based metric for ranking the transient stability of generators,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 3, pp. 1450–1459, Mar. 2019.
  * [18] L. Strenge, H. Kirchhoff, G. L. Ndow, and F. Hellmann, “Stability of meshed DC microgrids using Probabilistic Analysis,” _2017 IEEE Second International Conference on DC Microgrids (ICDCM)_ , pp. 175–180, Jun. 2017.
  * [19] J. F. Wienand, D. Eidmann, J. Kremers, J. Heitzig, F. Hellmann, and J. Kurths, “Impact of network topology on the stability of DC microgrids,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 29, no. 11, p. 113109, Nov. 2019.
  * [20] Potsdam Institute for Climate Impact Research, T. Kittel, S. Auer, W. Verschueren, and J. Liße, “JuliaEnergy/PowerDynamics.jl: v2.3.2,” Aug. 2020. [Online]. Available: https://github.com/JuliaEnergy/PowerDynamics.jl
  * [21] K. Strunz, “Benchmark systems for network integration of renewable and distributed energy resources,” _Task Force C6.04.02_ , 2013.
  * [22] Forum Netztechnik/Netzbetrieb im VDE, “VDE-AR-N 4110 - Technical requirements for the connection and operation of costumer installations to the medium voltage network (TAR medium voltage),” 2017.
  * [23] P. K. Mogensen, “JuliaNLSolvers/NLsolve.jl: v4.4.1,” Aug. 2020. [Online]. Available: https://doi.org/10.5281/zenodo.3986014
  * [24] C. Rackauckas and Q. Nie, “Differentialequations.jl – a performant and feature-rich ecosystem for solving differential equations in julia,” _Journal of Open Research Software_ , vol. 5, 5 2017.
  * [25] DIgSILENT GmbH, “PowerFactory User Manual Version 2020,” _Online Edition_ , February 2020.
  * [26] B. R. Oswald, _Berechnung von Drehstromnetzen_. Springer Vieweg, 2017, vol. 3.
  * [27] A. Agresti and B. A. Coull, “Approximate is better than “exact” for interval estimation of binomial proportions,” _The American Statistician_ , vol. 52, no. 2, pp. 119–126, 1998.
  * [28] S. Leng, W. Lin, and J. Kurths, “Basin stability in delayed dynamics,” _Scientific Reports_ , vol. 6, no. 1, p. 21449, Aug. 2016.
  * [29] M. Lindner and F. Hellmann, “Stochastic basins of attraction and generalized committor functions,” _arXiv Nonlinear Sciences_ , pp. 1–20, Mar. 2018.
  * [30] P. Schultz, F. Hellmann, K. N. Webster, and J. Kurths, “Bounding the first exit from the basin: Independence times and finite-time basin stability,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 28, no. 4, p. 043102, Apr. 2018.

  *[DERs]: distributed energy resource
  *[ADNs]: active distribution network
  *[ADN]: active distribution network
  *[FRT]: fault ride through
  *[DG]: distributed generator


# Locality-constrained autoregressive cum conditional normalizing flow for
lattice field theory simulations

Dinesh P. R Project Associate, Indian Institute of Science Education and
Research, Pune. Email<EMAIL_ADDRESS>
###### Abstract

Normalizing flow-based sampling methods have been successful in tackling
computational challenges traditionally associated with simulating lattice
quantum field theories. Further works have incorporated gauge and
translational invariance of the action integral in the underlying neural
networks, which have led to efficient training and inference in those models.
In this paper, we incorporate locality of the action integral which leads to
simplifications to the input domain of conditional normalizing flows that
sample constant time sub-lattices in an autoregressive process, dubbed local-
Autoregressive Conditional Normalizing Flow (l-ACNF). We find that the
autocorrelation times of l-ACNF models outperform an equivalent normalizing
flow model on the full lattice by orders of magnitude when sampling $\phi^{4}$
theory on a 2 dimensional lattice.

##  1 Introduction

Solving path integrals in quantum field theories for theories with large
couplings involves discretization of the underlying spacetime as lattice and
numerically sampling the fields using Markov Chain Monte Carlo (MCMC)
algorithms- referred to as lattice quantum field theory[9]. For large lattice
sizes and choices of action parameters that lead to small lattice spacing and
large correlation lengths, MCMC methods tend to suffer from long correlation
times leading to exponentially diverging computational costs- a phenomenon
known as critical slowing down (CSD)[17]. While a few non-local update
algorithms have been developed for specific models to address CSD [13, 16],
they cannot be applied for many key theories including quantum chromodynamics
(QCD).

In recent times, machine learning-based methods [19, 18] have been explored
for building generative models of statistical and field theories on a lattice.
In particular, a technique based on training normalizing flows followed by
independent Metropolis-Hastings sampling [2] was shown to nearly eliminate CSD
while modelling the $\phi^{4}$ theory in 2 dimensions. Subsequent works [7, 3]
proposed ways of incorporate gauge invariance of the action integral when
modelling gauge field theories using normalizing flows. This eliminated
unnecessary degrees of freedom in the variational distribution modeled by the
flow when local $SU(N)$ symmetries are present and reducing computational
costs of training and inference of the model. This addresses a well-known
problem in machine learning traditionally known as the curse of
dimensionality[8] where training time and model complexity are expected to
increase exponentially with the size (dimensionality) of input/output data.
Even though deep neural networks are shown to break this curse in many
cases[11], this is only valid under approximation that the target function
depends on low-dimensional projections of the input[4]. This is a pretty bad
approximation for generative models of lattice field theories when correlation
lengths are proportional to the lattice size, hence eliminating unnecessary
degrees of freedom is of vital importance. Another study [6] explores several
ways of tackling poor convergence in normalizing flows when the underlying
distribution is multimodal- like in broken symmetry phases. It also introduces
masked convolutional networks as conditional maps of coupling layers that
enforce (broken) translational invariance and a notion of locality with
stronger local dependence within the distribution.

In this work, we propose an autoregressive process of conditional normalizing
flows which enforces locality of the action integral more strictly. By doing
so, we reduce the effective input size of the conditional normalizing flow to
$O(L^{d-1})$ from the upper bound of $L^{d}$, hence further eliminating
unnecessary degrees of freedom in the spirit of tackling the curse of
dimensionality. We dub this _local-Autoregressive Conditional Normalizing
Flow_ (l-ACNF). After briefly introducing the system of interest and notations
in sec. 2.1, we analyse the dependencies of autoregressive conditional
distrbutions (introduced in sec. 2.2) in sec. 2.3 in fair detail. Sec. 2.4
discusses the mathematical structure of the l-ACNF model along with the
accompanying Metropolis-Hamilton process. In sec. 3, we describe the numerical
implementation for sampling the $\phi^{4}$ theory in a 2D lattice using
l-ACNF.

##  2 Mathematical background

###  2.1 Scalar lattice field theory

The system of interest consists of a hypercubic lattice of length $L$ in $d$
dimensions where every position is labeled using a $d$-dimensional vector
$\bm{r}\in[1,L]^{d}$. For open boundary conditions, the lattice positions
$\bm{r}\notin[1,L]^{d}$ do not exist. In case of periodic boundary conditions,
we map each component of $\bm{r}\in\mathbb{Z}^{d}$ to
$r_{i}\rightarrow(r_{i}-1)\mod L+1$, . State/configuration of the system is
described using scalar field values $\phi(\bm{r})$. The configurations obey
the Boltzmann distribution:

$p\left(\phi(\bm{r})\right)=e^{-S[\phi]}/Z$ (1)

where the _action_ $S[\phi]$ is a functional of the field $\phi$ and $Z$ is a
normalizing constant known as partition function. For the scalar field theory
with $\phi^{4}$ interactions, the action is given by:

$S[\phi]={\displaystyle\sum_{\bm{r}\in[1,L]^{d}}}\left[\phi(\bm{r}){\displaystyle\sum_{\bar{\bm{r}}}}\boxempty(\bm{r},\bar{\bm{r}})\phi(\bar{\bm{r}})+m^{2}\phi(\bm{r})^{2}+\lambda\phi(\bm{r})^{4}\right]$
(2)

where $m$, $\lambda$ are the bare mass and coupling. The d’Alembertian in the
lattice approximation is given by:

$\phi(\bm{r}){\displaystyle\sum_{\bar{\bm{r}}}}\boxempty(\bm{r},\bar{\bm{r}})\phi(\bar{\bm{r}})=\sum_{\mu=1}^{d}2\phi(\bm{r})^{2}-\phi(\bm{r})\phi(\bm{r}-\hat{\mu})-\phi(\bm{r})\phi(\bm{r}+\hat{\mu})$

where $\phi(\bm{r})$ can take any real value.

The action $S[\phi]$ in 2 has many interesting mathematical properties. It
contains only nearest neighbour product/interaction terms
$\phi(\bm{r})\phi(\bm{r}-\hat{\mu})$ and $\phi(\bm{r})\phi(\bm{r}+\hat{\mu})$,
besides powers of $\phi(\bm{r})$. This makes the action _local_ , a common
feature in many fundamental physical theories. The action is also
symmetric/invariant to the transformation
$\phi(\bm{r})\rightarrow-\phi(\bm{r})$, along with various discrete rotations
and reflections of the lattice (hypercubic symmetry group).

###  2.2 Autoregressive process

The Boltzmann distribution in 1 depends on a large number of random variables
$L^{d}=N$ and it’s usually difficult to sample from it or evaluate the PDF
directly. One way to circumvent is to model it as a product of conditional
distribtutions over individual random variables. Let’s first replace the
labels of lattice positions from vectors $\bm{r}\in[1,L]^{d}$ to an integer
_ordering_ $k\in[1,N]$. We’ll use a particular ordering that maps from vector
labels as:

$k(\bm{r})=\left({\displaystyle\sum_{i=1}^{d}}(r_{i}-1)L^{i-1}\right)+1$ (3)

Using a chain rule obtained from repeated application of Bayes theorem, we can
decompose the probability into a product of $N$ conditional probabilities, in
the above ordering:

$\displaystyle p(\\{\phi_{k}|k\in[1,N]\\})$
$\displaystyle=p(\phi_{1},\phi_{2}\dots\phi_{N})=p(\phi_{1})p(\phi_{2}|\phi_{1})\dots
p(\phi_{N}|\phi_{N-1}\dots,\phi_{1})$ $\displaystyle\implies\log
p(\\{\phi_{k}|k\in[1,N]\\})$ $\displaystyle=\log p(\phi_{1})+\log
p(\phi_{2}|\phi_{1})+\dots\log p(\phi_{N}|\phi_{N-1}\dots,\phi_{1})$ (4)
$\displaystyle=\sum_{k\in[1,N]}\log p(\phi_{k}|\phi_{<k})$ (5)

This is referred to as an _autoregressive process_. If we can evaluate and
sample the conditional probabilities sequentially, this would yield a sample
and log-probability of the entire Boltzmann distribution. Exact forms of
$p(\phi_{k}|\phi_{<k})$ are still analytically intractable starting from 1. So
approaches using this principle usually optimize the KL-divergence between a
machine learning-based variational distribution and the unnormalized Boltzmann
distribution[19].

###  2.3 $d-1$ dimensional dependency sets for local interactions

##### Ising model and open boundary conditions

Examining the $k$th conditional probability $p(\phi_{k}|\phi_{<k})$ in 5, its
distribution in general depends on $k-1$ values in
$\phi_{<k}=\\{\phi_{k-1},\dots\phi_{1}\\}$. This means the complexity of these
distributions can explode if the number of lattice points $N$ is large, which
is typically the case of interest. However for systems with nearest neighbour
interactions, the _dependency set_ is significantly smaller [12]. It’s easier
to show this (without loss of generality) for the nearest neighbour Ising
model whose action is given by:

$S[\phi]=-\beta
J\sum_{\mu=1}^{d}\sum_{\bm{r}}\phi(\bm{r}-\hat{\mu})\phi(\bm{r})$ (6)

where $\phi(\bm{x})$ takes values $\pm 1$. We’ll assume _open_ boundary
conditions for now and relax it later on. Restating the Boltzmann distribution
for the Ising model as an autoregressive process111$k-\hat{\mu}$ should be
understood as the lattice position $\bm{x}-\hat{\mu}$ where $\bm{x}$ maps to
$k$ according to the given ordering:

$\prod_{k=1}^{N}p(\phi_{k}|\phi_{<k})=p(\phi)=\exp\left(-\beta
J\sum_{\mu=1}^{d}\sum_{k}\phi_{k}\phi_{k-\hat{\mu}}\right)/Z$ (7)

From Bayes theorem, we can relate this conditional probability to the
unconditional joint probabilities of the first $k$ and $k-1$ spins, which can
in turn be written as integrated forms of the Boltzmann distribution:

$p(\phi_{k}|\phi_{<k})=\frac{p(\phi_{1},\dots\phi_{k})}{p(\phi_{1},\dots\phi_{k-1})}=\frac{{\displaystyle\sum_{\phi_{k+1}\dots\phi_{N}}}p(\phi)}{{\displaystyle\sum_{\phi_{k}\dots\phi_{N}}}p(\phi)}$

Expanding the $p(\phi)$ for the Ising model:

$\displaystyle p(\phi_{k}|\phi_{<k})$
$\displaystyle=\frac{{\displaystyle\sum_{\phi_{N},\dots\phi_{k+1}}}\exp\left(-\beta
J{\displaystyle\sum_{l=k}^{N}}\left(\phi_{l}{\displaystyle\sum_{\mu}}\phi_{l-\hat{\mu}}\right)+\delta(\phi_{<k})\right)}{{\displaystyle\sum_{\phi_{N},\dots\phi_{k}}}\exp\left(-\beta
J{\displaystyle\sum_{l=k}^{N}}\left(\phi_{l}{\displaystyle\sum_{\mu}}\phi_{l-\hat{\mu}}\right)+\delta(\phi_{<k})\right)}$

Since the values in $\phi_{<k}$ are fixed and not summed over, the terms
$\delta(\phi_{<k})$ containing only them cancel from both the numerator and
denominator, leaving us with:

$\displaystyle p(\phi_{k}|\phi_{<k})$
$\displaystyle=\frac{{\displaystyle\sum_{\phi_{N},\dots\phi_{k+1}}}\exp\left(-\beta
J{\displaystyle\sum_{l=k}^{N}}\phi_{l}{\displaystyle\sum_{\mu}}\phi_{l-\hat{\mu}}\right)}{{\displaystyle\sum_{\phi_{N},\dots\phi_{k}}}\exp\left(-\beta
J{\displaystyle\sum_{l=k}^{N}}\phi_{l}{\displaystyle\sum_{\mu}}\phi_{l-\hat{\mu}}\right)}$
(8)

Even though $\phi_{<k}$ contains $k-1$ values, the conditional probability
$p(\phi_{k}|\phi_{<k})$ depends only on those positions within $\phi_{<k}$
that are nearest neighbours of the positions in $\phi_{\geq k}$. We can draw
the same conclusion for scalar lattice field theory by replacing the sums with
integrals and including terms like $\phi_{l}^{2}$ and $\phi_{l}^{4}$ in the
above expression. The number of elements in the dependency set is bounded
above by $L^{d-1}$ or $N/L$ for our choice of ordering (see figure 1a for an
illustration on a $10\times 10$ 2D lattice) which is orders of magnitude
smaller than the original upper bound $N$. In fact, we can join the 2 strips
of blue spins in figure 1a into a single 1D line of length $L^{2-1}=10$, and
the conditional distribution on $\phi_{k}$ simply depends on the values along
this line.

(a) The conditional probability $p(\phi_{k}|\phi_{<k})$ of the field at red
position depends only on the nearest neighbours of the positions in
$\phi_{\geq k}$ (coloured red/white), within $\phi_{<k}$. The dependency set
consists only of blue positions while the black positions within $\phi_{<k}$
can be ignored. (b) The conditional probability of the field at red positions
at a constant time $p(\\{\phi(x,t)|t=5\\}|\\{\phi(x,t)|t<5\\})$ depends only
on the positions with $t=4$ marked blue $\\{\phi(x,t)|t=4\\}$ and not the
remaining ones marked black.

Figure 1: Dependency sets for the conditional probabilities of field value at
single position (1a) and joint conditional probability of positions in a
constant time manifold (1b), assuming open boundary conditions.

We can draw the same conclusion for scalar lattice field theory by replacing
the sums with integrals and including terms like $\phi_{l}^{2}$ and
$\phi_{l}^{4}$ in 8. Let’s split the lattice vector $\bm{r}$ into a $d-1$
dimensional _spatial_ vector $\bm{x}=[x_{1},\dots x_{d-1}]$ and _time_ $t$ so
that $\bm{r}=[r_{1}=x_{1},\dots r_{d-1}=x_{d-1},r_{d}=t]$. In figure 1a for a
$10\times 10$ lattice, the dependency set of blue positions is them simply a
line of length $L$ along the spatial dimension $x$. More generally, the
dependency set of $\phi_{k(\bm{r})}$ can be cast as a $d-1$ dimensional
spatial sub-lattice that’s defined parametrically using:

$B_{\bm{r}}^{o}=\left\\{\bm{r}^{\prime}|\bm{r}^{\prime}=\begin{cases}[\bm{x}^{\prime},t]&\text{if
}k([\bm{x}^{\prime},t])<k(\bm{r})\\\ {}[\bm{x}^{\prime},t-1]&\text{if
}k([\bm{x}^{\prime},t])>k(\bm{r})\end{cases}\right\\}$ (9)

We’ll call this the _dependency surface_ at $\bm{r}$. From 8, we can write
down the joint conditional probabilities for more than one variable. For
example, we can write down this $m$ variable joint distribution for the Ising
model:

$\displaystyle
p\left(\left\\{\phi_{k},\dots\phi_{k+m-1}\right\\}|\phi_{<k}\right)$
$\displaystyle=p(\phi_{k}|\phi_{<k})\dots p(\phi_{k+m-1}|\phi_{<k+m-1})$
$\displaystyle=\frac{{\displaystyle\sum_{\phi_{N},\dots\phi_{k+m}}}\exp\left(-\beta
J{\displaystyle\sum_{l=k}^{N}}\left(\phi_{l}{\displaystyle\sum_{\mu}}\phi_{l-\hat{\mu}}\right)\right)}{{\displaystyle\sum_{\phi_{N},\dots\phi_{k}}}\exp\left(-\beta
J{\displaystyle\sum_{l=k}^{N}}\left(\phi_{l}{\displaystyle\sum_{\mu}}\phi_{l-\hat{\mu}}\right)\right)}$
(10)

For the choice of $\bm{r}=[1,\dots 1,t]$ and $m=L^{d-1}$, the set
$\left\\{\phi_{k},\dots\phi_{k+m-1}\right\\}$ is the _constant time sub-
lattice_ $\\{\bm{r}^{\prime}|t^{\prime}=t\\}$, the dependency surface is the
constant time manifold at the previous time step:

$B_{t}^{o}=\left\\{\bm{r}^{\prime}|t^{\prime}=t-1\right\\}$ (11)

See figure 1b for an example in a $10\times 10$ lattice. This essentially
reduces the overall autoregressive process to a first order Markov process in
the time dimension.

##### Periodic boundary conditions

The general observation that the dependency surface of $\phi_{k(\bm{r})}$
contains only the nearest neighbours of $\phi_{\geq k}$ within $\phi_{<k}$ is
true regardless of boundary conditions as long as the action is local.
However, the dependency surface consists of two $n-1$ dimensional sub-lattices
in case of periodic boundary conditions:

$B_{t}^{p}=B_{t}^{o}\cup\left\\{\bm{r}|r_{d}=1\right\\}$ (12)

which includes the sub-lattice corresponding to the _initial time_ sub-
lattice. Note that periodic boundaries topologically renders the lattice as
$d$ dimensional torus, hence $B_{t}^{o}=\\{r_{d}=t-1\\}$ would be a $d-1$
torus along the spatial dimensions. This torus effectively “blocks” the
influence of all equal time $r_{d}<t$ tori except for the initial one
$r_{d}=1$ which is connected to $r_{d}>t$ on the other side. This is
illustrated for $d=2$ in the figure 2. In other words, only the sub-lattices
at $r_{d}=t-1$ and $r_{d}=1$ are _connected_ to the sub-lattice $r_{d}=t$
either directly or via unknown/future positions.

Figure 2: In case of a periodic boundary, distributions of field values at a
single red position (top) or positions in a constant time sub-lattice (bottom)
are conditioned on the dependency sets (marked blue) that are similar to that
of an open boundary, but including the initial time sub-lattice at $t=1$.

We’ll assume periodic boundary conditions in our models since we can take
advantage of the translational invariance of the action222Circular permutation
invariance is the accurate description for finite lattices with periodic
boundary condition, while true translational invariance strictly occurs only
for an infinite lattice..

###  2.4 Autoregressive Conditional Normalizing Flows and Metropolis sampling

####  2.4.1 Conditional masked normalizing flow

Conditional normalizing flow[15] models a bijective map $f_{\theta}$ between
samples of a standard prior distribution $z\sim r(z)$ to samples of a
variational distribution conditioned on $\chi$,
$\tilde{\phi}\sim\tilde{q}_{f,\theta}(\tilde{\phi}|\chi)$, so that
$\tilde{\phi}=f_{\theta}(z)$. From the change of variables formula for
probability distributions, we can obtain the log PDF of the variational
distribution:

$\displaystyle\log r(z)$
$\displaystyle=\log\tilde{q}_{f,\theta}(f_{\theta}(z)|\chi)+\log\left|\det\frac{\partial
f_{\theta}(z)}{\partial z}\right|$
$\displaystyle\implies\log\tilde{q}_{f,\theta}(f_{\theta}(z)|\chi)$
$\displaystyle=\log r(z)-\log\left|\det\frac{\partial f_{\theta}(z)}{\partial
z}\right|$ (13)

The map $f_{\theta}$ is constructed using a sequence of $n$ coupling layers
$g^{(i)}$ whose inputs are acted alternatively by complementary binary masks
$m$ and $1-m$333All products are element-wise unless stated otherwise:

$\displaystyle f_{\theta}(z)$ $\displaystyle=\left(g_{1}\circ\dots\circ
g_{n}\right)(z)$ $\displaystyle g_{i}(z_{i},\chi)$
$\displaystyle=m^{(i)}(z_{i})+(1-m^{(i)})\left(\tilde{\phi}e^{s^{(i)}(m^{(i)}(z_{i}),\chi)}+t^{(i)}(m^{(i)}(z_{i}),\chi)\right)$
$\displaystyle m^{(i)}$ $\displaystyle=\begin{cases}m&i\text{ even}\\\
1-m&i\text{ odd}\end{cases}$ (14)

where $z_{1}=z$ and $z_{i}=g_{i-1}(z_{i-1},\chi)$ for $i>1$, so
$z_{n+1}=\tilde{\phi}$. The functions $s^{(i)}$ and $t^{(i)}$ are in turn
modeled using feed-forward neural networks. The log determinant in 13 then
becomes:

$\log\left|\det\frac{\partial f_{\theta}(z)}{\partial
z}\right|=\sum_{i=1}^{n}\log\left|\det\frac{\partial g_{i}(z_{i})}{\partial
z_{i}}\right|=\sum_{i=1}^{n}\text{sum}(s^{(i)}(m^{(i)}(z_{i}),\chi))$

where the sum function adds up an input array of values. The couplings can be
inverted and composed together for the inverse flow:

$\displaystyle g_{i}^{-1}(z_{i+1},\chi)$
$\displaystyle=m^{(i)}(z_{i+1})+(1-m^{(i)})\left(z_{i+1}+t^{(i)}(m^{(i)}(z_{i+1}),\chi)\right)e^{-s^{(i)}(m^{(i)}(z_{i+1}),\chi)}$
$\displaystyle f_{\theta}^{-1}(\tilde{\phi})$
$\displaystyle=\left(g_{n}^{-1}\circ\dots\circ
g_{1}^{-1}\right)(\tilde{\phi})$ (15)

which is useful to evaluate the log PDF of a known sample $\tilde{\phi}$ of
the variational distribution:

$\log\tilde{q}_{f,\theta}(\tilde{\phi}|\chi)=\log
r(f_{\theta}^{-1}(\tilde{\phi}))-\sum_{i=1}^{n}\text{sum}(s^{(i)}(m^{(i)}(z_{i+1}),\chi))$
(16)

####  2.4.2 Autoregressive process

We split the Boltzmann distribution as a product of conditional distributions
of the field at single positions in 5. We then generalized this and evaluated
the dependency sets for joint conditional distribution of the field values in
a constant time sub-lattice in 11. Let’s rewrite this joint conditional
distribution for periodic boundary conditions, but this time conditioned on
the coupling $\lambda$ in 2 as well:

$\displaystyle p\left(\\{\phi(\bm{r})|t\\}|\phi_{<[1,\dots t]},\lambda\right)$
$\displaystyle=p\left(\\{\phi(\bm{r})|t\\}|\\{\phi(\bm{r}^{\prime})|t^{\prime}=t-1,1\\}\cup\lambda\right)$
$\displaystyle\implies\log p(\phi|\lambda)$ $\displaystyle=\sum_{t=1}^{L}\log
p\left(\\{\phi(\bm{r})|x_{d}=t\\}|\\{\phi(\bm{r})|x_{d}=t-1,1\\}\cup\lambda\right)$

where the second line is the full Boltzmann distribution in logarithmic form.
We use the output variational distribution from a conditional normalizing flow
$\tilde{q}_{\theta}(\tilde{\phi}|\chi,t)$ to model the individual conditional
distributions in the RHS, where
$\chi=\\{\phi(\bm{r})|x_{d}=t-1,1\\}\cup\lambda$. We sample one constant time
manifold of $\tilde{\phi}(t)=\\{\phi(\bm{r})|x_{d}=t\\}$ at a time and supply
them as inputs to subsequent calls of the normalizing flow:

$\displaystyle\tilde{\phi}(1)$
$\displaystyle\sim\tilde{q}_{\theta}(\tilde{\phi}|\Theta,\Theta,\lambda,1)=\tilde{q}_{\theta}(1)$
$\displaystyle\vdots$ $\displaystyle\tilde{\phi}(k)$
$\displaystyle\sim\tilde{q}_{\theta}(\tilde{\phi}|\tilde{\phi}(k-1),\tilde{\phi}(1),\lambda,k)=\tilde{q}_{\theta}(k)$
$\displaystyle\vdots$ $\displaystyle\tilde{\phi}(L)$
$\displaystyle\sim\tilde{q}_{\theta}(\tilde{\phi}|\tilde{\phi}(L-1),\tilde{\phi}(1),\lambda,L)=\tilde{q}_{\theta}(L)$
(17)

where $\theta$ are neural network weights and $\Theta$ are junk/missing
values444We supply zero fields along with a boolean flag for the initial
sampling. See code implementation for details.. Note that we explicitly suppy
time as an input since the conditional distributions aren’t necessarily time
independent/invariant. Finally, we can concatenate the samples along the time
direction and add the log PDFs to obtain the sample and corresponding log PDF
of the final model distribution $q_{\theta}(\phi)$ of l-ACNF:

$\displaystyle\phi$ $\displaystyle=\text{cat}_{t=1}^{L}\tilde{\phi}(t)$
$\displaystyle\log q_{\theta}(\phi)$
$\displaystyle=\sum_{t=1}^{L}\log\tilde{q}_{\theta}(t)$

We optimize the weights $\theta$ such that the KL-divergence between the
distributions $q_{\theta}$ and $p$ is minimized:

$\displaystyle\text{KL}(p\bigparallel q_{\theta})$
$\displaystyle=\mathbb{E}_{q_{\theta}(\phi)}\left[\log q_{\theta}(\phi)-\log
p(\phi)\right]$ $\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\left[\log
q_{\theta}(\phi^{(i)})+S[\phi^{(i)}]\right]\qquad\phi^{(i)}\in\\{\phi\sim
q_{\theta}\\}^{M}$ (18)

where we sample mini-batches of $\phi$ (of size $M$), evaluate the mean of the
quantity in square brackets and use the gradient of this quantity wrt $\theta$
to update the weights at every step. We ignore the log of partition function
$\log Z$ in the expansion of $\log p(\phi)$ since it’s a constant independent
of $\phi$.

####  2.4.3 Checkerboard masks and convolutional networks

Particular choices of neural network architecture for $s^{(i)}$and $t^{(i)}$,
and the binary mask $m$ used in 14 have algorithmic benefits in the context of
generative modeling of the scalar lattice field theory. Neural networks
containing only $d-1$ dimensional convolutional layers with circular padding
(along with activations and residual connections) and checkerboard mask for
$m$ enable a broken translational invariance along the $d-1$ spatial
dimensions of the lattice[6] (as long as it’s a symmetry of the prior
distribution $r(z)$ in 13). No such invariances exist explicitly along the
time dimension since we sequentially sample along this direction, though a
well-trained model is expected to approximately satisfy them. Fully
convolutional networks also enable sampling configurations independent of the
lattice size $L$ along the spatial dimensions. Autoregressive sampling in 17
features dependency sets independent of $L$ along the time dimension, hence we
can sample lattices of different sizes compared to the $L$ that was used to
train the model- making it _scalable_. $L$ is still constrained to be an even
integer in order to avoid “leaks” in the normalizing flow due to checkerboard
masks.

For benchmarking our model, we also consider a masked normalizing flow with
$d$ dimensional convolutional layers and circular padding where the NF
generates entire lattice configurations (instead of constant time sub-
lattices) conditioned only on the coupling $\lambda$, similar to the
implementation described in [6]. This admits translational invariance along
all $d$ directions and scalable as well. l-ACNF has the (theoretical)
advantage of better tackling the _curse of dimensionality_ with smaller
effective input size ($O(L^{d-1})$) compared to the latter ($L^{d}$).

####  2.4.4 Metropolis-Hastings sampling

Any bias in the variational distribution after training to minimize the KL-
divergence in 18 results in bias for various observables/estimators evaluated
after sampling from it. To get unbiased samples of $p$ from $q_{\theta}$, we
run its samples through a Metropolis-Hastings procedure- a Markov Chain Monte
Carlo (MCMC) technique. Starting from an initial sample $\phi^{(1)}$ and
corresponding log PDF $\log q_{\theta}(\phi^{(1)})$, we update the chain in
the $i$th step using current sample $\phi^{(i)}$ and the proposed sample
$\phi^{\prime}$ using:

$\displaystyle A(\phi^{(i-1)},\phi^{\prime})$
$\displaystyle=\min\left(1,\exp\left(\log q_{\theta}(\phi^{(i-1)})-\log
p(\phi^{(i-1)})-\log q_{\theta}(\phi^{\prime})+\log
p(\phi^{\prime})\right)\right)$ $\displaystyle=\min\left(1,\exp\left(\log
q_{\theta}(\phi^{(i-1)})-S[\phi^{(i-1)}]-\log
q_{\theta}(\phi^{\prime})+S[\phi^{\prime}]\right)\right)$
$\displaystyle\phi^{(i)}$ $\displaystyle=\begin{cases}\phi^{\prime}&A<\pi\\\
\phi&\text{else}\end{cases}\qquad\text{where }\pi\sim\text{Uniform}(0,1)$ (19)

A necessary condition for the asymptotic convergence of Metropolis-Hastings
samples to $p$ is _ergodicity_ which requires:

$q_{\theta}(\phi)>0\quad\text{or}\quad\log
q_{\theta}(\phi)>-\infty\quad\forall\text{$\phi$}\in\mathbb{R}^{L^{d}}$

The normalizing flow in $q_{\theta}$ satisfies this condition if the prior
distribution of the flow satisfies $r(z)>0$ $\forall z\in\mathbb{R}$ which is
true in case of a Gaussian prior, for example.

####  2.4.5 Symmetrizing and adiabatic retraining

In the Metropolis-Hastings update step 19, it’s useful to ensure various
symmetry transformations of the proposed configuration (that render the action
$S$ invariant) leave the acceptance probability $A$ unchanged. This includes
translations and reflections of the lattice, translations along every
dimension and the negation $\phi\rightarrow-\phi$. This amounts to averaging
over the PDF $\text{Avg}_{\mathfrak{g}}q_{\theta}(\mathfrak{g}(\phi))$ after
applying the transformations $\mathfrak{g}$ and choosing a randomly
transformed field at every step $\mathfrak{g}(\phi)$:

$\displaystyle\log q_{\theta}(\phi)$
$\displaystyle\rightarrow\text{logsumexp}_{\mathfrak{g}\in G}\left(\log
q_{\theta}\left(\mathfrak{g}(\phi)\right)\right)-\log\left|G\right|$
$\displaystyle\phi$
$\displaystyle\rightarrow\mathfrak{g}(\phi)\qquad\mathfrak{g}\sim G$ (20)

where $G$ is the set of all symmetry transformations we consider. This is one
way of tackling poor convergence when the distribution $p(\phi)$ has well-
separated modes as studied in [6]. In our case, rotations between the time and
spatial dimensions also address the inherent asymmetrical treatment between
spatial and time dimensions as well as lack of invariances along time
dimension in l-ACNF, as mentioned in 2.4.3. The scalability of the model to
sample lattice sizes different from the value it trained on allows us to
perform adiabatic retraining- where we initially train with smaller $L$ and
use this network as an initial state for training larger $L$, similar to [3].

##  3 Numerical experiments

###  3.1 Network architecture and training

The python deep learning library Pytorch[10] was used to build all models used
in this study. We model the $\phi^{4}$ theory in a 2 dimensional lattice
($d=2$) with length 16 ($L=16$). For Model 1, the conditional normalizing flow
network contains $n=8$ coupling layers described in 14. As mentioned in 2.4.3,
we use 1D checkerboard binary masks $m$ and 1D gated convolutional
networks[14] for the functions $s^{(i)}$ and $t^{(i)}$, where the conditional
variables $\chi$ are processed by a common chain of 4 gated-convolutional
networks and passed as a skip connection to another gated-convolutional
network with input $m^{(i)}z_{i}$. The model is trained using the KL-
divergence loss in 18 after setting $m^{2}=-4$ and
$\lambda\sim$Uniform$(4.8,5.4)$ is sampled at every step- for 9000 steps using
Adam optimizer and learning rate $\eta=10^{-3}$ with Step schedule for
$\approx 90$ minutes. Model 2 was constructed by adiabatic-retraining Model 1,
with $L=20$ for $\approx 40$ minutes using $\eta=10^{-5}$ with the same
optimizer and schedule. Model 3 is constructed using a 2D masked normalizing
flow as described in the second paragraph of 2.4.3 for benchmarking Model 1.
It uses 2D checkerboard masks and 2D gated convolutional networks with
circular padding for $n=12$ coupling layers. This ensures equal number of
convolutional operations and non-linear activations across all models, with
slightly higher number of weights due to 2D convolutional layers. Model 3 was
trained with $L=16$ for $\approx 90$ minutes using $\eta=10^{-3}$ to maintain
equivalence with Model 1.

All training and inference was done using an Nvidia RTX 2060 mobile GPU.
Source code is available at https://github.com/dinesh110598/l-ACNF-scalar. The
data and code for all plots is available in an interactive Julia notebook
here: https://github.com/dinesh110598/l-ACNF-
scalar/blob/main/Figures/notebook.jl.

###  3.2 Parameters and observables

The parameters $m^{2}$ and $\lambda$ in 2 were chosen close to the
$m_{p}L\approx 4$ critical line at the symmetrical phase of the system,
similar to the study in [2]. The difference is we use a single trained model
to evaluate all observables at different lattice lengths $L$ and couplings
$\lambda$. Observables of interest include the 2-point Green’s function:

$G_{c}(\bm{r})=\frac{1}{N}\sum_{\bm{r}^{\prime}}\left(\langle\phi(\bm{r}^{\prime})\phi(\bm{r}^{\prime}+\bm{r})\rangle-\langle\phi(\bm{r}^{\prime})\rangle\langle\phi(\bm{r}^{\prime}+\bm{r})\rangle\right)$

its spatial fourier transform

$\tilde{G}_{c}(\bm{p},t)=\frac{1}{L^{d-1}}\sum_{\bm{x}}e^{i\bm{p}.\bm{x}}G_{c}(\bm{x},t)$

an estimator of pole mass

$m_{p}(t)=\text{arccosh}\left(\frac{\tilde{G}_{c}(0,t-1)+\tilde{G}_{c}(0,t+1)}{2\tilde{G}_{c}(0,t)}\right)$
(21)

the 2-point susceptibility

$\chi_{2}=\sum_{\bm{r}}G_{c}(\bm{r})$

and the average Ising energy density

$E=\frac{1}{d}\sum_{\mu=1}^{d}G_{c}(\hat{\mu})$

We also study an autocorrelation function based on the accept/reject
statistics of the Metropolis-Hastings chain of length $T$:

$\widehat{\rho(\tau)/\rho(0)}_{acc}=\frac{1}{T-\tau}\sum_{j=1}^{T-\tau}\prod_{i=1}^{\tau}\mathbf{1}_{rej}(i+j)$
(22)

where $\bm{1}_{rej}$ is the identity map of a rejected MH step. Another set of
estimators of the autocorrelation is based on the 2-point observables we
defined above

$\widehat{\rho(\tau)/\rho(0)}_{\mathcal{O}}=\frac{\frac{1}{T-\tau}\sum_{i=1}^{T-\tau}(\mathcal{O}_{i}-\bar{\mathcal{O}})(\mathcal{O}_{i+\tau}-\bar{\mathcal{O}})}{\frac{1}{T}\sum_{i=1}^{T}(\mathcal{O}_{i}-\bar{\mathcal{O}})^{2}}$
(23)

The integrated autocorrelation time wrt acceptance statistics is given by

$\tau_{acc}^{int}=\frac{1}{2}+\lim_{\tau_{\text{max}}\rightarrow\infty}\sum_{\tau=1}^{\tau_{\text{max}}}\frac{\rho_{acc}(\tau)}{\rho_{acc}(0)}\approx\frac{1}{2}+\lim_{\tau_{\text{max}}\rightarrow\infty}\sum_{\tau=1}^{\tau_{\text{max}}}\widehat{\rho(\tau)/\rho(0)}_{acc}$
(24)

and the dynamical critical exponent $z$ for the above quantity is fit using

$\tau_{acc}^{int}(L)\approx cL^{z}$ (25)

We set $T=10^{5}$ and $L\in[8,10,12,14,16,18,20]$. In every Metropolis-
Hastings step given by 19, symmetrization in 20 was done to transform
$\phi^{(i)}$ and $\log q_{\theta}(\phi^{(i)})$ after every ($i$th) step. While
this lead to a slower MCMC evaluation, it significantly improved rates from
between $50$-$55\%$ to $65$-$70\%$. In table 1, we list the values of
$\lambda$ used in [2] for $L\leq 14$ and the remaining fit by trial and error
using Model 1 from 21.

$L$ | 8 | 10 | 12 | 14 | 16 | 18 | 20
---|---|---|---|---|---|---|---
$\lambda$ | 6.008 | 5.550 | 5.276 | 5.113 | 4.99 | 4.89 | 4.82
Table 1: Couplings $\lambda$ corresponding to $L\leq 14$ were used in a
previous study so that the system lies along the $m_{p}L\approx 4$ critical
line. Values for $L\geq 16$ were determined by trial and error evaluations of
$m_{p}L$ for different $\lambda$ using our model. See figure 3 for a plot of
$m_{p}L$ vs $L$ at these values of $\lambda$.

Figure 3 has a plot of $m_{p}L$ vs $L$ at the specified $\lambda$ and the
points are close to the red-dashed critical line at these points.

Figure 3: Product of pole mass and lattice length $m_{p}L$ plotted against $L$
evaluated using Model 1. Error bars indicate 68% confidence intervals obtained
from moving block bootstrap resampling with block size 100

Autocorrelation functions wrt acceptance/rejection statistics (Model 1) and
the observables $G_{c}(0),\chi_{2},E$ are plotted based on 22 and 23 in figure
4.

Figure 4:  Plots of the autocorrelation functions wrt acceptance statistics
and different observables at various lattice lengths evaluated using Model 1

The orange and red curves corresponding to $L=18,20$ seem to decrease slowly
in comparison to others. In the plot integrated autocorrelation times in
figure 5, we evaluate $\tau_{int}^{acc}$ from 24 against $L$ using all of
models 1, 2 and 3. The most striking feature in this plot is that the green
line representing the autocorrelation times of model 3 has values that are
orders of magnitude larger than that of models 1,2 which autoregressively
sample constant time sub-lattices according to our proposal. The values for
the adiabatically-retrained model 2 in the blue line is flatter and has a
smaller dynamical critical exponent compared to that model 1 represented by
the red line.

Figure 5: Integrated autocorrelation times for acceptance statistics plotted
against $L$ in the logarithm scale using Model 1, 2 and 3

The autocorrelation times generally remain roughly constant between $L=8$ and
$L=14$ and increase between $L=16$ and $L=20$, suggesting that the quality (in
terms of closeness to the actual Boltzmann distribution $p$) of $q_{\theta}$
drops when sampling lattice sizes larger than the value $L=16$ it was trained
on. This could also be the reason for growing error bars above $L=16$ of the
susceptibility plot (Model 1) in figure 6 as well as for the pole masses in
figure 3. Still, the dynamical critical exponents obtained by fitting with the
curve in 25 are close to zero, which is usually not the case with traditional
methods555The scaling behavior of flow-based models usually depends on
specific details of the implementation and curve-fitting to an exponential
function can be inaccurate which is evident from the lack of straight-line
behavior in our log-scale plots as well. See the study in [1] for more
details..

Figure 6: Plot of observables $\chi_{2}$ and $E$ against $L$ using Model 1.
The error bars indicate 95% confidence intervals obtained from moving block
bootstrap resampling with block size 100

##  4 Concluding remarks

In this work, we proposed a generative model (dubbed l-ACNF) that
autoregressively samples “equal-time” sub-lattices of a scalar lattice field
using masked conditional normalizing flows. From the locality of action $S$
for the $\phi^{4}$ theory in 2, we determined that these normalizing flows are
conditioned on a smaller sub-lattice of size $O(L^{d-1})$ according to 12
instead of the generic upper bound, $L^{d}$. In other words, locality of the
action reduces the autoregressive expansion of the Boltzmann distribution to a
Markov process in time. Models constructed according to our proposal have
reported orders of magnitude smaller autocorrelation times compared to an
equivalently-constructed existing model used in a previous study in [6]. This
is a consequence of smaller input/output space for the conditional normalizing
flows in our model, essentially addressing the curse of dimensionality in
existing approaches- avoiding exponentially higher number of training epochs
and/or using bulkier neural networks. Constraints in computing resources
limited our study to models that train and converge within hours on relatively
modest hardware. Yet, our models achieve small autocorrelation times and near-
zero dynamic critical exponents similar to large-scale studies like [2]. If
we’re working with the true space-time specification of $d=4$, l-ACNF uses 3D
convolutions for sampling constant time sub-lattices instead of 4D
convolutions. When running on GPU devices, the CUDNN library[5] (used by any
deep learning framework to access GPU-optimized functions/methods) contains
optimized convolution kernels for only upto 3 dimensions due to their relative
practical utility, which offers a slight edge to l-ANF over alternatives.

Since locality is a generic property of all lattice quantum field theories,
our approach can be extended to more complex theories as well. For example,
the Wilson action for lattice gauge fields

$S_{W}=-\sum_{p}\frac{2}{g^{2}}\text{Re}(\text{Tr}(U(p)))$

is a sum over plaquettes $p$ which are loops of links $U_{\mu}(\bm{r})$
through nearest neighbours at every position $\bm{r}=[\bm{x},t]$ and along the
direction $\mu$. We can extend the analysis in section 2.3 and determine that
the conditional distribution of links on the equal-time sub-lattice
$p\left(\\{U_{\mu}(\bm{r})|r_{d}=t\\}|\\{U_{\mu}(\bm{r})|r_{d}<t\\}\right)$
has the dependency set given by (assuming periodic boundary conditions):

$B(\bm{r}|r_{d}=t)=\\{U_{\mu}(\bm{r})|r_{d}=t-1,\mu<d\\}\cup\\{U_{\mu}(\bm{r})|r_{d}=1\\}$

which contains links inside the $r_{d}=t-1$ sub-lattice (except along the time
dimension since they are not part of/connected to plaquettes containing links
along $t=t_{o}$) and the $r_{d}=1$ sub-lattice, which again reduces the input
size of the corresponding flow networks to $O(L^{d-1})$. Combining our
approach with gauge equivariant flows proposed in [7] and [3] is an
interesting direction to investigate from here, but is unfortunately beyond
the scope of this study.

## References

  * [1] Ryan Abbott et al.  “Aspects of scaling and scalability for flow-based sampling of lattice QCD”  In _arXiv preprint arXiv:2211.07541_ , 2022
  * [2] Michael S Albergo, Gurtej Kanwar and Phiala E Shanahan  “Flow-based generative models for Markov chain Monte Carlo in lattice field theory”  In _Physical Review D_ 100.3 APS, 2019, pp. 034515
  * [3] Denis Boyda et al.  “Sampling using SU (N) gauge equivariant flows”  In _Physical Review D_ 103.7 APS, 2021, pp. 074504
  * [4] Michael M Bronstein, Joan Bruna, Taco Cohen and Petar Veličković  “Geometric deep learning: Grids, groups, graphs, geodesics, and gauges”  In _arXiv preprint arXiv:2104.13478_ , 2021
  * [5] Sharan Chetlur et al.  “cudnn: Efficient primitives for deep learning”  In _arXiv preprint arXiv:1410.0759_ , 2014
  * [6] Daniel C Hackett et al.  “Flow-based sampling for multimodal distributions in lattice field theory”  In _arXiv preprint arXiv:2107.00734_ , 2021
  * [7] Gurtej Kanwar et al.  “Equivariant flow-based sampling for lattice gauge theory”  In _Physical Review Letters_ 125.12 APS, 2020, pp. 121601
  * [8] Mario Köppen  “The curse of dimensionality”  In _5th online world conference on soft computing in industrial applications (WSC5)_ 1, 2000, pp. 4–8
  * [9] Colin Morningstar  “The Monte Carlo method in quantum field theory”  In _arXiv preprint hep-lat/0702020_ , 2007
  * [10] Adam Paszke et al.  “Pytorch: An imperative style, high-performance deep learning library”  In _Advances in neural information processing systems_ 32, 2019
  * [11] Tomaso Poggio et al.  “Why and when can deep-but not shallow-networks avoid the curse of dimensionality: a review”  In _International Journal of Automation and Computing_ 14.5 Springer, 2017, pp. 503–519
  * [12] Dinesh PR  “Analysis of Ising model using neural networks”  Indian Institute of Science EducationResearch digital repository, 2021  URL: http://dr.iiserpune.ac.in:8080/xmlui/handle/123456789/6014
  * [13] Robert H Swendsen and Jian-Sheng Wang  “Nonuniversal critical dynamics in Monte Carlo simulations”  In _Physical review letters_ 58.2 APS, 1987, pp. 86
  * [14] Aaron Van den Oord et al.  “Conditional image generation with pixelcnn decoders”  In _Advances in neural information processing systems_ 29, 2016
  * [15] Christina Winkler, Daniel Worrall, Emiel Hoogeboom and Max Welling  “Learning likelihoods with conditional normalizing flows”  In _arXiv preprint arXiv:1912.00042_ , 2019
  * [16] Ulli Wolff  “Collective Monte Carlo updating for spin systems”  In _Physical Review Letters_ 62.4 APS, 1989, pp. 361
  * [17] Ulli Wolff  “Critical slowing down”  In _Nuclear Physics B-Proceedings Supplements_ 17 Elsevier, 1990, pp. 93–102
  * [18] Dian Wu, Riccardo Rossi and Giuseppe Carleo  “Unbiased Monte Carlo cluster updates with autoregressive neural networks”  In _Physical Review Research_ 3.4 APS, 2021, pp. L042024
  * [19] Dian Wu, Lei Wang and Pan Zhang  “Solving statistical mechanics using variational autoregressive networks”  In _Physical review letters_ 122.8 APS, 2019, pp. 080602