# Interplay of electron correlations, spin-orbit couplings, and structural
effects for Cu centers in the quasi-two-dimensional magnet InCu2/3V1/3O3

R. Murugesan  Leibniz IFW Dresden, Helmholtzstr. 20, 01069 Dresden, Germany
M. S. Eldeeb  Leibniz IFW Dresden, Helmholtzstr. 20, 01069 Dresden, Germany
M. Yehia  Leibniz IFW Dresden, Helmholtzstr. 20, 01069 Dresden, Germany
Reactor Physics Department, Nuclear Research Center, Atomic Energy Authority,
Cairo, Egypt  B. Büchner  Leibniz IFW Dresden, Helmholtzstr. 20, 01069
Dresden, Germany  Institute for Solid State and Materials Physics and
Würzburg-Dresden Cluster of Excellence ct.qmat, TU Dresden, D-01062 Dresden,
Germany  V. Kataev  Leibniz IFW Dresden, Helmholtzstr. 20, 01069 Dresden,
Germany  O. Janson  Leibniz IFW Dresden, Helmholtzstr. 20, 01069 Dresden,
Germany  L. Hozoi  Leibniz IFW Dresden, Helmholtzstr. 20, 01069 Dresden,
Germany

###### Abstract

Less common ligand coordination of transition-metal centers is often
associated with peculiar valence-shell electron configurations and outstanding
physical properties. One example is the Fe+ ion with linear coordination,
actively investigated in the research area of single-molecule magnetism. Here
we address the nature of 3$d^{9}$ states for Cu2+ ions sitting in the center
of trigonal bipyramidal ligand cages in the quasi-two-dimensional honeycomb
compound InCu2/3V1/3O3, whose unusual magnetic properties were intensively
studied in the recent past. In particular, we discuss the interplay of
structural effects, electron correlations, and spin-orbit couplings in this
material. A relevant computational finding is a different sequence of the Cu
($xz$, $yz$) and ($xy$, $x^{2}\\!-\\!y^{2}$) levels as compared to existing
electronic-structure models, which has implications for the interpretation of
various excitation spectra. Spin-orbit interactions, both first- and second-
order, turn out to be stronger than previously assumed, suggesting that rather
rich single-ion magnetic properties can be in principle achieved also for the
3$d^{9}$ configuration by properly adjusting the sequence of crystal-field
states for such less usual ligand coordination.

##  I Introduction

Transition-metal (TM) ions with atypical ligand (L) coordination in TM
pnictides, chalcogenides, and halides may display quite unexpected features as
concerns their electronic structures and physical properties. A remarkable
example is monovalent Fe with linear L-Fe-L bonds, which due to its
impressively strong magnetic anisotropy is being presently investigated in the
context of single-molecule magnetism [1; 2]. Also rare is for instance the
square coordination realized for nickel (formally 3$d^{9}$) in NdNiO2.
Interestingly, superconductivity with $T_{\mathrm{c}}$’s of up to 15 K was
recently reported in this system [3], the available data suggesting physics
significantly different as compared to the Cu-oxide superconductors [4; 5].

The title compound of the present study, honeycomb InCu2/3V1/3O3, has
attracted significant attention due to its peculiar magnetic properties [6; 7;
8; 9; 10; 11; 12]. In particular, it has been found recently [13] that a
pronounced two-dimensional (2D) $XY$ anisotropy of the Cu-based honeycomb
spin-1/2 planes and strong inter-layer frustration yield an extended low-
temperature regime where a quasi-2D static magnetic state is stabilized,
showing some indications of a Berezinsky-Kosterlitz-Thouless transition to the
topologically ordered state of vortex-antivortex pairs predicted in the 2D XY
model [14; 15; 16; 17]. An interesting aspect is the uncommon type of ligand
environment encountered in this compound, with trigonal ligand bipyramids
around each copper ion. For $D_{3h}$ Cu-site symmetry as realized in this
situation, there are two sets of doubly degenerate orbitals ($xz$, $yz$ and
$xy$, $x^{2}\\!-\\!y^{2}$), as for magnetically anisotropic L-TM-L units with
either $D_{6h}$ [2; 18; 19; 20] or approximate $D_{3h}$ [1] symmetry. Finding
that spin-orbit (SO) interactions play an important role for such hole
configurations in various L-TM-L complexes motivates then a detailed
investigation of SO interactions for Cu 3$d^{9}$ hole states in InCu2/3V1/3O3.
Even if the $z^{2}$ ground-state character inferred from analysis of electron
spin resonance (ESR) spectra [6] and by band-structure calculations [21; 13]
indicates that to first order only the excited states would be affected by SO
couplings in this particular cuprate, such a study should still provide
insights into the prospect of strong single-ion magnetic anisotropy for
3$d^{9}$ centers with suitably engineered environments in related systems. In
addition to first-order SO interactions, here we also analyze in detail
second-order effects, in particular, their impact on the ground-state $g$
factors and on the structure of the on-site excitation spectrum.

##  II Ground-state $\mathbf{g}$ factors, experiment vs ab initio computations

Within the family of layered Cu-oxide compounds, the type of $g$-factor
anisotropy featured by Cu centers in InCu2/3V1/3O3, $g_{ab}\\!>\\!g_{c}$ , is
unusual since for the large majority of cuprates $g_{ab}\\!<\\!g_{c}$ is
commonly found ($g_{ab}$ and $g_{c}$ are here the in-plane and out-of-plane
components of the ${\bf g}$ tensor, respectively). This aspect was first
noticed in Ref. [6], confirmed in Ref. [9], and attributed to a less common
ground-state configuration having one hole in the $z^{2}$ Cu 3$d$ orbital. To
better substantiate these earlier findings and interpretations, we perform in
this work detailed electronic-structure calculations for this material. To
illustrate and further document the experimental basis and to facilitate the
discussion of the computational results we first present additional ESR data
measured on a powder sample of InCu2/3V1/3O3, synthesized and thoroughly
characterized in Ref. [7]. The powder was mixed with an epoxy resin and the
mixture was hardened in a magnetic field of several Tesla. As a result, the
so-prepared sample acquired a well-defined anisotropy axis (referred to
hereafter as orientation $o$ axis) [6] which according to x-ray diffraction
measurements is perpendicular to the crystallographic $c$ axis.

Figure 1:  Angular dependence of the $g$ factor (a) and of the peak-to-peak
ESR linewidth $\Delta H_{\rm pp}$ (b) of the oriented powder sample of
InCu2/3V1/3O3 measured at room temperature at the X-band frequency
$\nu\\!=\\!9.58$ GHz (symbols). The inset shows a typical ESR signal (field
derivative of the microwave absorption). Dashed and solid curves are the
results of numerical modeling (see main text).

The dependence of the $g$ factor $g=h\nu/\mu_{\rm B}\mu_{0}H_{\rm res}$ and of
the ESR peak-to-peak linewidth $\Delta H_{\rm pp}$ on the angle $\alpha$ which
the applied magnetic field ${\bf H}$ makes with the $o$ axis are shown in Fig.
1, as measured with a commercial X-band Bruker spectrometer at a frequency
$\nu=9.59$ GHz and room temperature. $h$, $\mu_{\rm B}$, and $\mu_{0}$ are
here the Planck constant, the Bohr magneton, and the vacuum permeability,
respectively. The parameters of the ESR signal $H_{\rm res}$ (resonance field)
and $\Delta H_{\rm pp}$ are defined in the inset of Fig. 1. The $g$ factor for
$\alpha=0$ corresponds to the single-crystalline in-plane $g$ factor $g_{ab}$,
whereas for $\alpha=\pm 90^{\circ}$ it yields the powder averaged $g$ factor
$g_{\rm pow}$. To model the experimental dependence of $g$ and $\Delta H_{\rm
pp}$ on $\alpha$ we performed numerical simulations of the ESR signal which
account for the $\alpha$-dependent averaging of the contributions from the
individual crystallites in the sample. The $g(\alpha)$ dependence can be very
well reproduced, yielding the single-crystalline $g$ values $g_{c}=2.020\pm
0.015$ and $g_{ab}=2.230\pm 0.015$ [dashed line in Fig. 1(a)]. The $g$-factor
anisotropy should give rise to a dependence of the measured peak-to-peak
linewidth $\Delta H_{\rm pp}$ on the angle $\alpha$. For $\alpha\\!=\\!0$, the
ESR signal has the width corresponding to the single-crystalline line $\Delta
H$, $\Delta H_{\rm pp}(\alpha\\!=\\!0)\\!=\\!\Delta H$. For $\alpha\\!=\\!\pm
90^{\circ}$, the ESR linewidth $\Delta H_{\rm pp}(\alpha\\!=\\!\pm
90^{\circ})$ corresponds to that of the powder sample and should be larger
because of the distribution of the resonance fields of the individual
crystallites arising from the $g$-factor anisotropy, $\Delta H_{\rm
pp}(\alpha\\!=\\!\pm 90^{\circ})\\!>\\!\Delta H$. However, the respective
$\Delta H_{\rm pp}(\alpha)$ modeled curve noticeably deviates from experiment
[dashed line in Fig. 1(b)], indicating that the $\Delta H_{\rm pp}(\alpha)$
dependence is not entirely due to averaging signals with anisotropic $g$
factors. Therefore, in addition, one could also consider a possible angular
dependence of the width of individual single-crystalline ESR lines of the form
$\Delta H=A+B\cdot(1+\cos^{2}\theta)$, where $\theta$ denotes the angle
between ${\bf H}$ and the $c$ axis. Here, the first term is the linewidth in
the absence of spin correlations at $T\rightarrow\infty$ and the second term
is the angular dependent contribution due to the three-dimensional (3D)
antiferromagnetic correlations at finite temperature [22; 23; 24; 25; 26].
Indeed, accounting for this additional contribution in the simulation, with
$A=571$ Oe and $B=45$ Oe, significantly improves the agreement with the
experiment [solid line in Fig. 1(b)]. The 3D-like behavior of InCu2/3V1/3O3 at
elevated temperature is not surprising since the in-plane spin correlations
turning this compound into a 2D-$XY$ magnet begin to grow at temperatures
$\lesssim 50$ K [13].

Finally, a more accurate value can be obtained for $g_{ab}$ from the frequency
$\nu$ versus field $H$ dependence of the ESR signal, which we measured using
the home-made high-field ESR setup described in [8] at a lower temperature
$T\\!=\\!60$ K (Fig. 2). A fit of the data to the resonance condition
$\nu=g_{\rm ab}\mu_{\rm B}\mu_{0}H/h$ yields $g_{ab}=2.25\pm 0.01$. With
$g_{\rm pow}=2.23$, one obtains from the relation $g_{\rm pow}=(1/3)g_{\rm
c}+(2/3)g_{ab}$ [6] the $c$-axis $g$ factor $g_{c}=1.99$. The resulting ${\bf
g}$ tensor $(g_{ab},g_{c})=(2.25,1.99)$ is, within error bars of $\pm$5–10%,
fully consistent with values obtained for other samples of InCu2/3V1/3O3 from
various sources and for batches earlier studied in Refs. [6; 9], providing
altogether a solid reference for theoretical analysis.

Figure 2:  Frequency $\nu$ versus magnetic field $H$ dependence of the ESR
signal of the oriented powder sample of InCu2/3V1/3O3 measured at $T=60$ K
with $\mathbf{H}\parallel o$ axis (open circles). A corresponding fit is shown
as dashed line. Exemplary ESR signals are shown for frequencies of 160 and 345
GHz. The asterisk labels a small additional peak presumably due to Cu2+ spins
at defect sites.

