# Faster estimation of the Knorr-Held Type IV space-time model

Fredrik Lohne Aanes
Norwegian Computing Center Email<EMAIL_ADDRESS>Geir Storvik
University of Oslo and Norwegian Computing Center Email<EMAIL_ADDRESS>corresponding author

###### Abstract

In this paper we study the type IV Knorr Held space time models. Such models
typically apply intrinsic Markov random fields and constraints are imposed for
identifiability. INLA is an efficient inference tool for such models where
constraints are dealt with through a conditioning by kriging approach. When
the number of spatial and/or temporal time points become large, it becomes
computationally expensive to fit such models, partly due to the number of
constraints involved. We propose a new approach, HyMiK, dividing constraints
into two separate sets where one part is treated through a mixed effect
approach while the other one is approached by the standard conditioning by
kriging method, resulting in a more efficient procedure for dealing with
constraints. The new approach is easy to apply based on existing
implementations of INLA.

We run the model on simulated data, on a real data set containing dengue fever
cases in Brazil and another real data set of confirmed positive test cases of
Covid-19 in the counties of Norway. For all cases we get very similar results
when comparing the new approach with the tradition one while at the same time
obtaining a significant increase in computational speed, varying on a factor
from 2 to 4, depending on the sizes of the data sets.

Keywords: Big data, INLA, linear constraints, Space-time interaction

##  1 Introduction

The availability of temporal areal data has greatly improved in recent years
and Bayesian analysis have gained increasing interest for such data. Temporal
and spatial dependence which are not accounted for by available explanatory
variables are typically dealt with through random effects with specific
structures. Knorr-Held (2000) introduced the Knorr-Held separable space-time
models for intrinsic Gaussian Markov random fields. These (and similar) models
have become highly popular within Bayesian disease mapping (Lawson, 2018),
analysis of mortality rates (Khana et al., 2018), analysis of environmental
systems (Wikle, 2003), election forecasts (Pavía et al., 2008) and social
sciences (Shoesmith, 2013; Vicente et al., 2020; Williams et al., 2019) as
well as in many other applications.

For the (Bayesian) analysis of spatio-temporal areal data, there is a huge
computational challenge. This is particular the case for the Knorr-Held
models. Methods and generic packages such as CARBayesST (Lee et al., 2018) for
Markov chain Monte Carlo (MCMC) are available, but can be inherently slow for
large data sets. Integrated nested Laplace approximation (INLA, Rue et al.,
2009) is very often used to fit such Bayesian models. In Nazia et al. (2022),
a review over different techniques used to model Covid-19, 14 out of the 18
Bayesian studies used the INLA approach. INLA uses a combination of analytical
approximations and numerical algorithms to approximate the marginal posterior
distributions. One can then avoid using the traditional MCMC sampling method
to obtain inference. However, handling a large number constraints, typically
needed for the Knorr-Held models, is a challenge also within the INLA
framework.

There are four different Knorr-Held models, where the complexity of the model
depends on the structure for the interaction term. We focus on the type IV-
model. It is the most complex model as it assumes the precision matrix of an
interaction between space and time term is both spatially and temporally
structured. An additional complexity is due to that the random effects
typically are modeled through intrinsic Gaussian Markov random fields which
are improper Gaussian models (singular precision matrices). Typically this
also leads to non-identifiable models. Additional identifiability constraints
are needed to make the model identifiable.

To make the model estimable, Schrödle and Held (2011) suggest to make the
Knorr-Held term proper by imposing several constraints on the prior. As an
alternative, Goicoa et al. (2018), suggest to make the posterior proper, only
including enough constraints so that the model is identifiable. Several
constraints are needed when using either sets of constraints, although the
latter will have somewhat fewer constraints. The high number of constraints
needed make the models hard to compute for large studies (where either
$n_{T}$, the number of time steps, and/or $n_{S}$, the number of spatial
regions, is very large). The latest version of INLA has been made much more
efficient by better use of parallelism and state-of-the-art sparse linear
solvers (Gaedke-Merzhäuser et al., 2023; Van Niekerk et al., 2023), but still
require substantial computational resources.

Incorporation of constraints can be performed in different ways. Conditioning
by kriging (Cressie, 1993) is a post-processing procedure which ”corrects”
unconstrained random processes into constrained ones through a linear
transformation. This is done in each iteration of the estimation procedure.
The computational cost of constraints when using this approach is
$\mathcal{O}(nk^{2})$ where $k$ the number of constraints and $n$ is the
dimension of the field (see e.g Bakka et al., 2018). Conditioning by kriging
is the default procedure within INLA.

An alternative approach is to formulate the model as a mixed effect model with
the inclusion of an appropriate design matrix that accounts for the
constraints. No post-processing is then needed, but the resulting precision
matrices are more dense. This approach can be implemented through the use of
the ”z-model” in INLA. In this paper we show that by _combining_ conditioning
by kriging with the mixed effect model in a suitable way, significant
computational efficiency can be gained for the Knorr-Held model type IV term.
A subset of the constraints will be dealt with through the mixed effect model,
while the remaining ones use the conditioning by kriging step approach. By
suitable constructions of such subsets, the increasing computational effort
due to the mixed effect model will be much less than the decreasing
computational effort due to the number of constraints that needs to be dealt
with through the conditioning by kriging step. This holds irrespective of the
computational resources.

Recently Fattah and Rue (2022) proposed a new way of estimating constrained
intrinsic Gaussian Markov random fields. They utilize the correspondence
between the constraints suggested by Schrödle and Held (2011) and the Moore-
Penrose inverse of the precision matrix. The procedure is still
computationally demanding, but 40 fold faster inference can be obtained
paralleling the computational efforts on a large server. They run it on 25
nodes Cascade Lake nodes, each with 40 cores, 2.50 GHz, 384 GB/usable 350 GB.
We do not have the same computational resources available, and cannot run the
models equally fast. We therefore do not try to compare our approach with the
Fattah and Rue (2022) approach. However, as we will see, we can obtain similar
speed-ups with our method without the need for numerous cores.

The paper is structured as follows. In Section 2 we formally introduce space-
time models and discuss how constrained fields are estimated in INLA. In
Section 3 we present the new method of estimating constrained fields. In
Section 4 we study three datasets, one simulated and two real data sets. We
also discuss the results. Finally, Section 5, provide a summary and a
discussion.

##  2 Model formulation

We consider a modeling framework where we have data
$\\{y_{t,s},t=1,...,n_{T},s=1,...,n_{S}\\}$ corresponding to observations at
time $t$ within area/region $s$. The observations are assumed conditionally
independent with

$\displaystyle y_{t,s}|\eta_{t,s}\sim$ $\displaystyle
p_{y}(y_{t,s}|\eta_{t,s},\phi)$ (1a) $\displaystyle\eta_{t,s}=$
$\displaystyle\mu+\bm{u}_{t,s}^{T}\bm{\beta}+\alpha_{t}+\gamma_{s}+\delta_{t,s}=(\mathcal{A}\bm{x})_{t,s}$
(1b)

for $t=1,...,{n_{T}},s=1,...,{n_{S}}$. Here, $p_{y}(\cdot)$ is a distribution
within the exponential family with $E[y_{t,s}]=g^{-1}(\eta_{t,s})$ for some
appropriate link function $g(\cdot)$ while $\phi$ is a dispersion parameter.
The $\mathcal{A}$-matrix is a sparse design matrix, while
$\bm{x}=\\{\mu,\bm{\beta},\bm{\alpha},\bm{\gamma},\bm{\delta}\\}$ is the
latent field. The intercept $\mu$ and the regression coefficients $\bm{\beta}$
are treated as fixed effects, the latter related to covariates
$\bm{u}_{t,s}$111We have used the somewhat unusual notation $\bm{u}$ for the
covariates here in order to follow the standard notation of $\bm{x}$ for
latent variables later on.. The two main random effects $\\{\alpha_{t}\\}$
(temporal component) and $\\{\gamma_{s}\\}$ (spatial component) are both
assumed to be zero mean Gaussian Markov random fields with precision matrices
$\bm{Q}_{\alpha}$ and $\bm{Q}_{\gamma}$, respectively. Also the interaction
term $\\{\delta_{t,s}\\}$ is assumed to following a zero mean Gaussian
distribution, with precision matrix $\bm{Q}_{\delta}$. Typically, the
precision matrices for the random effects are singular, corresponding to
_intrinsic_ Gaussian Markov random field models. In such settings, the model
can become non-identifiable and constraints are in such cases required on the
random effects.

For the temporal component, typical model choices are the random walk of order
1 (RW(1)), in which case $\alpha_{t}-\alpha_{n_{T}-1}\sim
N(0,\tau_{\alpha}^{-1})$ or the random walk of order 2 (RW(2)), in which case
$(\alpha_{t}-\alpha_{n_{T}-1})-(\alpha_{n_{T}-1}-\alpha_{t-2})\sim
N(0,\tau_{\alpha}^{-1})$. For the RW(1) model, a sum to zero constraint
$\sum_{t=1}^{{n_{T}}}\alpha_{t}=0$ will make the model proper. For the RW(2)
model, the additional constraint $\sum_{t=1}^{{n_{T}}}t\alpha_{t}=0$ results
in a proper model. Note however that depending on which fixed effects
included, the $\\{\alpha_{t}\\}$ process might be _identifiable_ without the
latter constraint. If a linear temporal effect is included as a fixed effect,
the additional constraint needs to be included.

The _intrinsic conditional autoregressive model_ (ICAR) was introduced by
Besag (1974) as a two-dimensional generalization of the discrete penalty
smoother by Whittaker (1922) and has been a popular choice for the spatial
component. Here we assume

$\displaystyle\gamma_{s}|\bm{\gamma}_{-s}\sim
N\left(\gamma_{s};{|\mathcal{N}_{s}|}^{-1}\sum_{s^{\prime}\in\mathcal{N}_{s}}\gamma_{s^{\prime}},(\tau_{\gamma}|\mathcal{N}_{s}|)^{-1}\right),$
(2)

where $\mathcal{N}_{s}$ is the set of nodes that are neighbors to node $s$ and
$|\mathcal{N}_{s}|$ is the number of neighbors for node $s$. In this case,
$\bm{Q}_{\gamma}$ will have $\tau_{\gamma}|\mathcal{N}_{s}|$ on the diagonals,
$-\tau_{\gamma}$ for all off-diagonals corresponding to neighbor nodes and
zero otherwise. Assuming all nodes are connected through a path of neighbor
nodes, $\bm{Q}_{\gamma}$ has rank deficiency 1 and a typical constraint
imposed is $\sum_{s=1}^{n_{S}}\gamma_{s}=0$.

Note that for the models discussed above, we have
$\bm{Q}_{\alpha}=\tau_{\alpha}\bm{R}_{\alpha}$ and
$\bm{Q}_{\gamma}=\tau_{\gamma}\bm{R}_{\gamma}$ where $\bm{R}_{\alpha}$ and
$\bm{R}_{\gamma}$ are known (singular) matrices, typically denoted as the
structural matrices. In more general settings, some unknown parameters might
be involved in the structural matrices as well, although we will not consider
such cases in this paper.

###  2.1 The space time interaction and corresponding constraints

