# On the damping of power systems during inverter-based blackstart ††thanks:
A.Anta and D.Cifelli are with AIT Austrian Institute of Technology GmbH,
Vienna, Austria, email: {Adolfo.Anta<EMAIL_ADDRESS>
Adolfo Anta, Diego Cifelli

###### Abstract

The massive deployment of inverter-based generation poses several challenges
to system operators but also offers new opportunities. In the context of grid
recovery, inverter-based generation is expected to take over new
responsibilities, as conventional generation that are blackstart capable are
being decommissioned. Inverter-based facilities, being very flexible and
controllable, have the ability to contribute during grid restoration, although
it is still unclear how to deal with the limited ratings of batteries, PV,
etc. and how to address nonlinear effects such as inrush currents or
harmonics. In this direction, the requirements of grid-forming inverters to
comply with black start capabilities are still to be defined. During the first
steps of a grid restoration process the dynamics of the grid largely vary from
the grid behaviour under normal operation. In particular, the grid is nearly
unloaded and therefore poorly damped. Indeed, overvoltages and resonances are
a common concern during the energization of long transmission lines. In this
work we analyse the spectral properties of unloaded transmission lines, and
describe the impact of voltage controllers on the damping of the system. As a
byproduct, this analysis can be leveraged to propose high level requirements
for grid-forming devices, in order to provide adequate damping during the
early stages of a grid restoration process.

###### Index Terms:

Grid restoration; grid-forming control; resonances; eigenvalue analysis;
Toeplitz matrices.

##  I Introduction and motivation

The introduction of inverter-based generation at large scale is causing issues
in the planning and operation of power grids. At the same time, it raises the
opportunity to upgrade and optimize standard approaches using inverter-based
resources, as inverters are much more flexible and controllable than
conventional generation. This transition is also affecting ongoing strategies
for grid restoration. Inverter-based generation have the ability to contribute
during the black start [1, 2, 3], although it is still unclear how these
inverters should behave, and how system-level strategies may need to be
adapted.

There has been considerable recent work at the unit level to assess the black
start capabilities of inverter-based generation. However, when it comes to
system-level analysis, grid restoration has been widely studied from a static
perspective, with a focus on designing topologies for each step of the
blackstart. Even in the few cases where dynamic aspects are considered, the
focus lies on grids dominated by conventional power plants, thus modelled
using swing equations and ignoring line dynamics. In this work we focus
instead on the dynamic aspects during the early stages of a grid restoration
under inverter-based generation, and in particular we analyse the damping (or
lack thereof) during the energization of transmission lines by an inverter-
based unit.

In terms of strategies, grid restoration typically follows a divide-and-
conquer approach: a grid is divided in several regions that are independently
restored by black-start capable power plants. This approach creates many
regional constraints and faces the issue of resynchronization among the
regions. A faster and simpler strategy consists in defining a backbone
structure, as extensive as possible, that is energized at once. This approach
is being tested already by system operators with promising results [4]. The
backbone energization is performed by means of a single unit, typically using
a so-called soft start (a ramp in voltage magnitude), that minimizes certain
adverse effects such as inrush currents [5]. During a backbone energization,
the lack of loads connected in the grid implies a very low level of damping,
and resonances and overvoltages are common and can jeopardize an adequate
black start. This is especially the case under the presence of transformers,
whose energization can lead to a significant harmonic content. In traditional
black start approaches, loads may be connected in order to avoid these issues,
given that conventional power plants are characterized by large ratings. Due
to the limited ratings of inverters, it is likely that all loads need to be
disconnected in order to allow the inverters to energize the grid. Moreover,
the expected resonances are different from the ones at normal operation [6]
(and partially unknown due to the uncertainties during black start), hence
designing damping controls targeting a particular frequency may be
ineffective. In particular, the first resonant frequency can be lower than
during normal operation, and therefore more detrimental [7].

This paper addresses these issues by analysing the dynamics during the early
stages of a black start, considering transmission lines using $\pi$-sections
and simple models of grid-forming units. Unlike in previous works, we do not
rely on numerical studies but rather develop an analytical framework to study
the impact of the grid and different control parameters on the system damping.
We first derive the dynamics of a transmission line under no load, and its
similarities with tridiagonal 2-Toeplitz matrices are pointed out. These
matrices have very special characteristics, and expressions for the
eigenvalues can be analytically computed. Using these expressions, we later
study the influence of certain control structures on the energizing units,
where a clear trade-off between reference tracking and damping can be
observed. Moreover, we show how adding virtual resistance loops can contribute
to the overall damping of the system. We hope this analysis contributes
towards the definition of high level specifications for black start
functionality from grid-forming units. Indeed, customized black start control
modes may be required during the energization. Most of the focus on grid-
forming research has centered around frequency and angle dynamics, while for
this problem the voltage magnitude is highly pertinent.

The rest of the paper is structured as follows: first, we model a transmission
line using $\pi$-sections, to be energized by controllable voltage sources,
and construct a state space representation of the line. Next, we establish the
connection between our dynamic model and the well-studied field of Toeplitz
matrices and Chebyshev polynomials. Using this parallelism, we derive a set of
nonlinear equations implicitly defining the eigenvalues of the system for
different cases: in open loop, under a proportional voltage controller, and
using a virtual resistance. Finally, an example shows the impact of the
controller parameters on the damping of the system. We conclude the paper
enumerating a list of possible future steps to better understand the grid
dynamics during this early stage of the grid restoration, and thereby defining
possible control structures for black-start capable grid-forming converters.

Figure 1: Transmission line model with cascaded $\pi$ networks

##  II Modelling

We consider a transmission line (see Figure 1) representing the grid backbone,
to be modelled by means of a cascaded connection of $\pi$-sections. Unlike
distributed parameter line models, these sections allow us to represent the
dynamics through a linear state space model and a finite number of states [8].
For simplicity, a single phase system is selected. At both ends of the line
voltage sources, representing grid-forming devices, are responsible for the
energization of the line. Writing one differential equation per dynamic
element (capacitors and inductances) leads us to a standard linear
representation of the type $\dot{x}=Ax+Bu$ such as in (27), with $R$, $G$, $L$
and $C$ representing the resistance, conductance inductance and capacitance
per $\pi$-section, $i$ and $v$ denoting the current and the voltage for each
branch and node, and $v_{a}$, $v_{b}$ the voltage sources. Non-mentioned
entries in the $A$ matrix in (27) are equal to zero.

The $A$ matrix in the state space representation is not only a tridiagonal
matrix, but in particular it is consistent with the structure of a 2-Toeplitz
matrix, for which eigenvalues can be analytically computed, as explained in
the next section. Grid dynamics have been widely studied, but in general with
the focus on analyzing stability properties [9]. At the same time, damping in
power grids has also been previously explored [10, 11], but ignoring line
dynamics and focusing on the dynamics of conventional generation. It has been
conjectured as well that resonant frequencies and amplitudes depend on the
short circuit power, the degree of meshing and amount of connected loads.
However, a theoretical understanding of these design choices is missing.
Moreover, it is unclear whether this intuition still holds when we replace
synchronous machines with power converters.

$\displaystyle\frac{d}{dt}\left[\begin{array}[]{c}i_{1}\\\ v_{1}\\\ i_{2}\\\
\vdots\\\ v_{s-1}\\\ i_{s}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{rrrrrr}-R/L&-1/L&&&&\\\
1/C&-G/C&-1/C&&&\\\ &1/L&-R/L&-1/L&&\\\ &&\ddots&\ddots&\ddots&\\\
&&&1/C&-G/C&-1/C\\\
&&&&1/L&-R/L\end{array}\right]\left[\begin{array}[]{c}i_{1}\\\ v_{1}\\\
i_{2}\\\ \vdots\\\ v_{s-1}\\\
i_{s}\end{array}\right]+\left[\begin{array}[]{cccccc}1/L&0\\\ 0&0\\\ 0&0\\\
\vdots&\vdots\\\ 0&0\\\
0&-1/L\end{array}\right]\left[\begin{array}[]{c}v_{a}\\\
v_{b}\end{array}\right]$ (27)

The controllable voltage sources $v_{a},v_{b}$ track a certain voltage
reference, either a step or a ramp. We start with a simple proportional
controller, for instance:

$v_{a}(t)=k_{V}(v_{ref_{A}}-v_{1}),\>\>v_{b}(t)=k_{V}(v_{ref_{B}}-v_{s-1})$
(28)

where $k_{V}>0$ is the proportional gain and $v_{1},v_{s-1}$ are the voltages
at their local node, which is the typical measurement of the inverter at the
point of common coupling. The reference voltage $v_{ref}$ for the energizing
unit(s) during black start is typically selected to be a ramp, in order to
minimize the impact of resonances, overvoltages and inrush phenomena [5]. A
proportional-integral (PI) voltage controller may not be adequate since it may
demand unnecessary reactive power, and perfect tracking is not necessary.
Likewise, there is typically no clear reactive power reference to be used for
the voltage regulation scheme in (28), since the the reactive power demand is
unclear during black start. High control gains would be prefered to guarantee
adequate tracking, but this may reduce overall damping, as we will see in the
next section. Moreover, should this soft energization be carried out by
several units simultaneously, it is necessary to avoid having a very
aggressive regulation at each unit that may lead to large power flows between
the energizing units.

While in some cases the ramp may take in the order of seconds [4], it is also
common to slow down the ramps to take up to several minutes due to concerns of
possible resonances. Slowing down has two drawbacks: first, protections may be
triggered as the voltage rises very slowly. Second, it slows down the whole
grid restoration process. It is also unclear how slow the ramp needs to be in
order to avoid resonances. Ideally, we would find improved methods to tackle
possible resonances, especially given the flexibility and controllability of
power converters.

######  Remark 1.

In this paper we model the energizing unit as an ideal, controllable voltage
source, since the focus of this work is not on the internal controls of grid-
forming units but rather to describe the expected high level behaviour of the
grid-forming unit during grid restoration, with the goal of defining a
possible high level specification for black start capable inverters. This is
typically a valid assumption if the AC and DC sides of the converter are
fairly decoupled [12].

##  III Spectral analysis and analytical characterization

In order to analyze the eigenvalues of the system, we draw in this section
parallelisms between the state space representation in (27) and Toeplitz
matrices. This helps us deriving analytical expressions for the spectrum and
its dependence on the voltage source and its controller gains.

###  III-A An overview to 2-Toeplitz matrices and Chebyshev polynomials

A matrix $A$ is tridiagonal if $a_{ij}=0$ whenever $|i-j|>1$. Moreover, $A$ is
r-Toeplitz of order $n$ if $a_{i+r,j+r}=a_{ij}$, for $i,j=1,2,...n-r$ [13].
With this definition, it can be seen that the state space matrix $A$ in (27)
is a 2-Toeplitz matrix. The inverse and spectra of these matrices have been
widely studied in mathematics, and has many applications in different topics
such as wave dispersion models or fiber optic design [14, 15]. In this section
we briefly summarize existing results on $2$-Toeplitz matrices that will be
helpful in computing analytical expressions for the eigenvalues of the $A$
matrix in (27).

The characteristic polynomial of these matrices can be written as a function
of Chebyshev polynomials. We first define the Chebyshev polynomials of second
kind, that satisfy the following 3-point recurrent relationship:

$U_{n+1}(\lambda):=\lambda U_{n}(\lambda)-U_{n-1}(\lambda),\qquad\forall
n=1,2,\ldots\>\>\lambda\in\mathbb{C}$ (29)

with initial conditions $U_{0}(\lambda)=1$ and $U_{1}(\lambda)=2\lambda$.
Moreover, $U_{n}$ satisfies:

$U_{n}(\lambda)=\frac{\sin((n+1)\theta)}{\sin(\theta)},\text{with
}\lambda=\cos(\theta),\>0\leq\theta<\pi$ (30)

For convenience we now define the following functions:

$\displaystyle g(\lambda)$ $\displaystyle:=$
$\displaystyle\left(\lambda+\frac{R}{L}\right)\left(\lambda+\frac{G}{C}\right)$
(31) $\displaystyle P_{n}(\lambda)$ $\displaystyle:=$
$\displaystyle\frac{1}{\left(LC\right)^{n}}U_{n}\left(\frac{1}{2}LC\lambda+1\right)$
(32)

With these definitions, we can already state that the eigenvalues of the $A$
matrix in (27) correspond to the solution of the characteristic polynomial

$\Delta_{2n+1}(\lambda):=\left(\lambda+\frac{R}{L}\right)P_{n}(g(\lambda))=0$
(33)

if the dimensions of the $A$ matrix are odd, and

$\Delta_{2n}(\lambda):=P_{n}(g(\lambda))-\frac{1}{LC}P_{n-1}(g(\lambda))=0$
(34)

if the dimensions of $A$ are even, see [13, 16] for a detailed
explanation111Notice that, unlike in [16], in our case the upper diagonal
elements are negative. Furthermore, in our $A$ matrix the upper diagonal is
identical to the lower diagonal with opposite signs, although this does not
seem to lead to a major simplification in the eigenvalue computation..

###  III-B Open Loop

In the case of our transmission line in (1), the dimensions of the state space
representation are $2n+1$, where $n$ is the number of $\pi$-sections. The
equation in (33) lead us to two nonlinear equations, using the definition in
(30) and identifying the argument inside the Chebyshev polynomial $U_{n}$ in
(32) to $\cos(\theta)$:

$\displaystyle\left(\lambda+\frac{R}{L}\right)\frac{1}{(LC)^{n}}\frac{\sin((n+1)\theta)}{\sin(\theta)}=0$
(35)
$\displaystyle\cos(\theta)=\frac{1}{2}LC\left(\lambda+\frac{R}{L}\right)\left(\lambda+\frac{G}{C}\right)+1$
(36)

Equation (36) is quadratic in $\lambda$, so it can be solved as a function of
$\cos(\theta)$. At the same time, (35) defines the conditions for $\theta$
plus another solution for $\lambda$, leading to the following eigenvalues:

$\displaystyle\lambda_{0}=\frac{-R}{L},$ (37) $\displaystyle\lambda_{\pm
k}=-\frac{1}{2}\left(\frac{R}{L}+\frac{G}{C}\right)\pm\sqrt{\frac{1}{4}\left(\frac{R}{L}-\frac{G}{C}\right)^{2}-\frac{2}{LC}\left(1-\cos(\theta)\right)}$

where $\theta=\frac{k\pi}{n+1},k=1,2\ldots n$. Notice that (35) is well
defined for $\theta=0$, and in particular $\lambda=\lambda_{0}$ corresponds to
$\theta=0$ in (36).

A few insights can be drawn from this formula. First, in the case of a single
$\pi$-section and negligible values of $R$ and $G$, the convoluted equation in
(LABEL:lambdaComplex) simplifies to the well known expression of $\lambda_{\pm
1}=1/\sqrt{LC}$ for the complex eigenvalues, as expected. Moreover, it can be
easily seen how, for $n>1$, this term $1/\sqrt{LC}$ can be fairly inaccurate
in describing the natural frequency of the eigenvalues [17]. For the typical
values of $R$, $L$, $G$, $C$ and $n$ the eigenvalues $\lambda_{\pm k}$ are
always complex, with the imaginary part being much larger than the real part,
hence the poor damping of the unloaded system. The imaginary part is always
bounded by $2/\sqrt{LC}$, given that the term $2/LC$ inside the square root is
much larger than the other terms. However, this bound is not very meaningful
since we are mainly interested in the first resonances, as the higher modes
are rarely excited. Moreover, if we increase the number of $\pi$-sections to
model a line (but correspondingly reducing the values of $R$,$L$, $G$ and $C$
per section), the real part stays constant but the imaginary part of the
eigenvalues increases proportional to the number of sections, thereby
modelling a higher frequency range in a more accurate way. This is expected
and consistent with the understanding of $\pi$-section modelling.

###  III-C The effect of voltage control

We consider now the impact of the voltage regulator on the eigenvalues and the
damping of the grid dynamics. As mentioned before, relatively high gains may
be needed in order to properly track the voltage reference. For simplicity, we
consider only the voltage source $v_{b}$, although by symmetry the effect of
the other unit will be identical. Closing the loop by means of the voltage
regulator in (28) entails modifying the $(2n+1,2n)$ element of the state space
matrix $A$, from $1/L$ to $1/L+k_{V}/L$. The state space matrix is no longer
Toeplitz but nonetheless we can derive the characteristic polynomial by
computing the determinant of $\lambda I-A$ along the last row using Laplace
expansion:

$\displaystyle\Delta_{2n+1}(\lambda)=$
$\displaystyle\left(\lambda+\frac{R}{L}\right)\Delta_{2n}(\lambda)+\left(\frac{1}{L}+\frac{k_{V}}{L}\right)\frac{1}{C}\Delta_{2n-1}(\lambda)$
$\displaystyle=$
$\displaystyle\left(\lambda+\frac{R}{L}\right)(P_{n}(g(\lambda))-\frac{1}{LC}P_{n-1}(g(\lambda)))$
$\displaystyle+\left(\frac{1}{L}+\frac{k_{V}}{L}\right)\frac{1}{C}\left(\lambda+\frac{R}{L}\right)P_{n-1}(g(\lambda))$
$\displaystyle=$
$\displaystyle\left(\lambda+\frac{R}{L}\right)\left(P_{n}(g(\lambda))+\left(\frac{k_{V}}{LC}\right)P_{n-1}(g(\lambda))\right)$
(39)

where we have used the expressions in (34) and (33) for odd and even order. As
in the open loop case in Section III-B, these leads to the following two
nonlinear equations:

$\displaystyle\left(\lambda+\frac{R}{L}\right)\frac{\sin((n+1)\theta)+k_{V}\sin(n\theta)}{\sin(\theta)}=0$
(40)
$\displaystyle\cos(\theta)=\frac{1}{2}LC\left(\lambda+\frac{R}{L}\right)\left(\lambda+\frac{G}{C}\right)+1$
(41)

where (41) is identical to (36). As before, the eigenvalues can be computed as
a function of the angle $\theta$ using (41), reaching to the same expressions
as in (37), (LABEL:lambdaComplex), but where now the angle satisfies the
following trigonometric equation:

$\sin((n+1)\theta)+k_{V}\sin(n\theta)=0$ (42)

Notice that again $\theta=0$ does not necessarily satisfy (40), and thus
$k_{V}$ affects all complex eigenvalues. However, the unique real eigenvalue
is not affected by the control gain and satisfies (40) for $\lambda=-R/L$ and
(41) with $\theta=0$. We can clearly see that $k_{V}$ only affects the angle
$\theta$, and therefore, under typical values of $L$, $C$, $G$, $R$ and $n$,
it does not influence the real component but only the imaginary part of the
eigenvalues in (LABEL:lambdaComplex). Hence, it can be already hinted that the
proportional control influences negatively the damping of the system. Even
though the relationship between the controller gain and $\theta$ is not
trivial, we can actually prove that the term $1-\cos(\theta)$ is monotonically
non-decreasing with respect to $k_{V}$, as stated in the following
proposition:

######  Proposition 1.

Consider $\theta(k_{V})$ as implicitly defined by (42). Then,

$\frac{d}{dk_{V}}(1-\cos(\theta))\geq 0,\quad\forall\>\theta\in[0,\pi[$ (43)

###### Proof.

We take as starting point the implicit relationship between $\theta$ and the
gain $k_{V}$ in (42), and differentiate with respect to $k_{V}$

$\displaystyle\frac{d}{dk_{V}}\left(\sin((n+1)\theta)+k_{V}\sin(n\theta)\right)=$
$\displaystyle\cos((n+1)\theta)(n+1)\frac{d\theta}{dk_{V}}+\sin(n\theta)+k_{v}\cos(n\theta)n\frac{d\theta}{dk_{V}}=0$

By replacing $k_{V}$ using again (42) we obtain a closed expression for
$d\theta/dk_{V}$:

$\displaystyle\frac{d\theta}{dk_{V}}=$ (44)
$\displaystyle\frac{-\sin(n\theta)^{2}}{(n+1)\cos((n+1)\theta)\sin(n\theta)-n\cos(n\theta)\sin((n+1)\theta)}$

The numerator is clearly non-positive, and the denominator can be shown to be
always non-positive using trigonometric identities and the fact that
$\sin(n\theta)\leq n\sin(\theta)$ (omitted here for conciseness). Moreover,
using L’Hôpital’s rule it can be shown that (44) is well defined for
$\theta\in[0,\pi[$, since the the denominator is only 0 at $\theta=0$.
Therefore, $d\theta/dk_{V}\geq 0$. Finally, we compute the dependence of
$1-\cos(\theta)$ as a function of the controller gain:

$\frac{d}{dk_{V}}(1-\cos(\theta))=\sin(\theta)\frac{d\theta}{dk_{V}}$

which is also greater or equal than 0 in the domain of $\theta$, see (30). ∎

######  Remark 2.

The zeros of $\frac{d}{dk_{V}}(1-\cos(\theta))$ correspond to $\theta=0$ and
$\sin(n\theta)=0$. In both cases, it can be seen that these values of $\theta$
do not satisfy (40) and (41) for the complex eigenvalues, and therefore for
all values of $\theta$ of interest the derivative is strictly positive.

Given that the square root is a monotonic function, it can be stated that
increasing $k_{V}$ exclusively increases the imaginary part of the complex
eigenvalues, that is:

$\displaystyle\frac{\partial}{\partial
k_{V}}\mathfrak{Im}\\{\lambda\\}>0,\qquad\frac{\partial}{\partial
k_{V}}\mathfrak{Re}\\{\lambda\\}=0$ (45)

This result clearly shows the existing trade-off between damping and tracking
error.

######  Remark 3.

An explicit solution for (42) can be easily found for $k_{V}=1$. In that case,
using simple trigonometric identities it is straightforward to see that:

$\sin((n+1)\theta)+\sin(n\theta)=2\sin\left(\left(n+\frac{1}{2}\right)\theta\right)\cos\left(\frac{\theta}{2}\right)=0$
(46)

Hence, $\theta=\frac{k\pi}{n+1/2}$, for $k=1,2,\ldots n$ is a solution, which
is indeed a larger angle than in the open loop, as expected from the previous
discussion. For the general case, it has been numerically observed that
$\theta=\frac{k\pi}{n+1/(1+k_{V})}$ represents a fairly accurate
approximation. Interestingly, this implies that, for large values of $k_{V}$,
the eigenvalues converge to the open loop eigenvalues of a model with $n-1$
$\pi$-sections. An explicit expression for the angle of the form
$\theta=\frac{k\pi}{n+1/g(k_{V})}$ (for some monotonically increasing function
$g$) may not exist but will be the focus of future research.

###  III-D The effect of virtual resistance and other variables

The approach in the previous subsection can be applied to other grid
parameters and other control loops. For instance, given that the proportional
controller decreases the (already low) damping of the system, a black start
unit could somehow incorporate another loop that compensates the effect of the
voltage feedback loop, using available local measurements. In particular, the
voltage control can be expanded as:

$v_{b}=k_{V}(v_{ref}-v_{s-1})+k_{I}i_{s}$ (47)

where the converter also measures the current at its terminal222The virtual
resistance appears as a positive feedback loop in our mathematical formulation
because of the sign convention used for currents in our modelling. and
modifies the voltage accordingly. Intuitively speaking, this can be seen as a
virtual resistance loop, which has been already suggested for black start
purposes in [18], although the focus was on reducing inrush and sympathetic
currents. From an algebraic perspective, this increases the magnitude of one
of the diagonal elements and therefore it is expected to improve the overall
damping. In this case, the virtual resistance modifies the last diagonal term
in the $A$ matrix (that is, the $(2n+1,2n+1)$ term), from $-R/L$ to
$-R/L-k_{I}/L$.

For simplicity, we analyze here just the effect of $k_{I}$. The combined
effect of $k_{V}$ and $k_{I}$ can also be computed using similar strategies.
Following a similar approach as in Section III-C, the characteristic
polynomial under the presence of a virtual resistance can also be derived as
before:

$\displaystyle\Delta_{2n+1}(\lambda)=$
$\displaystyle\left(\lambda+\frac{R}{L}+\frac{k_{I}}{L}\right)\Delta_{2n}(\lambda)+\left(\frac{1}{LC}\right)\Delta_{2n-1}(\lambda)$
$\displaystyle=$
$\displaystyle\left(\lambda+\frac{R}{L}+\frac{k_{I}}{L}\right)\left(P_{n}(g(\lambda))-\frac{1}{LC}P_{n-1}(g(\lambda))\right)$
$\displaystyle+\left(\frac{1}{LC}\right)\left(\lambda+\frac{R}{L}\right)P_{n-1}(g(\lambda))$
$\displaystyle=$
$\displaystyle\left(\lambda+\frac{R+k_{I}}{L}\right)P_{n}(g(\lambda))-\frac{k_{I}}{L^{2}C}P_{n-1}(g(\lambda))$
(48)

Similar expressions can be found in[16]. From this characteristic polynomial,
we can again derive one nonlinear equations of the form:

$\left(\lambda+\frac{R}{L}+\frac{k_{I}}{L}\right)\sin((n+1)\theta)-\frac{k_{I}}{L}\sin(n\theta)=0$
(49)

plus the equation for $\cos(\theta)$ as in (41). In this case the equations
are not partially decoupled as before, since the control gain $k_{I}$ does not
affect exclusively the angle $\theta$, so the derivation becomes more
cumbersome. We focus on the complex eigenvalues, which are the ones of
interest in order to determine the damping in the system. It can be seen in
(49) that the corresponding $\theta$ needs to be also complex for (49) to
hold. The fact that the angle is complex already hints that changes in the
angle will now affect the real part of the eigenvalues.

To proceed, we can carry out a series expansion for (49), given that the
function is also analytic in the complex plane. It can be seen as a function
of $3$ variables, namely $\theta$, $k_{I}$ and $\lambda$ evaluated at the
operating point corresponding to $k_{I_{*}}=0$, with the corresponding values
for $\theta_{*}$ and $\lambda_{*}$ defined in Section III-B:

$\displaystyle
F_{1}=\left(\lambda+\frac{R}{L}+\frac{k_{I}}{L}\right)\sin((n+1)\theta)-\frac{k_{I}}{L}\sin(n\theta)$
$\displaystyle\approx F_{1}(\lambda_{*},k_{I_{*}},\theta_{*})+\frac{\partial
F_{1}}{\partial\theta}\Bigr{|}_{\begin{subarray}{c}\theta_{*}\\\
\lambda_{*}\\\ k_{I_{*}}\end{subarray}}\Delta\theta+\frac{\partial
F_{1}}{\partial k_{I}}\Bigr{|}_{\begin{subarray}{c}\theta_{*}\\\
\lambda_{*}\\\ k_{I_{*}}\end{subarray}}\Delta k_{I}+\frac{\partial
F_{1}}{\partial\lambda}\Bigr{|}_{\begin{subarray}{c}\theta_{*}\\\
\lambda_{*}\\\ k_{I_{*}}\end{subarray}}\Delta\lambda=$
$\displaystyle\left(\left(\lambda_{*}+\frac{R}{L}+\frac{k_{I_{*}}}{L}\right)(n+1)\cos((n+1)\theta_{*})-\frac{k_{I_{*}}}{L}n\cos(n\theta_{*})\right)\Delta\theta$
$\displaystyle+\frac{1}{L}\left(\sin((n+1)\theta_{*})-\sin(n\theta_{*})\right)\Delta
k_{I}+\sin((n+1)\theta_{*})\Delta\lambda=0$ (50)

Taking into account that the operating point satisfies $F_{1}=0$ and therefore
$\sin((n+1)\theta_{*}=0$, we can reach a simple relationship between
variations in the control gain $k_{I}$ and the angle $\theta$:

$\Delta\theta=-\frac{1}{L}\frac{\sin\left(\frac{n}{n+1}\pi
k\right)}{(n+1)\left(\lambda_{*}+\frac{R}{L}\right)}\Delta k_{I}$ (51)

for $k=1,2,\ldots n$. Given that $\lambda_{*}$ is mainly imaginary, and the
rest of the terms are real, it can be seen that a variation in the control
gain leads mainly to a change in the imaginary part of $\theta$. Moreover, the
sensitivity is always positive, that is, the imaginary magnitude of the angle
always increases with $\Delta k_{I}$. Likewise, we can now look at the
relationship between the angle and the eigenvalues from (36):

$\displaystyle
F_{2}=\cos(\theta)-\frac{1}{2}LC\left(\lambda+\frac{R}{L}\right)\left(\lambda+\frac{G}{C}\right)-1$
$\displaystyle\approx F_{2}(\lambda_{*},\theta_{*})+\frac{\partial
F_{2}}{\partial\theta}\Bigr{|}_{\begin{subarray}{c}\theta_{*}\\\
\lambda_{*}\end{subarray}}\Delta\theta+\frac{\partial
F_{2}}{\partial\lambda}\Bigr{|}_{\begin{subarray}{c}\theta_{*}\\\
\lambda_{*}\end{subarray}}\Delta\lambda$
$\displaystyle=-\sin(\theta_{*})\Delta\theta-
LC\left(\lambda_{*}+\frac{1}{2}\left(\frac{G}{C}+\frac{R}{L}\right)\right)\Delta\lambda=0$
(52)

Plugging the value for $\lambda_{*}$ it can be easily seen that an increase of
the imaginary value in $\Delta\theta$ leads to an increase in the real
magnitude of $\Delta\lambda$, and therefore the control gain positively
affects the damping of the system. Unlike in Section III-C, this is a local
result that does not hold for extremely large values of $k_{I}$. This local
result is nonetheless of interest, given the expected range of control gain
values.

The derived results for $k_{V}$ and $k_{I}$ hint at two clear statements for
the behaviour of the grid-forming unit energizing the grid: large values of
$k_{V}$ improve the tracking accuracy of the voltage reference at the cost of
reducing the damping of the resonances, but this effect can be compensated by
increasing the values of $k_{I}$.

Before finalizing this section, we would like to emphasize the potential of
this framework to analyze the effect of other design choices or parameters.
For instance, we can consider a control strategy of the voltage source $v_{b}$
that can also read the current at the other end of the line $i_{1}$. This can
be seen as a communication between two voltage sources attempting at
energizing the same transmission line, in order to improve coordination. In
this case, the element $(2n+1,1)$ is no longer zero but instead $k_{X}/L$, for
some control gain $k_{X}\geq 0$. Using the same formulation as before, the
characteristic polynomial of the system under this controller turns out to be:

$\displaystyle\Delta_{2n+1}(\lambda):=\left(\lambda+\frac{R}{L}\right)P_{n}(g(\lambda))+\frac{k_{X}}{L}\frac{1}{(LC)^{n}}$
$\displaystyle=\left(\lambda+\frac{R}{L}\right)\frac{1}{(LC)^{n}}\frac{\sin((n+1)\theta)}{\sin(\theta)}+\frac{k_{X}}{L}\frac{1}{(LC)^{n}}=0$

This expression can be further analyze the derive the influence of $k_{X}$ on
the real part of the eigenvalues. For instance, it is straightforward to see
that $-R/L$ is no longer an eigenvalue of the system.

Likewise, using this framework the impact of a virtual reactance can be
analysed. The effect of equipping the controller with a virtual reactance is
somewhat similar to the combined effect of a voltage control and a current
control, as expected from the closed loop matrix structure, hence increasing
the real and the imaginary components of the eigenvalues. Other actions that
imply minor changes in the matrix structure have also been analysed using this
framework, for instance topological changes. Preliminary results show that
adding loops to a radial topology does not decrease the damping of the first
resonance in the grid.

##  IV An academic example

To illustrate our results, we consider a $100\text{\,}\mathrm{k}\mathrm{m}$
transmission line with R = $0.02$ $\Omega\mathrm{/}\mathrm{k}\mathrm{m}$, L =
$0.5\cdot 10^{-3}$ $\mathrm{H}\mathrm{/}\mathrm{k}\mathrm{m}$, C = $0.4\cdot
10^{-6}$ $\mathrm{F}\mathrm{/}\mathrm{k}\mathrm{m}$ and G = $0$
$\mathrm{S}\mathrm{/}\mathrm{k}\mathrm{m}$. These values refers to a typical
$110\text{\,}\mathrm{k}\mathrm{V}$ transmission line. The line is divided in 6
$\pi$-section of equal length. Therefore, the system has $13$ states
variables.

We first evaluate the root locus of the system with two voltage sources at
each extreme of the line. For ease of explanation, only $v_{b}$ is equipped
first with a voltage controller and then with a current loop. Fig. 2 and 3
show the evolution of the eigenvalues as a function of the proportional gain
$k_{V}$ and the virtual resistance $k_{I}$ respectively. As expected from the
theoretical results above, independently of the values of $k_{V}$ and $k_{I}$,
the system has a pure real eigenvalue and 5 complex eigenvalues pairs. The
increase of $k_{V}$ increases the imaginary part, thereby reducing the
damping. On the other hand, the increase of $k_{I}$ shifts all the eigenvalue
to the left increasing the damping. It is important to emphasize the different
sensitivities of each eigenvalue: in the case of $k_{V}$, the smallest
eigenvalues (the most relevant ones for our application) do not vary vastly as
a function of the control gain. This could also be deducted from the
expressions in Remark 3 and its conjectured explicit dependence on $k_{V}$:
for $n=6$ and $k=1$ (the first complex eigenvalue), the effect of $k_{V}$ on
the angle $\theta$ is minimal. In the case of $k_{I}$, the intermediate
eigenvalues exhibit larger sensitivity than the other ones. This is however
dependent on the line parameters.

Figure 2:  Root locus of (27) as a function of $k_{V}$ Figure 3:  Root locus
of (27) as a function of $k_{I}$

Similar effects can be seen looking at the frequency response of the system.
Figure 4 displays the frequency response between the voltage source and the
voltage in the intermediate node $v_{3}$. While $k_{V}$ does not reduce the
damping, it does shift the resonances towards higher frequencies. On the other
hand, the effect of $k_{I}$ is to reduce the resonances in the system.

We also show in Fig. 5 the time evolution of the voltage at the third node
$v_{3}$ under different control gains. The reference for both voltage sources
is a $50\text{\,}\mathrm{H}\mathrm{z}$ sine wave333The voltage sources $v_{a}$
and $v_{b}$ are assumed to be perfectly synchronized, given that our modelling
does not include a synchronization algorith.. The controlled voltage source
$v_{b}$ injects a $250\text{\,}\mathrm{H}\mathrm{z}$ component, of a magnitude
equal to $2\%$ of the main $50\text{\,}\mathrm{H}\mathrm{z}$ component. In the
open loop case, the $250\text{\,}\mathrm{H}\mathrm{z}$ component is clearly
amplified, leading to a worse harmonic distortion. The presence of the voltage
controller leads to a lower distortion, but not because of improved damping,
but rather because the resonances are moved up. Nonetheless, this strategy of
shifting the frequencies cannot reliably avoid resonances, given the general
uncertainties in the location of resonances and harmonics injection. On the
other hand, $k_{I}$ reduces the distortion by overall increasing the damping,
leading to a more robust strategy.

Finally, we analyze the impact of the virtual resistance in the energization
of a grid. We consider the well-known 9-bus system, where one inverter
energizes the unloaded grid. Each line is $100$km long and is modelled using
$3$ $\pi$-sections, with the same line parameters as before. Fig. 6 shows the
current at one node when the energization is carried out at $t=0.1s$ in open
loop (that is, the inverter generates a fixed voltage profile without
feedback) and under the presence of a voltage controller and a virtual
resistance as in (47), with $k_{V}=10$ and $k_{I}=100$. The inverter not only
injects the fundamental 50Hz component but also odd harmonics up to the
$9^{\text{th}}$ component, with an amplitude of 1% of the main component. Not
only we see a large overshoot in current at the beginning, but also in steady
state the undamped resonances lead to larger steady state current values. This
implies either an unnecessary overdimensioning of the inverter, or the need to
reduce the initial grid topology to be energized by the inverter. Similar
behaviour has been observed in other nodes across the grid and in the voltage
waveforms, leading to dangerous overvoltages that jeopardize the black start
of the system. The time response of the system with only the voltage
controller active is very similar to the open loop response and thus has been
omitted.

Figure 4: Effect of the controller gains on the frequency response Figure 5:
Time evolution of the voltage under harmonic injection Figure 6: Energization
of the 9-bus system with and without virtual resistance

##  V Outlook

In this paper we have explored the energization of a long transmission line
during black start. In particular, the damping properties and eigenvalues of
the system under different controllers have been analyzed. We hope that the
employed framework, drawing parallelisms with Toeplitz matrices and Chebyshev
polynomials, can be further utilized to understand the impact of other
parameters and design choices, such as the output filters of the energizing
filters (which are known to create resonances at low frequencies), the network
topology, the impact of a virtual reactance, the location of loads [19] or the
presence of transformers at the end of the lines, for instance. The end goal
of this work is twofold: first, help defining a set of requirements for black-
start capable grid-forming inverters; second, construct a set of design rules
for a successful black start strategy relying on a backbone energization.

## References

  * [1] L. Noris, J. Rueda, E. Rakhshani, and A. Korai, “Power system black-start and restoration with high share of power-electronic converters,” in _2019 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2019, pp. 1–5.
  * [2] National Grid, “Demonstration of black start from ders,” Technical report[online]. Accessed: 2024/02/14. Available at: https://www.nationalgrideso.com/future-energy/projects/distributed-restart, 2021.
  * [3] H. Jain, G.-S. Seo, E. Lockhart, V. Gevorgian, and B. Kroposki, “Blackstart of power grids with inverter-based resources,” in _IEEE Power & Energy Society General Meeting (PESGM)_, 2020.
  * [4] J. M. Aniceto, J. Bosch, and A. Anta, “Towards a full automatizated ultrafast black start through inverter-based generation,” in _2023 IEEE Belgrade PowerTech_. IEEE, 2023, pp. 1–6.
  * [5] A. Jain, O. Saborío-Romano, J. N. Sakamuri, and N. A. Cutululis, “Blackstart from hvdc-connected offshore wind: Hard versus soft energization,” _IET Renewable Power Generation_ , vol. 15, no. 1, pp. 127–138, 2021.
  * [6] C. Saldana and G. Calzolari, “Methodology utilized in black-start studies on ehv power networks,” in _IV International Conference on Power System Transients, IPST_ , 2001.
  * [7] Z. Emin, M. Martinez-Duro, M. V. Escudero, R. Adams, H. S. Bronzeado, B. Caillault, N. Chiesa, D. Jacobson, L. Kocis, T. Martinich _et al._ , “Transformer energization in power systems: A study guide,” _CIGRE WG C_ , vol. 4, 2014.
  * [8] J. R. Macias, A. G. Exposito, and A. B. Soler, “A comparison of techniques for state-space transient analysis of transmission lines,” _IEEE Transactions on Power Delivery_ , vol. 20, no. 2, pp. 894–903, 2005.
  * [9] F. Dörfler, J. W. Simpson-Porco, and F. Bullo, “Electrical networks and algebraic graph theory: Models, properties, and applications,” _Proceedings of the IEEE_ , vol. 106, no. 5, pp. 977–1005, 2018.
  * [10] E. Mallada and A. Tang, “Improving damping of power networks: Power scheduling and impedance adaptation,” in _2011 50th IEEE Conference on Decision and Control and European Control Conference_. IEEE, 2011, pp. 7729–7734.
  * [11] T. S. Borsche, T. Liu, and D. J. Hill, “Effects of rotational inertia on power system damping and frequency transients,” in _2015 54th IEEE conference on decision and control (CDC)_. IEEE, 2015, pp. 5940–5946.
  * [12] A. Tayyebi, A. Magdaleno, D. Vettoretti, M. Chen, E. Prieto-Araujo, A. Anta, and F. Dörfler, “System-level performance and robustness of the grid-forming hybrid angle control,” _Electric Power Systems Research_ , vol. 212, p. 108503, 2022.
  * [13] M. J. Gover, “The eigenproblem of a tridiagonal 2-toeplitz matrix,” _Linear Algebra and its Applications_ , vol. 197, pp. 63–78, 1994.
  * [14] H. Al Ba’ba’a, M. Nouh, and T. Singh, “Pole distribution in finite phononic crystals: Understanding bragg-effects through closed-form system dynamics,” _The Journal of the Acoustical Society of America_ , vol. 142, no. 3, pp. 1399–1412, 2017.
  * [15] M. V. Bastawrous and M. I. Hussein, “Closed-form existence conditions for bandgap resonances in a finite periodic chain under general boundary conditions,” _The Journal of the Acoustical Society of America_ , vol. 151, no. 1, pp. 286–298, 2022.
  * [16] C. da Fonseca, “The characteristic polynomial of some perturbed tridiagonal k-toeplitz matrices,” _Appl. Math. Sci_ , vol. 1, no. 2, pp. 59–67, 2007.
  * [17] W. Leterme, E. Heylen, and D. Van Hertem, “On the use of the frequency domain in assessing resonant overvoltages during transformer energization.” _Electric Power Systems Research_ , vol. 196, p. 107164, 2021.
  * [18] A. Jain, O. Saborío-Romano, J. N. Sakamuri, and N. A. Cutululis, “Virtual resistance control for sequential green-start of offshore wind power plants,” _IEEE Transactions on Sustainable Energy_ , vol. 13, no. 3, pp. 1420–1429, 2022.
  * [19] A. Anta and D. Cifelli, “The impact of load placement on grid resonances during grid restoration,” _12th IFAC Symposium on Control of Power & Energy Systems_, 2024.
  * [20] W.-C. Yueh and S. S. Cheng, “Explicit eigenvalues and inverses of tridiagonal toeplitz matrices with four perturbed corners,” _the ANZIAM Journal_ , vol. 49, no. 3, pp. 361–387, 2008.
  * [21] C. M. da Fonseca, “Comments on the spectrum of a tridiagonal k-toeplitz matrix,” _Journal of Computational and Applied Mathematics_ , vol. 375, p. 112793, 2020.
  * [22] R. Álvarez-Nodarse, J. Petronilho, and N. R. Quintero, “On some tridiagonal k-toeplitz matrices: Algebraic and analytical aspects. applications,” _Journal of computational and applied mathematics_ , vol. 184, no. 2, pp. 518–537, 2005.
  * [23] A.-J. Nikkilä, A. Kuusela, M. Weixelbraun, L. Haarla, M. Laasonen, and A. Pahkin, “Fast restoration of a critical remote load area using a gradual voltage build-up procedure,” _IET Generation, Transmission & Distribution_, vol. 14, no. 7, pp. 1320–1328, 2020.