Two complementary computational techniques were employed for addressing the Cu
3$d$-shell electronic structure in InCu2/3V1/3O3 and peculiar features such as
the $g_{ab}\\!>\\!g_{c}$ anisotropy of the $g$ factor. In a first step we
determined the optimal atomic positions in the frame of density-functional
calculations with periodic boundary conditions. We used to this end the full-
potential local-orbital code fplo, version 18 [27]. For nonmagnetic scalar-
relativistic as well as fully relativistic computations, we utilized the
generalized gradient approximation (GGA) [28]. A $k$-mesh of
14$\times$14$\times$7 points was employed in either case. Typical for undoped
cuprates, GGA yields a spurious metallic state in InCu2/3V1/3O3. This well-
known artifact stems from underestimating electronic correlations with
conventional functionals such as GGA. The insulating state can be restored in
DFT++ approaches, where interactions are considered explicitly on the level of
(extended) multi-orbital Hubbard models [29]. Nevertheless, the non-
interacting part of such models is adopted from conventional DFT calculations.
Hence, the respective model parameters, transfer integrals and crystal-field
splittings, are not affected by shortcomings of the GGA.

InCu2/3V1/3O3 is a quasi-2D material whose honeycomb planes are only weakly
coupled to each other. All samples synthesized so far feature a sizable number
of stacking faults, impeding a complete characterization of the crystal
structure [7]. Hence, the first step of our computational analysis was the
construction of a 3D structural model, based on the experimentally available
information: lattice constants and the internal coordinates of a single
honeycomb plane with Cu and V occupying the same crystallographic site [7]. As
explained in Ref. [13], a minimal 3D structural model for InCu2/3V1/3O3
employs the orthorhombic space group $Cmcm$ (63), which is compatible with a
fully ordered honeycomb lattice of Cu atoms and at the same time supports the
simplest type of stacking of the honeycomb planes. The respective atomic
coordinates are provided in Table 1. The Cu-site point group symmetry is in
this space group $C_{s}$ instead of $D_{3h}$.

Table 1:  Atomic sites, Wyckoff positions (Wyck. pos.), and internal coordinates of the two different crystal structures used in this study. The space group is $Cmcm$ (63), with $a$ = 10.0527, $b$ = 5.80393, $c$ = 11.9012 Å. In the computationally optimized structure, the O positions were relaxed within the antiferromagnetic (Néel) structure to yield the lowest total energy at the GGA+$U$ level. Site | Wyck. |  | Identical O5 cages |  | Optimized O positions
---|---|---|---|---|---
pos. |  | $x/a$ | $y/b$ | $z/c$ |  | $x/a$ | $y/b$ | $z/c$
Cu | $8g$ |  | $\nicefrac{{1}}{{3}}$ | $\nicefrac{{1}}{{3}}$ | 0.25 |  |  |  |
V | $4c$ |  | 0 | $\nicefrac{{1}}{{3}}$ | 0.25 |  |  |  |
In1 | $4a$ |  | 0 | 0 | 0 |  |  |  |
In2 | $8e$ |  | $\nicefrac{{2}}{{3}}$ | 0 | 0 |  |  |  |
O1 | $4c$ |  | 0 | 0 | 0.25 |  | 0 | 0.03400 | 0.25
O2 | $8g$ |  | $\nicefrac{{1}}{{3}}$ | 0 | 0.25 |  | 0.35045 | 0.98309 | 0.25
O3 | $8f$ |  | 0 | $\nicefrac{{2}}{{3}}$ | 0.5864 |  | 0 | 0.66708 | 0.58500
O4 | $16h$ |  | $\nicefrac{{1}}{{3}}$ | $\nicefrac{{1}}{{3}}$ | 0.0864 |  | 0.33313 | 0.33359 | 0.09067

The parent structural model [7] does not distinguish between Cu and V,
featuring identical O5 cages around Cu and V atoms. However, the honeycomb
planes in InCu2/3V1/3O3 may undergo a breathing distortion, giving rise to
unequal Cu-O and V-O bond lengths. To address this possibility, we performed a
structural relaxation of O positions within the space group $P2_{1}/m$,
supporting Néel antiferromagnetic order. To avoid inaccuracies stemming from
the spurious metallic ground state, the optimization was carried out with
respect to the GGA+$U$ energy. Following earlier studies on low-dimensional
cuprates [30; 31], we use a Coulomb repulsion parameter $U$ = 8.5 eV, a Hund’s
exchange $J$ = 1 eV, and the fully localized limit for the double-counting
correction. We find this way large differences between the in-plane Cu-O and
V-O bond lengths, 2.04 vs 1.74 Å, consistent with considerations concerning
the difference between the formal oxidation states of copper and vanadium in
this system [7]. Also, the computationally derived lattice parameters imply
significantly longer in-plane Cu-O links as compared to the apical Cu-O bonds
(2.04 vs 1.90 Å).

Using crystallographic data as listed in Table 1, we further performed many-
body quantum chemical calculations [32] on a finite atomic fragment having a
CuO5 unit as central region (all relevant computational details are provided
in Appendix). To describe configurational mixing effects mediated by SO
interactions and related $g$-factor anisotropies we relied on spin-orbit
multiconfiguration and multireference numerical schemes as implemented in the
quantum chemical package molpro [33] by Berning et al. [34]. An active space
defined by the Cu 3$d$ orbitals (nine electrons in five orbitals) was used to
this end in preliminary multiconfiguration computations [32]. All possible
states associated with this 3$d$-shell filling were considered, which in
quantum chemical terminology is referred to as complete-active-space self-
consistent-field (CASSCF) methodology 111 However, for the 3$d^{9}$ manifold
and sufficiently high symmetry, this essentially corresponds to single-
configuration restricted open-shell Hartree-Fock (ROHF) [32]. . The $z^{2}$
character of the Cu 3$d$ hole proposed on the basis of the measured
$g_{ab}\\!>\\!g_{c}$ structure of the $g$ factors [6] is also confirmed in the
CASSCF calculation. For making direct connection with the experimental data,
$g$ factors were also computed: using CASSCF wave-functions, related angular-
momentum and SO matrix elements [34], and the methodology described in Ref.
[36], we arrive to $g_{ab}\\!=\\!2.23$ and $g_{c}\\!=\\!1.96$; on the basis of
multireference configuration-interaction (MRCI) wave-functions built by
additionally considering single and double excitations [37] out of the Cu 3$d$
and O 2$p$ orbitals of the CuO5 polyhedron, we obtain $g_{ab}\\!=\\!2.21$ and
$g_{c}\\!=\\!1.97$. These CASSCF and MRCI values are within 98–99$\%$ of the
estimates derived from the ESR data and presented above. For a more detailed
picture on the Cu 3$d$-shell electronic structure in this compound, we analyze
in the following the Cu $d$-$d$ excitations since the $g$-factor values depend
on the relative energies of the excited states.

Table 2:  Cu2+ 3$d^{9}$ multiplet structure (relative energies in eV) using atomic positions as obtained by density-functional optimization. The MRCI values include Davidson corrections [32]. The character of the SO wave-functions is also specified.  $3d^{9}$ states | CASSCF | MRCI | SO-MRCI
---|---|---|---
${}^{2}\\!A_{1}^{\prime}$ | 0 | 0 | 0 (99$\%$ $A_{1}^{\prime}$ character)
${}^{2}\\!E^{\prime\prime}$ | 1.03 | 1.20 | 1.12 (89$\%$ $E^{\prime\prime}$, 11$\%$ $E^{\prime}$)
|  |  | 1.21 (99$\%$ $E^{\prime}$)
${}^{2}\\!E^{\prime}$ | 1.12 | 1.30 | 1.28 (99$\%$ $E^{\prime\prime}$)
|  |  | 1.46 (89$\%$ $E^{\prime}$, 11$\%$ $E^{\prime\prime}$)

##  III Excited states

CASSCF, MRCI, and spin-orbit MRCI (SO-MRCI) results for the Cu2+ multiplet
structure using atomic positions as determined by density-functional lattice
optimization are listed in Table 2. Notations corresponding to $D_{3h}$ point-
group symmetry are employed, where the Cu $z^{2}$, ($xy$,
$x^{2}\\!-\\!y^{2}$), and ($xz$, $yz$) orbitals belong to the irreducible
representations $A_{1}^{\prime}$, $E^{\prime}$, and $E^{\prime\prime}$,
respectively. It is found that the two sets of doubly degenerate levels lie
rather close in energy and that the MRCI treatment brings corrections of
0.1–0.2 eV to the CASSCF crystal-field splittings. Rectifications in this
range were previously computed at the MRCI level in other copper oxide
compounds [38; 39] and mainly originate from L 2$p$ – TM 3$d$ charge-transfer-
type correlation effects. SO interactions lift the degeneracy of the
$E^{\prime}$ and $E^{\prime\prime}$ states and additionally lead to some
degree of $E^{\prime}\\!-\\!E^{\prime\prime}$ mixing. This is an aspect that
deserves attention: since $A_{1}^{\prime}$ to $E^{\prime\prime}$ dipole
transitions are not allowed, the three-peak structure of the optical
absorption spectrum [7] can be then qualitatively understood on the basis of
the fact that only three of the four SO excited states have significant
$E^{\prime}$ character (see the right column in Table 2).

The $E^{\prime}\\!-\\!E^{\prime\prime}$ mixing is obviously the effect of
second-order SO couplings. Such physics was not addressed in detail in earlier
analysis of the Cu $d$–$d$ excitation spectrum of InCu2/3V1/3O3 but given the
near-degeneracy of the $E^{\prime}$ and $E^{\prime\prime}$ crystal-field-like
states it plays a significant role. To better illustrate this point, we
performed SO calculations for each of the $E^{\prime}$ and $E^{\prime\prime}$
terms separately. The splittings induced by first-order SO interactions came
out as 200 meV and 101 meV, respectively. Taking as reference the MRCI
relative energies provided in Table 2 (1.20 and 1.30 eV) and neglecting
second-order SO couplings, these first-order splittings translate into
excitation energies of 1.15, 1.25 eV ($E^{\prime\prime}$ terms) and 1.20, 1.40
eV ($E^{\prime}$ terms) for the SO excited states, different from the results
of the full SO-MRCI computation (1.12, 1.21, 1.28, and 1.46 eV). Compared to
earlier calculations based on the angular overlap model (AOM) [7], our ab
initio quantum chemical results differ therefore in two aspects: a different
sequence of the $E^{\prime}$ and $E^{\prime\prime}$ terms in the absence of
first-order SO interactions
($E(^{2}\\!E^{\prime})\\!>\\!E(^{2}\\!E^{\prime\prime})$ in Table 2 while
$E(^{2}\\!E^{\prime})\\!<\\!E(^{2}\\!E^{\prime\prime})$ in the AOM model [7])
and significantly stronger first-order SO interactions for the
$E^{\prime\prime}$ states (split in first order by 101 meV here, see above,
and by 62 meV (500 cm-1) in the AOM model [7]).

As concerns a direct comparison between the relative energies of the SO states
with significant $E^{\prime}$ character listed in Table 2 and the positions of
the peaks in the optical absorption spectrum [7], it is seen that the
agreement is not perfect: 1.12, 1.21, and 1.46 eV vs 1.18, 1.35, and 1.60 eV
(9500, 10900, and 12900 cm-1 in Ref. [7]). Deviations of $\sim$0.1 eV between
MRCI excitation energies and experimental peak positions, with MRCI constantly
underestimating experimental values, are however common for Cu2+ 3$d^{9}$
oxide compounds [38; 39].

