# Bus Admittance Matrix Revisited: Is It Outdated on Modern Computers?

Hantao Cui This work is supported by U.S. National Science Foundation, Award
#2226826.

###### Abstract

Bus admittance matrix is widely used in power engineering for modeling
networks. Being highly sparse, it requires fewer CPU operations when used for
calculations. Meanwhile, sparse matrix calculations involve numerous indexing
and scalar operations, which are unfavorable to modern processors. Without
using the admittance matrix, nodal power injections and the corresponding
sparse Jacobian can be computed by an element-wise method, which consists of a
highly regular, vectorized evaluation step and a reduction step. This paper
revisits the admittance matrix from the computational performance perspective
by comparing it with the element-wise method. Case studies show that the
admittance matrix method is generally slower than the element-wise method for
grid test cases with thousands to hundreds of thousands of buses, especially
on CPUs with support for wide vector instructions. This paper also analyzes
the impact of the width of vector instructions and memory speed to predict the
trend for future computers.

###### Index Terms:

Bus admittance matrix, sparse matrix, high-performance computing,
vectorization, single-instruction multiple data (SIMD).

##  I Introduction

Bus admittance matrix (also known as admittance matrix, $\mathbf{Y}$ matrix,
or $\mathbf{Y}_{\text{bus}}$) is ubiquitous in power engineering for network
modeling [1]. The bus admittance matrix reduces a power grid with $n_{b}$
buses and $n_{l}$ lines into an $n_{b}$-port network by relating bus current
injections to bus voltages. In realistic systems, a bus is only connected to a
few other buses, thus the admittance matrix is highly sparse. By adopting
sparse matrix techniques, calculations based on $\mathbf{Y}_{\text{bus}}$
requires significantly less memory and arithmetic operations, making it
superior to the then-popular impedance matrix method [2].

In many power systems studies such as power flow and stability simulations,
$\mathbf{Y}_{\text{bus}}$ is used to calculate bus power injections, as well
as the Jacobian matrix that contains the partial derivatives of injections
concerning voltages and angles. The computational performance to obtain
network injections and Jacobians is highly relevant because their computation
time is only second to solving sparse linear equations. Over the past decades,
the central processing units (CPUs) in computers have seen significant
upgrades and shifted growth paradigms [3], but not all computations can
benefit equally. Calculations using $\mathbf{Y}_{\text{bus}}$ require sparse
matrix-vector products, where both the sparse matrix and vector are composed
of complex numbers. Such calculations are non-atomic to processors, thus the
performance gain depends on multiple factors and remains to be investigated.

Various techniques have been developed for modern CPUs with the primary goal
to reduce latency [4, 5]. In the early 1980s, arithmetic calculations were
already faster than accessing data in memory, thus caches were added to CPU
chips. Instruction pipelining was later added to execute multiple, similar,
and independent instructions in different stages of completion. In the 90s,
branch prediction and out-of-order execution were added for relaxing dependent
calculation. With a probability of success, dependent calculations can be
executed in parallel. Another relevant technique is vector instruction known
as single-instruction multiple data (SIMD), which allows the processing of
data in packs. To benefit from these techniques and achieve high performance,
the rule of thumb is to design algorithms and data structures that exhibit
regularity [6], which translates to predictability for compilers and the
hardware.

The sparse $\mathbf{Y}_{\text{bus}}$ method is not unique for network
injection calculation. An alternative method is termed the “element-wise
method” (or the vectorized method), which calculates the element-wise power
injections to form branch-power vectors. Power injections into the line
terminal are then summed up at the connected buses to obtain the bus-wise
power injections [7, 8]. The equations for the element-wise method are not
unfamiliar, but the method is often dismissed because typical power grids have
more lines than buses. Nevertheless, calculations in the element-wise method
exhibit repetitiveness, namely, the same calculations on different data, which
may be favored by modern CPUs.

This paper aims to provide a high-performance computing perspective on the
$\mathbf{Y}_{\text{bus}}$ method and the element-wise method for a) bus power
injection calculation and b) sparse Jacobian matrix formulation. The method
formulations, optimized implementations, data access patterns, and
computational complexity are presented. While this work does not study the
solution of equations, the formulated equations and sparse Jacobians can be
readily interfaced with efficient solvers. Rather, the objective is to
identify the method with higher performance and understand its scalability
limitations on modern CPUs.

The main contribution of this paper is the finding that the
$\mathbf{Y}_{\text{bus}}$ method is generally slower than the element-wise
method for large-scale power systems. This finding is based on rigorous
implementation and benchmarks of highly optimized implementations. All
benchmarks performed on recent Intel and Apple platforms support the
observation, using test systems ranging from 14 to 82,000 buses and more.
Also, the reasons why the $\mathbf{Y}_{\text{bus}}$ method is trailing are
discussed. Further, the impact of CPU features and memory performance are
analyzed to assess the scalability promise of the two methods.

This paper is organized as follows: Section II discusses the fundamental
formulations of the $\mathbf{Y}_{\text{bus}}$ method and element-wise method.
Section III briefly introduces CPU features and computer architecture relevant
for result interpretation. Section IV discusses data structures, computation
steps, and implementation techniques for high performance. Section V presents
the benchmark data to compare different CPU generations across different
architectures. Section VI concludes the finding and predicts the scalability
of two methods on future hardware.

##  II Method Formulations

This section discusses the formulations to calculate network power injections
and Jacobians.

###  II-A $\mathbf{Y}_{\text{bus}}$ Method

The complex power injections are calculated by

$\bm{S}=\bm{V}\cdot\left(\mathbf{Y}_{\text{b}us}\bm{V}\right)^{*}$ (1)

where $\bm{S}$ and $\bm{V}$ are two complex-number vectors, respectively, for
power injections and voltage phasors.

The corresponding Jacobian matrix is given by

$\bm{J}=\begin{bmatrix}\bm{J}_{11}&\bm{J}_{12}\\\
\bm{J}_{21}&\bm{J}_{22}\end{bmatrix}=\begin{bmatrix}\real({\partial\bm{S}}/{\partial\bm{\theta}})&\real({\partial\bm{S}}/{\partial\bm{V}})\\\
\imaginary({\partial\bm{S}}/{\partial\bm{\theta}})&\imaginary({\partial\bm{S}}/{\partial\bm{V}})\end{bmatrix}$
(2)

where the derivative sub-matrices are calculated by [9]

$\frac{\partial\bm{S}}{\partial\bm{\theta}}=j\cdot\text{diag}(\bm{V})\left[\text{diag}(\bm{I}_{c}^{*})-(\mathbf{Y}_{\text{bus}}\text{diag}(\bm{V}))^{*}\right]$
(3) $\frac{\partial\bm{S}}{\partial\bm{V}_{m}}=\text{diag}(\bm{I}_{c}^{*}\cdot
e^{j\bm{\theta}})+\text{diag}(\bm{V})(\mathbf{Y}_{\text{bus}}\text{diag}(e^{j\bm{\theta}}))^{*}$
(4)

Although (3) and (4) are compact, they are not of top efficiency as both
require multiple sparse-sparse multiplications. Section III will discuss an
efficient two-pass algorithm to simultaneously compute the two derivative
matrices.

###  II-B Element-wise Method

Consider a transmission line or a two-winding transformer whose off-nominal
tap ratio and phase shift $me^{j\phi}$ is applied at the primary side. The
complex power injections $S$ into the two terminals, $h$ and $k$, are given by

$S_{hk}=v_{h}^{2}\frac{(y_{h}+y_{hk})^{*}}{m^{2}}-v_{h}v_{k}\frac{y_{hk}^{*}}{m}e^{j(\theta_{h}-\theta_{k}-\phi)}$
(5)
$S_{kh}=-v_{h}v_{k}\frac{y_{hk}^{*}}{m}e^{-j(\theta_{h}-\theta_{k}-\phi)}+v_{k}^{2}(y_{k}+y_{hk})^{*}$
(6)

where $y_{h}$ are the shunt admittance (including transformer core losses and
magnetizing reactance) at the $h$ terminal, $y_{k}$ is the shunt admittance at
the $k$ terminal, $y_{hk}$ is the series admittance, $v_{h}$ and $v_{k}$ are
the voltage magnitudes, and $\theta_{h}$ and $\theta_{k}$ are the voltage
phase angles of buses $h$ and $k$.

The power injection into bus $i$ can be obtained by adding up all the
injections from connected lines:

$S_{i}=\sum_{t=1}^{n_{b}}S_{it}\,,(i,t)\in L$ (7)

where $t$ is the other terminal of an existing line that connects buses $i$
and $t$, and $L$ is the set of lines denoted by bus pairs. Each of (5) and (6)
needs to be evaluated for $n_{l}$ times, and (7) requires up to $2n_{l}$
additions.

To obtain the sparse Jacobian matrix, one first calculates the non-zero
Jacobian entries and then assemble them into a sparse matrix. The four power
injection equations ($P_{h}$, $Q_{h}$, $P_{k}$ and $Q_{k}$) can be obtained
from (5) and (6). Taking the partial derivative of the four equations with
respect to four variables ($v_{h}$, $\theta_{h}$, $v_{k}$, and $\theta_{k}$)
yields 16 expressions. All 16 expressions need to be evaluated for all $n_{l}$
lines to obtain the complete Jacobian elements. Given a large number of
Jacobian expressions and the number of lines, this method is perceived as
having high computational complexity.

##  III Modern CPU Features

This section will walk through the features of modern CPUs to provide a
fundamental understanding of comparing the two methods. These features include
a) cache and prefetch, b) single-instruction multiple-data (SIMD)
vectorization, c) pipelining and branch prediction. All these features serve
one goal: to reduce the overall latency.

Modern CPUs have multiple levels of on-chip cache. These caches are made of
static random-access memory (SRAM) so is faster to access than memory, which
is made of dynamic random-access memory (DRAM). For example, an Intel Core
i9-12900K has three levels of cache, namely, L1, L2, and L3, with respective
sizes of 80 KB (per core), 1.25 MB (per core), and 30 MB (total). The time to
access the three levels of cache is around 1 ns, 4 ns, and 40 ns, compared
with a typical memory access time of 80 ns. If the requested data is in the
cache, the access is a “cache hit”. Otherwise, it is a “cache miss” and will
incur time to access the memory.

To utilize the on-chip cache and improve the cache hit rate, “prefetch” is
used to preemptively load data from memory to cache. Prefetch is accomplished
by both the CPU hardware and the compiler, which may insert prefetch
instructions into the machine code. Generally, data prefetch for an algorithm
will have a higher success rate if the algorithm access data in a more regular
manner, such as visiting a memory block linearly. Conversely, if the data
locations can only be determined at run time, prefetch success will likely
decrease.