We will in this section discuss the constraints needed for the spatio-temporal
interaction terms $\\{\delta_{t,s}\\}$. Knorr-Held (2000) considered four
types of interaction models corresponding to unstructured/structured precision
matrices for the temporal/spatial components of $\bm{Q}_{\delta}$ where
unstructured means independence while structured corresponds to the RW(1/2)
and to the ICAR models discussed in the previous section. We will here focus
on Type IV interactions where both the temporal and the spatial parts are
structured with $\bm{Q}_{\delta}=\tau_{\delta}\bm{R}_{\delta}$ and

$\bm{R}_{\delta}=\bm{R}_{\alpha}\otimes\bm{R}_{\gamma}.$ (3)

Denote the rank deficiencies for $\bm{R}_{\alpha}$ and $\bm{R}_{\gamma}$ by
$k_{\alpha}$ and $k_{\gamma}$, respectively. Then the rank deficiency for
$\bm{R}_{\delta}$ is
$k_{\delta}=k_{\alpha}*{n_{S}}+k_{\gamma}*{n_{T}}-k_{\alpha}*k_{\gamma}$. In
order to make the distribution for
$\bm{\delta}=(\delta_{1,1},...,\delta_{1,{n_{S}}},\delta_{2,1},...,\delta_{{n_{T}},{n_{S}}})$
a proper distribution, $k_{\delta}$ constraints are needed. However, in order
to make the model _identifiable_ , fewer constraints might be needed. For
instance, if the model is

$\eta_{t,s}=\delta_{t,s}$

so no fixed effects and no other random effects, all the $\delta_{t,s}$’s will
be identifiable (given that there is one observation for each $(t,s)$
combination). In practice, at least an intercept, a temporal $\alpha_{t}$ as
well as a spatial $\gamma_{s}$ effect are typically included for ease of
interpretation. In these situations, constraints are needed.

Several sets of constraints have been suggested in the literature. In the
following, we will always assume that the intercept $\mu$ is included as a
fixed effect, in which case the sum constraints

$\bm{1}_{n_{S}}^{T}\bm{\gamma}=0,\quad\bm{1}_{n_{T}}^{T}\bm{\alpha}=0$

are included in order to avoid identifiability issues between the intercept
and the main effects. Here, $\bm{1}_{n}$ is a vector of length $n$ with all
components equal to 1.

For the interaction term, more care is needed with respect to design of
constraints, due to that it can both have identifiability issues related to
the fixed effects and to the main random effects. The constraints suggested by
Goicoa et al. (2018) for the interaction term, named the GC-constraints in the
following, are given in the first row of Table 1. These constraints are based
on that main random effects always are included and focus on constraining
effects that are non-identifiable regardless of data available. Due to missing
data (e.g no observations for a specific region) some other effects might also
be non-indentifiable. We will however only consider cases where the remaining
effects _are_ identifiable, typically based on that there is at least one
observation for each combination $(t,s)$. We note that the GC-constraint
system is overdetermined as the rank of the constraint matrix is less than the
number of constraints. In practice, a set of constraints of full rank can be
obtained by removing one row in either
$\bm{I}_{n_{T}}\otimes\bm{1}_{n_{S}}^{T}$ or
$\bm{1}_{n_{T}}^{T}\otimes\bm{I}_{n_{S}}$.

Schrödle and Held (2011) suggest to use constraints that give a prior that is
_proper_ , a stronger requirement than the identifiability constraint. For
some settings a proper prior might be desirable, although, as noted by Goicoa
et al. (2018), the extra constraints in the SC sets are not needed from an
identifiability point of view unless a common linear trend and area specific
linear trends are included in the model as fixed effects. The SC constraints
are given in the second row of Table 1. Also this set of constraints are
overdetermined and deleting one row from either
$\bm{I}_{n_{T}}\otimes\bm{1}_{n_{S}}^{T}$ or from
$\bm{1}_{n_{T}}^{T}\otimes\bm{I}_{n_{S}}$ in addition to removing one row from
$\bm{d}^{T}\otimes\bm{I}_{n_{S}}$ will give a full rank set. Assuming an RW(2)
model for the temporal structure and proper priors also are wanted for the
main effects, the constraint $\bm{d}^{T}\bm{\alpha}=0$ should additionally be
imposed.

Depending on what fixed effects are included, additional constraints might be
needed. For instance, if a linear effect in time is part of the fixed effects,
the constraint $\bm{d}^{T}\bm{\alpha}=0$ should always be included. Note that
for both constraint sets, a sum to zero condition is implied, although not
explicitly listed.

The SC constraints in Table 1 can be seen as defined within the space spanned
by the eigenvectors corresponding to the eigenvalues equal to zero for
$\bm{Q}_{\delta}$. Schrödle and Held (2011) more generally suggest defining
the constraint matrix $\bm{A}$ by those eigenvectors of $\bm{R}_{\delta}$
which span the null space. Fattah and Rue (2022) considered this approach,
using the Moore-Penrose pseudo-inverse of the prior precision matrix $\bm{Q}$
to obtain such constraints. Typically, however, the resulting _covariance
matrix_ will be dense resulting in high computational complexity. Fattah and
Rue (2022) took advantage of multi-core architectures and abundance of memory
to deal with the computational burden, though still resulting in a
computational demanding procedure.

Constraint name | Constraint
---|---
Goicoa (GC) | $[\bm{I}_{n_{T}}\otimes\bm{1}_{n_{S}}^{T}]\bm{\delta}=\bm{0}_{n_{T}},[\bm{1}_{n_{T}}^{T}\otimes\bm{I}_{n_{S}}]\bm{\delta}=\bm{0}_{n_{S}},$
Schrödle (SC) |  $[\bm{I}_{n_{T}}\otimes\bm{1}_{n_{S}}^{T}]\bm{\delta}=\bm{0}_{n_{T}},[\bm{1}_{n_{T}}^{T}\otimes\bm{I}_{n_{S}}]\bm{\delta}=\bm{0}_{n_{S}},$ $[\bm{d}^{T}\otimes\bm{I}_{n_{S}}]\bm{\delta}=\bm{0}_{n_{S}}$
Table 1: Constraints on the interaction terms in the space-time model with
smooth main spatial effect $\bm{\gamma}$ smooth main temporal effect
$\bm{\alpha}$ and smooth interaction effect $\bm{\delta}$. Here
$\bm{d}=(1,2,...,n_{T})^{T}$.

In general, precision matrices involved will be sparse. This can be exploited
by estimation algorithms to obtain faster inference. However, including
constraints will affect the sparseness of the precision matrix for $\bm{x}$
and algorithms have to correct for the constraints ways such that fast
inference is still achievable. We will consider a Bayesian approach where we
assume there exist priors for the parameters involved. The form of these
priors will be specified within the numerical studies in Section 4, but a main
restriction is that the parameters related to the fixed effects follow a
(possible noninformative) Gaussian distribution. The constraining approach to
be presented in Section 3 should in principle work for any type of prior for
the remaining parameters (mainly scaling factors for the precision matrices
involved).

###  2.2 Integrated nested Laplace approximation

Integrated nested Laplace approximation (INLA, Rue et al., 2009) has been a
popular method for performing inference within latent Gaussian models such as
model (1). We provide here a rough overview of the method. Let $\bm{\theta}$
be the hyperparameters involved (typically precision and dispersion
parameters) and $\bm{x}$ be the set of latent variables, including both random
effects, the intercept and regression coefficients. INLA starts by making an
approximation

$\tilde{p}(\bm{\theta}|\bm{y})\propto\frac{p(\bm{x},\bm{\theta},\bm{y})}{\tilde{p}_{G}(\bm{x}|\bm{\theta},\bm{y})}\bigg{|}_{\bm{x}=\bm{x}^{*}(\bm{\theta})}$
(4)

where $\tilde{p}_{G}(\bm{x}|\bm{\theta},\bm{y})$ is a Gaussian approximation
to $p(\bm{x}|\bm{\theta},\bm{y})$ and $\bm{x}^{*}(\bm{\theta})$ is the mode of
$p(\bm{x}|\bm{\theta},\bm{y})$ for a given $\bm{\theta}$. Note that the
Gaussian prior on $\bm{x}$ make such an approximation reasonable. The marginal
distributions $p(x_{i}|\bm{y})$ are then approximated by

$\tilde{p}(x_{i}|\bm{y})=\sum_{k}\tilde{p}(x_{i}|\bm{\theta}_{k},\bm{y})\tilde{p}(\bm{\theta}_{k}|\bm{y})\Delta_{k}$

where $\tilde{p}(x_{i}|\bm{\theta}_{k},\bm{y})$ is an approximation of
$p(x_{i}|\bm{\theta}_{k},\bm{y})$ (Rue et al., 2009, considered three
alternatives with a Gaussian approximation being the simplest one). The
$\bm{\theta}_{k}$ with their associated weights $\Delta_{k}$ are carefully
distributed around the mode of $\tilde{p}(\bm{\theta}|\bm{y})$, see Rue et al.
(2009) for details.

The main computational cost is related to the Gaussian approximation
$\tilde{p}_{G}(\bm{x}|\bm{\theta},\bm{y})$. Due to the Gaussian assumption on
$\bm{x}$ and the conditional independence assumption on the observations, we
have

$\displaystyle
p(\bm{x}|\bm{y},\bm{\theta})\propto\exp\left(-\frac{1}{2}\bm{x}^{T}\bm{Q}\bm{x}+\sum_{t,s}g_{t,s}(\mathcal{A}\bm{x})\right),$
(5)

where $\mathcal{A}\bm{x}$ is specified through (1b). Further, $\bm{Q}$ is the
full precision matrix for all the elements in $\bm{x}$. The sparsity
structures of the precision matrices for the different random and fixed
effects in $\bm{x}$ will be inherited into $\bm{Q}$. The Gaussian
approximations involved in INLA is based on Taylor expansions of the
$g_{t,s}(\bm{(}\mathcal{A}\bm{x})_{t,s})$ terms. The Gaussian approximation is
in particular obtained by the Taylor approximation

$g_{t,s}((\mathcal{A}\bm{x}^{m})_{t,s})\approx
g_{t,s}(\mathcal{A}\bm{x}^{m})_{t,s}+b_{t,s}(\mathcal{A}\bm{x}^{m})_{t,s}-\tfrac{1}{2}c_{t,s}(\mathcal{A}\bm{x}^{m})_{t,s}^{2}$
(6)

based on the mode $\bm{x}^{m}$ of $p(\bm{x}|\bm{\theta},\bm{y})$. The mode is
computed iteratively using a Newton-Raphson algorithm. The computational steps
utilize the Cholesky decomposition of $\bm{Q}$ which is sparse if $\bm{Q}$ is
sparse. See Van Niekerk et al. (2023) for further details on the
implementation and computational cost related to INLA.

###  2.3 INLA and constraints

As discussed in Section 2.1, (linear) constraints are usually needed in order
to make intrinsic models identifiable. There are different ways in which such
constraints can be incorporated into the inferential procedure. In case of
constraints $\bm{A}\bm{x}=\bm{0}$, the density to be considered is
$p(\bm{x}|\bm{A}\bm{x}=\bm{0},\bm{y},\bm{\theta})$.