For completeness, we additionally provide in Table 3 on-site excitation
energies computed on the basis of the more idealized crystal structure with
identical MO5 cages. The crystal-field $a_{1}^{\prime}\\!-\\!e^{\prime}$ and
$a_{1}^{\prime}\\!-\\!e^{\prime\prime}$ splittings (which we denote as
$\delta^{\prime}$ and $\delta^{\prime\prime}$, respectively) are different in
this case, with $\delta^{\prime}\\!<\\!\delta^{\prime\prime}$. A reversed
sequence of the $e^{\prime}$ and $e^{\prime\prime}$ levels in the two
different crystal structures is also found by density-functional theory, by
Wannier projection of the GGA bands onto the five 3$d$ orbitals at a Cu site,
as discussed in Appendix. Most importantly, with SO couplings accounted for,
the calculated excitation energies are far from the peak positions found in
optics (see Table 3). In particular, ‘pairs’ of nearly-degenerate SO states
are obtained computationally for the more idealized crystal structure,
different from the three-peak structure observed experimentally [7]. A more
consistent picture in modeling experimental data was also found in the AOM
frame [7] when elongated in-plane Cu-O bonds are employed.

Table 3:  Results for the Cu2+ 3$d^{9}$ multiplet structure (relative energies in eV) using the idealized lattice configuration proposed in Ref. [7]. The MRCI values include Davidson corrections [32]. The character of the SO wave-functions is also provided.  $3d^{9}$ states | CASSCF | MRCI | SO-MRCI
---|---|---|---
${}^{2}\\!A_{1}^{\prime}$ | 0 | 0 | 0 (99$\%$ $A_{1}^{\prime}$ character)
${}^{2}\\!E^{\prime}$ | 0.98 | 1.17 | 1.08 (99$\%$ $E^{\prime}$)
|  |  | 1.12 (72$\%$ $E^{\prime\prime}$, 28$\%$ $E^{\prime}$)
${}^{2}\\!E^{\prime\prime}$ | 1.03 | 1.22 | 1.30 (99$\%$ $E^{\prime\prime}$)
|  |  | 1.34 (72$\%$ $E^{\prime}$, 28$\%$ $E^{\prime\prime}$)

##  IV Conclusions

In sum, a detailed analysis of the Cu $d$-shell electronic structure in the
quasi-2D honeycomb compound InCu2/3V1/3O3 is performed, with focus on the
interplay of electron correlations, spin-orbit couplings, and structural
effects in this material. Ab initio computational data are compared to the
outcome of ESR experiments and of previous optical absorption measurements
[7]. The three-peak structure of the optical spectra is reproduced by many-
body quantum chemical calculations only for atomic positions obtained from a
prior lattice optimization relying on density-functional theory and not for a
more idealized crystal structure with identical MO5 cages. It seems to be
related with the fact that $E^{\prime}\\!-\\!E^{\prime\prime}$ mixing occurs
for only two of the spin-orbit excited states. To describe electron
correlation effects, we relied on multireference configuration-interaction
calculations with single and double excitations out of the copper 3$d$ and
oxygen 2$p$ shells within a given CuO5 unit. Although the correlation-induced
corrections obtained this way to the $d$–$d$ excitation energies are sizable,
0.1–0.2 eV, deviations in the range of 0.1 eV still remain as compared to
experimental peak positions. On the other hand, the computed ground-state $g$
factors are in excellent agreement with values derived from all available ESR
data, within 98–99$\%$. The peculiar coordination generates strong axial
anisotropy at copper sites in InCu2/3V1/3O3, reminiscent to some exent of, e.
g., linearly coordinated transition-metal ions [1; 2; 18; 19; 20]. Adjusting
the sequence of crystal-field levels through adequate design of the
chemical/electrostatic environment [36; 40] such that the $z^{2}$ orbital is
filled leaves room for interesting single-ion magnetic properties even for the
3$d^{9}$ electron configuration, as pointed out in Refs. [19; 20].

## Acknowledgements

We thank A. Jesche for useful discussions, U. Nitzsche for technical
assistance, and A. Möller for providing a sample of InCu2/3V1/3O3 for the ESR
experiment and for important comments to the manuscript. L. H. acknowledges
financial support from the German Research Foundation (Deutsche
Forschungsgemeinschaft, DFG), grant HO-4427/3. O. J. was supported by the
Leibniz Association through the Leibniz Competition.

##  Appendix A Computational details

To evaluate crystal-field parameters and the SO coupling constant $\zeta$
within the GGA, we performed Wannier projections of the GGA bands onto the
five 3$d$ orbitals at a Cu site. As these orbitals hybridize with oxygen $2p$
and vanadium $3d$ orbitals, the resulting hoppings sensibly depend on the
energy window used for projections. Moreover, the Fourier-transformed Wannier
functions show deviations from the GGA bands at the bottom of the valence
band, especially for the experimental structure. However, we are interested
exclusively in the local terms, and these show only a small (on the order of
several meV) dependence on the energy window. Using parameters obtained from
such a procedure, we construct the local Hamiltonian $H_{0}$. Corresponding
eigenvalues are provided for both crystal structures in Table 4.

For both lattice configurations, $H_{0}$ features two pairs of nearly
degenerate eigenstates: $\lambda_{1}\simeq\lambda_{2}$ and
$\lambda_{3}\simeq\lambda_{4}$. But the gaps between them,
$\Delta_{1}\\!=\\!\lambda_{3}\\!-\\!\lambda_{2}$ and
$\Delta_{2}\\!=\\!\lambda_{5}\\!-\\!\lambda_{4}$, differ significantly in
InCu2/3V1/3O3: $\Delta_{1}\\!\simeq\\!0.30$ (0.19) eV is substantially smaller
than $\Delta_{2}\\!\simeq\\!1.06$ (1.10) eV in the idealized (optimized)
structure. As concerns the sequence of the two groups of nearly degenerate
states, we found that in the more idealized crystal structure $\lambda_{1}$
has dominant (84%) $xz$ character while $\lambda_{2}$ is dominated (84%) by
the $yz$ orbital. Similarly, $\lambda_{3}$ and $\lambda_{4}$ pertain to the
$xy$ and $x^{2}\\!-\\!y^{2}$ orbitals, respectively. But in the optimized
structure, the order is reversed: the lower-lying eigenstates in each group,
$\lambda_{1}$ and $\lambda_{3}$, are dominated by $yz$ (94%) and
$x^{2}\\!-\\!y^{2}$ (61%) orbitals.

Table 4:  Eigenvalues $\lambda_{n}$ of the local Hamiltonian $H_{0}$ parameterized using Wannier projection. In the fully relativistic calculations (GGA+SOC), all eigenvalues are doubly degenerate. All values are given in eV with respect to the Fermi level.  Eigenvalues | Identical O5 cages | Optimized O positions
---|---|---
of $H_{0}$ | GGA | GGA+SOC | GGA | GGA+SOC
$\lambda_{1}$ | $-1.61$ | $-1.65$ | $-1.44$ | $-1.48$
$\lambda_{2}$ | $-1.60$ | $-1.60$ | $-1.43$ | $-1.46$
$\lambda_{3}$ | $-1.31$ | $-1.37$ | $-1.25$ | $-1.29$
$\lambda_{4}$ | $-1.31$ | $-1.24$ | $-1.25$ | $-1.17$
$\lambda_{5}$ | $-0.26$ | $-0.26$ | $-0.15$ | $-0.15$

Finally, we applied a fully relativistic GGA treatment and performed Wannier
projections to estimate the SO coupling constant $\zeta$. The respective
eigenvalues are provided in Table 4 next to those for the scalar-relativistic
case 222Note that the fully relativistic calculations are nonmagnetic and all
eigenvalues are doubly degenerate.. The SO coupling $\zeta$ can be then
estimated by averaging over two offdiagonal elements of $H_{0}$:

$\zeta=-\frac{1}{2\sqrt{3}}\left(i\langle{}\phi_{z^{2}}^{\uparrow}|H_{0}|\phi_{yz}^{\downarrow}\rangle+\langle{}\phi_{z^{2}}^{\uparrow}|H_{0}|\phi_{xz}^{\downarrow}\rangle\right),$
(1)

where $\phi_{\alpha}^{\sigma}$ is the basis (Wannier) function of orbital
$\alpha$ and spin $\sigma$. We obtain this way $\zeta$ = 68(60) meV for the
idealized (optimized) structure.

As concerns the quantum chemical calculations, the actual cluster considered
in the computations consists of a central CuO5 unit along with three adjacent
VO5 polyhedra, three Cu nearest neighbors, and six adjacent In ions. The
remaining part of the extended crystalline surroundings was modeled as an
effective electrostatic field [42]. For the central CuO5 unit, all-electron
Douglas-Kroll basis sets of triple-$\zeta$ quality with polarization functions
were used. For farther oxygen ligands in the cluster and for the nearby V
ions, we applied basis sets of double-$\zeta$ quality. The results discussed
in the main text were obtained by modeling the Cu2+ nearest neighbors as
closed-shell Zn2+ total-ion potentials (TIP’s), an approximation also employed
in earlier quantum chemical investigations [38; 39; 43]. Test CASSCF
calculations in which the three adjacent Cu species are represented as 2+
$S\\!=\\!1/2$ 3$d^{9}$ ions provide central-site crystal-field splittings that
agree within 0.01 eV with those computed with Zn2+ TIP’s. We used the latter
for all subsequent computations since this makes the spin-orbit calculations
and related analysis much less complex. TIP’s were also applied for the nearby
In3+ ions.