Figure 1: Illustration of pipelining six instructions in a five-stage
pipeline.

To reduce the number of instruction cycles for the same computations, a
single-instruction multiple-data (SIMD) instruction set is introduced. SIMD
instructions in the x86 platform include the Advanced Vector eXtension (AVX)
instructions, which can perform simultaneous arithmetic operations on four to
eight double-precision floating points, using AVX2 and AVX512, respectively.
To utilize SIMD, data operands need to be in the cache and ready to be loaded
into vector registers. On the other hand, using SIMD will add overhead to
calculations that are done on scalars with mixed memory access. In such cases,
compilers will avoid generating SIMD instructions but instead use scalar
operations.

Another relevant technique is instruction pipelining. Multiple stages are
needed to execute an instruction, including fetch, decode, execute, memory
access, and register write-back. Each stage is performed by a part of the CPU.
Pipelining aims to keep all parts busy while executing multiple instructions.
Figure 1 illustrates the execution of six instructions in a five-stage
pipeline. The first instruction takes T1 to T5 to complete. Starting from T6,
each instruction will only take one time unit.

Although pipelines are efficient once warmed up, there is a downside. If one
instruction needs to be changed within the pipeline, the pipeline will stall
and needs to be flushed. Programs with if-else statements, for example, will
likely introduce instruction changes due to jumping. Modern x86 processors
have long pipelines, so for remediation, a branch prediction mechanism is
introduced. If the prediction is correct, the pipeline can continue, and only
a few cycles will be wasted. Otherwise, the time penalty will be much greater.

In summary, improving the cache hit rate, utilizing SIMD as much as possible,
and reducing instruction changes will improve the performance of a program. As
alluded to before, algorithms with a regular data storage and access pattern
will improve prefetch and improve the cache hit rate. Also, programs with
fewer conditional jumps will have better pipeline performance. These metrics
will be utilized to qualitatively assess the efficiency of the
$\mathbf{Y}_{\text{bus}}$ and element-wise method.

##  IV Data Structure, Steps, and Implementations

This section presents sufficient details on the data structure, computation
steps, and implementations that are crucial for a rigorous assessment of the
two methods.

###  IV-A Data Structure

####  IV-A1 $\mathbf{Y}_{\text{bus}}$ method

The prominent data structure for storing $\mathbf{Y}_{\text{bus}}$ is the
sparse matrix technique, which is comprehensively covered in the literature
[10, 11]. The most common sparse formats to store $\mathbf{Y}_{\text{bus}}$
are Compressed Sparse Column (CSC) and Compressed Sparse Row (CSR), and a
trade-off exists between them. The CSR format is faster than CSC for sparse
matrix-vector product operation but slower for factorization. If
$\mathbf{Y}_{\text{bus}}$ is in CSR, the resultant $\bm{J}$ will be in CSR by
default. Then, $\bm{J}$ will need to be converted to CSC every time to
interface with fast solvers like KLU [12]. As a result, it is more common to
store $\mathbf{Y}_{\text{bus}}$ and $\bm{J}$ in the CSC format.

####  IV-A2 Element-wise method

The element-wise method does not require special data structures besides the
array. To facilitate SIMD and pipelining, the same parameter or variable of
all device instances is stored in a one-dimensional array. For example, an
array $\bm{v}_{h}$ will hold the sending-end voltage magnitudes of all lines.
The data structure guarantees that all elements of $\bm{v}_{h}$ are in
contiguous memory.

The contiguous-memory requirement is also applied to the arrays of complex
numbers. That is, the real and imaginary parts of complex numbers are stored
in two separate arrays with an element type of floating point numbers. Also,
one can separate the real power from reactive power in (5) and (6) so that no
complex-valued arrays are needed.

###  IV-B $\mathbf{Y}_{\text{bus}}$ Method: Computation Steps and Algorithms

####  IV-B1 Power injection computation

The following steps are performed in order to calculate the power injections:

  1. 1.

Update the complex arrays $\bm{U}=e^{j\bm{\theta}}$ and
$\bm{V}=\bm{V}_{m}\cdot\bm{U}$

  2. 2.

Calculate sparse matrix-vector product
$\bm{I}_{\text{b}us}=\bm{Y}_{\text{b}us}\bm{V}$

  3. 3.

Calculate two vector product $\bm{V}\cdot\bm{I}^{*}_{\text{b}us}$

where $\bm{\theta}$ is the array of voltage phase angles, $\bm{V}_{m}$ is the
array of voltage magnitudes, and $\bm{I}_{bus}$ is the array of current
injections.

Steps 2 and 3 favor different data structures for $\bm{V}$ for performance.
For Step 2, the access into $\bm{V}$ is random, which depends on the current
column index of the non-zero element in $\mathbf{Y}_{\text{bus}}$. Placing
together the real and imaginary parts of each complex element in $\bm{V}$ will
enable the retrieval of two values in one access. For Step 3, storing the real
and imaginary parts of $\bm{V}$ and $\mathbf{I}_{\text{bus}}$ separately and
linearly will enable SIMD support with improved pipeline efficiency. This
conflict is fundamental and can only be alleviated by copying memory, which
will incur overhead.

####  IV-B2 Jacobian calculation

A two-pass algorithm [13] is adopted to compute
$\partial\bm{S}/\partial\bm{\theta}$ and $\partial\bm{S}/\partial\bm{V_{m}}$
simultaneously and reuse intermediate values. This method observes the
sparsity patterns of $\bm{J}$ and $\mathbf{Y}_{\text{bus}}$. That is,
calculations in (3) and (4) only involves the product of
$\mathbf{Y}_{\text{bus}}$ by diagonal matrices and the addition (or
subtraction) of diagonal matrices. Since these operations will not introduce
non-zero elements at locations where the corresponding
$\mathbf{Y}_{\text{bus}}$ element is zero, the Jacobians
$d\bm{S}/d\bm{\theta}$ and $d\bm{S}/d\bm{V}_{m}$ will have the same sparsity
pattern as $\mathbf{Y}_{\text{bus}}$. More precisely, the locations of
nonzeros in the Jacobians are subsets of the nonzeros in
$\mathbf{Y}_{\text{bus}}$ because some values may become zeros. Nonetheless,
these small numbers of zeros can be stored to reuse the sparsity pattern of
$\mathbf{Y}_{\text{bus}}$.

Algorithm 1 adopted from [14] presents the pseudo-code for the two-pass method
for CSC sparse matrix. In the inputs, $\bm{Y}_{p}$, $\bm{Y}_{i}$ and
$\bm{Y}_{v}$ are the column pointers, row indices, and nonzeros of the CSC
storage of $\mathbf{Y}_{\text{bus}}$, $\mathbf{I}_{\text{bus}}$ is the array
for bus current injections. In the initialization phase, copying the non-zero
values of $\mathbf{Y}_{\text{bus}}$ to that of $d\bm{S}/d\bm{\theta}$ and
$d\bm{S}/d\bm{V}_{m}$ only involve linear memory access and is thus rapid. The
first pass computes the $\mathbf{Y}_{\text{bus}}$-vector product in (3) and
(4), and the second pass computes the remainder.

Algorithm 1 Two-pass method for Jacobian using CSC

1: Input: $\bm{Y}_{p}$, $\bm{Y}_{i}$, $\bm{Y}_{v}$, $d\bm{S}_{\theta}$,
$d\bm{S}_{V}$, $\mathbf{I}_{\text{bus}}$, $\bm{V}$, $\bm{U}$

2: Initialize:

3: $\bm{I}_{bus}=0$, $d\bm{S}_{\theta}.{\text{nzval}}\mathrel{.}=\bm{Y}_{v}$,
$d\bm{S}_{V}.{\text{nzval}}\mathrel{.}=\bm{Y}_{v}$

4: Pass 1: for $j=1:n_{b}$ ($j$ is column index)

5: for $k=\bm{Y}_{p}[j]:(\bm{Y}_{p}[j+1]-1)$ ($k$ is $j$’s index range)

6: $\bm{I}_{bus}[\bm{Y}_{i}[k]]\mathrel{+}=\bm{Y}_{v}[k]*\bm{V}[j]$

7: $d\bm{S}_{\theta}.\text{nzval}[k]\mathrel{\times}=\bm{V}[j]$

8: $d\bm{S}_{V}.\text{nzval}[k]\mathrel{\times}=\bm{U}[j]$

9: end for $k$

10: end for $j$

11: Pass 2: for $j=1:n_{b}$

12: for $k=\bm{Y}_{p}[j]:(\bm{Y}_{p}[j+1]-1)$

13: $i=\bm{Y}_{i}[k]$ ($i$ is element $k$’s row number)

14: $d\bm{S}_{V}.\text{nzval}[k]=\bm{V}[i]\times(d\bm{S}_{V}[K])^{*}$

15: if $i==j$

16: $d\bm{S}_{\theta}.\text{nzval}[k]\mathrel{-}=\bm{I}_{\text{b}us}[j]$

17:
$d\bm{S}_{V}.\text{nzval}[k]\mathrel{+}=(\bm{I}_{bus}[j])^{*}\times\bm{U}[j]$

18: end if

19:
$d\bm{S}_{\theta}.\text{nzval}[k]=(1im)\times(d\bm{S}_{\theta}.\text{nzval}[k])^{*}\times\bm{V}[i]$

20: end for $k$

21: end for $j$

22: Return: $d\bm{S}_{\theta}$, $d\bm{S}_{V}$

The two-pass method is more efficient than directly performing sparse matrix
multiplications because

  1. 1.

Memory allocation for sparse matrices is avoided due to the reuse of sparsity
pattern.

  2. 2.

Memory access to the non-zero values of $d\bm{S}/d\bm{\theta}$ and
$d\bm{S}/d\bm{V}_{m}$ are linear, thus data are more likely to have been
prefetched to cache.

Also, note that a slightly more efficient version can be obtained for CSR-
stored $\mathbf{Y}_{\text{bus}}$. Since many sparse linear solvers only
support the CSC format, a CSR-to-CSC conversion will be required at every
step. The overhead of the conversion may out-weight the gain from using the
CSR format.

####  IV-B3 Jacobian matrix assembly

The $\bm{J}$ matrix with real values can be obtained by concatenation as given
in (2), but sparse matrix concatenation is slow due to index sorting and
memory allocation every time two sub-matrices are concatenated.

Algorithm 2 is proposed to facilitate rapid assembly of the real and imaginary
parts of the complex Jacobians into $\bm{J}$, which is of the size
$2n_{b}\times 2n_{b}$. With pre-allocated sparsity pattern for $\bm{J}$, the
assembly operation is essentially copying the derivative elements to proper
locations in $\bm{J}$. Since the sparsity patterns do not change, one can pre-
calculate the indexing arrays $\bm{P}_{xy}$ into $\bm{J}$’s non-zeros for
corresponding sub-matrices.