The most standard way of incorporating constraints within INLA is based on the
conditioning by kriging technique (Rue and Held, 2005, sec 2.3.3) using the
posterior distribution. First, denote by $\widetilde{\bm{Q}}$ the precision
matrix corresponding to the Gaussian approximation (5). Assume $\bm{x}^{*}\sim
N(\bm{\nu},\widetilde{\bm{Q}}_{\epsilon}^{-1})$, where

$\widetilde{\bm{Q}}_{\varepsilon}=\widetilde{\bm{Q}}+\bm{I}\varepsilon$ (7)

where $\varepsilon>0$ is a small number included for numerical reasons. Then

$\bm{x}=\bm{x}^{*}-\widetilde{\bm{Q}}_{\epsilon}^{-1}\bm{A}^{T}(\bm{A}\widetilde{\bm{Q}}_{\epsilon}^{-1}\bm{A}^{T})^{-1}\bm{A}\bm{x}^{*}=[\bm{I}-\widetilde{\bm{Q}}_{\epsilon}^{-1}\bm{A}^{T}(\bm{A}\widetilde{\bm{Q}}_{\epsilon}^{-1}\bm{A}^{T})^{-1}\bm{A}]\bm{x}^{*}$
(8)

will have the right conditional distribution (Rue and Held, 2005, sec 2.3.3).
This relationship is utilized in the optimization step by correcting the
solution at each iteration based on the relation using (8) based on

$p(\bm{x}|\bm{A}\bm{x}=\bm{0},\bm{y})=\frac{p(\bm{A}\bm{x}|\bm{x},\bm{y})p(\bm{x}|\bm{y})}{p(\bm{A}\bm{x}|\bm{y})}I(\bm{A}\bm{x}=\bm{0})$
(9)

which is evaluated at the $\bm{x}$ derived by (8). Due to the additional
$\varepsilon$ term, the conditioning by kriging correction will not give exact
results, although the effect should be small.

This correction step is done in each iteration of the INLA procedure. The
additional cost of this correction is $\mathcal{O}(nk^{2})$ where $k$ is the
number of constraints. For small $k$, conditioning by kriging is very
efficient. However, for spatio-temporal random effects $k$ is quite large.

When intrinsic models are used, it is common to add a small value
$\varepsilon$ to the precision matrices involved in order to make the models
proper. If the distribution is proper, adding the value $\varepsilon$ to the
diagonal should only give marginal modifications of the distribution. The
eigendecomposition of
$\bm{Q}_{\varepsilon}^{-1}=(\bm{Q}+\varepsilon\bm{I})^{-1}$ is given by

$\bm{Q}^{-1}_{\varepsilon}=\sum_{i=1}^{n}\frac{1}{\varepsilon+\lambda_{i}}\bm{a}\bm{a}_{i}^{T}$
(10)

where $\\{\bm{a}_{i}\\}$ are the eigenvectors of $\bm{Q}$ with corresponding
eigenvalues $\\{\lambda_{i}\\}$. Note that for eigenvalues $\lambda_{i}=0$,
the additional $\varepsilon$ terms remove the problems of calculating inverses
in (8) but can be numerically unstable due to the resulting
$\tfrac{1}{\varepsilon}$ terms in (10). We have a similar situation for the
posterior precision $\tilde{\bm{Q}}$.

##  3 The hybrid mixed effect and kriging method for constraining intrinsic
Gaussian models

In this section we propose a new way of dealing with constraints for the
spatio-temporal interaction terms $\bm{\delta}$. We will first demonstrate how
constraints can be dealt with through the mixed effect model approach. This
procedure will in itself not be very efficient, but we will then propose a new
method where we combine the mixed effect approach with the conditioning by
kriging approach in order to obtain a computationally efficient procedure.

###  3.1 The mixed effect constraining method for Gaussian models

The trick in the mixed effect model approach is to reformulate the model by
first assuming $\bm{x}^{*}$ has precision matrix $\bm{Q}_{\varepsilon}$ and
then define

$\bm{\eta}=\bm{U}\bm{\beta}+\bm{Z}\bm{x}^{*},$ (11)

for an appropriate design matrix $\bm{Z}$ such that the constraints are
directly incorporated into $\bm{Z}\bm{x}^{*}$. Assume that the constraints are
given by

$\bm{a}_{i}^{T}\bm{x}=0,i=n-k+1,...,n$

where now the $\bm{a}_{i}$’s are orthonormal eigenvectors of the prior
precision matrix $\bm{Q}$, corresponding to eigenvalues equal to zero. Then
the condition by kriging transformation (8) reduces to

$\bm{x}=\bm{x}^{*}-\sum_{i=n-k+1}^{n}\bm{a}_{i}\bm{a}_{i}^{T}\bm{x}^{*}=[\bm{I}-\sum_{i=n-k+1}^{n}\bm{a}_{i}\bm{a}_{i}^{T}]\bm{x}^{*}$

so by choosing $\bm{Z}=\bm{I}-\sum_{i=n-k+1}^{n}\bm{a}_{i}\bm{a}_{i}^{T}$ we
directly obtain the desired constraints. Note that $\bm{Z}^{2}=\bm{Z}$ showing
that $\bm{Z}$ is a projection matrix.

In INLA, the use of mixed effect models is implemented through the ”z-model”
(R-INLA, 2020), but has, to our knowledge, not been used for incorporating
constraints. The method is implemented by working with the full precision
matrix for both $\bm{x}^{*}$ and $\bm{x}=\bm{Z}\bm{x}^{*}$. Even though
$\bm{x}$ is a deterministic function of $\bm{x}^{*}$, for numerical reasons,
it is assumed that $\bm{x}|\bm{x}^{*}\sim
N(\bm{Z}\bm{x}^{*},\kappa^{-1}\bm{I})$ for a large precision value $\kappa$ in
addition to that $\bm{Q}$ is replaced by $\bm{Q}_{\varepsilon}$. This results
in the joint distribution for $(\bm{x},\bm{x}^{*})$ being a zero-mean GMRF
with joint precision matrix

$\bm{Q}_{J}=\begin{bmatrix}\kappa\bm{I}&-\kappa\bm{Z}\\\
-\kappa\bm{Z}&\kappa\bm{Z}+\bm{Q}_{\varepsilon}\end{bmatrix}.$ (12)

The model (11) is then related to this extended random vector through
$\bm{\eta}=\bm{U}\bm{\beta}+\bm{x}$. In contrast to the kriging approach
discussed in Section 2.3, the constraints are now corrected for in the prior
distribution so that it follows a constrained distribution _before_ performing
the Laplace approximation.

###  3.2 The hybrid mixed effect and kriging method

In general, the $\bm{Z}$ matrix can be quite dense, making the use of the
mixed effect model approach computationally inefficient. Our suggested
approach is based on _combining_ the mixed effect approach with the
conditioning by kriging approach by dividing the set of constraints into two
parts. A first set of constraints, $\bm{A}_{1}\bm{x}=\bm{0}$ with the $k_{1}$
rows being orthonormal are dealt with through the mixed effect model approach.
The remaining constraints $\bm{A}_{2}\bm{x}=\bm{0}$ are dealt with through the
conditioning by kriging procedure. For $\bm{A}_{2}$, the rows do not need to
be orthonormal. Note that through the mixed effect formulation (11), the
remaining constraints are put on $\bm{x}=\bm{Z}\bm{x}^{*}$ and not on
$\bm{x}^{*}$. This procedure, which we name _HYbrid MIxed effect and Kriging
method for constraining intrinsic Gaussian models_ (HyMiK), results in fewer
constraints needed to be corrected for through the conditioning by kriging
equation (8) at the cost of a larger precision matrix $\bm{Q}_{J}$. However,
if $\bm{Q}_{J}$ is sparse, the increased cost of the larger latent field can
be minor compared to the reduction in computing time of the conditioning by
kriging step.

In addition to computational gain, we also anticipate that the new method will
give more correct parameter estimates when using non-Gaussian observation
models. The reason is that $\bm{x}$ has been corrected for several constraints
before performing the approximation (4). Due to the intrinsic structures of
the precision matrices involved, some directions in the $\bm{x}$ space will
have large variances, giving sources for numerical instability. When some
constraints are accounted for prior to this approximation, less directions
will also have this problem with larger variances. Also, less constraints are
needed to be corrected for through the conditioning by kriging approach. Due
to that the Gaussian approximation within the HyMiK approach is performed on a
much more restricted space, we expect this approximation to be more accurate.
Since the approximation is more correct, the kriging step should also be more
accurate. For some models, the kriging step of HyMiK is not needed so that no
post-processing is performed. Therefore HyMiK is better than using the
traditional conditioning by kriging method.

The split of the constraints should be performed such that both the joint
precision matrix $\bm{Q}_{J}$ in (12) be reasonable sparse and the remaining
constraints used in the conditioning by kriging procedure is low in number.
From the GC constraints listed in Table 1, there are two obvious choices.
Either $\bm{A}_{1}=\bm{A}_{S,1}\equiv\bm{e}_{n_{T}}^{T}\otimes\bm{I}_{n_{S}}$
or $\bm{A}_{1}=\bm{A}_{T,1}\equiv\bm{I}_{n_{T}}\otimes\bm{e}_{n_{S}}^{T}$ can
be used. Here $\bm{e}_{n}=\bm{1}_{n}/\sqrt{n}$ is the normalized version of
$\bm{1}_{n}$. Which constraint matrix to choose will depend on the size of
$n_{T}$ compared to $n_{S}$. If $n_{S}>n_{T}$, the first choice is the best
one while the second one is better when $n_{T}>n_{S}$.

Considering the choice $\bm{A}_{1}=\bm{A}_{S,1}$, there are $n_{S}n_{T}^{2}$
non-zero elements in $\bm{Z}$, which means that a fraction $1/n_{S}$ of the
elements are non-zero. If $n_{T}=10,n_{S}=544$ (as in the simulation example
to be considered later), then 0.1838 % of the elements are non-zero. If we
however consider using the temporal constraints $\bm{A}_{1}=\bm{A}_{T,1}$,
then $1/n_{T}=10\%$ of the elements of the projector matrix $\bm{Z}$ is non-
zero. We observe a striking difference. This affects the run time.

For the SC-constraints, we have the additional constraints
$[\bm{d}\otimes\bm{I}_{n_{S}}]\bm{\delta}=\bm{0}$. These can be incorporated
to the projection part as well by making these orthogonal to the constraints
already included, as well as scaling them to have norm 1. Assuming
$\bm{A}_{1}=\bm{A}_{1,S}$, this can be achieved by standard scaling and
centering of $\bm{d}$, that is

$\tilde{d}_{t}=\frac{d_{t}-\bar{d}}{\sqrt{n_{T}^{-1}\sum_{i=1}^{n_{T}}(d_{i}-\bar{d})^{2}}}.$

We then can use $\bm{A}_{1}$ as the combination of $\bm{A}_{S,1}$ and
$\bm{[}\tilde{\bm{d}}\otimes\bm{I}_{n_{S}}]$.

###  3.3 Does the mixed effect and kriging methods give the same posterior?

The mixed effect approach is used to the correct the prior while the kriging
method is used to correct the distribution during the execution of the
algorithm. We now show that they provide the same inference. Assuming
$\bm{A}\bm{x}=\bm{0}$, we have