## References

  * Zadrozny  _et al._ [2013] J. M. Zadrozny, D. J. Xiao, M. Atanasov, G. J. Long, F. Grandjean, F. Neese, and J. R. Long, Nat. Chem. 5, 577 (2013).
  * Jesche  _et al._ [2014] A. Jesche, R. W. McCallumm, S. Thimmaiah, J. L. Jacobs, V. Taufour, A. Kreyssig, R. S. Houk, S. L. Bud’ko, and P. C. Canfield, Nat. Commun. 5, 3333 (2014).
  * Li  _et al._ [2019] D. Li, K. Lee, B. Y. Wang, M. Osada, S. Crossley, H. R. Lee, Y. Cui, Y. Hikita, and H. Y. Hwang, Nature 572, 624 (2019).
  * Hepting  _et al._ [2020] M. Hepting, D. Li, C. J. Jia, H. Lu, E. Paris, Y. Tseng, X. Feng, M. Osada, E. Been, Y. Hikita, Y.-D. Chuang, Z. Hussain, K. J. Zhou, A. Nag, M. Garcia-Fernandez, M. Rossi, H. Y. Huang, D. J. Huang, Z. X. Shen, T. Schmitt, H. Y. Hwang, B. Moritz, J. Zaanen, T. P. Devereaux, and W. S. Lee, Nat. Mater. 19, 381 (2020).
  * Jiang  _et al._ [2020] M. Jiang, M. Berciu, and G. A. Sawatzky, Phys. Rev. Lett. 124, 207004 (2020).
  * Kataev  _et al._ [2004] V. Kataev, A. Möller, U. Löw, W. Jung, N. Schittner, M. Kriener, and A. Freimuth, J. Magn. Magn. Mater. 290-291, 310 (2004).
  * Möller  _et al._ [2008] A. Möller, U. Löw, T. Taetz, M. Kriener, G. André, F. Damay, O. Heyer, M. Braden, and J. A. Mydosh, Phys. Rev. B 78, 024420 (2008).
  * Yehia  _et al._ [2010] M. Yehia, E. Vavilova, A. Möller, T. Taetz, U. Löw, R. Klingeler, V. Kataev, and B. Büchner, Phys. Rev. B 81, 060414 (2010).
  * Okubo  _et al._ [2011] S. Okubo, H. Wada, H. Ohta, T. Tomita, M. Fujisawa, T. Sakurai, E. Ohmichi, and H. Kikuchi, J. Phys. Soc. Jpn. 80, 023705 (2011).
  * Yan  _et al._ [2012] Y. J. Yan, Z. Y. Li, T. Zhang, X. G. Luo, G. J. Ye, Z. J. Xiang, P. Cheng, L. J. Zou, and X. H. Chen, Phys. Rev. B 85, 085102 (2012).
  * Jia  _et al._ [2017a] S.-Q. Jia, Q.-W. Wang, J. Liu, and L.-J. Zou, J. Phys.: Conf. Ser. 827, 012010 (2017a).
  * Jia  _et al._ [2017b] S.-Q. Jia, Q.-W. Wang, X.-L. Yu, and L.-J. Zou, AIP Advances 7, 055825 (2017b).
  * Iakovleva  _et al._ [2019] M. Iakovleva, O. Janson, H.-J. Grafe, A. P. Dioguardi, H. Maeter, N. Yeche, H.-H. Klauss, G. Pascua, H. Luetkens, A. Möller, B. Büchner, V. Kataev, and E. Vavilova, Phys. Rev. B 100, 144442 (2019).
  * Berezinskii [1971] V. L. Berezinskii, Sov. Phys JETP-USSR 32, 493 (1971).
  * Kosterlitz and Thouless [1972] J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys. 5, L124 (1972).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys. 6, 1181 (1973).
  * Kosterlitz [1974] J. M. Kosterlitz, J. Phys. C: Solid State Phys. 7, 1046 (1974).
  * Xu  _et al._ [2017] L. Xu, Z. Zangeneh, R. Yadav, S. Avdoshenko, J. van den Brink, A. Jesche, and L. Hozoi, Nanoscale 9, 10596 (2017).
  * Jesche  _et al._ [2015] A. Jesche, L. Ke, J. L. Jacobs, B. Harmon, R. S. Houk, and P. C. Canfield, Phys. Rev. B 91, 180403 (2015).
  * V. P. Antropov and V. N. Antonov [2014] V. P. Antropov and V. N. Antonov, Phys. Rev. B 90, 094406 (2014).
  * Lado and Pardo [2016] J. L. Lado and V. Pardo, Phys. Rev. B 94, 115141 (2016).
  * Kubo and Tomita [1954] R. Kubo and K. Tomita, J. Phys. Soc. Jpn. 9, 888 (1954).
  * Huber [1972] D. L. Huber, Phys. Rev. B 6, 3180 (1972).
  * Richards and Salamon [1974] P. M. Richards and M. B. Salamon, Phys. Rev. B 9, 32 (1974).
  * Benner and Boucher [1990] H. Benner and J. Boucher, in  _Magnetic Properties of Layered Transition Metal Compounds_, edited by L. de Jongh (Kluwer Academic Publishers, Dordrecht, 1990) Chap. Spin Dynamics in the Paramagnetic Regime: NMR and EPR in Two-Dimensional Magnets, pp. 323–378.
  * Zeisner  _et al._ [2019] J. Zeisner, A. Alfonsov, S. Selter, S. Aswartham, M. P. Ghimire, M. Richter, J. van den Brink, B. Büchner, and V. Kataev, Phys. Rev. B 99, 165109 (2019).
  * Koepernik and Eschrig [1999] K. Koepernik and H. Eschrig, Phys. Rev. B 59, 1743 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Lichtenstein and Katsnelson [1998] A. I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B 57, 6884 (1998).
  * Ahmed  _et al._ [2015] N. Ahmed, A. A. Tsirlin, and R. Nath, Phys. Rev. B 91, 214413 (2015).
  * Nawa  _et al._ [2017] K. Nawa, O. Janson, and Z. Hiroi, Phys. Rev. B 96, 104429 (2017).
  * Helgaker  _et al._ [2000] T. Helgaker, P. Jørgensen, and J. Olsen,  _Molecular Electronic-Structure Theory_ (Wiley, Chichester, 2000).
  * Werner  _et al._ [2012] H. J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, Wiley Rev: Comp. Mol. Sci. 2, 242 (2012).
  * Berning  _et al._ [2000] A. Berning, M. Schweizer, H.-J. Werner, P. J. Knowles, and P. Palmieri, Mol. Phys. 98, 1823 (2000).
  * Note [1] However, for the 3$d^{9}$ manifold and sufficiently high symmetry, this essentially corresponds to single-configuration restricted open-shell Hartree-Fock (ROHF) [32].
  * Bogdanov  _et al._ [2015] N. A. Bogdanov, V. M. Katukuri, J. Romhányi, V. Yushankhai, V. Kataev, B. Büchner, J. van den Brink, and L. Hozoi, Nat. Commun. 6, 7306 (2015).
  * Knowles and Werner [1992] P. J. Knowles and H.-J. Werner, Theor. Chim. Acta 84, 95 (1992).
  * Hozoi  _et al._ [2011] L. Hozoi, L. Siurakshina, P. Fulde, and J. van den Brink, Sci. Rep. 1, 65 (2011).
  * Huang  _et al._ [2011] H.-Y. Huang, N. A. Bogdanov, L. Siurakshina, P. Fulde, J. van den Brink, and L. Hozoi, Phys. Rev. B 84, 235125 (2011).
  * Zangeneh  _et al._ [2019] Z. Zangeneh, S. Avdoshenko, J. van den Brink, and L. Hozoi, Phys. Rev. B 100, 174436 (2019).
  * Note [2] Note that the fully relativistic calculations are nonmagnetic and all eigenvalues are doubly degenerate.
  * Klintenberg  _et al._ [2000] M. Klintenberg, S. E. Derenzo, and M. J. Weber, Comput. Phys. Commun. 131, 120 (2000).
  * Bogdanov  _et al._ [2017] N. A. Bogdanov, V. Bisogni, R. Kraus, C. Monney, K. Zhou, T. Schmitt, J. Geck, A. O. Mitrushchenkov, H. Stoll, J. van den Brink, and L. Hozoi, J. Phys. Condens. Matter 29, 035502 (2017).


# Semi-Supervised Multi-Task Learning Based Framework for Power System
Security Assessment

M. E. Za’ter1, A. Sajadi12, B. M. Hodge12

1University of Colorado Boulder, Boulder, CO 80303  2National Renewable Energy
Laboratory (NREL), Golden, CO 80401

###### Abstract

This paper develops a novel machine learning-based framework using Semi-
Supervised Multi-Task Learning (SS-MTL) for power system dynamic security
assessment that is accurate, reliable, and aware of topological changes. The
learning algorithm underlying the proposed framework integrates conditional
masked encoders and employs multi-task learning for classification-aware
feature representation, which improves the accuracy and scalability to larger
systems. Additionally, this framework incorporates a confidence measure for
its predictions, enhancing its reliability and interpretability. A topological
similarity index has also been incorporated to add topological awareness to
the framework. Various experiments on the IEEE 68-bus system were conducted to
validate the proposed method, employing two distinct database generation
techniques to generate the required data to train the machine learning
algorithm. The results demonstrate that our algorithm outperforms existing
state-of-the-art machine learning based techniques for security assessment in
terms of accuracy and robustness. Finally, our work underscores the value of
employing auto-encoders for security assessment, highlighting improvements in
accuracy, reliability, and robustness. All datasets and codes used have been
made publicly available to ensure reproducibility and transparency.

###### Index Terms:

machine learning, multi-tasking learning, power system stability, power system
security assessment

##  I Introduction

Power system security assessment requires high accuracy, robustness,
scalability, and interpretability [1], and the lack thereof could result in
interruptions of power delivery and potentially lead to disruption of service,
or blackouts. Recently, with the consistent increase in the periods of
operation dominated by variable renewable energy, it has become essential to
enhance the security assessment of the power system to accommodate higher
degrees of uncertainty, variability, and complexity, which are mainly
introduced by the nature of solar and wind generation [2, 3]. Such security
assessment system improvements are pivotal to assist system operators in
determining near real-time corrective actions to maintain system security
throughout the year [4].

Solutions for dynamic security assessment have evolved in two main directions;
model-based and model-free. Model-based techniques rely on mathematical models
of the power system to assess its security. They use models of power system
components to derive explicit nonlinear differential equations that describe a
power system, including its different components, and solve these equations in
closed form. Although they provide relatively accurate solutions depending on
the granularity of the model used, they have huge computational burdens, thus
offering limited applicability for large power systems [5]. Similarly,
electromagnetic transient (EMT) models are considered an accurate option for
dynamic security assessment, but they are also only practical for small-scale
system, due to their high computational complexity[6]. As an alternative to
model-based techniques model-free techniques have been developed which do not
rely on explicit mathematical models of the power system [7]. These algorithms
learn directly from historical and/or simulated data to assess the security of
the power system. In the model-free approach, the system is analyzed by
observing the patterns and relationships in the data, without requiring an
explicit understanding of the underlying system. These techniques are
typically faster and computationally cheaper, and hence scalable to larger
systems, but are often not as accurate as model-based techniques as they are
approximations of the actual mathematical representations [8]. Nonetheless,
they are practical for real-time operational purposes due to their
computational efficiency, and therefore their potential has fueled researchers
efforts into developing more accurate and reliable model-free methods [9, 10].

Over the past decade, model-free techniques for power system security
assessment have converged chiefly towards the use of machine learning (ML)
[11, 12, 13], mainly encouraged by the wider adoption of advanced measurement
instruments, e.g., especially Phasor Measurement Units (PMUs) [14]. This
sensing technology has enabled an unprecedented degree of visibility into the
power grid by generating an abundance of data, presenting immense
opportunities for data-driven methodologies within power system applications
[15]. A variety of ML-based techniques, such as decision trees (DT) and random
forests (RF), have been investigated for power system security assessment [16,
17, 18, 19, 20]. Although these methods offer high interpretability, which is
highly desirable for power system operators, they tend to fall short in terms
of achieving the desired accuracy, robustness, and scalability [21, 22].
Therefore, alternative solutions such as support vector machines (SVM)[23],
logistic regression [14], long short-term memory [24], and ensemble methods
[25, 23] have been explored to attain improved accuracy. Nevertheless, when
these algorithms are applied to larger power systems, such as those typically
encountered in real-world scenarios, their accuracy fails to meet the expected
performance standards [14, 16, 26, 24].

For applications such as power system security assessment that are safety-
critical and continuously changing (i.e different contingencies and
topologies), it is crucial to have high accuracy, stability against minor
disturbances and adversarial threats, awareness of system topology, and the
capacity to handle a broad range of possible scenarios from type of
contingency to different topologies. ML techniques for power system security
assessment strive to achieve a balance between interpretability and
scalability to effectively support the ever-evolving needs of the power
industry [27]. Interpretability is important to trace the decision-making
process [28] and recommend corrective actions,while scalability is important
because power systems are continuously subject to topological changes, arising
from various factors such as maintenance activities or environmental
conditions [29]. With all of these demanding requirements for power system
security assessment, conventional machine learning models have not been able
to satisfy the desired performance requirements [30, 31]. However, with the
huge surge of data and digitized measurements in power system, it became
feasible to apply more advanced, data-hungry algorithms of machine learning
and deep learning that have proven to be very effective in other fields. These
advanced methods, including the application of deep auto-encoders [14],
Bayesian methods [29], and tree-regularization techniques [32], have been
utilized to push the boundaries of what is achievable in power system security
assessment for larger systems. ML methods thus far have demonstrated promising
utility to tackle power system security assessment, exhibiting the capacity to
accommodate numerous contingencies and scenarios in offline computations,
surpassing the capabilities of deterministic and probabilistic approaches in
terms of the number of scenarios and contingencies considered [33]. However,
despite their high accuracy, these models have not yet fully satisfied the
stringent critical requirements inherent to the security assessment of power
systems [34].

In this paper we develop a new framework for dynamic security assessment in
power systems. At the core of this framework resides an ML-based technique
that leverages power-system-specific features for enhanced awareness and uses
a deep autoencoder for feature extraction. The proposed algorithm’s training
technique employs multi-task learning and semi-supervised learning approaches
to develop a classification-aware conditional deep autoencoder, customizing
the feature representation to enhance the classifier. The results demonstrate
the proposed framework outperforming the existing state-of-the-art
methodologies, achieving a 5% improvement in F1-score over the method in [29],
and a 7.2% increase in F1-score compared to [11] when evaluated on the IEEE
68-bus system. This performance was also maintained when tests included bad
data. All generated datasets and the code for the proposed algorithms have
been made available to the public at no cost for the sake of reproducibility
[35].