Algorithm 2 is generalized for cases where all non-zeros of $\bm{J}_{xy}$
exist in $\bm{J}$ but not necessarily the other way. For a basic power flow
calculation, each sub-matrices will have identical sparsity pattern to
$\mathbf{Y}_{\text{bus}}$. This may not be the case when series devices other
than lines and transformers (such as TCSC) exist. Algorithm 2 is thus designed
to obtain the indices based on matching row and column coordinates.

Algorithm 2 Obtain indices to write to sparsity pattern

1: Input: $\bm{J}_{p}$, $\bm{J}_{i}$, $\bm{Y}_{\text{row}}$,
$\bm{Y}_{\text{col}}$, $\bm{P}_{xy}$ where $x,y\in(1,2)$

2: Initialize:

3: $\bm{P}_{xy}\mathrel{.}=0$ where $x,y\in(1,2)$

4: Loop for $P_{xy}$: for $i=1:\text{length}(\bm{Y}_{\text{col}})$

5: $c_{ol}=Y_{col}[i]+(y-1)*n_{b}$

6: $r_{ow}=Y_{row}[i]+(x-1)*n_{b}$

7: for $k=\bm{J}_{p}[c_{ol}]:(\bm{J}_{p}[c_{ol}+1]-1)$

8: if $J_{i}[k]\mathrel{==}r_{ow}$

9: $P_{xy}[i]=k$; break;

10: end if

11: end for $k$

12: end for $i$

13: Run the loop with $(x,y)=(1,1),(1,2),(2,1),\text{and}~{}(2,2)$

14: Return: $\bm{P}_{11}$, $\bm{P}_{12}$, $\bm{P}_{21}$, $\bm{P}_{22}$

Using the pre-calculated indices, the Jacobian submatrices can be rapidly
copied by J.nzval[$P_{xy}$] .= $J_{xy}$.nzval. This method is efficient
because

  1. 1.

no random access to the sparse matrix $\bm{J}$ is needed. Time is thus saved
by avoiding index look-ups,

  2. 2.

the memory copy operations visit the sub-matrix non-zeros linearly, which
improves cache performance.

###  IV-C Element-wise Method Computation Steps

####  IV-C1 Vectorized calculation

The element-wise method does not require sophisticated algorithms for
calculation. The principle of the element-wise method is to separate CPU-bound
vectorized evaluations from memory-bound operations and maximize SIMD
utilization for the former. Since all parameters and variables are stored in
arrays of an equal length, the evaluation of all expressions for the same
device can be packed in the same loop, in which the same index is used to
access the arrays. An illustration of the calculation loop is given in
LABEL:lst:g-elemwise-loop,

⬇

1for m in 1:n_line

2 Ph[m] = vh[m]*(vh[m]*itap2ghghk[m] - ...)

3 Pk[m] = vk[m]*(vk[m]*ghk[m] - ...) ...

4end

Listing 1: Illustration of element-wise power injection calculation. Part of
the expressions is omitted for brevity.

where all constants, such as itap2ghghk=$(g_{h}+g_{hk})/T^{2}_{ap}$, are
calculated once ahead of time.

There are three benefits of element-wise evaluations:

  1. 1.

the same parameters and variables are used in multiple expressions. As a
result, cached values may be reused.

  2. 2.

since the expressions across devices are homogeneous, the compiler can
generate SIMD instructions to simultaneously evaluate one expression for
multiple instances.

  3. 3.

no branch exists in the calculation code, so the execution pipeline can be
sustained.

After calculating the line-wise power injections, the bus-wise power
injections are calculated using (7). This step is termed “reduction” as it
reduces the number of injections from $n_{l}$ to $n_{b}$. Given that multiple
lines can connect to the same bus, each reduction involves three operations:
load, add, and store. Also, for a given order of lines, their connected buses
can be random, thus random access is required for the reduction.

####  IV-C2 Jacobian matrix assembly

There are multiple ways to assemble the sparse Jacobian matrix from elements.
This step is also a “reduction” because multiple non-zero values for the same
location need to be reduced into one by summation. The quickest implementation
is to convert the coordinates-value triplets to a new sparse matrix. Due to
memory allocation and sorting, this method is usually the slowest.

Alternatively, one can pre-build the sparsity pattern matrix as in the
$\mathbf{Y}_{\text{bus}}$ method and copy-add elements to the non-zero array
of the sparsity pattern matrix. Figure 2 illustrates this process, where each
element in $\bm{J}_{\text{elem}}.{\text{nzval}}$ require four operations for
an element in $\bm{J}_{\text{shape}}.{\text{nzval}}$. These four operations
include two loads, one add, and one store.

Figure 2: Copy-adding non-zero elements to the sparsity pattern storage.
Figure 3: Two-step technique with separate copying and copy-adding.

We propose a two-step implementation technique to separate the initial copy-
only operations from subsequent copy-add operations. The two steps are
illustrated in Figure 3. The first step will copy elements from
$\bm{J}_{\text{elem}}.{\text{nzval}}$ only when the destinations in
$\bm{J}_{\text{shape}}.{\text{nzval}}$ are visited for the first time. It
reduces the operations in the first step to one load and one store. Subsequent
copy-add operations will cover the remaining elements, which will require four
operations. This method eliminates one load and one add operation for
length($\bm{J}_{\text{shape}}.{\text{nzval}}$) times. Given the reduced memory
access, the total time will be considerably reduced. Still, the reduction step
has no parallelism and can be more time-consuming than the equation
evaluation.

##  V Case Studies

This section presents comprehensive benchmark studies on the two methods for a
variety of power grid test cases, ranging from 14 to 82,000 buses and more
[15]. Micro-benchmarks are first shown to demonstrate the performance gains by
the techniques described in Section IV. Next, the $\mathbf{Y}_{\text{bus}}$
and the element-wise methods are benchmarked for grid test cases to identify
the top performer, analyze the time breakdown, and understand the limitations.
Further, multiple x86 and ARM-based computer systems are utilized to
understand the trend of scalability as computer systems evolve.

All the discussed algorithms are implemented and optimized in Julia. As a
compiled language, Julia provides superior performance for code that is
optimally implemented. Generated machine code is inspected to guarantee SIMD
vectorization whenever possible. All floating point numbers are stored in the
64-bit, IEEE 754-compliant format. For complex numbers, the real and imaginary
parts are stored using 64 bits.

###  V-A Micro-benchmarks of Computation Steps

This section uses the Synthetic USA system with 82,000 buses on a computer
with Intel Core i9-12900K and 3600 MHz DDR4 RAM. The tested 12900K CPU only
supports AVX2, which is a 256-bit wide SIMD instruction.

####  V-A1 $\mathbf{Y}_{\text{bus}}$ method

The time to compute power injections is first measured. Specifically, the
power injections computation contains two steps: a) phasor calculations and b)
sparse matrix and vector calculations, which include steps 2 and 3 in Section
IV.B. Cases are compared with the real and imaginary parts of complex phasors
stored interleaved or separately. As shown in Table I, separate storage
reduces the time to update the voltage vectors by using SIMD, although extra
time is required to change the data structure before the power calculations.
Note that the best time performance is given in the bold text in all
subsequent tables.

TABLE I: P/Q calculation time in $\mu s$ for $\mathbf{Y}_{\text{bus}}$ with different storage |  $\bm{V}$ and $\bm{U}$ Calc. | Power Calc. | Total Time
---|---|---|---
Interleave | 557 | 722 | 1279
Separate | 137 | 806 | 943

The Jacobian computation has two steps: a) $d\bm{S}/d\bm{\theta}$ and
$d\bm{S}/d\bm{V}_{m}$ calculation, and b) assembly into $\bm{J}$. For Step a),
matrix multiplication and the two-pass method are compared. For Step b), the
concatenation and in-place methods are compared. It can be seen that the two-
pass method is several times faster than the plain matrix multiplication, and
the in-place copy approach is nearly 25 times faster than concatenation. In
following comparisons with the element-wise method, the fastest methods will
be used.

TABLE II: Jacobian calculation and assembly time in $\mu s$ for $\mathbf{Y}_{\text{bus}}$ method Derivatives Calculation | Matrix Assembly
---|---
Matrix Multiplication | Two-Pass | Concatenation | In-Place Copy
13,263 | 2,733 | 15,350 | 602

####  V-A2 Element-wise method

Both the power injection and Jacobian calculations in the element-wise method
has two steps: evaluation and reduction. For power injection calculation, with
and without SIMD are tested. Table III shows that with the array-based data
structure, the SIMD-enabled case only takes about 1/4 of the time without
SIMD. This is expected because the processor supports AVX2, which packs four
64-bit floating-point calculations in one instruction.

TABLE III: Element-wise power injection computation time in $\mu s$ Evaluation Methods | Reduction
---|---
No SIMD | SIMD |
1,179 | 299 | 183

For the Jacobian computation, the evaluation time and the three reduction
methods are compared in Table IV. For the reduction part, creating a new
matrix is much more expensive than the copy-add and the two-step method. Also,
calculating all 1,665,936 Jacobian elements takes merely 848 $\mu s$, compared
with 1,830 $\mu s$ for the two-step reduction method. Therefore, the
bottleneck for Jacobian computation in the element-wise method is in the
reduction.

TABLE IV: Element-wise Jacobian computation time $\mu s$ Vectorized Evaluation | Reduction Methods
---|---
| New matrix | Copy-add | Two-step
848 | 13,873 | 2,025 | 1,830

####  V-A3 Comparing $\mathbf{Y}_{\text{bus}}$ and element-wise methods

The best computation times for both methods are given in Table V. It can be
seen that the $\mathbf{Y}_{\text{bus}}$ method is slower than the element-wise
method for both power injection and Jacobian calculations. Since the
theoretical computational complexity is in the same order of magnitude (with
the $\mathbf{Y}_{\text{bus}}$ method having slightly less calculations), the
time difference needs to be investigated at the execution level.

TABLE V: Comparison of best computation time in $\mu s$ |  $\mathbf{Y}_{\text{bus}}$ Method | Element-wise Method
---|---|---
Power Computation | 943 | 482
Jacobian Computation | 3,335 | 2,678

We use the CPU statistics by the Linux perf monitoring tool to compare the
element-wise Jacobian evaluation and the two-pass method. Note that the
element-wise method has more complexity on paper due to the number of
elements, but the two-pass method has multiple jumps. The two functions are
executed for 1,000 runs, and the CPU statistics are shown in Table VI. It can
be seen that the $\mathbf{Y}_{\text{bus}}$ method requires more cache access
and has nearly twice the cache miss rate than the element-wise method. Also,
the $\mathbf{Y}_{\text{bus}}$ method has significantly more branches and much
higher branch mispredictions. Therefore, the total number of instructions for
the $\mathbf{Y}_{\text{bus}}$ method is higher, resulting in more cycles and a
longer execution time.