$\displaystyle p(\bm{x}|\bm{A}\bm{x}=\bm{0},\bm{y})=$
$\displaystyle\frac{p(\bm{x})p(\bm{A}x=\bm{0}|\bm{x})p(\bm{y}|\bm{x},\bm{A}\bm{x}=\bm{0})}{p(\bm{y}|\bm{A}\bm{x}=\bm{0})p(\bm{A}\bm{x}=\bm{0})}$
$\displaystyle=$
$\displaystyle\frac{p(\bm{x})p(\bm{y}|\bm{x})}{p(\bm{y}|\bm{A}\bm{x}=\bm{0})p(\bm{A}\bm{x}=\bm{0})}$
$\displaystyle=$
$\displaystyle\frac{p(\bm{x}|\bm{A}\bm{x}=\bm{0})p(\bm{y}|\bm{x})}{p(\bm{y})},$

where we in the second equality have used that $p(\bm{A}x=\bm{0}|\bm{x})=1$
when $\bm{A}\bm{x}=\bm{0}$ and that given $\bm{x}$, the extra conditioning
$\bm{A}\bm{x}=\bm{0}$ does not provide extra information. In the last equality
we have used that $p(\bm{y}|\bm{A}\bm{x}=\bm{0})$ reduce to $p(\bm{y})$
assuming that the constraints imposed only are related to identifiability so
that the density of $\bm{y}$ do not depend on these constraints. From this we
see that conditioning in the posterior is equivalent to conditioning in the
prior.

##  4 Numerical studies

In this section we perform numerical studies in order to compare the HyMiK
method with the standard conditioning by kriging method. Three data sets are
considered. The first is a simulated data set where we use a neighbor
structure based on districts within Germany. The second example is a real data
set of dengue fewer cases in Brazil, previously studied in Lowe et al. (2021).
The last example is a real data set consisting of cases of positive tests of
COVID-19 in Norway.

The computational gain is measured as the fraction between computing time
using the standard conditioning by kriging method and the time for the HyMiK
method. We call this the _computational factor_.

All computations were performed on a linux server (kernel 5.4.0-1100.azure)
using INLA version 22.05.07 with the Pardiso-extension (Van Niekerk et al.,
2021). The Pardiso extension uses a modern sparse-matrix library, which gives
faster inference for complex models. A maximum number of 10 threads were used
in all calls to INLA.

All datasets and scripts are available at the github repository
https://github.com/geirstorvik/INLAconstraints. When calculating marginal
(log)likelihoods, care has to be taken due to that for some models the
normalizing constant, i.e. the term $0.5\log(|\bm{R}|)$ (with $\bm{R}$ being
the structural matrix taking constraints into account), is missing.
Calculation of this term was performed as a post-processing step. Further, the
interpretation of the marginal likelihoods for the GC constraints is somewhat
difficult due to that the priors are still improper under the GC constraints.
The prior is first made proper by considering $\bm{Q}_{\varepsilon}$, see (7).
The constraints only make the posterior proper and not the prior, which have
some directions with very large variances. Technically, the variance in some
directions is still $1/\varepsilon$ under the constraints, which gives what we
call an improper prior. For the SC constraints, the priors are proper under
the constraints (where the additional constraint $\sum t\alpha_{t}$ is needed
on the temporal term). We still consider $\bm{Q}_{\varepsilon}$ but the
constraints make the variance in the direction which originally had variance
$1/\varepsilon$ now 0. This distribution is singular, but proper.

###  4.1 Simulated data

We simulated data from model (1) with $p_{y}$ corresponding to the Poisson
model with log-link. The spatial structure was based on the neighbor structure
of the $n_{S}=544$ districts of Germany. This neighbor structure is available
in the INLA package through the Germany data, see e.g.
https://rdrr.io/github/inbo/INLA/src/demo/Graph.R. This neighbor structure was
used to define the precision matrix for the main spatial random effect through
the ICAR model (2). For the temporal effect, we used the RW(2) model and the
interaction term followed the type IV model of Knorr-Held (2000), that is
equation (3). When simulating data, the SC constraints were used to obtain a
proper prior of the interaction effect and also appropriate constraints where
used to obtain proper main spatial and temporal effects. A summary of the
procedure is given in the supplementary material, section A-1.

For inference, we use the true model but consider either the GC or the SC
constraints. Table 2 shows the estimates of the hyperparameters and the fixed
effect for the two sets of constraints while Figure 1 compare the posterior
mean of the interaction random effects for GC and SC, respectively. Figure A-1
in the appendix shows the posterior standard deviations. Exact correspondence
is not to be expected, which has been explained in Section 3. Even though the
new approach introduce an extra approximative $\kappa$ term (adding the value
to get a proper normal distribution), the new method is expected to be more
accurate. We do however see that there is a high correspondence between the
two methods: the posterior means of the random effects are very similar, while
the uncertainties (Figure A-1) are somewhat lower for the HyMiK method. The
most striking difference between the parameter estimates is the estimate of
$\tau_{\delta}$ which is somewhat larger for the standard method compared to
the HyMiK method. This is valid for GC constraints only. We discuss this and
computing times in Section 4.4.

| GC | SC
---|---|---
| Standard | HyMiK | Standard | HyMiK
| mean | sd | mean | sd | mean | sd | mean | sd
$\mu$ | 1.50 | 0.00 | 1.50 | 0.00 | 1.50 | 0.00 | 1.50 | 0.00
$\tau_{\alpha}$ | 61.29 | 27.36 | 58.45 | 26.21 | 57.08 | 25.10 | 57.46 | 25.44
$\tau_{\gamma}$ | 9.71 | 0.62 | 9.70 | 0.61 | 9.70 | 0.61 | 9.70 | 0.61
$\tau_{\delta}$ | 24.99 | 0.92 | 20.82 | 0.82 | 17.07 | 0.66 | 17.07 | 0.66
Av. marg.lik | -2.184 | -2.183 | -2.165 | -2.183
Comp. time | 37 sec | 26 sec | 93 sec | 24 sec
Table 2: Simulated data: Estimates of the intercept and the hyperparameters based on the GC and the SC constraints for both the standard (conditioning by kriging) and the new HyMiK approach. |
---|---
Figure 1: Simulated data: Comparison of the posterior means between the
standard method ($x$-axis) and the HyMiK approach ($y$-axis) for the
interaction effects, based on the GC (left) and the SC (right) constraints.
Standard deviations are shown in Figure A-1 in the appendix.

###  4.2 Dengue case study

In Lowe et al. (2021) a spatio temporal model is fitted on data recording
dengue fewer numbers between 2000 and 2019 given for the 558 microregions in
Brazil. Several covariates and models are tested. The authors use distributed
lag-non-linear models (Gasparrini, 2014). Essentially, the model is given by

$y_{t,s}|\bm{\eta}\stackrel{{\scriptstyle ind}}{{\sim}}\text{Neg-
Binom}(E_{s}\exp(\eta_{t,s}))$ (13)

with $\eta_{t,s}$ following the structure of equation (1b) where the
covariates $\bm{u}_{t,s}$ contain several derived features based on different
risk factors at different time-lags (we have chosen a total of 10 features,
somewhat smaller than used in the original work). The reader is referred to
Gasparrini (2014); Gasparrini et al. (2017) and the scripts within our GitHub
repository for details. Further, $E_{s}$ is the population in microregion $s$
divided by $100\,000$, so that response is equivalent to an incidence rate per
100 000 people.

From Figure 2 we observe very good correspondence for both the GC and SC
constraints. We observe, however, from Figure A-2 that the standard errors of
HyMiK parametrization are somewhat smaller than for the original
parametrization. The estimates of the regression parameters are given in Table
A-1, where we observe very good correspondence. The estimates of the
hyperparameters and the intercept are given in Table 3. We observe very good
correspondence, except for $\tau_{\delta}$ for the GC constraints, where there
is some deviation. There is also a relatively large deviation for the
$\tau_{\delta}$ parameter when estimating with SC constraints. We provide a
general discussion of this and also compare computing times in Section 4.4.

|
---|---
Figure 2: Dengue data: Comparison of the posterior means between the standard method and the new approach for the interaction effects, based on the GC (left) and the SC (right) constraints. Standard deviations are shown in Figure A-2 in the appendix. | GC | SC
---|---|---
| Standard | HyMiK | Standard | HyMiK
| mean | sd | mean | sd | mean | sd | mean | sd
$\mu$ | -1.732 | 0.708 | -1.844 | 0.708 | -0.130 | 9.589 | -0.013 | 9.924
1/Disper | 0.860 | 0.008 | 0.882 | 0.009 | 0.724 | 0.009 | 0.725 | 0.006
$\tau_{\alpha}$ | 5.957 | 1.129 | 5.916 | 0.997 | 11.592 | 1.550 | 2.637 | 0.630
$\tau_{\gamma}$ | 0.241 | 0.016 | 0.236 | 0.019 | 0.302 | 0.068 | 0.254 | 0.017
$\tau_{\delta}$ | 0.231 | 0.013 | 0.114 | 0.006 | 0.172 | 0.007 | 0.175 | 0.007
Av. marg.lik | -2.864 | -2.835 | -2.837 | -2.837
Comput.time | 187 sec | 103 sec | 386 sec | 101 sec
Table 3: Estimates of intercept and hyperparameters for the Dengue data based
on the GC and the SC constraints for both the standard and the new HyMiK
approach. Estimates of the remaining regression coefficients are given in
Table A-1.

###  4.3 Covid-19 data

Here we consider the daily numbers of positive tests within each county in
Norway during the Covid-19 pandemic. The data are available at
https://github.com/folkehelseinstituttet/surveillance_data/tree/master/covid19
and we have here considered data in the period 2020-09-02 to 2022-01-14 giving
$n_{T}=500$ days with counts available for each of the $n_{S}=11$ counties in
Norway. These data are given as cumulative numbers over time, but with some
inconsistencies in that in a few cases the cumulative numbers are lower than
the previous day. In those cases we have put the daily numbers to zero. We
assume the model

$\displaystyle y_{t,s}|\bm{\eta}\stackrel{{\scriptstyle
ind}}{{\sim}}\text{Poisson}(E_{s}\exp(\eta_{t,s}))$

with $\eta_{t,s}$ given in (1) without any covariates. Here $E_{s}$ is the
population size within county $s$, assumed to be constant within the time
period considered. Figure 3 indicates that the means of the latent variables
are very similar for both the GC and SC constraint. Also Figure A-3 indicates
very good correspondence between the standard and the HyMiK approaches for
making inference in the Knorr Held interaction models for both GC and SC
constraints. The estimates for the hyperparameters, Table 4, are also very
similar, except for the precision parameters for the interaction effects,
valid when estimating with GC constraints and also SC constraints. See Section
4.4 for a discussion.