The contributions of this paper are three-fold. First, we present an enhanced
classifier that outperforms traditional methods like the Fisher score [36],
Principal Component Analysis (PCA) [25], and deep learning-based classifiers
[29, 37], in addition to integrating a conditional deep autoencoder. This
integration enhances the model’s contingency awareness, improving power system
security assessment and leading to better accuracy and sensitivity [14].
Second, we shift from the Bayesian inference proposed in [29, 38] to a feed-
forward neural network (FFNN) classifier that embeds a novel distance measure
based on the Mahalanobis distance [7, 39] specifically tailored for power
systems. This twofold enhancement not only sidesteps the noise robustness and
scalability challenges inherent in Bayesian inference [29], but also ensures a
more reliable and robust classifier that deeply considers power system
characteristics. Finally, we address a substantial gap in the current
literature by introducing a method that amplifies topology awareness grounded
in power system attributes, offering a more discerning measure of topological
differences compared to prevalent empirical methods [29] using Singular Value
Sequence (SVS), a method known for effectively capturing the characteristics
of power systems [7].

##  II Proposed Framework and Associated Processes and Algorithms

In this section, we discuss the structure and processes of our novel ML-based
framework for dynamic security assessment. Our framework approaches this
complex nonlinear problem in a multi-faceted manner, addressing the challenges
in database generation, limitations of machine learning algorithms, in
addition to the complex and nonlinear nature of power systems, as illustrated
in Fig. 1.

Figure 1: The proposed Machine Learning framework for Dynamic Security
Assessment (DSA)

As illustrated in 1, we divide the proposed framework into two main stages;
offline implementation, and online security assessment. The offline stage is
conducted before using the model for classifying the state of the system and
consists of preparing and generating the training dataset, in addition to the
model training and evaluation. The online stage includes classifying the
current and new operating conditions into secure and insecure conditions.

The offline stage starts with data collection on diverse system
characteristics; a complex database is created using computer simulation to
manage the data collected from power system dynamic simulations, which are
crucial for training the machine learning algorithms (II-B). After the dataset
is generated, the machine learning models are trained on these datasets
(II-C), which are then evaluated across different operating conditions on a
testset to ensure that the trained algorithm performs as desired II-F).
Whenever a new power system topology is introduced the similarity model
component assesses whether the topological changes necessitate database
updates and re-training the machine learning model (II-D). This step helps in
improving the topological awareness of the model, especially improving its
ability to minimize missed alarms, as described in Section II-D.

In the online stage, the trained model is used to classify new operating
conditions into secure or insecure states. With every state prediction
generated from the model a distance measure component provides confidence
levels for the predictions to provide a level of interpretability for
operators (II-E). The remainder of this section explore each of the components
mentioned above in detail.

###  II-A Data Collection and Assumptions

Data collection is a critical phase in any machine learning-based algorithm
process, as it is essential for training and testing. The more comprehensive,
well-distributed, and representative of the various scenarios that the model
might encounter the training data is, the more reliable, robust, and accurate
the performance of the ML algorithms will be.

Here, we rely on commercially available power system simulation tools (ANDES
[40]) to generate our data. The simulations are tailored to cover a vast range
of system operating conditions and scenarios, ensuring that the dataset is
diverse and comprehensive as further discussed in II-B. We consider this data
in lieu of high-precision high-frequency measurement data that are used in
real-time operations, such as Phasor Measurement Units (PMUs). Therefore,
parameters such as voltage levels, phase shifts, and real and reactive power
are all computed and integrated into the dataset.

###  II-B Database Generation

Ideally, machine learning models are trained on historical data of the power
systems, as historical data offers exact information on previous operating
conditions. However, historical data does not always cover a wide range of the
potential operating points, and often contains few insecure points, which are
insufficient to train an accurate data-driven technique. Therefore,
constructing a database beyond the available historical data is crucial to
build a machine learning based system that is accurate and reliable for system
security assessment. In this work, we employ two database construction
methodologies as detailed below.

####  II-B1 Time-Domain-Simulation

Time-domain-Simulation (TDS) is the most common method for database generation
in the literature [24] as it accurately determines both secure and insecure
operating points, and can take into consideration many aspects of security
such as transient, voltage, static, and small-signal stability. The process of
generating the database using TDS largely depends on the nature and
availability of historical data and observations, as well as user choices.
However, the majority of the literature [13, 1, 41] adheres to the flow shown
in Fig 2 which what is followed in this work.

Figure 2: Process flow diagram of time domain simulations to create a database

Initially, power system configuration data and available historical
observations are collected. Subsequently, a statistical model is applied to
these historical observations to generate a larger number of samples to cover
more (and more extreme) potential operating points [42]. These samples are
then considered as pre-fault operating conditions used to perform time-domain
simulations after setting up specific contingencies of interest [43]. In this
work, we attempt to consider more than one aspect of power system and not to
be limited to one, such as transient or voltage security [44, 12, 45]. The
labeling of secure and insecure for each operating state is based on the
following four criteria:

  * •

Transient Stability: For each specific contingency, the system is considered
transiently insecure if the transient stability index (TSI) defined by 1 is
less than 10%. In equation 1, $\triangledown\delta_{max}$ is the maximum
angular separation between any two rotor angles in degrees. The TSI in
Equation 1 is based on the TSAT power swing-based algorithm [25].

$TSI=\frac{360-\triangledown\delta_{max}}{360+\triangledown\delta_{max}}*100$
(1)

  * •

Small-signal stability: Entergy requires a 3% damping ratio [16]. This
criterion is applied to the inter-area oscillation modes of generator rotor
angles, with a frequency range of 0.25-1.0 Hz and varying amplitudes.

  * •

Voltage Stability: A system is considered insecure if any bus voltage deviates
from the range of 0.8 pu to 1.1 pu for more than 0.5 seconds [20].

  * •

Static Security: The overload index, as calculated in equation 2[13], was
taken into account.

$f_{x}=\sum_{i=1}^{N_{l}}wf_{i}(\frac{S_{mean,i}}{S_{max,i}})^{p}$ (2)

Where $f_{x}$ represents the overload performance index for the operating
point x, $N_{l}$ denotes the total number of transmission lines, and
$S_{mean,i}$ and $S_{max,i}$ indicate the average and maximum apparent power
flows of the $i$-th line.

####  II-B2 Efficient Database Generation

While TDS is an accurate and effective method for database generation, it has
its own limitations. The computational complexity limits the number of
operating conditions it can consider, which hinders its scalability, and thus
limits the size of system to which it can be applied for security assessment.
As a result, there has been a need to develop an efficient, modular and
parallelizable approach for generating databases, that are more easily
scalable to larger power systems [46, 28].

In light of this limitation researchers have started employing infeasibility
certificates to eliminate large portions of the search space of power system
operating conditions, if they are practically infeasible. After significantly
reducing the search scope, the security criteria is determined (i.e transient
stability, voltage stability, etc) where a detailed security boundary is built
upon this criteria. Consequently, more samples will be drawn close to a
security boundary using historical observations. This technique allows for
high parallelization as multiple processes can run concurrently to sample at
different security boundaries. This approach was first proposed by [46] and
then extended to include small-signal stability in [28] by assigning the step
direction according to the sensitivity of the damping ratio.

###  II-C Machine Learning Algorithm Training

Once the necessary database for training has been constructed, several steps
are needed to ensure efficient and accurate training of a machine learning
algorithm. These steps encompass data cleaning and pre-processing, feature
engineering and extraction, after which the data is ready for training the
model, as depicted in Fig. 3:

Figure 3: Machine Learning model training stage

In Fig. 3, each database represents the generated samples for a single
topology $T$, along with their corresponding binary labels (secure/insecure)
based on the previously mentioned criteria. Each training $X$ consists of pre-
fault operating conditions, including active and reactive power, power flows,
phase angles and voltages for each bus. The total number of features is $X\in
R^{(C\times n)\times m}$, where C represents the number of contingencies, $n$
denotes the number of operating conditions from one topology, and $m$ is the
sum of the number of generators, loads, power flows, voltages, and phase
angles in the system. The output label $Y\in R^{n}$ where $Y$ is a vector in
the size of the data representing the labels for each operating condition.

####  II-C1 Proposed Machine Learning Algorithm

In this work, we propose a novel technique that integrates the capabilities of
conditional masked deep auto-encoders that is trained jointly with the neural
network based classifier, named Semi-supervised Multi-Task Learning (SS-MTL),
as shown in Fig. 4.

Figure 4: Multi-task Semi-supervised Conditional masked Deep Auto-encoder

This proposed method utilizes a multi-task learning approach, which involves
training two tasks concurrently and minimizing a joint loss function across
these tasks. The two tasks include feature representation through the
conditional deep auto-encoder and the neural network classifier. By utilizing
this technique, we create a classification-aware feature representation
specifically tailored to our classification objective. In brief, deep auto-
encoders convert input features into a dense vector that represents the entire
set of input features (power system operating condition). The conditional
masking of these auto-encoders improves performance across multiple
contingencies. The resulting dense vector representations from the conditional
auto-encoders are then fed into a classifier, which is trained to distinguish
between secure and insecure power system states. Algorithm 1 demonstrates how
the proposed framework works which will be explained in later sections in
detail.

Algorithm 1 SS-MTL Framework

1: Define Learning rate, dropout, data size, batch size, Adam optimizer, etc.

1: Each batch $B$ contains samples only from one topology.

2: Initialize all parameters Function(SS-
MTL)($X_{train}^{T_{0}}$,$Y_{train}^{T_{0}}$)

3: repeat

4:

5: Assign 0 to $\beta$ and 1 to $\alpha$

6: repeat

7: Mini-batch $M$ optimization

8: until 2 Epochs

9: until Convergence

10: repeat

11: Mini-batch $M$ optimization

12: Calculate derivation w.r.t. $\theta$ in (6)

13: Update $\theta$: $\theta=\theta+\lambda\delta\hat{\theta}$

14: until $\theta$ has been optimized

15:

16: return $f_{\hat{W}}$

17: Feed-forward inference through $f_{\hat{W}}^{SSMTL}$ using
$X_{test}^{T_{k}}$, $Y_{test}^{T_{k}}$ Function(Distance measure)($D$, $th$)

18: Calculate $D$

19:

20: return Confidence

21: if New topology then

22: Define Threshold $h$ Function(Similarity Index)(power system data)

23: Calculate SVS

24: Calculate RMSE

25: if RMSE $>$ $h$ then

26: Train model again

27: else

28: Keep Model

29: end if

30: end if

##### Conditional masking auto-encoder

Auto-encoders are trained with the objective of reconstructing the input data,
and are commonly used as a method of feature representation in a lower
dimension than the input dimension. It consists of encoder and a decoder
networks, where the encoder is responsible for converting the higher
dimensional input into a dense lower dimensional vector (called latent
representation), whereas the decoder generates an output from the latent
representation that attempts to match the input vector.

Encoder network The encoder network consists of multiple feed forward layers,
each applying a nonlinear activation function to a linear transformation of
its input:

$h_{i}=f(W_{i}\cdot h_{i-1}+b_{i})$ (3)

where $h_{i}$ is the output of layer $i$, $W_{i}$ and $b_{i}$ are the weight
matrix and bias vector for layer $i$, and $f$ is a nonlinear activation
function, such as the ReLU or sigmoid function [21]. The input to the first
layer is $h_{0}=x$ the input data. The latent representation, $z$, is the
output of the final encoder layer, $z=h_{L}$, where $L$ is the number of
layers in the encoder network.

Decoder network: The decoder network is the inverse of the encoder network,
where it also consists of multiple layers, each applying a nonlinear
activation function to a linear transformation of its input:

$g_{i}=f(V_{i}g_{i-1}+c_{i})$ (4)

where $g_{i}$ is the output of layer $i$, $V_{i}$ and $c_{i}$ are the weight
matrix and bias vector for layer $i$, and $f$ is a nonlinear activation
function. The input to the first layer is $g_{0}=z$, the latent
representation. The reconstruction, $\hat{x}$, is the output of the final
decoder layer, $\hat{x}=g_{M}$, where $M$ is the number of layers in the
decoder network.

The conditional part of the auto-encoder is added in order for the model to be
more aware of the different contingencies. Given input data $x$, we first
condition the auto-encoder on some auxiliary information $c$, where this
information represents the label of the contingency in the training, in the
form of a one-hot encoder. The encoder maps the input $x$ and the conditioning
information $c$ to a latent representation $z$ as
$z=f_{\text{enc}}(x,c;\theta_{\text{enc}})$ where $f_{\text{enc}}$ is the
encoder function and $\theta_{\text{enc}}$ represents its parameters. The
decoder then maps the latent representation $z$ and the conditioning
information $c$ back to the reconstructed input $\hat{x}$ as
$\hat{x}=f_{\text{dec}}(z,c;\theta_{\text{dec}})$ where $f_{\text{dec}}$ is
the decoder function and $\theta_{\text{dec}}$ represents its parameters.

$L_{recon}=||x-D(E(x,c),c)||^{2}$ (5)

where $x$ is the input data and $C$ represents the one-hot encoded
representation of the contingencies, whereas the encoder $E$ maps the input to
a latent representation $z$, and the decoder $D$ reconstructs the input from
the latent representation.

##### Supervised classifier

The neural network classifier is trained using the latent representation $z$
obtained from the encoder. The classifier aims to minimize the classification
loss $L_{class}$:

$L_{class}=-\sum_{i=1}^{N}y_{i}\log C(E(x_{i},c_{i}))$ (6)

Here, $N$ represents the number of labeled samples, $y_{i}$ is the true label
for the $i$-th sample, and $c_{i}$ is the conditioning information for the
$i$-th sample.

To train these two tasks simultaneously on the combined loss function, which
is a combination of the reconstruction loss $L_{recon}$ and the classification
loss $L_{class}$, hyper-parameters $\alpha$ and $\beta$ are used to balance
the importance of the two losses:

$L=\alpha L_{recon}+\beta L_{class}$ (7)

The goal of the semi-supervised architecture is to minimize the combined loss
$L$ by optimizing the parameters of the encoder, decoder, and classifier. This
results in a model that leverages both the unsupervised reconstruction task
and the supervised classification task to improve the overall performance. The
training process for this method proceeds as follows: (i) Auto-encoder warmup,
where initially the auto-encoder is trained exclusively by setting the $\beta$
in 7 to 0 and $\alpha$ to 1. This step is introduced to prevent the classifier
network from training on representations when they are closer to random values
than representative weights. (ii) The values of $\alpha$ and $\beta$ are
adjusted to 0.5, which prompts the optimizer to update both the classifier and
auto-encoder weights. (iii) In the final stage, the values of $\alpha$ and
$\beta$ are set to 0.75 and 0.25, respectively, allowing for the fine-tuning
of the classifier for secure/insecure predictions.

###  II-D Topological Similarity Index Model

Topological changes in a power system can arise due to various reasons, such
as planned maintenance, equipment failures, or the integration of renewable
energy sources. These changes can significantly affect the performance of a
machine learning-based power security assessment too, as it might alter the
behaviour of such systems. Therefore, we propose adding an additional layer
that detects topological similarity and evaluates the need to re-train the
algorithm on the newly introduced topology.

Multiple indices of topological similarity have been proposed and evaluated in
the literature [47, 48] for power systems applications, with the Singular
Value Sequence (SVS) being recognized as the most effective in terms of
accuracy and efficiency [48]. SVS is calculated using singular value
decomposition (SVD) [48, 21]. The singular value sequence (SVS) is the ordered
set of singular values. In the context of power systems, SVS offer insights
about the power system characteristics that have been utilized in previous
applications such as fault detection and diagnosis, and stability analysis.

For power system topological similarity between two different topologies, we
employ the Root Mean Square Error (RMSE) between the SVS of two different
topologies to quantify their similarity. The lower the RMSE value, the more
similar the topologies. Let SVS
$S_{1}=\\{\sigma_{11},\sigma_{12},...,\sigma_{1n}\\}$ and
$S_{2}=\\{\sigma_{21},\sigma_{22},...,\sigma_{2n}\\}$ be the singular value
sequences of the two different power system topologies. The RMSE between
$S_{1}$ and $S_{2}$ can be calculated as follows:

$\text{RMSE}(S_{1},S_{2})=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(\sigma_{1i}-\sigma_{2i})^{2}}$
(8)

where $\sigma$ are the singular values obtained from the SVD.

###  II-E Distance Measure

Finally, after the predictions are generated through the classifier, we employ
the Mahalanobis distance [39] to estimate the uncertainty in the prediction,
in order to provide the operator with insights into the prediction, as well as
establish confidence in their reliability. The Mahalanobis distance is a
metric used to measure the distance between two data points in a multivariate
space while taking into account the correlations between variables. It is
particularly useful for measuring distances between data points when the
variables have different scales or units. The Mahalanobis distance has been
utilized to measure the distance between two operating conditions in a power
system by considering the correlations between power system variables and the
different scales at which they operate [49]. The Mahalanobis distance between
two data points $x$ and $y$ is defined as [50]:

$D_{M}(x,y)=\sqrt{(x-y)^{T}\cdot S^{-1}\cdot(x-y)}$ (9)

where $S$ is the covariance matrix of the data, and $S^{-1}$ is its inverse.
The covariance matrix captures the correlations between the different
variables in the data, and the inverse of this matrix is used to normalize the
distance calculation.

In the context of power systems, operating conditions can be represented by a
set of variables such as voltage magnitudes, angles, active and reactive power
injections, and power flows. When measuring the distance between two operating
conditions, the Mahalanobis distance takes into account the correlations
between these variables and the different scales at which they operate,
providing a more meaningful measure of distance than the Euclidean distance
[39].

To provide confidence in the predictions of the network state, the above-
mentioned distance is used to measure the distance between the current
operating point (where the model generated the label) and all of the operating
points in the dabase on which the model was trained. This assumes that it is
more likely that the trained model will generate more accurate predictions if
the current operating point is close to the trained operating conditions.
Nevertheless, comparing the current operating point with all the operating
points in the database may not result in a significant difference when all the
distances are averaged. For this reason, we calculate the average of the
Mahalanobis distances between the current operating condition and its 50
nearest points. This average is then used to determine the confidence level
associated with the current condition.

###  II-F Model Evaluation:

Power system security assessment is a safety critical application where the
cost of missed alarms out-weighs the cost of false ones. Therefore we need to
establish a metric which penalizes missing alarms. Resilience against
imperfect data is also an important aspect of a strong framework, hence we
will evaluate the algorithm’s resilience to better understand how it would
perform under realistic operating conditions with measurement uncertainty,
communication delays, and missing data.

####  II-F1 F-beta Score:

The F-beta score is a measure of a classification model’s accuracy, which
considers both precision ($\phi$) and recall ($\rho$). It is a weighted
harmonic mean of precision and recall, with the beta parameter determining the
importance of recall relative to precision. The F-beta score is given by [21]:

$F_{\beta}=(1+\beta^{2})\frac{\mathrm{\phi}\cdot\mathrm{\rho}}{(\beta^{2}\cdot\mathrm{\phi})+\mathrm{\rho}}$
(10)

The selection of $\beta$ was chosen to be 2 since, in the context of power
system security assessment it emphasizes penalizing missed alarms greater than
false alarms. This is because missed alarms could result in significant
reliability repercussions, i.e. potential blackouts, while false alarms,
although not harmless, may merely cause operators to take unnecessary
preventive actions.

####  II-F2 Bad Data Resilience

In practical scenarios, it is frequent to encounter inaccurate or missing
data. Thus, it is crucial to assess the ability of algorithms to manage poor
data, whether it results from faulty measurements, or external malicious
attacks [51]. Bad data in general exploits the vulnerabilities of machine
learning models in terms of their resilience, particularly deep neural
networks [52], which might cause them to produce incorrect outputs.

In order to simulate such conditions, we will introduce perturbations to the
input data, which are designed to deceive the model while appearing virtually
identical to the original input. A common approach to generating adversarial
examples to simulate bad data is the Fast Gradient Sign Method (FGSM)[34, 53].
Given a model with loss function $J(\theta,x,y)$, where $\theta$ represents
the model’s parameters, $x$ is the input data, and $y$ is the true label, the
bad data example $x_{adv}$ is generated as follows [54]:

$x_{\text{adv}}=x+\epsilon\cdot\text{sign}(\nabla_{x}J(\theta,x,y))$ (11)

Here, $\epsilon$ is a small constant that controls the magnitude of the
perturbation, and $\nabla xJ(\theta,x,y)$ is the gradient of the loss function
with respect to the input data x. The sign function returns the element-wise
sign of the gradient, ensuring that the perturbation is small yet effective.
Opting for minor perturbations over significant ones is due to the fact that
large disturbances are easily identified by both the algorithm and the model,
making it less likely for them to mislead the model. In this approach, a
neural network was trained using a portion of the original training data to
produce undesirable data examples using equation 11. It is important to note
that this equation is frequently utilized to create adversarial examples [55],
which rank among the most prevalent instances of flawed data or efforts to
trick the model.

##  III Implementation of SS-MTL and Testing Methodology

###  III-A Implementation

####  III-A1 Database Construction