TABLE VI: CPU statistics for Jacobian calculation |  $\mathbf{Y}_{\text{bus}}$ J Calculation | Element-wise J Calculation
---|---|---
Cycles | 14,327,078,176 | 4,577,698,819
Cache Access | 828,580,410 | 581,952,974
Cache Misses | 156,647,072 | 60,248,850
Branches | 5,474,516,729 | 19,537,258
Branch Mispredict. | 173,767,459 | 2,538
Instructions | 50,308,631,851 | 3,690,449,376

###  V-B Benchmark for Multiple Test Cases

This section continues to use Intel 12900K for benchmarking the two methods
for nine power grid test cases. Results are shown in Table VII. For the
element-wise method,  a) the Jacobian reduction is the bottleneck of all its
computations.  b) the power equation evaluation scales nearly linearly to
system size.  For the $\mathbf{Y}_{\text{bus}}$ method, the Jacobian
calculation is the bottleneck.

Several observations are made by comparison:

  1. 1.

For small cases such as the 14-bus and the 118-bus systems, the
$\mathbf{Y}_{\text{bus}}$ method is faster in both power injection and
Jacobian calculations. This is because all data are likely in the L1 or L2
cache, and the computation demand is insufficient to sustain the pipeline.

  2. 2.

The effect of pipelining can be seen by comparing the power calculation time
for the 118-bus and the 300-bus systems. The two systems are of modest sizes
so data can mostly fit into L1 and L2 caches. As the number of buses increases
by 2.5x, the calculation time increases by a factor of 2.6 for the
$\mathbf{Y}_{\text{bus}}$ method. For the element-wise method, the calculation
time only increases by 2.1x due to improved pipelining for SIMD.

  3. 3.

The power computation is faster using the element-wise method than the
$\mathbf{Y}_{\text{bus}}$ method for cases with greater than 300 buses. This
is possibly due to cache misses in the $\mathbf{Y}_{\text{bus}}$ method when
performing sparse matrix-vector multiplication.

  4. 4.

The Jacobian computation is faster using the element-wise method for systems
of any size. This is mostly due to the branching and jumping that are
unfriendly to pipelining, as discussed in Section V.A.3.

Overall, the element-wise method is faster than $\mathbf{Y}_{\text{bus}}$
method for both power injection and Jacobian computations for systems ranging
from 300 to 82k buses.

TABLE VII: Computation time for the $\mathbf{Y}_{\text{bus}}$ and the element-
wise methods for systems of different sizes.
Time is given in nanoseconds for rows 1 and 2 and in microseconds for the remainder. |  $\mathbf{Y}_{\text{bus}}$ Method | Element-wise Method
---|---|---
Case | Phasor Calc. | Power Calc. | J Calc. | J Red. | Power Calc. | Power Red. | J Calc. | J Red.
case14 (ns) | 29.86 | 93.21 | 333.37 | 88.93 | 59.84 | 105.69 | 65.46 | 87.41
case118 (ns) | 207.65 | 519.15 | 2,819.56 | 565.71 | 537.21 | 300.17 | 1,044.30 | 869.42
case300 | 0.51 | 1.37 | 7.19 | 1.64 | 1.18 | 0.54 | 2.58 | 2.04
case1354pegase | 2.28 | 5.94 | 31.82 | 8.05 | 5.66 | 2.65 | 12.51 | 13.29
case2736sp | 4.75 | 12.42 | 65.87 | 16.70 | 9.92 | 4.56 | 15.49 | 23.61
case9241pegase | 15.49 | 79.08 | 341.35 | 75.46 | 45.89 | 26.72 | 124.80 | 248.02
case_ACTIVSg25k | 41.77 | 214.54 | 818.18 | 176.18 | 92.62 | 49.70 | 230.34 | 491.11
case_ACTIVSg70k | 117.00 | 665.77 | 2,291.57 | 504.21 | 255.39 | 142.50 | 686.66 | 1,536.49
case_SyntheticUSA | 137.09 | 806.48 | 2,733.15 | 602.55 | 299.40 | 183.85 | 847.68 | 1,829.67

###  V-C Computational Complexity Due to Grid Size

While the element-wise method leads in performance, it is observed that the
power and Jacobian calculations for the two methods slow down differently.
That is, as the system size increases, the average computational complexity
increases by more for some operations than others. We can measure a rough
metric coined as the “computation time per bus“ for each step. As the grid
size increases, cache misses and memory access increase, thus a longer
computation time per bus is expected.

To characterize the relative computation performance across systems, we
calculate the computation-time-per-bus ratios of all systems to that of the
118-bus system. This ratio is termed the “time-cost ratio“. A time-cost ratio
smaller than one means that an operation takes less time per bus to perform,
as compared to the 118-bus system. Ratios greater than 1 correspond to super-
linear complexity (namely, unfavorable scalability), indicating worsening
performance as system size increases. The larger the time-cost ratio, the
worse the scalability is.

Results are plotted in Figure 4, which shows that the
$\mathbf{Y}_{\text{bus}}$ method has linear scalability up to the 2736-bus
system because data are likely in the CPU caches. The element-wise method for
power calculation has sub-linear scaling and is thus optimal compared with
$\mathbf{Y}_{\text{bus}}$. Further, the Jacobian calculation in the element-
wise method has worse scalability than the $\mathbf{Y}_{\text{bus}}$ method.
If this trend sustains, the $\mathbf{Y}_{\text{bus}}$ method can be faster to
compute the Jacobian than the element-wise method.

Figure 4: Time complexity ratio relative to the 118-bus system on 12900K

###  V-D Impact of SIMD Width

As noted, the tested 12900K CPU only supports AVX2, which is 256-bit wide. It
is anticipated that CPUs with AVX512 support will be faster for the element-
wise evaluations. But is it the case across the board?

We benchmark the nine cases on an Intel 11900K CPU with 32GB DDR4 memory
clocking at 3200 MHz. Compared with 12900K, this CPU is one generation older
in architecture, runs at a lower frequency, and has smaller caches. Also, the
memory is slower than that used in the 12900K system. On the contrary, this
CPU supports 512-bit wide SIMD through AVX512. Comparing the results shown in
Table VIII with Table VII, the following observations are made:

  1. 1.

Phasor evaluation in the $\mathbf{Y}_{\text{bus}}$ method and the power
evaluation in the element-wise method is faster for all systems due to AVX512
support. Given the system size, the operands can fit into the L2 cache. These
evaluations are thus CPU-bound.

  2. 2.

In the $\mathbf{Y}_{\text{bus}}$ method, the power evaluation cannot utilize
AVX512. Its computation performance is thus slower and proportional to the CPU
clock frequency.

  3. 3.

In the element-wise method, the Jacobian evaluation is faster for small
systems and slower for larger ones. The slowdown is due to memory access for
large systems. Therefore, this evaluation switches from CPU-bound to memory-
bound as system size increases.

It can be concluded that AVX512 will reduce the computation time for
vectorized evaluations by a modest percentage, and it is most conducive to
small and medium-sized systems. For large systems, the whole computation
becomes memory-bound for both methods.

TABLE VIII: Benchmark on Intel i9-11900K with AVX512 | Ybus Method | Element-wise Method
---|---|---
Case | Phasor Calc. | Power Calc. | J Calc. | J Red. | Power Calc. | Power Red. | J Calc. | J Red.
case14 (ns) | 17.46 | 123.93 | 463.99 | 111.28 | 56.24 | 145.56 | 71.34 | 90.59
case118 (ns) | 122.71 | 605.35 | 3,911.88 | 819.61 | 400.09 | 443.47 | 824.30 | 1,027.40
case300 | 0.30 | 1.66 | 10.06 | 2.20 | 0.72 | 0.78 | 1.93 | 2.43
case1354pegase | 1.38 | 6.89 | 43.85 | 9.63 | 3.39 | 3.57 | 10.86 | 19.00
case2736sp | 2.72 | 14.90 | 91.93 | 21.79 | 6.18 | 6.31 | 18.64 | 32.99
case9241pegase | 9.41 | 88.23 | 430.85 | 89.32 | 29.63 | 34.76 | 97.57 | 255.11
case_ACTIVSg25k | 25.40 | 249.08 | 1,053.98 | 205.37 | 58.59 | 65.77 | 195.06 | 394.81
case_ACTIVSg70k | 69.62 | 757.26 | 3,274.72 | 712.65 | 168.35 | 185.77 | 755.40 | 2,041.78
case_SyntheticUSA | 81.58 | 969.50 | 3,956.52 | 875.01 | 205.78 | 218.37 | 959.30 | 2,444.29
TABLE IX: Benchmark on Apple M2 | Ybus Method | Element-wise Method
---|---|---
Case | Phasor Calc. | Power Calc. | J Calc. | J Red. | Power Calc. | Power Red. | J Calc. | J Red.
case14 (ns) | 42.17 | 131.38 | 360.44 | 100.81 | 92.32 | 184.30 | 93.86 | 140.80
case118 (ns) | 281.18 | 666.93 | 2,916.67 | 761.69 | 823.96 | 529.02 | 688.47 | 1,366.60
case300 | 0.69 | 1.49 | 7.33 | 1.86 | 1.95 | 0.98 | 1.75 | 3.47
case1354pegase | 3.12 | 7.38 | 34.88 | 8.58 | 9.33 | 5.11 | 14.29 | 23.33
case2736sp | 6.29 | 14.71 | 77.46 | 17.67 | 16.42 | 7.79 | 18.17 | 32.75
case9241pegase | 21.13 | 88.42 | 355.29 | 69.13 | 95.33 | 48.83 | 124.08 | 222.04
case_ACTIVSg25k | 57.13 | 227.29 | 883.04 | 157.33 | 190.25 | 89.17 | 229.67 | 332.88
case_ACTIVSg70k | 160.08 | 645.92 | 2,592.08 | 452.63 | 521.21 | 250.63 | 1,025.92 | 962.33
case_SyntheticUSA | 187.79 | 769.29 | 2,821.54 | 524.63 | 619.17 | 296.88 | 1,351.75 | 1,194.46

###  V-E Impact of Improved Memory Speed

At this point, it is understood that memory access is the bottleneck for
large-scale system calculations. Solutions from computer systems are to clock
up the memory and physically move the processor and memory closer. A notable
system is the Apple Silicon, which places the RAM beside the processor. This
design has the benefit of reduced memory latency, which may lead to
performance gains for large-scale systems. Understanding how this design
impacts the performance of the two methods will help us to predict their
scalability.