|
---|---
Figure 3: Covid data: Comparison of the posterior means between the standard method (x-axis) and the new HyMiK approach (y-axis) for the interaction effects, based on the GC (left) and the SC (right) constraints. Standard deviations are shown in Figure A-3 in the appendix. | GC | SC
---|---|---
| Standard | HyMiK | Standard | HyMiK
| mean | sd | mean | sd | mean | sd | mean | sd
$\mu$ | -9.421 | 0.004 | -9.421 | 0.004 | -9.421 | 0.004 | -9.421 | 0.004
$\tau_{\alpha}$ | 9.522 | 0.706 | 9.521 | 0.711 | 9.528 | 0.711 | 9.498 | 0.696
$\tau_{\gamma}$ | 5.210 | 2.039 | 5.408 | 2.185 | 5.334 | 1.989 | 5.299 | 2.065
$\tau_{\delta}$ | 12.487 | 0.577 | 11.512 | 0.558 | 12.437 | 0.594 | 11.568 | 0.538
Av. marg.lik | -4.338 | -4.329 | -4.338 | -4.329
Comput.time | 31 sec | 10 sec | 30 sec | 9 sec
Table 4: Estimates of intercept and hyperparameters for the Covid data based
on the GC and the SC constraints for both the standardS and the new HyMiK
approach.

###  4.4 Discussion on results

We first consider the simulated and dengue fever datasets. Here $n_{S}\gg
n_{T}$. With the HyMiK method we can handle $n_{S}$ of the $n_{S}+n_{T}-1$ GC
constraints through the mixed effect model. For the SC constraints, $2n_{S}$
out of $2n_{S}+n_{T}-1$ constraints can be considered within the mixed effect
model. The remaining constraints
$[\bm{I}_{n_{T}}\otimes\bm{1}_{s_{N}}^{T}]\bm{\delta}=\bm{0}_{n_{T}}$ need to
be enforced through conditioning by kriging. This means the number of time
steps should not be too large. When considering the Covid-19 dataset, where
$n_{T}\gg n_{S}$, we can remove $n_{T}$ constraints from either set of
constraints. Conditioning by kriging is then used on the last constraints, a
relative low number due to that $n_{S}=11$ in this case.

The main difference between the dengue and simulated data cases and the
Covid-19 data, is the number of spatial and temporal regions. In the former
cases, we have many more regions than temporal time points and in the latter
case the opposite is true. This means we should take different approaches when
estimating the models. When $n_{S}\gg n_{T}$ we should take the approach used
in the dengue and simulated data cases, and when $n_{T}\gg n_{S}$ we should
use same approach as in the Covid-19 data problem,

The main benefit of the HyMiK method is the reduced run time of the model. We
observe that the computational factor is 3 for both the GC constraints and the
SC constraints when we consider the COVID data problem. When considering the
dengue dataset the computational factors are 1.8 and 3.8, where the factor is
greatest when using the SC constraint set. For the simulated dataset the
factors are 1.4 and 3.9. We therefore in all examples obtain a significant
reduction in computing time. We also observe that the computing time is
reduced the most for the SC constraints when we compare HyMiK and the original
method, which is as expected since the SC constraints have the largest number
of constraints in the original formulation, so that using the traditional
kriging method in INLA should take the largest amount of time. The reduction
in computing time for the large latent field in the dengue dataset problem is
considerable and very fast inference can be obtained with the HyMik approach.
This is natural as several constraints are dropped in HyMiK.

Somewhat surprisingly, HyMiK ran a little slower with GC constraints on the
COVID dataset and also the simulated dataset. This is not as expected. The
reason is that the projection matrix is more dense for the SC constrains, but
the other constraints remain unchanged.

For all case studies, the precision parameters are the parameters with
estimates that differ the most, somewhat surpricing given that the interaction
terms themselves are very similar. We find it plausible that the HyMiK is more
precise (as explained in Section 3), so that the estimates of the precision
parameters for the new method should be more accurate.

One way to tell how much of the variance that is explained by each of the main
effects and interaction effects is by considering the precision parameters. We
have previously noted that the hyperparameters differ when comparing the HyMiK
and the standard conditional by kriging method. We then get different results
when quantifying the proportion of variance explained by the different
effects. As previously explained, we expect the HyMiK method to be more exact,
as the kriging step is done in a lower dimensional space, meaning we should
consider the proportions by using the new method.

In this paper we have chosen not to scale the precision matrices when
performing inference (for scaling see Sørbye and Rue, 2014), as the focus is
on demonstrating differences and similarities of HyMiK and the condition by
kriging methods. However, some experiments where scaling was included gave
similar correspondence between the two methods for dealing with constraints.

We have also computed the average marginal likelihood for the different
models. As INLA does not include the determinant of $\bm{Q}$ for the generic0
model, we had to calculate this as a post-processing step to obtain the
correct average marginal likelihood (a routine for this correction is provided
in the GitHub repository). From Tables 2, 3 and 4 we see that the values are
very similar.

We have used the new HyMiK approach for estimation of models within INLA,
which makes the model less dependent on the number of observations. The
inference is still dependent on the size of the latent field. This is
demonstrated in the simulation experiment, where the number of observations is
30 times the size of the linear predictor in the usual space time set up.

##  5 Summary, general discussion and conclusion

The idea behind the HyMiK method is simple: Instead of handling all the
constraints by the conditioning by kriging approach, some of the constraints
are built in to the model through a mixed model approach. The remaining
constraints are handled by the conditioning by kriging approach, but now on a
(much) lower number of constraints. INLA is used to estimate the posterior as
it has the conditioning by kriging step and also the mixed model implemented.

Both the conditioning by kriging approach and the mixed model approach add
small numbers to the diagonals of the precision matrices in order for the
intrinsic models to become proper, which is needed for the numerical
calculations to not break down. In addition, INLA is based on a Gaussian
approximation for handling non-Gaussian data. Both methods thereby introduce
some numerical errors. We have however argued that the use of the mixed model
through the HyMiK approach should be numerically more accurate due to that
corrections for some of the constraints are performed prior to the Gaussian
approximation.

It is crucial to use the correct set of constraints when estimating the
different types of Knorr Held models. Fattah and Rue (2022) follow Schrödle
and Held (2011) while Gocioa and coauthors have proposed the constraints as
stated in Goicoa et al. (2018). Both say that constraining the field with
wrong constraint sets will give wrong estimates of model parameters. We can
fit both set of constraints. The user has to decide which constraints to use.
We note that while GC constraints may have less constraints, SC constraints
will always lead to a proper space time prior and can thus be used with
missing values present. Given the set of constraints, there is a design choice
in which constraints that should be handled by the mixed effect approach and
which by conditioning by kriging. We have suggested two possibilities, which
one to choose depends on the size of the temporal component compared to the
spatial component.

The recent paper Orozco-Acosta et al. (2023) also considers estimating large
spatio-temporal fields. The authors suggest to use a divide and conquer
approach. The approach is approximate since the spatial domain is divided into
smaller domains, where spatio-temporal models are fitted on the smaller
spatial domains. This approach can benefit from our model reparametrization as
well, as each submodel can utilize our HyMiK parametrization to obtain faster
estimation. The user can choose to use our approach with INLA, which is more
exact, or combine our reparametrization with Orozco-Acosta et al. (2023).
Combining Orozco-Acosta et al. (2023) with our reparametrization is further
work. We have focused on making a new computational efficient procedure. As we
can see, it can easily be combined with several other estimation procedures as
well. As Orozco-Acosta et al. (2023) also uses INLA, it should be relatively
straightforward to combine the two methods.

When comparing our method to the method of Fattah and Rue (2022), called
INLAPLUS, our method can also be used with GC constraints. This makes our
method more flexible. Our method can also be used on severs and ordinary
computers without vast number of cores available. We have found similar
speedups, which means our method is the preferred one when there are limited
computational resources available. However, INLAPLUS can give further
reduction in computing time for SC constraints when computational resources
are not an issue.

HyMiK is faster than the traditional way of estimating constrained models in
INLA, both for GC constraints and for SC constraints. We observe very close
correspondence between estimated mean and standard deviations of the field,
but that the hyperparameters can be somewhat different. Our method is expected
to be more accurate, as previously explained. We have implemented procedures
that can be used to set up the models easily. The preparation is fast. In
total, we argue that HyMiK should become the standard procedure of estimating
the Knorr Held type IV interaction model in INLA.

## Acknowledgment

This research is partly funded by the Norwegian Research Council research-
based innovation center BigInsight, project no 237718. We thank Professor
Håvard Rue for clarifying some issues related to how INLA is implemented. All
code is available at the GitHub repository
https://github.com/geirstorvik/INLAconstraints.