The initial step is constructing the database, encompassing (i) pre-fault
operating conditions (OCs) data and (ii) corresponding post-fault labels that
identify whether the system is safe or not. The pre-fault OCs encompass active
and reactive power (either generation $G^{or}_{ac}\in R^{n\times g}$,
$G^{or}_{re}\in R^{n\times g}$ or load $L^{or}_{ac}\in R^{n\times l}$,
$L^{or}_{re}\in R^{n\times l}$, power flows $F^{or}_{ac}\in R^{n\times f}$,
$F^{or}_{re}\in R^{n\times f}$, voltages $V_{\text{or}}\in R^{n\times v}$, and
phase angles $\Theta_{\text{or}}\in R^{n\times\theta}$ for each bus. Combined,
these simulations form the $m$ dimensional original training features
$X_{\text{or}}\in R^{n\times m}$, where $n$ represents the size of the entire
dataset for one topology and $m=2\times(g+l+f)+v+\theta$. The corresponding
post-fault labels are denoted as $Y_{\text{or}}\in R^{n}$, where each element
is denoted as $y_{i}$, where $y_{i}$ is assigned 1 when secure and 0
otherwise, and is computed based on the aforementioned indices provided in
II-B. We use $X^{T0}_{\text{train}}$, $Y^{T0}_{\text{train}}$,
$X^{T0}_{\text{test}}$, and $Y^{T0}_{\text{test}}$ to denote the data from the
first topology $T0$.

###  III-B Testing Methodology

####  III-B1 Case Study: 68-bus Power System

The IEEE 68-bus power system is a well established benchmark that contains
16-machines across 5 regions which is a reduced equivalent of the
interconnected New England test system (NETS) and New York power system (NYPS)
[56, 57]. For the sake of simplicity, we assume that PMU devices are installed
and provide real-time measurements. Therefore, at each bus voltage magnitudes,
phase angles, active and reactive power are available. On the other hand, the
power flow data is calculated by the solver while generating the data.
Following [29], for generating and sampling observations from a set of
predefined operating conditions, we drew the active load power from a
multivariate Gaussian distribution. The active load power is scaled within
+-50% if the nominal values, while the reactive power was scaled by the
assumption of having a constant impedance for the buses. The power factor was
restricted to $[0.95,1]$.

The contingencies considered are selected following the rules stated in [56],
whereas the location of the faults were selected to be close to the
generators, while the clearing of the fault is coupled with line tripping.
Examples of the considered contingencies can be shown in Table III in [57]. To
generate operating conditions from variety of topologies from IEEE 68-bus
system [56], we have followed the works in [29, 11, 14] that generated 44
different topologies. Examples of different topologies are the double line
connection between bus NO.27 and NO.53 which suggests that any disconnection
between these two buses would likely have a minor impact on the rest of the
network. Bus NO.17 carries the highest load in the system, indicating that a
disconnection at this bus could result in significant changes to the power
flow patterns [29].

####  III-B2 Testing Procedure

The time-domain simulation method for sample generation was implemented in
Matlab and executed on an AWS instance with 84 CPUs and 16 GB of RAM. Machine
learning algorithms were developed in Python [35]. Training was conducted on
an instance with 16 CPUs for decision trees and support vector machines, while
deep learning algorithms utilized an Nvidia T4 GPU.

To ensure a thorough and fair assessment of machine learning algorithms for
power system security, the training and testing processes must be
comprehensive to include all possible scenarios such as unseen topologies.
Some of the generated topologies were included in both the training and
testing databases, allowing for an evaluation of the algorithms’ performance
on familiar topologies. Additionally, to assess performance on new, unseen
topologies, certain topologies were exclusively reserved for testing.
Furthermore, during the training phase, each batch provided to the model
contained observations from only one topology.

##  IV Numerical Results and Discussion

The algorithms were assessed across 22 distinct topologies for the IEEE 68 bus
system on 4000 overall test samples. It is worth noting that two topologies
were intentionally selected for all the aforementioned systems to
significantly deviate from the default topology. One of these topologies was
included in both the training and testing operating conditions, while the
other was only incorporated in the testing phase to evaluate the effectiveness
of the similarity model. The error bars in the results figure illustrates the
variance of the scores across the 22 different topologies.

Fig. 5 demonstrates that the proposed algorithm outperforms all other
algorithms in terms of F2 score on both database generation techniques. The
average F2 score on all topologies for the proposed algorithm is around 0.95,
while the BDAC algorithm proposed in [29] has a 0.925 average F2 score.
Machine learning based classifiers such as DT and RF had average F2 scores of
0.74 and 0.72, respectively. SVM had the best performance of this category of
algorithms with a 0.85 F2 score. Deep learning based methods without auto-
encoders outperformed DT and RF and had similar performance to SVM methods
with 0.84 and 0.81 for FFNN and Recurrent Neural Networks (RNNs) respectively.

From Fig. 5 several observations can be made. First of all, the performance of
neural networks and their variants surpasses DT, RF, and SVM, presenting a
distinct trade-off between accuracy and interpretability. Techniques that draw
a decision boundary, i.e. DT, RF, and SVM, are significantly more affected by
the unseen topology (lowest score for each bar in the figure). Algorithms that
employ auto-encoders (deep , convolutional, variational auto-encoders, etc)
for feature extraction demonstrate significantly better performance compared
to classifiers implemented on raw data. On the database generation side, the
efficient database generation technique produced improved results across
almost all algorithms, which may be due to the larger database or the reduced
complexity in the data, as transient stability is not considered in the
criteria. However, the performance difference between the algorithms follows
the same pattern for both database generation techniques, which underscores
the importance of the database generation method in enriching the number and
quality of samples. Finally, the ensemble method (XGBoost) displays
exceptional performance on both dataset generation techniques, despite being
trained on raw data. This suggests that future work exploring this technique
combined with improved feature selection could potentially achieve a balance
between interpretability and accuracy.

Figure 5: Performance of the algorithms on the IEEE 68 bus system

###  IV-A Reliability and Vulnerability Assessment

To assess the vulnerability of the trained algorithms to bad data, we
evaluated the algorithms’ performance with and without the injection of bad
data in the test-set. Fig. 6 displays the performance of the various
algorithms on the IEEE 68 bus system for both database generation algorithms.
Small perturbations in the input significantly impacted the performance of
classifiers such as DT, RF, SVM, and XGBoost. This is expected since they
establish hard decision boundaries, thus making them susceptible to minor
perturbations in the input. Feed-forward Neural Networks and Recurrent Neural
Networks displayed better robustness compared to other classifiers. However,
autoencoder-based techniques demonstrated superior resilience in the presence
of bad data, aligning with the literature [34, 58]. Methods that draw a hard
boundary across the decision boundary (i.e. DT, RF, SVM, etc.) were more
affected by bad data as their F2 Score tended to drop by 0.2 points between
the cases, whereas neural network based methods such as RNN and FFNN showed
some drop in F2 score performance (0.14 and 0.17, respectively), but less
relative to the hard boundary methods. Finally, auto-encoder methods [29, 59],
which are known for their resiliency [60], showed a drop of only 0.03 in
F2-score between the cases.

Figure 6: Performance of the algorithms on the IEEE 68 bus system, with and
without bad data

The performance comparison of the implemented algorithms shows that DT and RF
provide moderate accuracy and varying levels of interpretability, robustness,
and scalability, with low training time and parameter tuning complexity. SVM
have moderate accuracy and scalability but low interpretability, with higher
training time demands. XGBoost delivers high accuracy and scalability, with
moderate robustness and training time, and requires moderate parameter tuning
complexity. FFNN and RNN also offer high accuracy and scalability, with
moderate robustness and similar training time, though they require high
parameter tuning complexity. DAC and CBDAC, along with the proposed SS-MTL
(which achieves the highest accuracy), stand out with very high accuracy,
robustness, and scalability but are complex to tune and demand more data and
time for training. These models are highly effective but computationally
intensive.

###  IV-B Similarity Index Integration

As previously mentinoed, we included 22 different topologies in the test, 21
included in the training and one that was not included in the training. In
both Figs. 5 and 6 it is evident that for most of the algorithms, there is one
outlier value with performance vastly worse than the others. In these cases it
is the out-of-training topology, indicating a need to update the model in case
of a significant topology change. In this study, we have added topological
awareness to machine learning frameworks by incorporating a similarity index
model (II-D. This model, which depends on calculating the SVS between the
newly introduced topology and those in the training set, helps determine
whether the machine learning algorithm’s weightings need to be updated to
accommodate the new topology. Therefore, it is necessary to establish a
similarity threshold, determined empirically by analyzing the relationship
between the Root Mean Square (RMS) difference (equation 8) and its effect on
the F-beta score. To compute the mean RMSE for any new topology, we average
the RMSE values between the new topology and all existing topologies.

Based on empirical results, the similarity threshold for the IEEE 68 bus
system was determined to be 16; when higher, the model weights need to be
updated, when lower it is safe to continue with the same model.

Figure 7: Performance of the algorithms with the similarity model

As seen in Fig. 7, incorporating the topological similarity model, improved
the model’s awareness of unseen topologies and removed the topology with the
lowest score from the training data for the updated model. It should be noted,
however, that the scores for other topologies changed only slightly, remaining
within a similar range. The consistent behavior demonstrated by all algorithms
attests to the efficacy of the similarity model. Finally, it is important to
note that these experiments were conducted exclusively with TDS.

###  IV-C Speed Comparison

Since the motivation behind using machine learning for dynamic security
assessment is driven by the reduced computation speed that would allow real-
time predictions, it is essential to compare the speed of prediction between
TDS (the most accurate method for dynamic security assessment) and machine
learning techniques. Table II shows the time required for TDS and SS-MTL for
different systems.

TABLE I: Speed Comparison between TDS and SS-MTL | System
---
| Time (ms) for TDS
---
| Time (ms) for SS-MTL
---
IEEE 14 Bus | 2450 | 45
IEEE 39 Bus | 8640 | 67
IEEE 68 Bus | 19560 | 87
Nesta 162 Bus | 25780 | 122

As can be seen, SS-MTL provides approximately 200 times speed up in security
assessment. It is important to note that these TDS results are per
contingency, whereas machine learning results are for 22 contingencies per
operating condition. It is worth noting that the time complexity of TDS with
respect to the number of buses is $O(N^{3})$ (could vary based on the solver
used) where N is the number of buses in the system, suggesting challenges in
scaling to interconnection-scale systems. The complexity of a trained neural
network is proportional to number of layers and neurons, not the size of the
power system, and therefore it will increase very little when scaled to larger
power systems.

##  V Conclusion

This paper explores dynamic security assessment in power systems using machine
learning techniques. It introduces a Semi-Supervised Multi-Task Learning (SS-
MTL) framework to improve performance, topological awareness, and resilience
to bad data. Results show the proposed algorithm outperforms existing methods,
particularly in reducing missed alarms. This development serves as an initial
screening tool for assessing power system security, providing operators
insights into grid state under credible contingencies. Future work will focus
on reducing missed alarms and testing the algorithm’s effectiveness on larger
systems to encourage industrial adoption.

## References

  * [1] O. A. Alimi, K. Ouahada, and A. M. Abu-Mahfouz, “A review of machine learning approaches to power system security and stability,” _IEEE Access_ , vol. 8, pp. 113 512–113 531, 2020.
  * [2] P. Panciatici, Y. Hassaine, S. Fliscounakis, L. Platbrood, M. Ortega-Vazquez, J. Martinez-Ramos, and L. Wehenkel, “Security management under uncertainty: From day-ahead planning to intraday operation,” in _2010 IREP Symposium Bulk Power System Dynamics and Control-VIII (IREP)_. IEEE, 2010, pp. 1–8.
  * [3] A. Sajadi, R. Preece, and J. Milanović, “Identification of transient stability boundaries for power systems with multidimensional uncertainties using index-specific parametric space,” _International Journal of Electrical Power & Energy Systems_, vol. 123, p. 106152, 2020.
  * [4] J. L. Cremer and G. Strbac, “A machine-learning based probabilistic perspective on dynamic security assessment,” _International Journal of Electrical Power & Energy Systems_, vol. 128, p. 106571, 2021.
  * [5] H.-D. Chiang, _Direct methods for stability analysis of electric power systems: theoretical foundation, BCU methodologies, and applications_. John Wiley & Sons, 2011.
  * [6] R. W. Kenyon, B. Wang, A. Hoke, J. Tan, and B.-M. Hodge, “Comparison of electromagnetic transient and phasor dynamic simulations: Implications for inverter dominated systems,” in _2023 IEEE Power & Energy Society Innovative Smart Grid Technologies Conference (ISGT)_. IEEE, 2023, pp. 1–5.
  * [7] H. Fan, Y. Chen, S. Huang, X. Zhang, H. Guan, and D. Shi, “Post-fault transient stability assessment based on k-nearest neighbor algorithm with mahalanobis distance,” in _2018 International Conference on Power System Technology (POWERCON)_. IEEE, 2018, pp. 4417–4423.
  * [8] L. Duchesne, E. Karangelos, and L. Wehenkel, “Using machine learning to enable probabilistic reliability assessment in operation planning,” in _2018 Power Systems Computation Conference (PSCC)_. IEEE, 2018, pp. 1–8.
  * [9] D. R. Gurusinghe and A. D. Rajapakse, “Post-disturbance transient stability status prediction using synchrophasor measurements,” _IEEE Trans. on Power Systems_ , vol. 31, no. 5, pp. 3656–3664, 2015.
  * [10] H. Deng, J. Zhao, Y. Liu, and X. Wu, “A real-time generator-angle prediction method based on the modified grey verhulst model,” in _2011 4th International Conference on Electric Utility Deregulation and Restructuring and Power Technologies (DRPT)_. IEEE, 2011, pp. 179–184.
  * [11] J.-M. H. Arteaga, F. Hancharou, F. Thams, and S. Chatzivasileiadis, “Deep learning for power system security assessment,” in _2019 IEEE Milan PowerTech_. IEEE, 2019, pp. 1–6.
  * [12] M. Gholami, M. J. Sanjari, M. Safari, M. Akbari, and M. R. Kamali, “Static security assessment of power systems: A review,” _International Trans. on Electrical Energy Systems_ , vol. 30, no. 9, p. e12432, 2020.
  * [13] F. R. S. Sevilla and L. Vanfretti, “Static stability indexes for classification of power system time-domain simulations,” in _2015 IEEE Power & Energy Society Innovative Smart Grid Technologies Conference (ISGT)_. IEEE, 2015, pp. 1–5.
  * [14] Z. Zhang, D. Zhang, and R. C. Qiu, “Deep reinforcement learning for power system applications: An overview,” _CSEE Journal of Power and Energy Systems_ , vol. 6, no. 1, pp. 213–225, 2019.
  * [15] X. Wang, “Machine learning applications in power systems,” 2020.
  * [16] I. Genc, R. Diao, V. Vittal, S. Kolluri, and S. Mandal, “Decision tree-based preventive and corrective control applications for dynamic security enhancement in power systems,” _IEEE Trans. on Power Systems_ , vol. 25, no. 3, pp. 1611–1619, 2010.
  * [17] L. Wehenkel and M. Pavella, “Decision tree approach to power systems security assessment,” _International Journal of Electrical Power & Energy Systems_, vol. 15, no. 1, pp. 13–36, 1993.
  * [18] R. Sunitha, R. S. Kumar, and A. T. Mathew, “Online static security assessment module using artificial neural networks,” _IEEE Trans. on Power Systems_ , vol. 28, no. 4, pp. 4328–4335, 2013.
  * [19] M. Vasconcelos, L. M. Carvalho, J. Meirinhos, N. Omont, P. Gambier-Morel, G. Jamgotchian, D. Cirio, E. Ciapessoni, A. Pitto, I. Konstantelos _et al._ , “Online security assessment with load and renewable generation uncertainty: The itesla project approach,” in _2016 International Conference on Probabilistic Methods Applied to Power Systems (PMAPS)_. IEEE, 2016, pp. 1–8.
  * [20] C. Liu, K. Sun, Z. H. Rather, Z. Chen, C. L. Bak, P. Thøgersen, and P. Lund, “A systematic approach for dynamic security assessment and the corresponding preventive control scheme based on decision trees,” _IEEE Trans. on Power Systems_ , vol. 29, no. 2, pp. 717–730, 2013.
  * [21] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _nature_ , vol. 521, no. 7553, pp. 436–444, 2015.
  * [22] J. R. Aguero, E. Takayesu, D. Novosel, and R. Masiello, “Modernizing the grid: Challenges and opportunities for a sustainable future,” _IEEE Power and Energy Magazine_ , vol. 15, no. 3, pp. 74–83, 2017.
  * [23] F. R. Gomez, A. D. Rajapakse, U. D. Annakkage, and I. T. Fernando, “Support vector machine-based algorithm for post-fault transient stability status prediction using synchronized measurements,” _IEEE Trans. on Power Systems_ , vol. 26, no. 3, pp. 1474–1483, 2010.
  * [24] D. A. Copp, F. Wilches-Bernal, I. Gravagne, and D. A. Schoenwald, “Time-domain analysis of power system stability with damping control and asymmetric feedback delays,” in _2017 North American Power Symposium (NAPS)_. IEEE, 2017, pp. 1–6.
  * [25] C. Liu, F. Tang, and C. Leth Bak, “An accurate online dynamic security assessment scheme based on random forest,” _Energies_ , vol. 11, no. 7, p. 1914, 2018.
  * [26] I. Konstantelos, M. Sun, S. H. Tindemans, S. Issad, P. Panciatici, and G. Strbac, “Using vine copulas to generate representative system states for machine learning,” _IEEE Trans. on Power Systems_ , vol. 34, no. 1, pp. 225–235, 2018.
  * [27] F. Xiao, J. D. McCalley, Y. Ou, J. Adams, and S. Myers, “Contingency probability estimation using weather and geographical data for on-line security assessment,” in _2006 international conference on probabilistic methods applied to power systems_. IEEE, 2006, pp. 1–7.
  * [28] A. Venzke, D. K. Molzahn, and S. Chatzivasileiadis, “Efficient creation of datasets for data-driven power system applications,” _Electric Power Systems Research_ , vol. 190, p. 106614, 2021.
  * [29] T. Zhang, M. Sun, J. L. Cremer, N. Zhang, G. Strbac, and C. Kang, “A confidence-aware machine learning framework for dynamic security assessment,” _IEEE Trans. on Power Systems_ , vol. 36, no. 5, pp. 3907–3920, 2021.
  * [30] H. Ahmadi and H. Ghasemi, “Maximum penetration level of wind generation considering power system security limits,” _IET Generation, Transmission & Distribution_, vol. 6, no. 11, pp. 1164–1170, 2012.
  * [31] G. Díaz and J. R. Bermejo, “Static analysis of source code security: Assessment of tools against samate tests,” _Information and software technology_ , vol. 55, no. 8, pp. 1462–1476, 2013.
  * [32] C. Ren, Y. Xu, and R. Zhang, “An interpretable deep learning method for power system transient stability assessment via tree regularization,” _IEEE Trans. on Power Systems_ , vol. 37, no. 5, pp. 3359–3369, 2021.
  * [33] V. J. Mishra and M. D. Khardenvis, “Contingency analysis of power system,” in _2012 IEEE Students’ Conference on Electrical, Electronics and Computer Science_. IEEE, 2012, pp. 1–4.
  * [34] C. Ren, X. Du, Y. Xu, Q. Song, Y. Liu, and R. Tan, “Vulnerability analysis, robustness verification, and mitigation strategy for machine learning-based power system stability assessment model under adversarial examples,” _IEEE Trans. on Smart Grid_ , vol. 13, no. 2, pp. 1622–1632, 2021.
  * [35] “GitHub - muhi-zatar/Machine-learning-for-power-system-security-assessment — github.com,” https://github.com/muhi-zatar/Machine-learning-for-power-system-security-assessment, [Accessed 17-06-2024].
  * [36] S. Fliscounakis, P. Panciatici, F. Capitanescu, and L. Wehenkel, “Contingency ranking with respect to overloads in very large power systems taking into account uncertainty, preventive, and corrective actions,” _IEEE Trans. on Power Systems_ , vol. 28, no. 4, pp. 4909–4917, 2013.
  * [37] L. Meegahapola, P. Mancarella, D. Flynn, and R. Moreno, “Power system stability in the transition to a low carbon grid: A techno-economic perspective on challenges and opportunities,” _Wiley Interdisciplinary Reviews: Energy and Environment_ , vol. 10, no. 5, p. e399, 2021.
  * [38] T. Liu, Y. Liu, J. Liu, L. Wang, L. Xu, G. Qiu, and H. Gao, “A bayesian learning based scheme for online dynamic security assessment and preventive control,” _IEEE Trans. on Power Systems_ , vol. 35, no. 5, pp. 4088–4099, 2020.
  * [39] L. Mili, N. Vichare, M. Cheniae, and P. Rousseeuw, “Robust mahalanobis distances in power system state estimation,” in _Proceedings of 36th Midwest Symposium on Circuits and Systems_. IEEE, 1993, pp. 1014–1018.
  * [40] H. Cui and F. Li, “Andes: A python-based cyber-physical power system simulation tool,” in _2018 North American Power Symposium (NAPS)_. IEEE, 2018, pp. 1–6.
  * [41] O. A. Alimi, K. Ouahada, and A. M. Abu-Mahfouz, “Real time security assessment of the power system using a hybrid support vector machine and multilayer perceptron neural network algorithms,” _Sustainability_ , vol. 11, no. 13, p. 3586, 2019.
  * [42] N. V. Tomin, V. G. Kurbatsky, D. N. Sidorov, and A. V. Zhukov, “Machine learning techniques for power system security assessment,” _IFAC-PapersOnLine_ , vol. 49, no. 27, pp. 445–450, 2016.
  * [43] J. D. Glover, M. S. Sarma, and T. Overbye, _Power system analysis & design, SI version_. Cengage Learning, 2012.
  * [44] X. Zhou, X. Guan, D. Sun, H. Jiang, J. Peng, Y. Jin, and Y. Zhang, “Transient stability assessment based on gated graph neural network with imbalanced data in internet of energy,” _IEEE Internet of Things Journal_ , vol. 9, no. 12, pp. 9320–9331, 2021.
  * [45] H. Hagmar, L. Tong, R. Eriksson _et al._ , “Voltage instability prediction using a deep recurrent neural network,” _IEEE Trans. on Power Systems_ , vol. 36, no. 1, pp. 17–27, 2020.
  * [46] F. Thams, A. Venzke, R. Eriksson, and S. Chatzivasileiadis, “Efficient database generation for data-driven security assessment of power systems,” _IEEE Trans. on Power Systems_ , vol. 35, no. 1, pp. 30–41, 2019.
  * [47] B. Bush, Y. Chen, D. Ofori-Boateng, and Y. R. Gel, “Topological machine learning methods for power system responses to contingencies,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 17, 2021, pp. 15 262–15 269.
  * [48] H. Cao, Y. Xu, J. Tang, T. Yu, and J. Zhang, “A study of reliability assessment and topology similarity: Methodology, relationship, and applications,” _International Trans. on Electrical Energy Systems_ , vol. 30, no. 1, p. e12168, 2020.
  * [49] L. Pagnier and M. Chertkov, “Physics-informed graphical neural network for parameter & state estimations in power systems,” _arXiv preprint arXiv:2102.06349_ , 2021.
  * [50] M. I. Jordan and T. M. Mitchell, “Machine learning: Trends, perspectives, and prospects,” _Science_ , vol. 349, no. 6245, pp. 255–260, 2015.
  * [51] A. Abedi, L. Gaudard, and F. Romerio, “Review of major approaches to analyze vulnerability in power system,” _Reliability Engineering & System Safety_, vol. 183, pp. 153–172, 2019.
  * [52] N. Narodytska and S. P. Kasiviswanathan, “Simple black-box adversarial attacks on deep neural networks.” in _CVPR Workshops_ , vol. 2, 2017, p. 2.
  * [53] Y. Zheng, Z. Yan, K. Chen, J. Sun, Y. Xu, and Y. Liu, “Vulnerability assessment of deep reinforcement learning models for power system topology optimization,” _IEEE Trans. on Smart Grid_ , vol. 12, no. 4, pp. 3613–3623, 2021.
  * [54] H. Abdi and L. J. Williams, “Principal component analysis,” _Wiley interdisciplinary reviews: computational statistics_ , vol. 2, no. 4, pp. 433–459, 2010.
  * [55] A. Pan, Y. Lee, H. Zhang, Y. Chen, and Y. Shi, “Improving robustness of reinforcement learning for power system control with adversarial training,” _arXiv preprint arXiv:2110.08956_ , 2021.
  * [56] B. Pal and B. Chaudhuri, _Robust control in power systems_. Springer Science & Business Media, 2006.
  * [57] G. Anagnostou and B. C. Pal, “Impact of overexcitation limiters on the power system stability margin under stressed conditions,” _IEEE Trans. on Power Systems_ , vol. 31, no. 3, pp. 2327–2337, 2015.
  * [58] Y. Chen, Y. Tan, and D. Deka, “Is machine learning in power systems vulnerable?” in _2018 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm)_. IEEE, 2018, pp. 1–6.
  * [59] M. Zhang, J. Li, Y. Li, and R. Xu, “Deep learning for short-term voltage stability assessment of power systems,” _IEEE Access_ , vol. 9, pp. 29 711–29 718, 2021.
  * [60] Y. Wang, H. Yao, and S. Zhao, “Auto-encoder based dimensionality reduction,” _Neurocomputing_ , vol. 184, pp. 232–242, 2016.