Table IX shows the benchmark results on an Apple M2 system with 24 GB LPDDR5
RAM running at 6400 MHz. Apple M2 has support for 128-bit SIMD, 192 KB L1
cache , and 16 MB L2 cache. Comparing Table IX with Table VII provides
interesting observations:

  1. 1.

The phasor calculation using $\mathbf{Y}_{\text{bus}}$ and the power and
Jacobian evaluation using the element-wise method on Apple M2 is consistently
within 2x the time of Intel 12900K. This is because these computations are
CPU-bound, in which the SIMD width plays a central role.

  2. 2.

The power evaluation in $\mathbf{Y}_{\text{bus}}$ is initially slower and
later faster on M2 than on 12900K. This is because, for small systems, the
larger cache in 12900K can help reduce cache misses. For large systems, the
problem will become memory-bound, thus M2 becomes faster.

  3. 3.

The Jacobian calculation in $\mathbf{Y}_{\text{bus}}$ is slower on M2 than on
12900K by a stable factor. This step cannot utilize SIMD and is impacted by
both the CPU frequency and memory speed.

  4. 4.

The Jacobian reduction is initially CPU-bound but quickly becomes memory
bound. Faster memory will bring major improvements to this step.

Figure 5: Computation time ratio relative to the 118-bus system on Apple M2

In addition, Figure 5 shows the time-cost ratio relative to the 118-bus system
on M2. It can be seen that for the 82k-bus Synthetic system, the element-wise
method scales slightly worse than the $\mathbf{Y}_{\text{bus}}$ method for
Jacobian calculation. This observation matches that on the Intel platform,
although the gap becomes smaller since faster memory is more helpful to the
element-wise method.

TABLE X: Computation time in microseconds using ultra-large grid test cases on Intel 12900K | Ybus Method | Element-wise Method
---|---|---
Case | Phasor Calc. | Power Calc. | J Calc. | J Red. | Power Calc. | Power Red. | J Calc. | J Red.
164k | 275 | 1,721 | 5,806 | 1,478 | 608 | 373 | 2,084 | 4,192
246k | 412 | 2,796 | 9,043 | 2,580 | 937 | 583 | 3,344 | 7,369
328k | 549 | 4,060 | 12,511 | 3,938 | 1,346 | 775 | 4,577 | 11,067
410k | 689 | 5,208 | 15,948 | 5,300 | 1,915 | 965 | 6,008 | 16,740
TABLE XI: Computation time in microseconds using ultra-large grid test cases on Apple M2 | Ybus Method | Element-wise Method
---|---|---
Case | Phasor Calc. | Power Calc. | J Calc. | J Red. | Power Calc. | Power Red. | J Calc. | J Red.
164k | 367 | 1,698 | 5,911 | 1,158 | 1,256 | 594 | 2,906 | 3,442
246k | 566 | 2,577 | 8,883 | 1,776 | 1,930 | 912 | 4,771 | 6,414
328k | 757 | 3,435 | 12,079 | 2,399 | 2,619 | 1,226 | 6,953 | 9,590
410k | 950 | 4,350 | 15,096 | 3,006 | 3,226 | 1,538 | 7,982 | 12,798

###  V-F Outlook for Ultra-Large Systems

The tested cases with up to 82k buses declare the element-wise method as a
winner for large-scale systems. Still, worse scalability is observed for the
Jacobian computation in the element-wise method. This indicates that, as the
system size further increases, the $\mathbf{Y}_{\text{bus}}$ method may take
over. This section investigates how the two methods will continue to scale for
ultra-large cases to understand the performance turning points.

The 82k-bus Synthetic USA system is repeated multiple times to create four
ultra-large systems. All buses and lines are duplicated by two to five times,
respectively, resulting in four cases. The performance of the two methods on
12900K and Apple M2 are reported in Table X and Table XI. It can be seen that
the Jacobian reduction step in the element-wise method continues to scale
poorly. On the Intel platform, the $\mathbf{Y}_{\text{bus}}$ method and the
element-wise methods will be similarly slow when the number of buses reaches
about 328k. On the M2 platform, the tie happens around 246k buses.

Regarding the question that which of the two methods will be faster on future
systems, we argue that the element-wise method will continue to outperform the
$\mathbf{Y}_{\text{bus}}$ method for large-scale systems ranging from 1k to
410k buses or more. If we take the best columns from Table X and Table XI, the
element-wise method is a clear winner for power calculation. For the Jacobian
calculation, the best columns of the two methods will take roughly the same
time for the 410k bus system. Nevertheless, faster memory and AVX2/AVX512 will
moderately benefit the element-wise Jacobian calculation but not as much for
the $\mathbf{Y}_{\text{bus}}$ method. In addition, the Jacobian reduction
method presented in Section IV still has room for optimization, but that is
beyond the scope of this paper.

##  VI Conclusions

In this paper, we investigated the computational performance of the
$\mathbf{Y}_{\text{bus}}$ method for power injection calculation and sparse
Jacobian formulation. The main findings are as follows:

  1. 1.

The major issue with the $\mathbf{Y}_{\text{bus}}$ method is the sparse
matrix-centered computations, which are not capable of utilizing SIMD and
incurs many jumps due to indexing.

  2. 2.

Compared with the element-wise method, the $\mathbf{Y}_{\text{bus}}$ method is
slower for practical systems with bus numbers ranging from thousands to a few
hundred thousand on the x86 CPUs.

  3. 3.

The bottleneck for the element-wise method is reducing the Jacobian elements
into a sparse matrix. The step scales poorly as the number of buses enters the
realm of hundreds of thousands. This step is primarily constrained by memory
speed.

We conclude that the $\mathbf{Y}_{\text{bus}}$ method trails the element-wise
method for practical power grids on modern computers, especially those with
wide SIMD support and fast memory. Bus admittance matrix, as a network
reduction method for increasing computation performance, is missing the
support for new processor features due to its computational irregularity.
Future work needs to investigate algorithms for reducing the time complexity
of the Jacobian reduction step in the element-wise method.

## References

  * [1] William.. Tinney and Clifford E. Hart  “Power Flow Solution by Newton’s Method”  In _IEEE Transactions on Power Apparatus and Systems_ PAS-86.11, 1967, pp. 1449–1460  DOI: 10.1109/TPAS.1967.291823
  * [2] William F. Tinney  “Compensation Methods for Network Solutions by Optimally Ordered Triangular Factorization”  In _IEEE Transactions on Power Apparatus and Systems_ PAS-91.1, 1972, pp. 123–127  DOI: 10.1109/TPAS.1972.293321
  * [3] Philip E. Ross  “Why CPU Frequency Stalled”, 2008  IEEE Spectrum  URL: https://spectrum.ieee.org/why-cpu-frequency-stalled
  * [4] David Kirk and Wen-mei Hwu  “Programming Massively Parallel Processors: A Hands-on Approach”  Amsterdam ; Cambridge, MA: Elsevier, 2017
  * [5] Daniel Etiemble  “45-Year CPU Evolution: One Law and Two Equations”  arXiv, 2018  DOI: 10.48550/arXiv.1803.00254
  * [6] Michael J Flynn and Kevin W Rudd  “Parallel Architectures”  In _ACM Computing Surveys_ 28.1, 1996, pp. 4
  * [7] Hantao Cui, Fangxing Li and Kevin Tomsovic  “Hybrid Symbolic-Numeric Framework for Power System Modeling and Analysis”  In _IEEE Transactions on Power Systems_ 36.2, 2021  DOI: 10.1109/TPWRS.2020.3017019
  * [8] Federico Milano  “Power System Modelling and Scripting”  In _Power System Modelling and Scripting_ , Power Systems  Berlin, Heidelberg: Springer, 2010, pp. 3–17
  * [9] F. Milano  “An Open Source Power System Analysis Toolbox”  In _IEEE Transactions on Power Systems_ 20.3, 2005, pp. 1199–1206  DOI: 10.1109/TPWRS.2005.851911
  * [10] Joe H. Chow and Juan J. Sanchez-Gasca  “Power System Modeling, Computation, and Control”  John Wiley & Sons, 2020  GOOGLEBOOKS: M8G˙DwAAQBAJ
  * [11] J. Glover, Thomas Jeffrey Overbye, Mulukutla S. Sarma and Adam Birchfield  “Power System Analysis and Design”  Cengage Learning, 2022
  * [12] Timothy A Davis and Eka Palamadai Natarajan  “User Guide for KLU and BTF”, 2022
  * [13] Florian Schäfer and Martin Braun  “An Efficient Open-Source Implementation to Compute the Jacobian Matrix for the Newton-Raphson Power Flow Algorithm”  In _2018 IEEE PES Innovative Smart Grid Technologies Conference Europe (ISGT-Europe)_ , 2018, pp. 1–6  DOI: 10.1109/ISGTEurope.2018.8571471
  * [14] Santiago Peñate Vera  “SanPen/GridCal:”, 2022  Zenodo  DOI: 10.5281/zenodo.6896466
  * [15] Ray Daniel Zimmerman, Carlos Edmundo Murillo-Sánchez and Robert John Thomas  “MATPOWER: Steady-State Operations, Planning, and Analysis Tools for Power Systems Research and Education”  In _IEEE Transactions on Power Systems_ 26.1, 2011, pp. 12–19  DOI: 10.1109/TPWRS.2010.2051168


the role of $\mathfrak{D}(m)$ in its proof, one verifies that it is of the
above form.

######  Lemma 4.5.

Let $\alpha\in(0,1)$ and $f\in L_{p,w_{2}}$ with $p\in(1,\infty)$, $w_{2}\in
A_{p}$. Then for any $\kappa\in(0,\alpha)$, we have

$\|\left(P_{\alpha}\ast\psi_{j}\right)\ast f(t,\cdot)\|_{L_{p,w_{2}}}\leq
N(d,\alpha,\kappa,[w_{2}]_{A_{p}})\left(2^{-2\kappa j/\alpha}t^{-\kappa}\wedge
1\right)\|f\|_{L_{p,w_{2}}}$ (4.14)

for $j=0,\pm 1,\pm 2,\ldots$, and

$\|\left(P_{\alpha}\ast\Psi\right)\ast f(t,\cdot)\|_{L_{p,w_{2}}}\leq
N(d,\alpha,[w_{2}]_{A_{p}})\|f\|_{L_{p,w_{2}}}.$ (4.15)

###### Proof.