## References

  * Bakka et al. (2018) Bakka, H., H. Rue, G.-A. Fuglstad, A. Riebler, D. Bolin, J. Illian, E. Krainski, D. Simpson, and F. Lindgren (2018).  Spatial modeling with R-INLA: A review.  Wiley Interdisciplinary Reviews: Computational Statistics 10(6), e1443.
  * Besag (1974) Besag, J. (1974).  Spatial interaction and the statistical analysis of lattice systems.  Journal of the Royal Statistical Society: Series B (Methodological) 36(2), 192–225.
  * Cressie (1993) Cressie, N. (1993).  Statistics for spatial data (Revised edition ed.).  Wiley.
  * Fattah and Rue (2022) Fattah, E. A. and H. Rue (2022).  Approximate Bayesian Inference for the Interaction Types 1, 2, 3 and 4 with Application in Disease Mapping.  arXiv preprint arXiv:2206.09287.
  * Gaedke-Merzhäuser et al. (2023) Gaedke-Merzhäuser, L., J. van Niekerk, O. Schenk, and H. Rue (2023).  Parallelized integrated nested Laplace approximations for fast Bayesian inference.  Statistics and Computing 33(1), 25.
  * Gasparrini (2014) Gasparrini, A. (2014).  Modeling exposure–lag–response associations with distributed lag non-linear models.  Statistics in medicine 33(5), 881–899.
  * Gasparrini et al. (2017) Gasparrini, A., F. Scheipl, B. Armstrong, and M. G. Kenward (2017).  A penalized framework for distributed lag non-linear models.  Biometrics 73(3), 938–948.
  * Goicoa et al. (2018) Goicoa, T., A. Adin, M. Ugarte, and J. Hodges (2018).  In spatio-temporal disease mapping models, identifiability constraints affect PQL and INLA results.  Stochastic Environmental Research and Risk Assessment 32(3), 749–770.
  * Khana et al. (2018) Khana, D., L. M. Rossen, H. Hedegaard, and M. Warner (2018).  A Bayesian spatial and temporal modeling approach to mapping geographic variation in mortality rates for subnational areas with R-INLA.  Journal of data science: JDS 16(1), 147.
  * Knorr-Held (2000) Knorr-Held, L. (2000).  Bayesian modelling of inseparable space-time variation in disease risk.  Statistics in medicine 19(17-18), 2555–2567.
  * Lawson (2018) Lawson, A. B. (2018).  Bayesian disease mapping: hierarchical modeling in spatial epidemiology.  Chapman and Hall/CRC.
  * Lee et al. (2018) Lee, D., A. Rushworth, and G. Napier (2018).  Spatio-temporal areal unit modeling in R with conditional autoregressive priors using the CARBayesST package.  Journal of Statistical Software 84, 1–39.
  * Lowe et al. (2021) Lowe, R., S. A. Lee, K. M. O’Reilly, O. J. Brady, L. Bastos, G. Carrasco-Escobar, R. de Castro Catão, F. J. Colón-González, C. Barcellos, M. S. Carvalho, et al. (2021).  Combined effects of hydrometeorological hazards and urbanisation on dengue risk in Brazil: a spatiotemporal modelling study.  The Lancet Planetary Health 5(4), e209–e219.
  * Nazia et al. (2022) Nazia, N., Z. A. Butt, M. L. Bedard, W.-C. Tang, H. Sehar, and J. Law (2022).  Methods Used in the Spatial and Spatiotemporal Analysis of COVID-19 Epidemiology: A Systematic Review.  International Journal of Environmental Research and Public Health 19(14), 8267.
  * Orozco-Acosta et al. (2023) Orozco-Acosta, E., A. Adin, and M. D. Ugarte (2023).  Big problems in spatio-temporal disease mapping: Methods and software.  Computer Methods and Programs in Biomedicine 231, 107403\.
  * Pavía et al. (2008) Pavía, J. M., B. Larraz, and J. M. Montero (2008).  Election forecasts using spatiotemporal models.  Journal of the American Statistical Association 103(483), 1050–1059.
  * R-INLA (2020) R-INLA (2020).  The z-model.  https://inla.r-inla-download.org/r-inla.org/doc/latent/z.pdf.
  * Riebler et al. (2016) Riebler, A., S. H. Sørbye, D. Simpson, and H. Rue (2016).  An intuitive bayesian spatial model for disease mapping that accounts for scaling.  Statistical methods in medical research 25(4), 1145–1165.
  * Rue and Held (2005) Rue, H. and L. Held (2005).  Gaussian Markov random fields: theory and applications.  CRC press.
  * Rue et al. (2009) Rue, H., S. Martino, and N. Chopin (2009).  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71(2), 319–392.
  * Schrödle and Held (2011) Schrödle, B. and L. Held (2011).  Spatio-temporal disease mapping using INLA.  Environmetrics 22(6), 725–734.
  * Shoesmith (2013) Shoesmith, G. L. (2013).  Space–time autoregressive models and forecasting national, regional and state crime rates.  International journal of forecasting 29(1), 191–201.
  * Sørbye and Rue (2014) Sørbye, S. H. and H. Rue (2014).  Scaling intrinsic Gaussian Markov random field priors in spatial modelling.  Spatial Statistics 8, 39–51.
  * Van Niekerk et al. (2021) Van Niekerk, J., H. Bakka, H. Rue, and O. Schenk (2021).  New Frontiers in Bayesian Modeling Using the INLA Package in R.  Journal of Statistical Software 100, 1–28.
  * Van Niekerk et al. (2023) Van Niekerk, J., E. Krainski, D. Rustand, and H. Rue (2023).  A new avenue for Bayesian inference with INLA.  Computational Statistics & Data Analysis 181, 107692.
  * Vicente et al. (2020) Vicente, G., T. Goicoa, and M. Ugarte (2020).  Bayesian inference in multivariate spatio-temporal areal models using INLA: analysis of gender-based violence in small areas.  Stochastic Environmental Research and Risk Assessment 34, 1421–1440.
  * Whittaker (1922) Whittaker, E. T. (1922).  On a new method of graduation.  Proceedings of the Edinburgh Mathematical Society 41, 63–75.
  * Wikle (2003) Wikle, C. K. (2003).  Hierarchical models in environmental science.  International Statistical Review 71(2), 181–199.
  * Williams et al. (2019) Williams, D., J. Haworth, M. Blangiardo, and T. Cheng (2019).  A spatiotemporal Bayesian hierarchical approach to investigating patterns of confidence in the police at the neighborhood level.  Geographical analysis 51(1), 90–110.

SUPPLEMENTARY MATERIAL

##  A-1 Description of procedure for constructing simulated data

The simulated data in section 4.1 are generated by the following procedure:

  * •

Use the Germany graph to make $\bm{Q}_{\alpha}$.

  * •

Let the temporal component be a random walk of order 2 with 10 time points

  * •

Scale the precision matrices $\bm{Q}_{\alpha}$ and $\bm{Q}_{\gamma}$ to obtain
scaled versions $\bm{Q}^{*}_{\alpha}$ and $\bm{Q}_{\gamma}^{*}$ and set
$\bm{Q}_{\delta}^{*}=\bm{Q}^{*}_{\alpha}\otimes\bm{Q}^{*}_{\gamma}$. The
function inla.scale.model is used to scale the matrices.

  * •

Generate one sample from each of the IGMRF models using the SC constraints on
the interaction term and appropriate sum to zero constraints on the two
remaining main effects. This means proper priors are used. The function
inla.qsample is used to generate the samples.

  * •

Simulate 30 response variables from Poisson($\bm{\eta}_{t,s}$)for each linear
predictor $\eta_{t,s}$.

  * •

Fit the models using INLA.

Scaling of a precision matrix means it is divided by a constant. This constant
is found by calculating a mean of the diagonal elements of the inverse of the
precision matrix under suitable constraints. Then the “typical” marginal
variance of the prior is 1. See e.g Sørbye and Rue (2014) or Riebler et al.
(2016) for details. The script for simulating the data is available at our
GitHub repository under the name SimulateModel.R.

##  A-2 Additional results

Here, figures showing the similarities between the posterior standard errors
of the interaction terms for the simulated data, the dengue data and the
Covid-19 data are shown in figures A-1, A-2 and A-3, respectively. Further,
Table A-1 contain estimates of the regression coefficients for the dengue data
set.

|
---|---
Figure A-1: Simulated data: Comparison of the posterior standard deviation between the standard method and the new approach for the interaction effects, based on the GC (left) and the SC (right) constraints. |
---|---
Figure A-2: Dengue data: Comparison of the posterior standard deviation between the standard method and the new approach for the interaction effects, based on the GC (left) and the SC (right) constraints. | GC | SC
---|---|---
| Standard | HyMiK | Standard | HyMiK
| mean | sd | mean | sd | mean | sd | mean | sd
$\beta_{1}$ | 0.029 | 0.003 | 0.029 | 0.003 | 0.027 | 0.003 | 0.027 | 0.003
$\beta_{2}$ | -0.121 | 0.182 | -0.091 | 0.181 | -0.254 | 0.185 | -0.239 | 0.185
$\beta_{3}$ | 1.048 | 0.147 | 1.033 | 0.147 | 1.031 | 0.150 | 1.019 | 0.151
$\beta_{4}$ | -0.568 | 0.115 | -0.557 | 0.115 | -0.412 | 0.118 | -0.424 | 0.119
$\beta_{5}$ | 1.201 | 0.571 | 1.277 | 0.570 | 0.776 | 0.583 | 0.793 | 0.584
$\beta_{6}$ | -1.476 | 0.466 | -1.352 | 0.466 | -1.293 | 0.477 | -1.307 | 0.478
$\beta_{7}$ | -0.321 | 0.372 | -0.294 | 0.372 | -0.036 | 0.385 | -0.062 | 0.386
$\beta_{8}$ | 1.953 | 0.135 | 1.993 | 0.135 | 1.974 | 0.139 | 2.001 | 0.139
$\beta_{9}$ | 2.465 | 0.109 | 2.419 | 0.108 | 2.346 | 0.111 | 2.321 | 0.111
$\beta_{10}$ | -1.769 | 0.081 | -1.767 | 0.081 | -1.843 | 0.085 | -1.858 | 0.085
Table A-1: Estimates of regression coefficients for the Dengue data. |
---|---
Figure A-3: Covid data: Comparison of the posterior standard deviation between
the standard method and the new approach for the interaction effects, based on
the GC (left) and the SC (right) constraints.

-


# Optimal Load Shedding for Public Safety Power Shutoffs

Aniruddha Rajendra Rao
Industrial AI Lab
Hitachi America, Ltd., R&D
Santa Clara, USA
<EMAIL_ADDRESS>
Chandrasekar Venkatraman
Industrial AI Lab
Hitachi America, Ltd., R&D
Santa Clara, USA
<EMAIL_ADDRESS>
Robert Ellis
CDS-DSL
Hitachi Limited
Kokubunji, Tokyo, Japan
<EMAIL_ADDRESS>
Chetan Gupta
Industrial AI Lab
Hitachi America, Ltd., R&D
Santa Clara, USA
<EMAIL_ADDRESS>

###### Abstract

Public utilities are faced with situations where high winds can bring trees
and debris into contact with energized power lines and other equipments, which
could ignite wildfires. As a result, they need to turn off power during severe
weather to help prevent wildfires. This is called Public Safety Power Shutoff
(PSPS). We present a method for load reduction using a multi-step genetic
algorithm for Public Safety Power Shutoff events. The proposed method
optimizes load shedding using partial load shedding based on load importance
(critical loads like hospitals, fire stations, etc). The multi-step genetic
algorithm optimizes load shedding while minimizing the impact on important
loads and preserving grid stability. The effectiveness of the method is
demonstrated through network examples. The results show that the proposed
method achieves minimal load shedding while maintaining the critical loads at
acceptable levels. This approach will help utilities to effectively manage
PSPS events and reduce the risk of wildfires caused by the power lines.

Keywords: Optimization, Genetic Algorithm, Load Distribution, PSPS, Load
Shedding.

##  1 Introduction

Strong winds can cause tree branches and debris to come in contact with active
power lines, resulting in equipment damage and an increase in the risk for
wildfires. Hence, utility companies use Public Safety Power Shutoff (PSPS)
during severe weather to prevent wildfires. The current methodology involves
forecasting of severe weather up to a week before power is shutoff and
notifying customers two days in advance. Efficiently determining areas for
power shutoff is challenging. Early shutoff causes revenue loss, while
delaying increases wildfire risk. Recent tragic wildfires highlight the need
for improved strategies due to electrical faults causing significant damage.
For example, in the year 2011, CBC [CBC, 2011] reported two wildfires in
Bastrop county that got started by trees coming in contact with nearby power
lines and that became the most destructive wildfires in Texas’s history. It
resulted in four deaths and caused over $300 million in damage. In California,
the 2018 Camp Fire, which was ignited by a power line, led to the death of 84
people and caused approximately $9.3 billion in residential property damage
alone [T. et al., 2019]. The severity of this and other fires that occurred
during the 2017 and 2018 California fire seasons led Pacific Gas & Electricity
(PG&E), the responsible utility company, to file for bankruptcy and face
charges for involuntary manslaughter [PG&E, 2019b, PG&E, 2020]. In 2021, it
was confirmed that California’s Dixie Fire – the second-largest wildfire in
California’s history - was sparked when power lines owned by Pacific Gas and
Electric (PG&E) came into contact with a tree. The Dixie fire went on to burn
over 950 thousand acres across five counties and lasted almost 4 months. The
more recent wildfires on going in different parts of North America and Europe
are causing catatrophic losses. There is a need to effectively mitigate such
situations.

The likelihood of wildfires caused by power infrastructure is higher during
windy conditions, which can cause power failures and make wildfires more
difficult to contain. Studies show that wildfires ignited by power lines are
usually larger and more damaging than others [Miller et al., 2017]. Power
infrastructure can cause ignitions in various ways, with the most common being
contact between vegetation and conductors [Jazebi et al., 2020]. Although
efforts have been made to reduce the probability of ignitions, day-to-day
operations leave utilities with limited options to reduce the wildfire risk.
The only way to completely avoid the risk of ignitions is to de-energize the
power line, this intentional black outs are referred to as Public Safety Power
Shutoffs (PSPS) [PG&E, 2019a, CEC, 2019]. But this impacts the power system’s
ability to provide reliable electricity and can result in intentional
blackouts affecting millions of people [Anderson and Bell, 2012]. Our goal is
to mitigate such a situation with optimal load shedding.

Load shedding is a crucial aspect of power system management during extreme
conditions, such as wildfires, storms or any other calamity, where there is a
risk of power failures and potential fire ignitions. Over the years,
researchers have focused on optimizing different approaches for load shedding
and energy generation to reduce the impact of power failures [Zanocco et al.,
2021, Kody et al., 2022, Rhodes et al., 2021]. However, there is still a gap
in addressing the problem of load shedding for public safety power shut-offs
(PSPS) with partial load shedding while also considering load importance.
Despite the availability of various optimization approaches, none of them
tackles these specific challenges associated with PSPS, which can have
significant economic and societal impacts. Hence, it is crucial to address
these challenges to ensure that load shedding is conducted in a way that
minimizes the adverse impact on the public while still maintaining safety
during extreme conditions.

Our study addresses utility companies’ challenges in managing PSPS while
minimizing load shedding. It lays the basis for a tool aiding incident
managers to identify optimal load shedding while maintaining grid stability.
Our novel contributions include an optimization framework for minimizing
wildfire and cascading failure risks during transmission line failures,
partial load shedding considering load importance, and demonstrating the
approach’s efficacy using the RTS-GMLC test case under N-1 contingency. The
paper’s structure comprises problem formulation and our load shedding approach
details in Section 2, RTS-GMLC case study in Section 3, and concluding remarks
with future research directions in Section 4.

##  2 Optimal Load shedding

One of the most accurate way to solve the load shedding problem would be
through discrete/binary combinatorial optimization. In this, we generate all
possible solutions and then select the best one according to a predefined
fitness function (objective function). For problems with ”N” variables, the
total number of solutions that need to be generated is $2^{N}$ for binary
optimization. However, generating so many solutions takes time, and the
limitations of computational devices further restricts solutions generation
[Wosley and Nemhauser, 2014]. Traditional deterministic approaches are rigid
and don’t factor uncertainty. To overcome this, recent studies have proposed
stochastic methods such as simulated annealing, ant colony, particle swarm
optimization (PSO), and genetic algorithm (GA) to reduce the time required to
find the best solution. Among these methods, GA [Rajappa, 2012] is effective
in solving discrete as well as binary combinatorial optimization problems. The
GA is capable of real-time applications as it converges rapidly.

Therefore, our research is based on applying GA for optimal load shedding for
PSPS and validating it with grid networks, primarily for transmission lines
[Fernando et al., 2019, Rao et al., 2013]. We further enhanced the research
output to address partial load shed, i.e. instead of completely turning off a
load in a bus, to shed partial load, thus achieving minimal load shed. The
second enhancement is related to ‘load importance’ – while minimal shedding is
the desired state, it is also important to consider where the power is being
consumed – such as hospitals, critical infrastructure, call centers, fire
stations, etc. Optimal load shedding given load importance is crucial to
minimise the adverse economic and societal impact of the load shedding on
consumers [Hassan et al., 2019, Jabian et al., 2018]. Our research addresses
how to provide recommendations to incident commanders once they have an idea
of ignition risk to optimally choose which loads to shed, as shown in Figure
1.

Figure 1: Overview of our Solution

In Genetic Algorithm (GA) we have to define a fitness function that we want to
maximize. We would like to define a fitness function such that the system is
stable and load shedding is minimal and practical. The algorithm execution
involves several steps.

In Step 1, a sample of the load system is taken from the grid network using a
power flow analysis tool to act as the setup for the GA. The GA first
initializes a population with N random solutions for optimization. Each bus
with load (called as gene) in the network is initialized with a random load
value from unit range. These random solutions are called chromosomes. The
length of the chromosome is equal to the number of loads in the system.

In Step 2, the fitness function is evaluated for each chromosome and ranked
accordingly. The most fit chromosomes (parents) are selected from the
population to generate the offspring for the next generation. Number of
parents selected will be M. Typically, this is done using some selection
strategy, such as tournament selection or roulette wheel selection.

In Step 3, offspring solutions (chromosomes) are created through crossover
operations, such as single-point or multi-point crossover. In single-point
crossover, genetic material between two parents is swapped at a random point,
yielding two offspring solutions with mixed genetic material. In multi-point
crossover, multiple points are chosen for swapping, resulting in even more
diverse offspring solutions compared to single-point crossover.

In Step 4 we perform mutation. A small random change is introduced into the
offspring solutions to explore new regions of the search space. This is
typically done using a mutation operator, which randomly alters a small
portion of the offspring’s genetic material.

In Step 5, the fitness function is computed for the new chromosomes. Only the
top N fittest chromosomes remain in the population as the new generation. Step
3 and Step 4 is repeated till specified number of iterations are completed or
stopping criteria is met. The chromosome with the lowest load shedding in the
final generation is selected as the best solution.

For our approach, we set the values of most genes in each chromosome as 1,
since ideally we would like to have all the loads operating at full capacity.
This helps us to reach the solution faster. We use single-point crossover and
a mutation factor of 0.1. A load with value 0 means that is it completely
shut, a load with value 1 implies it is at 100% capacity and a load with 0.8
value means it is operating at 80% capacity.

We incorporate load importance using multi-step genetic algorithm approach.
Multi-step GA is a type of optimization algorithm that uses a combination of
GA and multiple stages of optimization to solve complex problems. In the
context of optimal load shedding, multi-step GA can be used to optimize load
shedding for multiple loads with different priorities (importance) as seen in
Figure 2. We generate load importance values randomly using a standard
distribution but in practice, the incident managers will use domain knowledge
and experience to specify load importance.

Figure 2: Our multi-step partial load genetic algorithm approach for optimal
load shedding

The first step of the multi-step GA is to utilize the load importance values
to define different stages or levels of optimization, with each stage focusing
on a specific subset of loads based on their importance.

For example, in a power network with multiple loads, some loads may be more
critical than others in terms of their importance to the overall functioning
of the grid. By assigning importance values to each load, a multi-step
approach can be designed to optimize the loads in stages, with the least
important loads being optimized for shedding first, followed by the next set
of loads with increasing importance if necessary.

In the first stage, the least important loads are optimized using a GA, while
the other important loads are assumed to be 1 (functioning at full capacity),
with the fitness function being defined to make sure the network is stable and
has maximum load capacity (or minimum load shedding). Once the optimal
solution for this stage is found, the loads are adjusted accordingly in the
network to make it stable otherwise we move on to the next stage.

In the second stage, the next set of important loads are optimized along with
the least important loads, meaning the number of genes (loads) in a chromosome
has increased slightly, while keeping the more important loads set to 1. This
process is repeated for each subsequent stage until a solution is reached or
till the last stage where all the loads are included to find a solution. Each
stage is optimizing a higher subset of loads based on their importance.

By breaking down the optimization process into stages based on load
importance, the GA can be run sequentially to optimize the loads in a more
efficient manner. This approach also ensures that the most critical loads are
given priority in the optimization process, leading to a more reliable and
stable power grid. The main advantage of a multi-step GA is that it can handle
complex optimization problems that are difficult to solve using traditional
optimization methods. However, the algorithm requires careful design and
tuning to ensure that it converges to a good solution in a reasonable amount
of time.

##  3 Case Study

We assess our approach by conducting tests using a power system located in an
area that faces a high risk of wildfires. We consider the RTS-GMLC [Barrows et
al., 2020] (Reliability Test System - Grid Modernization Laboratory
Consortium) network, which is a widely used benchmark power system test case
for research and development in the field of power system analysis and
optimization. The test case was developed by the Grid Modernization Laboratory
Consortium of the US Department of Energy, and it includes detailed modeling
of the power grid, including generators, transmission lines, and loads.

The RTS-GMLC test case is based on the IEEE Reliability Test System (RTS) and
includes additional features such as renewable generation, demand response,
and energy storage systems. The RTS-GMLC test case includes 73 buses, 104
transmission lines (with the exception of the HVDC line being removed), 51
loads, 3 shunts and 18 generators as seen in Figure 3. The RTS-GMLC system is
geographically located in southern California, Nevada, and western Arizona.
These areas cover a lot of high risk wildfire areas.

The fitness function was designed to optimize for system safety and the
remaining load in the system, with the maximum fitness value being 18726. It
is defined as $FitnessFunction=I*(10000)+L$, where $I$ is an indicator
function for system safety and the $L$ is measuring the sum of the remaining
loads (addition of Active and Reactive power at each load). A system is
considered safe when none of the lines have a load percent more than 100 and
the load bus voltage is between 0.95 to 1.05 of the nominal voltage. The
safety criterion can be modified according to the problem in hand. We check
for system safety using PandaPower [Thurner et al., 2018]. It is a power
systems analysis toolbox which is used for automation of analysis and
optimization in power systems. The results obtained from this analysis can
provide insights into the optimal load shedding strategies for maintaining
system safety during line outages that can happen due to multiple reasons like
wildfires, repairs, etc. We loop through the cases where each of the 104 lines
in the network are down one at a time.

Figure 3: RTS-GMLC network with no power shutoff applied

Table 1 shows the result of optimal load shedding under different approaches
(Binary optimal load shedding and Partial optimal load shedding using GA) and
conditions (different saturation values). We can use the saturate
functionalities to stop the GA when the fitness value does not change for a
defined number of iterations. We ran the GA for 500 iterations under no
conditions and observed that for both approaches we are able to find solutions
in 17 cases when a line is down, there is no feasible solution in 2 cases when
a line is down, and for 85 cases there is no instability in the network when a
line goes down. No feasible solution might be because of the limitation in
pandapower to run contingency analysis for these 2 cases. This could be
overcome by using Linear Sensitivity analysis [Jobanputra et al., 2017] but
this is not considered in this work. Our results are insightful, because
mitigation measures, like adding tie lines, can be built for the 2 cases (Line
9 and Line 10) when there is a potential for cascading failure. This
demonstrates how our approach can be used to determine the robustness of a
power grid network. We get similar level of performance with saturation value
of 25 slight increase in number of no feasible solution. The saturation value
is a important hyper parameter which helps in increasing the speed but also
decreases efficacy if the value is too low. From the last column in Table 1,
we see that the total time to execute our approach on the 104 cases is minimal
compared to using brute force. The performance of Binary vs Partial is similar
though partial GA takes a bit longer to find the solution as we have increased
the search space for each load from 0 or 1 to 0, 0.1, 0.2,….,1 (a step
increase of 10%). We can reduce the step increment even further to decrease
the amount of load shed. Also, we our approach is more reliable than the
incident manager, who takes actions based on just experience without any
mathematical backing.

The fitness value is also higher for partial approach since we are considering
partial load shedding. This is visible in Table 2, where we compare the load
shed for each case when a line goes down and we must find an optimal load
shedding solution. We can see from the table that the partial load shedding
has a lower load shed amount (active power, p_mw) every time compared to
binary load shedding. In some cases, the gain in performance is huge. This
means that we can avoid unnecessary load shedding with this mechanism. The
solution of partial optimal load shedding can be improved upon by running for
large iterations over a broader search space.