For simplicity, set $P_{\alpha}=P$ and $P_{\alpha}\ast\psi_{j}=P_{j}$ for
$j=0,\pm 1,\pm 2,\ldots$. By using the integral representation (3.8) in Remark
3.10 with the fact that $\hat{P}_{j}=\hat{P}\,\hat{\psi_{j}}$ and
$\operatorname{supp}\hat{\psi}_{j}\subset\\{2^{j-1}\leq|\xi|\leq 2^{j}\\}$, we
verify that

$\left|D_{\xi}^{\gamma}\hat{P}_{j}(t,\xi)\right|\leq
N(d,\alpha,\kappa,\gamma)2^{-j|\gamma|}\left(2^{-2j/\alpha}t^{-1}\right)^{\kappa}1_{2^{j-1}\leq|\xi|\leq
2^{j+1}}$ (4.16)

for any multi-index $\gamma$ and for any $\kappa\in(0,\alpha)$. Indeed, for
$\gamma=0$ and for any $\kappa\in(0,\alpha)$, we have

$\left|\hat{P}_{j}(t,\xi)\right|\leq N1_{2^{j-1}\leq|\xi|\leq
2^{j+1}}\int_{0}^{\infty}\frac{r^{\alpha-1}}{r^{2\alpha}+2r^{\alpha}\cos{\alpha\pi}+1}e^{-rt|\xi|^{2/\alpha}}\,\mathrm{d}r$
$\leq N1_{2^{j-1}\leq|\xi|\leq
2^{j+1}}\int_{0}^{1}r^{\alpha-1}e^{-rt|\xi|^{2/\alpha}}\left(rt|\xi|^{2/\alpha}\right)^{\kappa}\left(rt|\xi|^{2/\alpha}\right)^{-\kappa}\,\mathrm{d}r$
$+N1_{2^{j-1}\leq|\xi|\leq
2^{j+1}}\int_{1}^{\infty}r^{-\alpha-1}e^{-rt|\xi|^{2/\alpha}}\left(rt|\xi|^{2/\alpha}\right)^{\kappa}\left(rt|\xi|^{2/\alpha}\right)^{-\kappa}\,\mathrm{d}r$
$\leq N1_{2^{j-1}\leq|\xi|\leq
2^{j+1}}\left(2^{-2j/\alpha}t^{-1}\right)^{\kappa}\left(\int_{0}^{1}r^{\alpha-\kappa-1}\,\mathrm{d}r+\int_{1}^{\infty}r^{-\alpha-\kappa-1}\,\mathrm{d}r\right)$
$\leq N(d,\alpha,\kappa)1_{2^{j-1}\leq|\xi|\leq
2^{j+1}}\left(2^{-2j/\alpha}t^{-1}\right)^{\kappa},$

where we use the fact that
$e^{-rt|\xi|^{2/\alpha}}\left(rt|\xi|^{2/\alpha}\right)^{\kappa+m}$ is
uniformly bounded in $(r,t,\xi)$ for any $m\geq 0$ and the last inequality is
due to $\kappa\in(0,\alpha)$. The same argument holds for $\gamma$ such that
$|\gamma|\neq 0$, with the help of $|D_{\xi}^{\gamma_{1}}\hat{\psi}_{j}|\leq
N(d)2^{-j|\gamma_{1}|}$ for any multi-index $\gamma_{1}$. Then

$\sup_{r>0}\left(r^{2|\gamma|-d}\int_{r<|\xi|<2r}\left|D_{\xi}^{\gamma}\hat{P}_{j}(t,\xi)\right|^{2}\,\mathrm{d}\xi\right)^{1/2}$
$=N2^{-j|\gamma|-2\kappa
j/\alpha}t^{-\kappa}\sup_{2^{j-2}<r<2^{j+1}}\left(r^{2|\gamma|-d}\int_{r<|\xi|<2r}1_{2^{j-1}\leq|\xi|\leq
2^{j+1}}\,\mathrm{d}\xi\right)^{1/2}$ $\leq
N\left(2^{-2j/\alpha}t^{-1}\right)^{\kappa}$

for $N=N(d,\alpha,\kappa,\gamma)>0$.

On the other hand, since $|D_{\xi}^{\gamma}\hat{P}(t,\xi)|\leq
N\||\cdot|^{|\gamma|}P(t,\cdot)\|_{L_{1}}$ and $P(1,x)$ satisfies (3.6) in
Remark 3.10, by direct calculation we get

$\left|D_{\xi}^{\gamma}\hat{P}_{j}(t,\xi)\right|\leq
N\sum_{|\gamma_{1}|\leq|\gamma|}1_{2^{j-1}\leq|\xi|\leq
2^{j+1}}2^{-j|\gamma_{1}|}t^{\left|\left(\gamma\right|-\left|\gamma_{1}\right|\right)\alpha/2}$

for any multi-index $\gamma$, and then,

$\sup_{r>0}\left(r^{2|\gamma|-d}\int_{r<|\xi|<2r}\left|D_{\xi}^{\gamma}\hat{P}_{j}(t,\xi)\right|^{2}\,\mathrm{d}\xi\right)^{1/2}\leq
N\sum_{|\gamma_{1}|\leq|\gamma|}t^{\left(\left|\gamma\right|-\left|\gamma_{1}\right|\right)\alpha/2}2^{j\left(\left|\gamma\right|-\left|\gamma_{1}\right|\right)},$

where $N=N(d,\alpha,\gamma)>0$. Observe that if $t\leq 2^{-2j/\alpha}$, the
last summation in the above inequality is bounded by
$\sum_{|\gamma_{1}|\leq|\gamma|}1=N(\gamma)$. Therefore, for any multi-index
$\gamma$ with $|\gamma|\leq d$, we have (4.12) with $\hat{P}_{j}$ in place of
$m$, where

$\mathfrak{D}(\hat{P}_{j})\leq
N(d,\alpha,\kappa)\left(2^{-2j/\alpha}t^{-1}\right)^{\kappa}\wedge 1.$

This along with Lemma 4.4 proves (4.14).

To prove (4.15), recall that $\hat{\Psi}\in C_{0}^{\infty}(\mathbb{R}^{d})$,
$\operatorname{supp}\hat{\Psi}\subset\\{|\xi|\leq 2\\}$ and $\hat{\Psi}=1$ on
$\\{|\xi|\leq 1\\}$. It is easily seen that $\hat{\Psi}$ satisfies (4.12) with
$\mathfrak{D}(\hat{\Psi})=N(d)$ and then by Lemma 4.4,

$\|\left(P_{\alpha}\ast\Psi\right)\ast
f(t,\cdot)\|_{L_{p,w_{2}}}=\|\Psi\ast\left(P_{\alpha}\ast
f\right)(t,\cdot)\|_{L_{p,w_{2}}}\leq N(d)\|P_{\alpha}\ast
f(t,\cdot)\|_{L_{p,w_{2}}}.$

Also, by a suitable dyadic decomposition, we have $P_{\alpha}\ast f\leq
N(d,\alpha)\mathcal{M}f$ where $\mathcal{M}$ is the Hardy-Littlewood maximal
operator. Then by Hardy-Littlewood maximal function theorem with
$A_{p}$-weights, we obtain (4.15). The lemma is proved. ∎

We prove the key estimate of this subsection in the following proposition.

######  Proposition 4.6.

Let $p,q\in(1,\infty)$, $w=t^{\mu}w_{2}(x)$ with $\mu\in(-1,q-1)$, $w_{2}\in
A_{p}$. Suppose $\alpha\in\left(\left(1+\mu\right)/q,1\right)$. Then, for
$f\in C_{0}^{\infty}(\mathbb{R}^{d})$, we have

$\|P_{\alpha}\ast f\|_{L_{p,q,w}(T)}\leq
N\|f\|_{B_{p,q,w_{2}}^{-\frac{2(1+\mu)}{q\alpha}}},$ (4.17)

where $N=N(d,\alpha,p,q,[w_{2}]_{A_{p}},\mu,T)>0$.

###### Proof.

We set $P_{\alpha}=P$ for notational simplicity. Also, for $g\in\\{f,P\\}$, we
set $g_{j}=g\ast\psi_{j}$, $j=1,2,\ldots$, and $g_{0}=g\ast\Psi$ (not
$g\ast\psi_{0}$). Observe that
$\hat{\Psi}=\hat{\Psi}(\hat{\Psi}+\hat{\psi}_{1})$,
$\hat{\psi}_{1}=\hat{\psi}_{1}(\hat{\Psi}+\hat{\psi}_{1}+\hat{\psi}_{2})$, and
$\hat{\psi}_{j}=\hat{\psi}_{j}(\hat{\psi}_{j-1}+\hat{\psi}_{j}+\hat{\psi}_{j+1})$
for $j=2,3,\ldots$, by the construction of $\Psi$ and $\psi_{j}$. Then

$P\ast f=P_{0}\ast f_{0}+P_{1}\ast
f_{0}+\sum_{j=1}^{\infty}\sum_{i=j-1}^{j+1}P_{i}\ast f_{j},$

and by Lemma 4.5, for some $\kappa\in(0,\alpha)$ ($\kappa$ is to be specified
below.)

$\|P(t,\cdot)\ast f(\cdot)\|_{L_{p,w_{2}}}\leq\|P_{0}\left(t,\cdot\right)\ast
f_{0}\left(\cdot\right)\|_{L_{p,w_{2}}}+\|P_{1}(t,\cdot)\ast
f_{0}\left(\cdot\right)\|_{L_{p,w_{2}}}$
$+\sum_{j=1}^{\infty}\sum_{i=j-1}^{j+1}\|P_{i}(t,\cdot)\ast
f_{j}(\cdot)\|_{L_{p,w_{2}}}$ $\leq
N\left(\||f_{0}|+|f_{1}|\|_{L_{p,w_{2}}}+\sum_{j=1}^{\infty}\left(2^{-2\kappa
j/\alpha}t^{-\kappa}\wedge 1\right)\|f_{j}\|_{L_{p,w_{2}}}\right),$

where $N=N(d,\alpha,\kappa,p,[w_{2}]_{A_{p}})>0$. From the above inequalities,
it suffices to estimate $\mathcal{I}$ and $\mathcal{J}$ where

$\mathcal{I}:=\int_{0}^{T}\||f_{0}|+|f_{1}|\|_{L_{p,w_{2}}}^{q}t^{\mu}\,\mathrm{d}t$

and

$\mathcal{J}:=\int_{0}^{T}\left(\sum_{j=1}^{\infty}\left(2^{-2\kappa
j/\alpha}t^{-\kappa}\wedge
1\right)\|f_{j}\|_{L_{p,w_{2}}}\right)^{q}t^{\mu}\,\mathrm{d}t.$

It is clear that

$\mathcal{I}\leq
N(\alpha,q,\mu,T)\left(\|f_{0}\|_{L_{p,w_{2}}}^{q}+\|f_{1}\|_{L_{p,w_{2}}}^{q}\right),$

where $s=-2(1+\mu)/(q\alpha)$. On the other hand, note that $\mathcal{J}$ is
bounded by a positive constant times

$\int_{0}^{T}\left(\sum_{j=1}^{\infty}\|f_{j}\|_{L_{p,w_{2}}}1_{t\leq
2^{-2j/\alpha}}\right)^{q}t^{\mu}\,\mathrm{d}t$
$+\int_{0}^{T}\left(\sum_{j=1}^{\infty}2^{-2\kappa
j/\alpha}t^{-\kappa}\|f_{j}\|_{L_{p,w_{2}}}1_{t>2^{-2j/\alpha}}\right)^{q}t^{\mu}\,\mathrm{d}t:=\mathcal{J}_{1}+\mathcal{J}_{2}.$

Estimate of $\mathcal{J}_{1}$: Observe that, for $t>2^{-2/\alpha}$, the
integrand of $\mathcal{J}_{1}$ is always zero. For each
$t\in(0,2^{-2/\alpha}]$, there exits the largest positive integer
$j_{1}=j_{1}(t)$ such that $t\leq 2^{-2j_{1}/\alpha}$. Then by Hölder’s
inequality, for some $r<0$,

$\left(\sum_{j=1}^{\infty}\|f_{j}\|_{L_{p,w_{2}}}1_{t\leq
2^{-2j/\alpha}}\right)^{q}=\left(\sum_{j=1}^{j_{1}(t)}\|f_{j}\|_{L_{p,w_{2}}}\right)^{q}1_{t\leq
2^{-2j/\alpha}}$
$\leq\left(\sum_{j=1}^{j_{1}(t)}\left(2^{-r\frac{2j}{\alpha}}\right)^{q^{\prime}}\right)^{q/q^{\prime}}\sum_{j=1}^{j_{1}(t)}\left(2^{r\frac{2j}{\alpha}}\right)^{q}\|f_{j}\|_{L_{p,w_{2}}}^{q}1_{t\leq
2^{-2j/\alpha}}$ $\leq
N\left(2^{-2j_{1}/\alpha}\right)^{qr}\sum_{j=1}^{j_{1}}\left(2^{r\frac{2j}{\alpha}}\|f_{j}\|_{L_{p,w_{2}}}\right)^{q}1_{t\leq
2^{-2j/\alpha}}$ $\leq
Nt^{qr}\sum_{j=1}^{j_{1}}\left(2^{r\frac{2j}{\alpha}}\right)^{q}\|f_{j}\|_{L_{p,w_{2}}}^{q}1_{t\leq
2^{-2j/\alpha}},$

where $N=N(r,\alpha,q)$ and the last inequality is due to $t\leq
2^{-2j_{1}/\alpha}$ and $r<0$. Then, by taking $r<0$ such that $qr+\mu>-1$,
i.e., $-(1+\mu)/q<r<0$ (recall that $-(1+\mu)/q<0$ since $\mu>-1$),

$\mathcal{J}_{1}\leq
N\sum_{j=1}^{\infty}\left(\int_{0}^{2^{-2j/\alpha}}t^{qr+\mu}\,\mathrm{d}t\right)2^{\frac{2j}{\alpha}qr}\|f_{j}\|_{L_{p,w_{2}}}^{q}\leq
N\sum_{j=1}^{\infty}2^{jqs}\|f_{j}\|_{L_{p,w_{2}}}^{q},$

where $s=-2(1+\mu)/(q\alpha)$ and $N=N(\alpha,q,\mu)$.

Estimate of $\mathcal{J}_{2}$: By repeating the above process with
$t>2^{-2j/\alpha}$, where $j_{2}=j_{2}(t)$ is the smallest positive integer
such that $t>2^{-2j_{2}/\alpha}$,

$\displaystyle\left(\sum_{j=1}^{\infty}2^{-2\kappa
j/\alpha}t^{-\kappa}\|f_{j}\|_{L_{p,w_{2}}}1_{t>2^{-2j/\alpha}}\right)^{q}$
$\displaystyle\leq t^{-\kappa
q}\left(\sum_{j=j_{2}}^{\infty}\left(2^{-\frac{2j}{\alpha}(\kappa+r^{\prime})}\right)^{q^{\prime}}\right)^{q/q^{\prime}}\sum_{j=j_{2}}^{\infty}2^{\frac{2j}{\alpha}qr^{\prime}}\|f_{j}\|_{L_{p,w_{2}}}^{q}1_{t>2^{-2j/\alpha}}$
$\displaystyle=Nt^{-\kappa
q}2^{-\frac{2j_{2}}{\alpha}q(\kappa+r^{\prime})}\sum_{j=j_{2}}^{\infty}2^{\frac{2j}{\alpha}qr^{\prime}}\|f_{j}\|_{L_{p,w_{2}}}^{q}1_{t>2^{-2j/\alpha}}$

for some $r^{\prime}>-\kappa$ such that $qr^{\prime}+\mu<-1$ and
$N=N(\alpha,q,\mu)$. Note that such $r^{\prime}$ always exists since
$\kappa<\alpha$ and $\alpha>(1+\mu)/q$. Then, from the fact that
$t>2^{-2j_{2}/\alpha}$,

$\mathcal{J}_{2}\leq
N\sum_{j=1}^{\infty}\left(\int_{2^{-2j/\alpha}}^{\infty}t^{qr+\mu}\,\mathrm{d}t\right)2^{\frac{2j}{\alpha}qr}\|f_{j}\|_{L_{p,w_{2}}}^{q}=N\sum_{j=1}^{\infty}2^{jqs}\|f_{j}\|_{L_{p,w_{2}}}^{q},$

where $s=-2(1+\mu)/(q\alpha)$ and $N=N(\alpha,q,\mu)$. Thus, by taking
$\kappa\in(0,\alpha)$ such that $\kappa>(1+\mu)/q$ (for example,
$2\kappa=\alpha+(1+\mu)/q$) and combining the estimates for $\mathcal{I}$ and
$\mathcal{J}$, we have

$\|P_{\alpha}\ast f\|_{L_{p,q,w}(T)}^{q}\leq
N\left(\|f_{0}\|_{L_{p,w_{2}}}^{q}+\sum_{j=1}^{\infty}2^{jqs}\|f_{j}\|_{L_{p,w_{2}}}^{q}\right),$

that is, (4.17) with $N=N(d,\alpha,p,q,[w_{2}]_{A_{p}},\mu,T)$. The
proposition is proved. ∎

With the help of Proposition 4.6, we prove Theorem 3.14, the extension theorem
for $\mathbb{H}_{p,q,w}^{\alpha,2}(T)$ and
$\mathcal{H}_{p,q,w}^{\alpha,1}(T)$.

###### Proof of Theorem 3.14.

($i$) It is enough to find $U(t,x)\in\mathbb{H}_{p,q,w}^{\alpha,2}(T)$
satisfying (3.12) for $u_{0}(x)\in C_{0}^{\infty}(\mathbb{R}^{d})$. Take $U$
as in (3.3), that is, $U(t,\cdot)=P_{\alpha}(t,\cdot)\ast u_{0}(\cdot)$. Then
by applying Proposition 4.6 to $(1-\Delta)u_{0}$, we have

$\||U|+|DU|+|D^{2}U|\|_{L_{p,q,w}(T)}\leq
N\|(1-\Delta)u_{0}\|_{B_{p,q,w_{2}}^{2\theta-2}}=N\|u_{0}\|_{B_{p,q,w_{2}}^{2\theta}},$

where $N=N(d,\alpha,p,q,[w_{2}]_{A_{p}},\mu,T)$. This implies (3.12) since $U$
satisfies $\partial_{t}^{\alpha}U=\Delta U$ in $\mathbb{R}^{d}_{T}$.

($ii$) Let $w_{0}=(1-\Delta)^{-1/2}v_{0}\in B_{p,q,w_{2}}^{2\theta}$. Then by
($i$), there is a solution $W\in\mathbb{H}_{p,q,w}^{\alpha,2}(T)$ to (3.4)
satisfying $W(0,\cdot)=w_{0}(\cdot)$ and

$\|W\|_{\mathbb{H}_{p,q,w}^{\alpha,2}(T)}\leq
N\|w_{0}\|_{B_{p,q,w_{2}}^{2\theta}}=N\|v_{0}\|_{B_{p,q,w_{2}}^{2\theta-1}},$

where $N=N(d,\alpha,p,q,[w_{2}]_{A_{p}},\mu,T)$. By taking
$V=(1-\Delta)^{1/2}W\in\mathcal{H}_{p,q,w}^{\alpha,1}(T)$, we obtain (3.13).
In particular, since $\alpha>(1+\mu)/q$, $U(0,\cdot)=u_{0}(\cdot)$ and
$V(0,\cdot)=v_{0}(\cdot)$ are well-defined in the trace sense by Theorems 3.11
and 3.12, respectively. The theorem is proved. ∎

######  Remark 4.7.

Recall that in order for the initial trace $u(0,\cdot)\in L_{p,w_{2}}(\Omega)$
to make sense for $u\in\mathbb{H}_{p,q,w}^{\alpha,0}(\Omega_{T})$, we need
$\alpha\geq(1+\mu)/q$ with $\mu\in(-1,q-1)$. On the other hand, if $u_{0}\in
L_{p,w_{2}}(\Omega)$, the initial value of the homogeneous solution
$U=P_{\alpha}(t,\cdot)\ast u_{0}(\cdot)$ is well-defined regardless of the
range of $\alpha\in(0,1)$. This is made possible by appropriately choosing the
time integrability $q$ and the power $\mu$ of weight $w_{1}(t)=t^{\mu}$. In
fact, this is also one of the advantages of the weighted $L_{p,q}$-theory. For
example, we simply consider an unweighted space, i.e., $w(t,x)=1$, and an
initial value $u_{0}\in L_{p}$. For a given $\alpha\in(0,1)$, we choose
$q\in(1,\infty)$ such that $1/\alpha<q<2/\alpha$, i.e., $1<q\alpha<2$. Then,
for $s=1-2/(q\alpha)<0$, we have

$\|u_{0}\|_{B_{p,q}^{s}}\leq N\|u_{0}\|_{L_{p}},$

which means that $u_{0}\in L_{p}\subset B_{p,q}^{1-2/(q\alpha)}$. Together
with the condition $\alpha>1/q$, this ensures that the solution
$U\in\mathcal{H}_{p,q}^{\alpha,1}(T)$ satisfies $U(0,\cdot)=u(\cdot)$ in the
trace sense by Theorem 3.12 and Theorem 3.14-($ii$). On the other hand, by
direct calculation, one can see that