Table 1: Results using Genetic Algorithm for N-1 contingency Approach and Condition |  Number of lines for which there is no feasible solution |  Number of lines for which there is solution |  Number of lines for which there is no instability |  Total Run Time (seconds)
---|---|---|---|---
Binary and no condition | 2 | 17 | 85 | 1209
Binary with Saturate=25 | 6 | 13 | 85 | 95
Partial and no condition | 2 | 17 | 85 | 1763
Partial with Saturate=25 | 3 | 16 | 85 | 17
Table 2: Comparison of Partial load shedding vs Binary load shedding using Genetic Algorithm. Line down (N-1 contingency) |  Total partial load shed (MW) |  Load number associated with partial load shedding (load percent) |  Total binary load shed (MW) |  Load number associated with binary load shedding
---|---|---|---|---
4 | 13.6 | 5 (0.9) | 136 | 5
8 | 40.8 | 5 (0.7) | 136 | 5
16 | 13.6 | 5 (0.9) | 74 | 3
18 | 34.2 | 7 (0.8) | 171 | 7
19 | 26.6 | 4 (0.9), 9 (0.9) | 97 | 1
23 | 54 | 2 (0.7) | 97 | 1
40 | 13.6 | 22 (0.9) | 272 | 5, 22
44 | 40.8 | 22 (0.7) | 244 | 0, 22
46 | 34.2 | 24 (0.8) | 261 | 5, 23
47 | 34.2 | 24 (0.8) | 331 | 9, 23
73 | 13.6 | 39 (0.9) | 136 | 39
75 | 32.4 | 34 (0.7) | 71 | 38
76 | 42.6 | 38 (0.6) | 74 | 37
77 | 40.8 | 39 (0.7) | 136 | 39
79 | 25 | 40 (0.8) | 171 | 41
80 | 34.2 | 41 (0.8) | 125 | 40
90 | 72 | 36 (0.6) | 180 | 36
Table 3: Comparison of Partial load shedding vs Partial load shedding with importance using Genetic Algorithm. Line down (N-1 contingency) |  Total partial load shed (MW) |  Load number associated with partial load shedding (load percent, importance) |  Total parital and important load shed (MW) |  Load number associated with partial and important load shedding (load percent, importance)
---|---|---|---|---
19 | 26.6 | 4 (0.9, 0.825), 9 (0.9, 0.959) | 27.2 | 5 (0.8, 0.711)
23 | 54 | 2 (0.7, 0.915) | 75.6 | 0 (0.3, 0.794)
80 | 34.2 | 41 (0.8, 0.900) | 37.5 | 40 (0.7, 0.734)
90 | 72 | 36 (0.6, 0.750) | 97.2 | 34 (0.1, 0.566)

Figure 4 shows the network information when line 77 is down. We see the
original network information on the left, followed by the network information
when line 77 is down. It can be seen that line 73 is overloaded (line load is
128%). The third sub-figure from the left shows the binary load shedding
solution and the last one is for the partial load shedding solution. In both
cases, load 39 at bus 53 is shut but in the partial case, load 39 is operating
at 70% capacity (in other words, 30% load is shed) for the system to be
stable, which is more optimal than completely shutting that load. Hence,
partial load shedding helps in minimizing the economic and societal impacts
while giving better service to the customers.

Figure 4: Solution comparison when line 77 is down

Next, we look into load importance. We assign random load importance to the 51
loads using beta distribution ($\alpha$=5 and $\beta$=1). In practice, these
information will be know by the incident managers. In our Multi-step GA
approach, in the first step, we only consider non-critical loads (load
importance $<$ 0.8) for shedding. If we don’t find an optimal load shedding
solution in the first step, then we consider all the loads for load shedding
in the second step. The performance efficacy is similar to Table 1 with slight
increase in execution time. Table 3 shows how partial optimal load shedding
performs with and without load importance. For example, we can see that when
line 23 is down, and if we don’t consider load importance, then the solution
is to set load 2 at 70% capacity resulting in 54 MW of load shedding but the
load importance of load 2 is 0.915 whereas when we consider load importance,
the solution is to set load 0 at 30% capacity resulting in 75.6 MW load
shedding and load 0 has lower importance (0.794) compared to load 2. Table 3
indicates that we are always successful in making sure no important load is
shed if an alternate solution exists. Also, we are able to achieve this with
marginal increase in load shedding. In the other cases, the load number shed
is the same with or without importance consideration and hence not shown in
the table. Also, we don’t consider the binary case for load importance as the
partial load shedding is superior.

Figure 5: Load convergence of GA.

The convergence of the load shedding is demonstrated in Figure 5. The GA runs
until the specified number of iterations is reached or a stopping criterion is
met to make the network safe and the total load shed is minimized. We can see
that the gains in load shedding after some iterations lead to minimal increase
in performance and independent runs (denoted using different colors) give a
similar amount of load shedding. The independent runs give different solutions
due to the stochastic nature of the GA, that can be useful for exploring the
solution space and identifying trade-offs between different objectives. The
incident managers can use this to their benefit along with domain knowledge to
select the best solution.

In addition, we also investigated optimal load shedding under the case of N-K
contingency, and the results showed similar success. Moreover, we analyzed a
larger network, CASE9241PEGASE [Fliscounakis et al., 2013], and were able to
successfully find solutions for it under N-1 contingency.

##  4 Conclusions

In this paper, our approach offers significant benefits, including optimal
load shedding recommendations within a quick time frame, flexibility in
choosing loads for shedding and the ability of partial load shedding to
minimize overall impact. Also, the approach assigns importance to critical
loads and enables the incident commander to make informed decisions quickly by
taking all relevant factors into account. We have shown the efficacy of our
proposed approach using the RTS-GMLC test case. Also, we have shown that our
approach gives multiple solutions, and along with external information, the
incident managers can take the best decision possible. Finally, along with
optimal load shedding, our approach help us to understand the electric grids
better, as well as determine crucial lines, critical loads, which can lead us
to improve the power system stability, reliability, and power quality.

There are several potential areas for further development and exploration. Our
solution can be applied at finer granularity for optimal partial load
distribution at bus level. Further development can involve testing on various
networks, extensive N-K contingency cases for robustness assessment, and
collaborating with utility companies for usability feedback.

## References

  * [Anderson and Bell, 2012] Anderson, B. and Bell, M. (2012).  Lights out impact of the august 2003 power outage on mortality in new york, ny.  Epidemiology (Cambridge, Mass.), 23:189–93.
  * [Barrows et al., 2020] Barrows, C., Bloom, A., Ehlen, A., Ikäheimo, J., Jorgenson, J., Krishnamurthy, D., Lau, J., McBennett, B., O’Connell, M., Preston, E., Staid, A., Stephen, G., and Watson, J.-P. (2020).  The ieee reliability test system: A proposed 2019 update.  IEEE Transactions on Power Systems, 35(1):119–127.
  * [CBC, 2011] CBC (2011).  Texas wildfire likely caused by power line sparks.
  * [CEC, 2019] CEC (2019).  Southern california edison company’s 2019 wildfire mitigation plan.
  * [Fernando et al., 2019] Fernando, W., Rankothge, W., Perera, A., Dissanayake, S., and De Silva, W. (2019).  Optimization of customer-friendly manual load shedding system.  In ICAC, pages 476–480.
  * [Fliscounakis et al., 2013] Fliscounakis, S., Panciatici, P., Capitanescu, F., and Wehenkel, L. (2013).  Contingency ranking with respect to overloads in very large power systems taking into account uncertainty, preventive, and corrective actions.  IEEE Transactions on Power Systems, 28(4):4909–4917.
  * [Hassan et al., 2019] Hassan, Y., Rashid, Y., and Tuamiah, F. (2019).  Demand priority in a power system with wind power contribution load shedding scheme based.  Journal of Engineering, 25:92–110.
  * [Jabian et al., 2018] Jabian, M. E., Funaki, R., and Murata, J. (2018).  Load shedding optimization considering consumer appliance prioritization using genetic algorithm for real-time application.  IFAC-PapersOnLine, 51(28):486–491.  10th IFAC Symposium on CPES 2018.
  * [Jazebi et al., 2020] Jazebi, S., de León, F., and Nelson, A. (2020).  Review of wildfire management techniques—part i: Causes, prevention, detection, suppression, and data analytics.  IEEE Transactions on Power Delivery, 35(1).
  * [Jobanputra et al., 2017] Jobanputra, J., Kotwal, C., and Thakkar, J. (2017).  Linear sensitivity analysis and series compensation for transmission congestion management.
  * [Kody et al., 2022] Kody, A., West, A., and Molzahn, D. K. (2022).  Sharing the load: Considering fairness in de-energization scheduling to mitigate wildfire ignition risk using rolling optimization.  In 2022 IEEE 61st Conference on Decision and Control (CDC), pages 5705–5712.
  * [Miller et al., 2017] Miller, C., Plucinski, M., Sullivan, A., Stephenson, A., Huston, C., Charman, K., Prakash, M., and Dunstall, S. (2017).  Electrically caused wildfires in victoria, australia are over-represented when fire danger is elevated.  Landscape and Urban Planning, 167:267–274.
  * [PG&E, 2019a] PG&E (2019a).  Pacific gas and electric company amended 2019 wildfire safety plan.
  * [PG&E, 2019b] PG&E (2019b).  Press release: Pg&e files for reorganization under chapter 11.
  * [PG&E, 2020] PG&E (2020).  Press release: Pg&e reaches plea agreement on state charges related to 2018 camp fire.
  * [Rajappa, 2012] Rajappa, G. P. (2012).  Solving combinatorial optimization problems using genetic algorithms and ant colony.
  * [Rao et al., 2013] Rao, K. U., Bhat, S. H., Ganeshprasad, G. G., Jayaprakash, G., and Pillappa, S. N. (2013).  A novel grading scheme for loads to optimize load shedding using genetic algorithm in a smart grid environment.  In 2013 IEEE Innovative Smart Grid Technologies-Asia (ISGT Asia), pages 1–6.
  * [Rhodes et al., 2021] Rhodes, N., Ntaimo, L., and Roald, L. (2021).  Balancing wildfire risk and power outages through optimized power shut-offs.  IEEE Transactions on Power Systems, 36(4):3118–3128.
  * [T. et al., 2019] T., J., S., Y., D., M., F., C., and R., T. (2019).  2019 wildfire risk report.
  * [Thurner et al., 2018] Thurner, L., Scheidler, A., Schafer, F., Menke, J. H., Dollichon, J., Meier, F., Meinecke, S., and Braun, M. (2018).  pandapower - an open source python tool for convenient modeling, analysis and optimization of electric power systems.  IEEE Transactions on Power Systems.
  * [Wosley and Nemhauser, 2014] Wosley, L. A. and Nemhauser, G. L. (2014).  Integer and Combinatorial Optimization.  Wiley-Intersience.
  * [Zanocco et al., 2021] Zanocco, C., Flora, J., Rajagopal, R., and Boudet, H. (2021).  When the lights go out: Californians’ experience with wildfire-related public safety power shutoffs increases intention to adopt solar and storage.  Energy Research & Social Science, 79:102183.