$\|U(t,\cdot)-u_{0}(\cdot)\|_{L_{p}(\Omega)}\to 0\quad\textrm{as}\quad t\to
0.$

For another perspective, we refer the reader to [21] and [13], which discuss
the $L_{2}$-theory.

The following is a version of Lemma 4.5 for the case of
$\beta=1+\alpha\in(0,1)$.

######  Lemma 4.8.

Let $\alpha\in(0,1)$, $\beta=1+\alpha$ and $f\in L_{p,w_{2}}$ with
$p\in(1,\infty)$, $w_{2}\in A_{p}$. Then for any $\kappa_{1}\in(0,\beta)$ and
$\kappa_{2}\in(0,\alpha)$, we have

$\|\left(P_{\beta}\ast\psi_{j}\right)\ast f(t,\cdot)\|_{L_{p,w_{2}}}\leq
N(d,\alpha,\kappa_{1},[w_{2}]_{A_{p}})\left(2^{-2\kappa_{1}j/\beta}t^{-\kappa_{1}}\wedge
1\right)\|f\|_{L_{p,w_{2}}},$ (4.18)
$\|\left(\widetilde{P}_{\beta}\ast\psi_{j}\right)\ast
f(t,\cdot)\|_{L_{p,w_{2}}}\leq
N(d,\alpha,\kappa_{2},[w_{2}]_{A_{p}})\left(2^{-2j/\beta}2^{-2\kappa_{2}j/\beta}t^{-\kappa_{2}}\wedge
t\right)\|f\|_{L_{p,w_{2}}}$ (4.19)

for $j=0,\pm 1,\pm 2,\cdots$, and

$\left\|\left(P_{\beta}\ast\Psi\right)\ast
f(t,\cdot)\right\|_{L_{p,w_{2}}}\leq
N(d,\alpha,[w_{2}]_{A_{p}})\|f\|_{L_{p,w_{2}}},$ (4.20)
$\left\|\left(\widetilde{P}_{\beta}\ast\Psi\right)\ast
f(t,\cdot)\right\|_{L_{p,w_{2}}}\leq
N(d,\alpha,[w_{2}]_{A_{p}})t\|f\|_{L_{p,w_{2}}}.$ (4.21)

###### Proof.

The inequality (4.18) can be obtained in the same way as in the proof of Lemma
4.5. The only difference is, since $\beta>1$, we need to use the
representation (3.7) in Remark 3.10 for two parametric Mittag-Leffler function
$E_{\beta,1}(=E_{\beta})$ instead of (3.8). Note that when we investigate
(3.7), the fact that $\cos(\delta/\beta)<0$ for any
$\delta\in(\frac{\pi\beta}{2},\pi]$ is useful. Similarly, we get (4.19) with
the help of the relation between the Fourier transform of
$\widetilde{P}_{\beta}$ and two parametric Mittag-Leffler function
$E_{\beta,2}$, that is, (3.9) in Remark 3.10. We also get (4.20) the exact
same way as in the proof of (4.15). Then (4.21) directly follows from (4.20),
since $\widetilde{P}_{\beta}(t,x)=\int_{0}^{t}P(s,x)\,\mathrm{d}s$. The lemma
is proved. ∎

By following the proof of Theorem 3.14 with the help of Lemma 4.8, we prove
Theorem 3.15.

###### Proof of Theorem 3.15.

We omit the proof since it is almost identical to the proof of Theorem 3.14,
except we use Lemma 4.8 instead of Lemma 4.5. We remark that we must have
$\alpha>(1+\mu)/q$ to be $\theta_{1}\in(0,1)$, but regarding $\theta_{0}$, we
always have $\theta_{0}\in(0,1)$ for any $\alpha\in(0,1)$ because $\mu<q-1$. ∎

## References

  * [1] Antonio Agresti, Nick Lindemulder, and Mark Veraar.  On the trace embedding and its applications to evolution equations.  Math. Nachr., 296(4):1319–1350, 2023.
  * [2] Mark Allen, Luis Caffarelli, and Alexis Vasseur.  A parabolic problem with a fractional time derivative.  Arch. Ration. Mech. Anal., 221(2):603–630, 2016.
  * [3] Mkhitar M. Djrbashian.  Harmonic analysis and boundary value problems in the complex domain, volume 65 of Operator Theory: Advances and Applications.  Birkhäuser Verlag, Basel, 1993.  Translated from the manuscript by H. M. Jerbashian and A. M. Jerbashian [A. M. Dzhrbashyan].
  * [4] Hongjie Dong and Doyoon Kim.  On $L_{p}$-estimates for elliptic and parabolic equations with $A_{p}$ weights.  Trans. Amer. Math. Soc., 370(7):5081–5130, 2018.
  * [5] Hongjie Dong and Doyoon Kim.  $L_{p}$-estimates for time fractional parabolic equations with coefficients measurable in time.  Adv. Math., 345:289–345, 2019.
  * [6] Hongjie Dong and Doyoon Kim.  An approach for weighted mixed-norm estimates for parabolic equations with local and non-local time derivatives.  Adv. Math., 377:107494, 44, 2021.
  * [7] Hongjie Dong and Doyoon Kim.  Time fractional parabolic equations with measurable coefficients and embeddings for fractional parabolic Sobolev spaces.  Int. Math. Res. Not. IMRN, (22):17563–17610, 2021.
  * [8] Hongjie Dong and Yanze Liu.  Weighted mixed norm estimates for fractional wave equations with VMO coefficients.  J. Differential Equations, 337:168–254, 2022.
  * [9] Samuil D. Eidelman and Anatoly N. Kochubei.  Cauchy problem for fractional diffusion equations.  J. Differential Equations, 199(2):211–255, 2004.
  * [10] Rudolf Gorenflo, Joulia Loutchko, and Yuri Luchko.  Computation of the Mittag-Leffler function $E_{\alpha,\beta}(z)$ and its derivative.  volume 5, pages 491–518. 2002.
  * [11] Gerd Grubb.  Pseudo-differential boundary problems in $L_{p}$ spaces.  Comm. Partial Differential Equations, 15(3):289–340, 1990.
  * [12] G. H. Hardy and J. E. Littlewood.  Some properties of fractional integrals. I.  Math. Z., 27(1):565–606, 1928.
  * [13] Yavar Kian and Masahiro Yamamoto.  Well-posedness for weak and strong solutions of non-homogeneous initial boundary value problems for fractional diffusion equations.  Fract. Calc. Appl. Anal., 24(1):168–201, 2021.
  * [14] Doyoon Kim, Kyeong-Hun Kim, and Kwan Woo.  Trace theorem and non-zero boundary value problem for parabolic equations in weighted Sobolev spaces.  Stoch PDE: Anal Comp, 2022.
  * [15] Doyoon Kim, Seungjin Ryu, and Kwan Woo.  Parabolic equations with unbounded lower-order coefficients in Sobolev spaces with mixed norms.  J. Evol. Equ., 22(1):Paper No. 9, 40, 2022.
  * [16] Ildoo Kim, Kyeong-Hun Kim, and Sungbin Lim.  An $L_{q}(L_{p})$-theory for the time fractional evolution equations with variable coefficients.  Adv. Math., 306:123–176, 2017.
  * [17] Ildoo Kim, Kyeong-Hun Kim, and Sungbin Lim.  A Sobolev space theory for stochastic partial differential equations with time-fractional derivatives.  Ann. Probab., 47(4):2087–2139, 2019.
  * [18] Kyeong-Hun Kim, Daehan Park, and Junhee Ryu.  An $L_{q}(L_{p})$-theory for diffusion equations with space-time nonlocal operators.  J. Differential Equations, 287:376–427, 2021.
  * [19] Anatoly N. Kochubei.  Cauchy problem for fractional diffusion-wave equations with variable coefficients.  Appl. Anal., 93(10):2211–2242, 2014.
  * [20] Adam Kubica, Katarzyna Ryszewska, and Masahiro Yamamoto.  Time-fractional differential equations—a theoretical introduction.  SpringerBriefs in Mathematics. Springer, Singapore, [2020] ©2020.
  * [21] Adam Kubica and Masahiro Yamamoto.  Initial-boundary value problems for fractional diffusion equations with time-dependent coefficients.  Fract. Calc. Appl. Anal., 21(2):276–311, 2018.
  * [22] Douglas S. Kurtz and Richard L. Wheeden.  Results on weighted norm inequalities for multipliers.  Trans. Amer. Math. Soc., 255:343–362, 1979.
  * [23] Ralf Metzler and Joseph Klafter.  The random walk’s guide to anomalous diffusion: a fractional dynamics approach.  Phys. Rep., 339(1):77, 2000.
  * [24] Norman G. Meyers and James Serrin.  $H=W$.  Proc. Nat. Acad. Sci. U.S.A., 51:1055–1056, 1964.
  * [25] Martin Meyries and Roland Schnaubelt.  Interpolation, embeddings and traces of anisotropic fractional Sobolev spaces with temporal weights.  J. Funct. Anal., 262(3):1200–1229, 2012.
  * [26] Martin Meyries and Mark C. Veraar.  Traces and embeddings of anisotropic function spaces.  Math. Ann., 360(3-4):571–606, 2014.
  * [27] Jan Prüss.  Evolutionary integral equations and applications. Monogr. Math., 87, Birkhäuser Verlag, Basel, 1993. xxvi+366 pp., ISBN:3-7643-2876-2.
  * [28] Bui Huy Qui.  Weighted Besov and Triebel spaces: interpolation by the real method.  Hiroshima Math. J., 12(3):581–605, 1982.
  * [29] Stefan G. Samko, Anatoly A. Kilbas, and Oleg I. Marichev.  Fractional integrals and derivatives.  Gordon and Breach Science Publishers, Yverdon, 1993.  Theory and applications, Edited and with a foreword by S. M. Nikol′skiĭ, Translated from the 1987 Russian original, Revised by the authors.
  * [30] Hans Triebel.  Interpolation theory, function spaces, differential operators. North-Holland Publishing Co., Amsterdam-New York, 1978.
  * [31] Hans Triebel.  Theory of function spaces. Monogr. Math., 78, Birkhäuser Verlag, Basel, 1983. 284 pp., ISBN:3-7643-1381-1.
  * [32] Rico Zacher.  Maximal regularity of type $L_{p}$ for abstract parabolic Volterra equations.  J. Evol. Equ., 5(1):79–103, 2005.
  * [33] Rico Zacher.  Quasilinear parabolic integro-differential equations with nonlinear boundary conditions.  Differential Integral Equations, 19(10):1129–1156, 2006.
  * [34] Rico Zacher.  Weak solutions of abstract evolutionary integro-differential equations in Hilbert spaces.  Funkcial. Ekvac., 52(1):1–18, 2009.