# Data-Driven Ground-Fault Location Method in Distribution Power System With
Distributed Generation

Mauro Caporuscio  Department of Computer Science and Media Technology,
Linnaeus University  Antoine Dupuis  Department of Electrical Engineering,
Uppsala University  Welf Löwe  Department of Computer Science and Media
Technology, Linnaeus University

###### Abstract

The recent increase in renewable energy penetration at the distribution level
introduces a multi-directional power flow that outdated traditional fault
location techniques. To this extent, the development of new methods is needed
to ensure fast and accurate fault localization and, hence, strengthen power
system reliability. This paper proposes a data-driven ground fault location
method for the power distribution system. An 11-bus 20 kV power system is
modeled in Matlab/Simulink to simulate ground faults. The faults are generated
at different locations and under various system operational states. Time-
domain faulted three-phase voltages at the system substation are then analyzed
with discrete wavelet transform. Statistical quantities of the processed data
are eventually used to train an Artificial Neural Network (ANN) to find a
mapping between computed voltage features and faults. Specifically, three ANNs
allow the prediction of faulted phase, faulted branch, and fault distance from
the system substation separately. According to the results, the method shows
good potential, with a total relative error of 0,4% for fault distance
prediction. The method is applied to datasets with unknown system states to
test robustness.

##  1 Introduction

Power reliability is a very important property for any power grid. Reliability
can be achieved through preemptive or fault location techniques. Preemptive
techniques aim at avoiding outages by leveraging pathways and equipment
redundancy, which requires large investments. Fault location techniques aim at
decreasing the time for fault clearance. European networks provide a high
level of reliability with a low number of faults and a short fault-clearance
time (between 15 to 400 min per customer/year) [4]. However, since about 80%
of the faults observed in power networks occur at the distribution level,
fault location in distribution networks has become an area of high interest
for academic and industrial research, previously reserved for transmission
networks [4].

Ground faults are very challenging to locate. Traditionally, the fault
location process is initiated by customers, which notify an operator about an
outage [6]. The trouble call only indicates one specific location compromised
by the outage, and additional calls are needed to approximate the affected
area. This information is then combined with knowledge about the network and
positions of fault-clearing devices (e.g., Circuit Breaker) to identify
possible fault locations. Relying on trouble calls as a fault indicator has
several shortcomings: (1) customers tend to postpone the fault report, (2)
reports are usually incomplete, (3) customers make mistakes and report false
outages, and (4) night-time faults are less likely to be reported. Once a
fault is verified and confirmed, a technician is usually dispatched to check,
classify, and recover the fault [2]. This can be a tedious and resource-
consuming task as the area in which the fault occurred might be very large,
and the whole process usually depends on acquired knowledge about the area,
such as previous experience and historical data [2] [6]. To this extent, data
analysis can help to improve the accuracy of the fault location, reduce the
time needed to fix the problem and optimize the resources [6].

The main methods for fault location through data analysis can be categorized
as [15]: Impedance and Other Fundamental Frequency Component-Based Methods
(IBFL) and High-Frequency Components and Traveling Wave Based Methods (TWBFL).
IBFLs are widely used in distribution systems because of their cost-
effectiveness. The method computes multiple estimations by identifying a
number of possible fault locations. In fact, the impedance is calculated
starting from the measuring point and identifies all possible points in the
network with an impedance equal to the assumed fault impedance. TWBFLs are
based on the reflection and transmission of the fault-generated traveling
waves on the faulted power network. Although the fault can be located with
high accuracy, the implementation of such a technique is complex and more
expensive than the implementation of impedance-based techniques.

Motivated by low cost and high potential, data-driven methods and associated
machine learning techniques have been applied to many different fields and
hence enjoy extensive literature and research nowadays. The field of fault
location is not an exception to that. Recently, data-driven fault location
(DDFL)methods have been proposed to reduce the number of real-time
calculations and, therefore, the computational load. In contrast to model-
based techniques, which base fault localization on physical models, DDFL finds
complex mappings based on observations, so-called training data. DDFL methods
have been based on different machine learning (ML) algorithms, including ANN
and Support Vector Machines [16, 17]. In this context, the main contributions
of this paper are:

  1. (i)

Present a DDFL method adaptable to specific distribution networks.

  2. (ii)

Compare the proposed method with traditional fault localization methods and
other DDFL methods and highlight their difference and respective advantages.

Our DDFL method (i) suggests the following steps, each detailed in a separate
section. For a new distribution network, define a simulation model to generate
(raw) training data as detailed in Section 3. Instead of training the ML
models with raw data, process them to extract useful information, which lets
training converge faster. Section 4 defines this data processing using
discrete wavelet transform (DWT). The transformed data is used for training ML
models for the localization of faults anywhere in the specific distribution
network. Section 5 introduces our ANN models for DDFL. They are validated in
Section 6. Comparing our DDFL method with traditional and other DDFL methods
(ii) is done in the background and related work Section 2; points in favor and
against our method are further discussed in Section 7. Finally, Section 8
concludes the paper and shows directions of future research.

##  2 Background and Related Work

Fast fault recovery is provided by several types of fault location methods
that have respective advantages, drawbacks, or limitations induced by their
nature but also by their field of application. Concretely, a method suitable
for one distribution power system (DPS) might be impractical for another one.
This section reviews the methods suggested in the literature and spotlights
their respective characteristics for application in distribution systems.
Before, we briefly set the background of electrical power systems.

###  2.1 Power Systems Background

Power systems consist of three parts: generation, transmission, and
distribution, where electrical power is first _generated_ at different power
plants from which the bulk power is _transmitted_ across relatively long
distances after passing through a step-up transformer. The electrical is then
_distributed_ to commercial, industrial, and residential loads at different
voltage levels.

Carrying a tremendous quantity of electrical power over long distances,
transmission lines’ voltage levels must be high to avoid power loss. The power
is then carried at a lower voltage through step-down transformers to be
distributed to and consumed by customers. Distribution voltage levels can be
divided into two groups: $4-35kV$ for the primary distribution system and
$120-240V$ for secondary distribution systems at which the residential loads
are connected. This occurs at primary distribution levels that DG units can
directly be connected to.

Differences in voltage level don’t impact fault location methods. However, the
difference in topology between transmission and distribution systems does.
While transmission systems are mainly purely radial, DPS consists of branches
and tapped laterals dispersed over different types of areas. The electrical
power is carried by the transmission system to the primary distribution
substation through a step-down transformer. It is then carried to secondary
distribution transformers via feeder conductors. DPS are also subject to
changes in line composition and phase number, often leading to unbalanced
operation. These characteristics make DPS more exposed to different fault
sources such as wind, vegetation, lighting, storm, birds contact, and
equipment failure. Also, it limits the use of methods that would be perfectly
adapted to radial systems like transmission, as explained in section 2.2.

In addition to higher vulnerability and different topologies, modern DPS
include dynamic load and DG from renewable sources. Adding a current source at
this level of the overall system reshapes the power flow traditionally going
from the transmission to load through distribution systems. As a result, power
tends to flow dynamically within distribution systems, again limiting the use
of classical methods for fault location.

###  2.2 Ground Fault Location Methods

Traveling wave-based fault location (TWBFL) and impedance-based fault location
(IBFL) fault location are commonly applied methods and give satisfying
results. However, they both present non-negligible disadvantages, and higher
performance must be reached to improve network stability further. To this end,
emerging data-driven fault location (DDFL)methods using machine learning
(ML)algorithms are investigated and promising results.

Some fault location algorithms require the determination of the faulted area
first to be applied successfully. For instance, several methods give the
distance between the distribution system substation and the fault as a result
which raises problems of multiple possible locations when the network presents
several branches, also called the multiple estimations problem. In this case,
solutions for finding the correct faulted branch, i.e., the faulted area, must
be designed and included in the algorithm. However, several independent means
of fault area location exist, such as trouble call, fault indicator, or
advanced metering infrastructure methods [2].

####  2.2.1 Traveling wave-based fault location

When a fault occurs, high-frequency traveling waves of current and voltage are
generated and propagate away from the fault point in both directions along the
line toward its ends. The end of the line constitutes a circuit discontinuity
in terms of wave propagation. Whenever a traveling wave reaches a transition
point, it is divided into two secondary waves, one reflected and the other
refracted. Their amplitudes or respective energy are attenuated due to medium
property every time they meet a transition point and divide again until the
subsequent refracted and reflected waves vanish, and the post-fault steady is
reached [1]. The foundation of TWBFL methods is based on exploiting these
traveling waves: by analyzing their speed, the distance they travel, and the
time they reach the relays located at the line terminals, one can determine
the location where they originated.

As transmission systems are purely radial, their topology is particularly well
suited for TWBFL methods. In fact, these methods were initially designed for
transmission line fault location and mostly applied at this level, for which
they are considered very mature and efficient techniques. On the contrary,
DPSs have lateral and sub-branches. This specific topology involves shorter
lines and more discontinuity, increasing the number and sampling rate of
measuring devices to be able to distinguish superimposed traveling waves,
making the implementation cost of this method higher at distribution level
[2]. Nonetheless, TWBFL can still successfully be applied at distribution
levels. The TWBFL adoption proposed in [10] was tested with simulation-
generated traveling wave data. Two different tests gave errors relative to the
total power line length of 0,45% and 1,36%, resp. The two-staged approach
proposed in [9] determines the fault area first before accurate fault
localization. For different fault locations, fault resistances, inception
angles, and fault types, the smallest relative error is 0,34%, and the largest
is 3,31% with an average of 1,69%. Finally, an original and robust TWBFL
method for single line to ground fault localization was proposed in [5].
Simulations involving different fault locations, fault resistances, inception
angles, power line section area, and disturbance are run on a 34-bus system
with and without two DG units to observe their impact on the method. Tests
show that the method performs well even with DG units included. The mean
relative error is 1,74%.

TWBFL methods have both great advantages and drawbacks, which also depend on
the applied method. Generally, the implementation and high cost of high
sampling rate measurement devices, synchronization, and communication devices,
as well as required knowledge about line configuration and parameters to
compute propagation velocity, constitute non-negligible requirements and
drawbacks. The limitation regarding the network topology is also a very
important characteristic of this method. Systems that are not purely radial,
just like DPS, present obstacles to the implementation of TWBFL methods. Also,
discontinuity in power line material and size implies a change in wave
velocity that increases errors. Moreover, it remains quite insensible to fault
impedance. Finally, this method is mostly independent of network data and
consequently insensitive to modeling errors.

####  2.2.2 Impedance-based fault location

Instead of transient signals, fundamental frequency root mean square (RMS)
values of current and voltage during fault are exploited in IBFL methods. In
most cases, phase voltages and currents are only available at the substation
of the distribution system and constitute the required measurements to perform
the various algorithms for fault location. The IBFL methods are based on
circuit analysis theory and particularly the use of Kirchoff’s voltage and
current law to estimate the fault distance. Starting from the first line
section of the DPS, the faulted conditions of the next sections are found
iteratively until the faulted section is reached, where the fault location can
be estimated. This fault location highly depends on the impedance modeling of
the overall system, such as self and mutual line impedance, load impedance as
well as equivalent laterals. Hence, the accuracy of the method relies on how
accurately the different system components are modeled and, therefore, highly
relies on system knowledge. The main challenges associated with modeling are
dynamic load profiles, distribution line heterogeneity, and distributed
generation. Moreover, DG units have an intermittent and unpredictable
generation that transforms the direction of power flow from radial to multi-
directional and its magnitude from quasi-constant to variable within the
distribution networks. Consequently, methods that claim to be suitable for
high penetration of DG must consider these changes. Hence, one IBFL method
very optimized for one type of DPS might not be adapted for a different one
[12].

In [13], the authors present a method to calculate ground fault location
adapted to unbalanced systems with DG. The problem is solved using RMS values
of voltage and current during the fault available at the system substation.
The adaptation of the method to the presence of DG shows good results with an
average relative error of 0,2%. An adaption of the IBFL method in [13] for
unbalanced DPS considering shunt admittance and dynamic load is explained in
[7]. This method is tested on a 13-buses system. The results show the effect
of the load profile matching algorithm, which allows a relative error of only
0,92% for a load variation of 50% compared to the initial load profile. When
the load doesn’t vary, this method gives a relative error of 0,12%. Another
related method was proposed on a 6-buses DPS [11]. In contrast to the methods
discussed earlier, the authors use RMS voltage and current measurements at
every bus, which simplifies the method formulation. Similar to [13], the fault
current is approximated to be able to calculate the fault location. A DG unit
was also added to test the impact of DG on the method. Results are promising,
with a maximal relative error of 0,124% showing the robustness of the method.
The average relative error is as low as 0,024% due to the larger availability
of real measurements within the system.

Overall, IBFL methods are simple to implement and relatively low cost, notably
without the need for high-frequency measurement devices. They constitute the
reference for fault location at distribution levels. In that sense,
significant progress has been made to overcome the main challenges brought by
the required accurate modeling of DPSs. Indeed, solutions and algorithms to
model dynamic load, laterals as well as penetration of distributed generation
have been designed. As an important drawback, system knowledge is mandatory
for IBFL methods to be applied since they are based on the equivalent
impedance of the system, such as line and load, and modeling of other system
components. Another requirement is the measurement of faulted RMS voltages and
currents at the substation.

####  2.2.3 Data-driven methods

are based on the use of a large amount of data containing knowledge about the
system under fault occurrence. The field of ML underwent major developments
during the past recent years, being applied to many different areas, including
the field of fault location. Techniques attempt to learn or train general
fault location functions (ML models) from examples of fault locations and
corresponding observations. ML models can predict either continuous or
discrete values, corresponding to classification and regression problems,
resp. The arguments are called features, and their corresponding outputs are
either called targets for regression and categories for classification
problems. A part of the dataset is used for training, and the remaining part
to testing the trained model. Since the model is fit the training dataset, the
testing dataset enables an unbiased evaluation of the trained model. While ML
techniques can find highly complex nonlinear functions, their interpretation
can be challenging, if not impossible, for human logic. Consequently, finding
the best ML model is often an error and trial process.

The selection of the dataset plays a key role in guaranteeing that the trained
model generalizes well to a new, unknown input. That is, the chosen dataset
and its features should be similar to the potential future input vector that
this model will be used on. When a fault occurs, the power system can operate
under various conditions, and the fault has different parameters. Operating
conditions mainly include the rate of DG penetration and load values; fault
parameters are the fault impedance, the faulted phases, and the angle of the
voltage phase when the fault occurs, called the inception angle. The unique
combination of different system operating conditions and fault parameters for
a given fault occurrence is called a fault scenario. To be generalizable, the
dataset should be associated with realistic and various fault scenarios. In
the case of the power system fault location problem, transient or steady-state
three-phase voltage and/or current measured at one or more locations within
the system during the fault are most often taken as the features for training
an ML model. Naturally, the location of the fault is the target of the model.
Measurements of voltage and current should be done under various fault
scenarios and at different locations within the system, enabling the ML model
to learn the mapping of these signals to the location of the fault.

The accuracy of DDFL methods highly depends on the degree of detail and
exactitude of the power system modeling. To perform well, the system model
used to generate the dataset should be as close as possible to the real system
on which the method will be applied in real life. In practice, when a real
fault occurs, measurements are done and provide the input to the ML model
trained on the training data generated by the system model. As a result, the
ML model predicts the fault location.

In [17], a method to locate ground fault uses support vector regression to
compute the fault location. Traveling wave data recorded at the substation are
used as raw data. It is then processed with modal transformation on which
discrete wavelet transform of scale from 1 to 6 is applied. The arrival time
delay between modal components and peak amplitude ratio of modal components
wavelet coefficients are used as input features to train the ML model. This
data is generated for a balanced long distribution system that includes two
laterals. Faults are simulated every 0,5 km. The number of simulations
performed is as low as 40 since the correlation between chosen features and
fault location is easily established and almost proportional. In terms of
results, this method predicts fault location with a mean relative error of
6,2% after testing the method with 17 samples.

The method proposed in [8] uses different variants of ANN s to classify and
locate faults within a transmission system. The data is recorded at the
faulted transmission line terminal and contains measurements for the half-
cycle duration of three-phase current, voltage, and power phasors. Samples of
each signal are arranged in individual matrices, and their respective
determinant is computed, serving as input features for the ANN s. In total,
3’000 fault scenarios are simulated at 27 different fault locations. Fault
type classification ANN uses 32 input features, whereas only 6 are used for
fault location ANN. Networks are trained with a training dataset of 400
simulations, and the remaining 80’600 simulations are used for testing. Even
with this small training dataset, faults are classified with an accuracy of
99,8% and location is predicted with a mean relative error of 6,3%.

The method described [14] exploits the relation between the energy spectrum of
voltage transient traveling waves recorded at the substation and fault
location. The energy spectrum is found by using the discrete wavelet transform
at scale 8 to 4’096, and input features are the energy content of each scale,
corresponding to 10 input features. To predict the fault location, an ANN is
chosen. The training is done with 3’474 samples of different fault scenarios
in a distribution system that contains 5 laterals. Applied to the test data,
the method has a mean relative error of 0,5%.

DDFL methods are flexible since they can be applied to any power system type,
architecture, and topology as long as this system can be modeled accurately
enough using EMTP s. The accuracy of the built model directly impacts the
precision of the method: the whole system and its numerical model should
behave similarly in faulted simulation. Hence, limitations might be found when
building a very detailed model, leading to a computation burden when the
faults are simulated for data generation. Moreover, in real life, various
disturbances of different natures that cannot be anticipated or modeled in the
numerical version could potentially lead to more significant errors in fault
location prediction. Since DDFL methods rely on accurate power system
modeling, for any change in the network topology or else, the numerical model
should be updated, data generated once again, and the ML model retrained.

##  3 Dataset Generation

The major obstacle to employing ML in practice is the lack of existing
datasets. Indeed, available data on faulted power systems are neither large
enough nor pertinent enough. To solve this problem, electromagnetic transient
program (EMTP)systems allow the modeling of power systems, the simulation of
ground faults as well as the measurement of voltage or current signals
everywhere in the system. Consequently, a large and relevant dataset of
faulted voltage and current can be generated with the help of simulation,
enabling the development of DDFL methods.

When a fault occurs, high-frequency traveling waves of current and voltage are
generated and propagate away along the line toward its ends. Therefore,
voltage and current waveforms undergo a discontinuity where high-frequency
transients are added to steady-state sinusoids until a new post-fault steady-
state is reached. These transients, their frequency, and their magnitude are
dependent on the system states and the fault locations. To this extent, the
transient three-phase voltages at the system substation are chosen to
constitute the raw dataset of the data-driven method.

The modeling of distribution power system (DPS)enables the on-demand
generation of large datasets required for the development of DDFL methods. To
generate three-phase faulted voltages at different locations and under
different system states, a DPS model shown in Figure 1 is developed in
SimulinkTM. It is a balanced 11-bus 20 kV DPS to which a distributed
generation unit modeled by a synchronous generator is added. Its topology and
network data are inspired from [14] and given in appendix A.

Figure 1: Distribution power system test model structure.

To be relevant, the dataset should map well the possible system states, that
is to have fault occurrence under various system conditions. Hence, a total of
168 different fault scenarios are chosen for fault occurrence. Two percentage
of DG penetration in the system (DG%: 10% and 50 seven faulted phase(s)
($F_{ph}$: a, b, c, ab, ac, bc, abc), four fault impedance ($Z_{f}$: 0,01,
0,1, 1 and 10$\Omega$), and three inception angles ($\theta_{i}$: 45°, 90° and
135°) are chosen. To reduce computation time, loads remain unchanged in all
scenarios. In addition, faults are simulated every 500 meters, giving a total
of 38 fault locations, and bringing the total number of fault simulations to
6384. Fault occurrence time is taken at $0,025$ s and the sampling frequency
is $0,67$ MHz;

Figure 2 shows four simulations of faulted three-phase voltages under
different fault scenarios. In particular, two-phase a-to-ground faults
occurring under the same system states but located at (c) 2000 m and (d) 8000m
on path 1 put into light the impact of the fault location on the three-phase
faulted voltages.

(a) $DG=10\%$ ; $F_{ph}=a$ ; $Z_{f}=0,01\Omega$ ; $\theta_{i}=45$ ; $x=500$m
on path 1

(b) $DG=50\%$ ; $F_{ph}=ab$ ; $Z_{f}=0,1\Omega$ ; $\theta_{i}=90$ ; $x=3000$m
on path 1

(c) $x=2000$m on path 1

(d) $x=8000$m on path 1

Figure 2: Three-phase voltages recorded at the system substation under various
fault scenarios and different fault locations.

In addition, the three time-domain modal voltages $V_{0}$, $V_{1}$ and $V_{2}$
are computed using modal transform. They form, together with the three-phase
voltages, the raw dataset of the method.

##  4 Data Processing

After simulation, the generated raw dataset must be processed to extract
useful information. The use of DWT is described in this section.

###  4.1 Discrete Wavelet Transform

Transient signals are located in time and so cannot be captured by classical
signal processing methods, such as the Fourier transform, which is used for
stationary signals, or short-time Fourier transform, which has fixed time and
frequency resolutions. These shortcomings are solved by wavelet transform
(WT)which allows simultaneous time-frequency analysis of signals at different
time-frequency resolutions.

The WT is based on a mother-wavelet function $\psi(t)$ that can be stretched
or compressed in time with the use of scaling parameter $a$. Each unique
scaled wavelet $\psi_{a,b}(t)$, called daughter-wavelets, is shifted in time
along the signal $s(t)$ via the shifting parameters $b$ the signal. The
mother-wavelet must be finite and well localized in time and satisfy specific
conditions, described in [3]. Daubechies-4 (db4) is chosen for this study. The
continuous wavelet transform produces coefficients given as follows:

$\begin{split}&C_{a,b}=\int_{-\infty}^{+\infty}s(t)\psi^{*}_{a,b}(t)dt\\\
&a>0,-\infty<b<+\infty.\end{split}$ (1)

Where $\psi^{*}_{a,b}(t)$ is the conjugate of the daughter wavelet defined by
(2)

$\psi_{a,b}(t)=\frac{1}{\sqrt{a}}\psi\left(\frac{t-b}{a}\right).$ (2)

The wavelets coefficients are obtained for a center pseudo-frequency of the
wavelet decomposition that can be expressed as:

$F=\frac{f_{c}f_{s}}{s}$ (3)

with $f_{c}$ the center frequency of the daughter wavelet and $f_{s}$ the
sampling frequency of the signal. Because the scaling parameter is inversely
proportional to the wavelet frequency $f_{c}$, the use of different scales $a$
allows the extraction of different frequencies from the original signal. This
way, the signal can be analyzed in different frequency bands. The energy
content of each scale $a$ or frequency, $E_{wave}$ is defined by (4) below:

$E_{wave}(a)=\sum_{n=0}^{N-1}(C(a,nT_{s}))^{2}.$ (4)

The numerical implementation of the continuous WT on a sampling signal $s(t)$
is as follows (according to [3]):

$\displaystyle
C_{a,iT_{s}}=\frac{T_{s}}{\sqrt{a}}\sum_{n=0}^{N-1}\psi^{*}\left(\frac{(n-i)T_{s}}{a}\right)s(nT_{s})$
(5) $\displaystyle i=0,1,2,...,N$

where $N$ is the number of signal samples and $T_{s}$ is the sampling period.

From (5), one can deduce that the number of wavelet coefficients is equal to
the number of samples times the number of scales which can lead to
computational burdens. The DWT allows quicker signal processing by using
dyadic scaling and shifting. This implies discrete scaling and shifting
parameters shown in (6).

$\begin{split}&a=2^{j}\\\ &b=2^{j}m\end{split}$ (6)

where $j$ and $m$ are integer numbers. Hence, the discrete daughter wavelets
are as follows

$\psi_{j,m}(t)=\frac{1}{\sqrt{2^{j}}}\psi\left(\frac{t-m2^{j}}{2^{j}}\right).$
(7)

The implementation of DWT decomposes the signal using high-pass and low-pass
filters, giving high-pass and low-pass sub-banks. The high-pass and low-pass
sub-banks are populated by detail coefficients $CD$ and approximation
coefficients $CA$, respectively. If $j>1$, a multi-level decomposition occurs
as described in Figure 3. After each level of decomposition $j$ with
associated scale $2^{j}$, the details coefficients become inputs of the next
level filtering with scale $2^{j+1}$ giving high-pass and low-pass sub-banks
of level $j+1$.

Figure 3: Discrete wavelet transform with multi-level decomposition where
$j=1,2,3$.

The covered frequency bands for each level decomposition $j$ are expressed
below (according to [9]):

$\begin{split}&CD_{j}:[2^{-(j+1)}f_{s},2^{-j}f_{s}]\\\
&CA_{j}:[0,2^{-(j+1)}f_{s}]\end{split}$ (8)

Corresponding values for $f_{s}=0.67$ MHz are presented in Table 1 for
$j=[1;8]$.

j | CD | Frequency band | CA | Frequency band
---|---|---|---|---
$1$ | $CD_{1}$ |  $167.5-335$ kHz | $CA_{1}$ |  $0-167.5$ MHz
$2$ | $CD_{2}$ |  $83.7-167.5$ kHz | $CA_{2}$ |  $0-83.7$ MHz
$3$ | $CD_{3}$ |  $41.9-83.7$ kHz | $CA_{3}$ |  $0-41.9$ kHz
$4$ | $CD_{4}$ |  $20.9-41.9$ kHz | $CA_{4}$ |  $0-20.9$ kHz
$5$ | $CD_{5}$ |  $10.5-20.9$ kHz | $CA_{5}$ |  $0-10.5$ kHz
$6$ | $CD_{6}$ |  $5.2-10.5$ kHz | $CA_{6}$ |  $0-5.5$ kHz
$7$ | $CD_{7}$ |  $2.6-5.2$ kHz | $CA_{7}$ |  $0-2.6$ kHz
$8$ | $CD_{8}$ |  $1.3-2.6$ kHz | $CA_{8}$ |  $0-1.3$ kHz
Table 1: Frequency band values of multi-level DWT analysis for levels
$j=[1;8]$

The higher the decomposition level is, the fewer the samples and hence,
wavelet coefficients. Therefore, high-frequency resolution but low time
resolution is obtained, and the other way around for small-level
decomposition.

The dataset is processed with DWT until level 8. The computation is carried
out on MATLABTM with the help of the Wavelet toolbox. The decomposition of a
three-phase voltages sample is shown in Figure 4.

Figure 4: Detail wavelet coefficients of phase A voltage for a phase A-to-
ground fault at 4500 m along path 1.

###  4.2 Relevant Features Extraction

After having processed voltage phases and modes signals, relevant features
should be extracted to constitute the final dataset. In this paper,
statistical quantities of the processed dataset are used as input to the
neural networks. These are standard deviation $(std)$, variance $(var)$, third
central moment $(mom3)$, skewness $(skn)$, and mode $(mode)$ and are applied
to the following data for every simulation:

  * •

Three-phase time-domain voltages: $V_{a}$, $V_{b}$ and $V_{c}$

  * •

Time-domain modal voltages: $V_{0}$, $V_{1}$ and $V_{2}$

  * •

Level 8 detail and approximation wavelet coefficients of each phase voltages:
$V_{a}CD_{8}$, $V_{b}CD_{8}$, $V_{c}CD_{8}$, $V_{a}CA_{8}$, $V_{b}CA_{8}$ and
$V_{c}CA_{8}$

  * •

Level 8 detail and approximation wavelet coefficients of each modal voltages:
$V_{0}CD_{8}$, $V_{1}CD_{8}$, $V_{2}CD_{8}$, $V_{0}CA_{8}$, $V_{1}CA_{8}$ and
$V_{2}CA_{8}$

In addition, frequency energy contents, $E_{wave}(V_{a}CD_{8})$,$...$,
$E_{wave}(V_{2}CA_{8})$, are also taken as input features.

Relation between these inputs and the networks’ targets can be directly
observed. Indeed, correlation with faulted phase(s) classification is shown in
Figure 5 where the standard deviation of $CD_{8}$ of each phase is plotted for
20 different simulations having the same faulted phases. One can observe that
there is a discernible difference between healthy (c) and faulty phases (ab).

Figure 5: Correlation between faulted phase and standard deviation.

Similarly, Figure 6 shows the effect of fault distance by giving the variance
of $CD_{8}$ for ab-to-ground fault simulation that differs only by the fault
distance along path 1. It is visible that the variance of $CD_{8}$ of the
faulted phases decreases with distance.

Figure 6: Correlation between fault distance and variance.

The faulted path is, however, harder to identify from the available features.
This is attested in Figure 7. It indicates the skewness of faulted phase
voltages from simulations with the same fault scenario but different locations
in distance and path. One can notice that faults occurring at the same
distance but on different paths aren’t well dissociated for paths 2 and 3.
Nevertheless, paths 1, 2, and 4 can be distinguished quite well.

Figure 7: Correlation between faulted path and skewness.

Only the statistical quantities having a noticeable correlation with the
networks’ targets are kept as input features for the fault location method.
They are referred to later in the next section.

##  5 Fault Location Models

After the determination of the relevant features, they are fed as input to
three different sets of ANN s. Each set has three distinct prediction targets,
namely the faulted phase(s), the fault distance, and the faulted path. The
fault location methodology using these different sets is described in this
section.

###  5.1 Artificial Neural Networks

ANN s are used for their strong predictive power as well as their high
versatility. Given enough training data, an ANN can find _any_ mapping between
input and output where other conventional methods can’t. The structure of an
ANN is shown in Figure 8.

Figure 8: Artificial neural network architecture.

Feed-forward neural networks are chosen to predict the fault distance, whereas
for faulted phase and path classification, pattern recognition type networks
are utilized. The Levenberg-Marquardt (LM)and the scaled conjugate gradient
(SCG)optimization algorithms are chosen to train the feed-forward and pattern
recognition neural network, respectively. Every individual ANN s deployed in
the method is described in Table 2, where $P$ is the number of input features,
$Q$ is the number of hidden neurons, and $T$ is the number of outputs. $f^{h}$
and $f^{o}$ are the hidden and output layer activation function, respectively.
Each ANN is trained with the different input features determined in the
previous Section 4 according to their relevancy for the respective ANN’s
prediction. These inputs are referred to in Appendix A. Optimum weights and
biases values aren’t presented here for clarity purposes (but can be made
available on demand). Finally, the structure, that is, the number of hidden
neurons of each individual ANN s, was based on a trial-and-error approach for
which $10\leq Q\leq 50$, and hence differs between ANN s. The available
datasets are divided as follows: 70% is used for training, and the remaining
30% is later used for testing and validation.

ANN | Data size | P | Q | $f^{h}$ | T | $f^{o}$
---|---|---|---|---|---|---
Faulted Phase Classification
Ph | 6384 | 15 | 10 | tansig | 7 | softmax
Fault Distance Regression
D-a | 912 | 23 | 10 | tansig | 1 | $f(x)=x$
D-b | 912 | 23 | 20 | tansig | 1 | $f(x)=x$
D-c | 912 | 23 | 30 | tansig | 1 | $f(x)=x$
D-ab | 912 | 38 | 10 | tansig | 1 | $f(x)=x$
D-ac | 912 | 38 | 30 | tansig | 1 | $f(x)=x$
D-bc | 912 | 38 | 10 | tansig | 1 | $f(x)=x$
D-abc | 912 | 45 | 30 | tansig | 1 | $f(x)=x$
Faulted Path Classification
Pa-a$h_{1}$ | 408 | 30 | 50 | tansig | 4 | softmax
Pa-b$h_{1}$ | 408 | 30 | 50 | tansig | 4 | softmax
Pa-c$h_{1}$ | 408 | 30 | 30 | tansig | 4 | softmax
Pa-ab$h_{1}$ | 408 | 37 | 50 | tansig | 4 | softmax
Pa-ac$h_{1}$ | 408 | 37 | 50 | tansig | 4 | softmax
Pa-bc$h_{1}$ | 408 | 37 | 30 | tansig | 4 | softmax
Pa-abc$h_{1}$ | 408 | 26 | 30 | tansig | 4 | softmax
Pa-a$h_{2}$ | 504 | 30 | 50 | tansig | 3 | softmax
Pa-b$h_{2}$ | 504 | 30 | 50 | tansig | 3 | softmax
Pa-c$h_{2}$ | 504 | 30 | 30 | tansig | 3 | softmax
Pa-ab$h_{2}$ | 504 | 37 | 50 | tansig | 3 | softmax
Pa-ac$h_{2}$ | 504 | 37 | 50 | tansig | 3 | softmax
Pa-bc$h_{2}$ | 504 | 37 | 30 | tansig | 3 | softmax
Pa-abc$h_{2}$ | 504 | 26 | 30 | tansig | 3 | softmax
Table 2: Artificial neural networks structure

###  5.2 Faulted Phase Prediction

The faulted phase is computed using a single ANN, namely Ph, that hence
predicts the fault type: $ag$, $bg$, $cg$, $abg$, $acg$, $bcg$ or $abcg$. It
has 15 input features and uses the 6384 samples generated dataset for training
and testing.

###  5.3 Fault Distance Prediction

The fault distance prediction strategy takes advantage of the faulted phase
classification results: depending on the fault type, one of the seven fault
distance ANN s returns the fault distance $x$ from the substation as output.
If the faulted phase is found to be phase a, D-a is used, and similarly for
other faulted phases. Each of them is trained exclusively with input data
related to the faulted phase(s) they are specialized on. Consequently, only a
ratio of the generated dataset can be used for their training, that is, 912
samples out of 6384 for each of the seven ANN s.

###  5.4 Faulted Path Prediction

To avoid multiple location estimation problems, the faulted path must be
determined too. To this extent, the method uses 14 different trained ANN s: in
addition to having one ANN per fault type, they are further divided into two
groups $h_{1}$ and $h_{2}$. The groups $h_{1}$ for fault distance occurring at
a distance of $x\leq 4500$ meters and $h_{2}$ for $x>4500$ meters. This is
illustrated in Figure 9. This allows us to reduce the number of categories
that can be predicted. The $h_{1}$ labeled ANN s can predict four different
paths, from 1 to 4, and only three for the others, paths 1, 5, and 6. Thus,
the sample set size related to each ANN is even smaller and goes down to 408
and 504 samples for $h_{1}$ and $h_{2}$ groups, respectively.

Figure 9: Faulted path classification ANNs groups $h_{1}$ and $h_{2}$.

The algorithm of the method, from faulted phase to faulted path
classification, including data processing and feature extraction steps, is
summarized in Figure 10.

Figure 10: Fault location algorithm using optimized workflow.

##  6 Method Validation

After training the different networks involved in the method algorithm, a
validation test is carried out of the remaining 30% of the different samples
set. This section presents the theoretical results of the method, summarized
in Table 3.

ANN | Classification |  Mean $e_{rel}$ | Total
---|---|---|---
Faulted Phase Classification
Tfph | 100% | n/a | 100%
Fault Distance Regression
Ofd-A | n/a | 0,25% | 0,40%
Ofd-B | n/a | 0,58%
Ofd-C | n/a | 0,43%
Ofd-AB | n/a | 0,3%
Ofd-AC | n/a | 0,55%
Ofd-BC | n/a | 0,32%
Ofd-ABC | n/a | 0,40%
Faulted Path Classification
Ofp-A$H_{1}$ | 73,77% | n/a | 75,15%
Ofp-B$H_{1}$ | 68,85% | n/a
Ofp-C$H_{1}$ | 67,21% | n/a
Ofp-AB$H_{1}$ | 77,6% | n/a
Ofp-AC$H_{1}$ | 65,5% | n/a
Ofp-BC$H_{1}$ | 65,5% | n/a
Ofp-ABC$H_{1}$ | 62,29% | n/a
Ofp-A$H_{2}$ | 98,68% | n/a
Ofp-B$H_{2}$ | 85,52% | n/a
Ofp-C$H_{2}$ | 76,31% | n/a
Ofp-AB$H_{2}$ | 62,29% | n/a
Ofp-AC$H_{2}$ | 80,26% | n/a
Ofp-BC$H_{2}$ | 93,42% | n/a
Ofp-ABC$H_{2}$ | 75% | n/a
Table 3: ANNs prediction results

###  6.1 Fault Location Results

####  6.1.1 Faulted Phase Classification

After testing, the faulted phase can be predicted with an accuracy of 100%,
showing the great potential of the method for faulted phase classification.

####  6.1.2 Fault Distance Regression

The performance of the method for fault distance prediction is evaluated with
the mean relative error computed with equation (9).

$e_{rel}=\frac{|D_{cal}-D_{real}|}{L}$ (9)

where $D_{cal}$ and $D_{real}$ are the calculated and real fault distance from
the substation respectively. For the distribution system used, the maximum
distance from the system substation is $L=11$ km.

As a result, the combined mean relative errors of the seven fault distance
prediction ANN s is as low as 0,40%, showing here again considerable
potential.

####  6.1.3 Faulted Path Classification

In terms of faulted path classification, the method shows an accuracy of
75,15%. In contrast to the two previous location types, this result remains
quite unsatisfying for fault location purposes.

###  6.2 Relevancy of Multiple Artificial Neural Networks

To justify and show the interest of dividing the prediction problem into
smaller groups, and hence, using several ANN s, another method, that only
differs at fault distance prediction and faulted path prediction stage where
it uses a single ANN for each of the prediction, is tested. The results are
compared with the multiple ANN developed in this paper in Table 4.

Prediction Type | Multiple ANN | Single ANN
---|---|---
Faulted Phase | 100%
Total Classification |
Fault Distance | 0,40% | 2,49%
Total Mean $e_{rel}$ |  |
Faulted Path | 75,15% | 62,12%
Total Classification |  |
Table 4: ANNs prediction results

According to Table 4, one can easily conclude that the use of several ANN s is
very beneficial, the multiple ANN method outstanding by far the single ANN
method. Indeed, the fault distance prediction average relative error increases
consequently to 2,49% whereas the faulted path prediction accuracy is
decreased down to 62,12%. Concerning faulted phase classification, since both
methods use a single and hence same ANN, the 100% accuracy holds for both.

##  7 Discussion

As a result, one can conclude that the use of a set of ANN increases the
overall accuracy of the method and lead to a competitive fault distance
location. However, the major drawbacks of the method should be mentioned.

###  7.1 Robustness of the Method

Since faults can occur under a quasi-unlimited number of scenarios as well as
load values, the ability of the method to show good performance for unknown
datasets, that is, its robustness, is a major concern. To study the robustness
of the method, the method is tested on two new datasets:

  * •

Dataset 1: of 532 samples with new fault scenario (%DG: 30%; $Z_{f}$: 0,5 and
5$\Omega$; $\theta_{i}$: 70° and 110°).

  * •

Dataset 2: of 532 samples with similar fault scenario (%DG: 10%; $Z_{f}$: 0,1
and 1$\Omega$; $\theta_{i}$: 45° and 90°) but with load values increased of 30
%.

In this paper, the choice has been made not to consider load variation for the
method implementation, but testing the effect of load variation on the method
would give an interesting insight into load generalization capacity. As a
result, the method performance is significantly decreased in both cases as
shown in Table 5, especially in terms of fault distance prediction going down
to 3,95% and 12,24% for dataset 2 and 1 respectively. The method appears to be
less affected by the load variation than by new fault scenarios. Nevertheless,
path classification isn’t as much impacted with 57,89% and 63,53% accuracy for
new fault scenario and changed system load dataset respectively. Concerning
faulted phase determination, the method shows great robustness, with results
still very close to 100%.

Results | Original | Dataset | Dataset
---|---|---|---
Dataset | 1 | 2
Number of Faults | 6384 | 532 | 532
Faulted Phase | 100% | 98,87% | 100%
Fault Distance | 0,40% | 12,24% | 3,95%
Faulted Path | 75,15% | 57,89% | 63,53%
Table 5: Results of fault location method application on different datasets

###  7.2 Prediction Error Correlation

The multiple ANN s strategies leads to potential correlation errors between
the different prediction steps. A too-large error in fault distance prediction
would lead to the use of the wrong ANN s for faulted path classification.
Similarly, a wrong faulted phase prediction would lead to the use of the wrong
input features for the rest of the workflow. For the sake of clarity, these
two phenomena are called path correlation error (PCE)and phase correlation
error (PHCE), respectively. To observe the impact of these errors, the
correlation errors impact on datasets 1 and 2 are shown in Table 6.

Results | Dataset 1 | Dataset 2
---|---|---
Dataset samples | 532 | 532
% of PHCE | 1,13% | 0%
% of PCE | 12,60% | 8,46%
% of PHCE leading to PCE / Number of cases | 83,33% / 6 | n/a
Mean $e_{rel}$ in case of PHCE / Number of cases | 31,02% / 6 | n/a
Potential mean $e_{rel}$ without PHCE / Number of cases | 12,02% / 526 | n/a
Faulted path accuracy in case of PHCE/ Number of cases | 66,67% / 6 | n/a
Faulted path accuracy in case of PCE / Number of cases | 25% / 67 | 44,44% / 45
Faulted path accuracy in case of both PHCE and PCE/ Number of cases | 60% / 5 | n/a
Potential faulted path accuracy without correlation errors / Number of cases | 62,5% / 464 | 65,30% / 487
Table 6: Impact of correlation errors on fault location for dataset 1 and 2

One can first observe that the direct impact on the distance prediction
accuracy is quite significant, with a mean $e_{rel}$ up to 31,02% in the case
of dataset 1. In the case of faulted path classification, the accuracy goes
down to 25% and 44,44% in datasets 1 and 2, respectively. Indeed, the overall
impact of correlation errors is largely reduced by the small occurrence rate
of these errors even for datasets 1 and 2. PHCE occurs for only 1,13% of the
cases for dataset 1, whereas PCE occurrence is higher with 12,60% and 8,46%
for dataset 1 and 2 respectively and directly related to $e_{rel}$. As a
result, the faulted path classification is decreased from 62,5% down to 57,89%
and from 65,30% down to 63,53% for dataset 1 and 2, respectively, whereas the
relative distance error is slightly increased from 12,02% to 12,24% for
dataset 1. It is also worth mentioning that most PHCE s lead to PCE s as well,
about 83% in the case of dataset 1.

Due to these observations, one can conclude that even if the method is not
affected by correlation errors, a better strategy on how to divide the faulted
path prediction problem could be reviewed for potentially better results.

###  7.3 Limitations

As the proposed method relies on simulation, it follows well-known systems
simulation best practices. Indeed, any time the real system changes, the
simulation model should be updated accordingly, the dataset generated once
again, and the ML model retrained.

Additionally, the root of this method relies on the distribution system
numerical model fidelity, in other words, how much the simulated faulted
three-phase voltages mirror the real faults situation. Unpredictable elements
such as interference, device measurement errors, and other disturbances cannot
be simulated and considered by the software, which can consequently impact the
method’s accuracy.

Finally, the dataset relevancy stands as a major drawback. Faults labeled with
path 1, due to its greater length, represent more than 53% of the dataset,
whereas path 3, which is very short, only represents 5% of the dataset. This
unbalanced training dataset leads the network to overfit on the most
represented category, that is, path 1, and hence inaccurate for prediction
purposes. Moreover, the use of specialized ANN s significantly reduces the
size of the dataset that they can use. Indeed, the faulted path classification
ANN s can only use as few as 408 or 504 of the 6384 samples generated, as
shown in Table 3. It is an important limitation knowing that the size of the
dataset is an important criterion for building robust and accurate ANN s.
However, the two former drawbacks aren’t directly associated with the method
itself since new datasets can be generated on demand.

###  7.4 Method Relevancy

The motivation for studying a new approach to fault location is to overcome or
reduce the limitations of classical methods as well as improve their accuracy.

At first sight, this method would need measurement devices such as transducers
only at the system substation. Unlike TWBFL methods, high-frequency
measurement isn’t necessarily required since the wavelet coefficients used
here cover a frequency band in the range of kHz, easing costs for measurement
devices.

Accurate knowledge about the network and associated components is also
required, similarly to IBFL methods. This method stands out in its
methodology, which remains the same for any kind of system and topology.
Moreover, the short execution time of the method, once the machine learning
model is trained, is appreciable, with the algorithm running in less than half
a second given the three-phase voltages raw data.

Finally, although at distribution levels IBFL methods are the most accurate
nowadays, the proposed method shows competitiveness in terms of accuracy, not
to mention the potential along with future improvements.

###  7.5 Towards Improvements

The main weakness of the method is its lack of robustness. A larger and more
diversified dataset appears as a natural solution to reduce this weakness. In
addition, the original dataset remains small when it comes to specialized ANN
s that only use a part of this dataset. Finally, unbalanced datasets also
negatively influence the overall accuracy of faulted path classification. To
this extent, the first step toward the improvement of the method accuracy
would be to retrain the related ANN s with a more complete, balanced, and
larger dataset as well as to improve the faulted path location ANN set
structure on of a more balanced (or a different) dataset.

Besides data quality, other factors can affect the accuracy of a given ANN,
such as their hyper-parameters or their number of layers. An ANN structure
remains basic compared to the power of deep learning. In that sense, trials
could be done on more complex ANN.

Finally, the accuracy of the distribution system modeling plays a key role in
the method, which is why further work should include frequency-dependent
component modeling, such as load and power lines, to gain model veracity and,
at the same time, observe the impact on simulation time, an important factor
for large dataset generation.

##  8 Conclusion

In this paper, a data-driven ground-fault location method in a distribution
power system is explained and tested. The transient faulted three-phase
voltage waveforms at the substation contain implicit information on fault
location. These signals are processed using discrete wavelet transform.
Mathematical statistic features of the wavelet coefficients, energy content,
and three-phase faulted voltages are then calculated.

We show that these features can help predict the fault location with a
relative error of 0,4% by using a workflow of several ANN s. Faulted phases
are classified with an accuracy of 100% whereas results for faulted path
classification are lower with 75,15% correct prediction.

Even though the obtained results are promising and show the potential of the
method, there is a need for future work to improve the method’s accuracy and
robustness. In particular, we plan to generate a larger and more diversified
dataset, as well as explore other faulted path classification ANN workflow.

Finally, leveraging the results obtained in this paper, future research
directions will focus on the development of digital twin solutions able to
deal with real-life issues, such as the potential interferences and
disturbances that can occur in the distribution system, and limit the method
implementation.

##  9 Acknowledgments

This research is partially funded by the Swedish Knowledge Foundation, Grant
No. 20150088 “Software Technology for Self-Adaptive Systems” and Grant No.
20200117 “Aligning Architectures for Digital twin of the Organization”.

## References

  * [1] M. A. Aftab, S. S. Hussain, I. Ali, and T. S. Ustun.  Dynamic protection of power systems with high penetration of renewables: A review of the traveling wave based fault location techniques.  International Journal of Electrical Power & Energy Systems, 114, 2020.
  * [2] A. Bahmanyar, S. Jamali, A. Estebsari, and E. Bompard.  A comparison framework for distribution system outage and fault location methods.  Electric Power Systems Research, 145:19–34, 2017.
  * [3] A. Borghetti, S. Corsi, C. A. Nucci, M. Paolone, L. Peretto, and R. Tinarelli.  On the use of continuous-wavelet transform for fault location in distribution power systems.  International Journal of Electrical Power & Energy Systems, 28(9), 2006.
  * [4] S. M. Brahma.  Fault location in power distribution system with penetration of distributed generation.  IEEE Transactions on Power Delivery, 26(3):1545–1553, 2011.
  * [5] X. Chen, X. Yin, X. Yin, J. Tang, and M. Wen.  A novel traveling wave based fault location scheme for power distribution grids with distributed generations.  In 2015 IEEE Power & Energy Society General Meeting, 2015.
  * [6] E. F. Ferreira and J. D. Barros.  Faults monitoring system in the electric power grid of medium voltage.  Procedia Computer Science, 130:696–703, 2018.  The 9th International Conference on Ambient Systems, Networks and Technologies (ANT 2018) / The 8th International Conference on Sustainable Energy Information Technology (SEIT-2018) / Affiliated Workshops.
  * [7] G. Ferreira, D. Gazzana, A. Bretas, and A. Netto.  A unified impedance-based fault location method for generalized distribution systems.  pages 1–8, 07 2012.
  * [8] A. Garg and B. Panigrahi.  Fault classification, location in a series compensated power transmission network using online sequential extreme learning machine.  In 2016 IEEE 1st International Conference on Power Electronics, Intelligent Control and Energy Systems (ICPEICES), pages 1–6. IEEE, 2016.
  * [9] M. Goudarzi, B. Vahidi, R. A. Naghizadeh, and S. H. Hosseinian.  Improved fault location algorithm for radial distribution systems with discrete and continuous wavelet analysis.  International Journal of Electrical Power & Energy Systems, 67, 2015.
  * [10] F. H. Magnago and A. Abur.  Fault location using wavelets.  IEEE Transactions on Power Delivery, 13(4):1475–1480, 1998.
  * [11] Y. Menchafou, H. El Markhi, M. Zahri, and M. Habibi.  Impact of distributed generation integration in electric power distribution systems on fault location methods.  In 2015 3rd International Renewable and Sustainable Energy Conference (IRSEC), 2015.
  * [12] J. Mora-Florez, J. Melendez, and G. Carrillo-Caicedo.  Comparison of impedance based fault location methods for power distribution systems.  Electric power systems research, 78(4), 2008.
  * [13] J. U. N. Nunes and A. S. Bretas.  Impedance-based fault location formulation for unbalanced primary distribution systems with distributed generation.  In 2010 International Conference on Power System Technology, pages 1–7, 2010.
  * [14] M. Pourahmadi-Nakhli and A. A. Safavi.  Path characteristic frequency-based fault locating in radial distribution systems using wavelets and neural networks.  IEEE Transactions on Power Delivery, 26(2):772–781, 2010.
  * [15] P. Stefanidou-Voziki, N. Sapountzoglou, B. Raison, and J. Dominguez-Garcia.  A review of fault location and classification methods in distribution grids.  Electric Power Systems Research, 209:108031, 2022.
  * [16] D. Thukaram, H. Khincha, and H. Vijaynarasimha.  Artificial neural network and support vector machine approach for locating faults in radial distribution systems.  IEEE Transactions on Power Delivery, 20(2):710–721, 2005.
  * [17] L. Ye, D. You, X. Yin, K. Wang, and J. Wu.  An improved fault-location method for distribution system using wavelets and support vector regression.  International Journal of Electrical Power & Energy Systems, 55:467–472, 2014.

##  Appendix A Distribution Power System Data

Table 7: Generation data of test distribution power system Infinite Grid
---
$V_{G}=20$ kV, $f=60$ Hz, $X/R$ ratio $=10$
Synchronous Generator
$V_{DG}=3,5$ kV, $f=60$ Hz, $P_{nominal}=20$ MVA
Table 8: Transformers data of test distribution power system Transformers | Nominal | Voltage | Connection
---|---|---|---
Power
(MVA) | (kV)
Infinite Grid-Bus 1 | 60 | 63/20 | Yg-Yg
Bus 8-DG | 10 | 3.5/20 | Yg-Yg
Bus 3-Load 3 | 10 | 20/0,4 | Yg-Yg
Bus 4-Load 4 | 10 | 20/0,4 | Yg-Yg
Bus 6-Load 6 | 10 | 20/0,4 | Yg-Yg
Bus 8-Load 8 | 10 | 20/0,4 | Yg-Yg
Bus 10-Load 10 | 10 | 20/0,4 | Yg-Yg
Bus 11-Load 11 | 10 | 20/0,4 | Yg-Yg
Table 9: Loads data of test distribution power system Load | Model | Connection | P | Q
---|---|---|---|---
Bus |  |  | (MW) | (kVar)
3 | Static, constant Z | Yg | 7,5 | 42,5
4 | Static, constant Z | Yg | 7,5 | 42,5
6 | Static, constant Z | Yg | 3,6 | 12,5
8 | Static, constant Z | Yg | 7,5 | 42,5
10 | Static, constant Z | Yg | 7,5 | 42,5
11 | Static, constant Z | Yg | 3,6 | 12,5
Total | 37.2 | 195
Table 10: Power lines data of test distribution power system Line | DC | Inner | Outer
---|---|---|---
Type |  Resistance $(\Gamma/km)$ | Radius | Radius
| $(\Gamma/km)$ | (cm) | (cm)
J Marti | 0,065 | 0,14 | 0,8
Path | Line between Buses | Length (km)
1 | 1-2 | 2
1 | 2-5 | 1,5
1 | 5-7 | 1,5
1 | 7-9 | 3
1 | 9-11 | 2
2 | 2-3 | 2
3 | 2-4 | 1
4 | 5-6 | 1
5 | 7-8 | 2
6 | 9-10 | 3
Figure 11: Power lines model of test distribution power system. Table 11: Relevant features for optimized workflow ANNs |  Relevant Features
---|---
Tfp |  $std(V_{a})$, $std(V_{b})$, $std(V_{c})$, $std(V_{a}CD_{8})$, $std(V_{b}CD_{8})$, $std(V_{c}CD_{8})$, $std(V_{a}CA_{8})$, $std(V_{b}CA_{8})$, $std(V_{c}CA_{8})$, $E_{wave}(V_{a}CD_{8})$, $E_{wave}(V_{b}CD_{8})$, $E_{wave}(V_{c}CD_{8})$, $E_{wave}(V_{a}CA_{8})$, $E_{wave}(V_{b}CA_{8})$, $E_{wave}(V_{c}CA_{8})$
Ofd-i |  $var(V_{i})$, $var(V_{i}CD_{8})$, $skn(V_{i})$, $skn(V_{i}CD_{8})$, $skn(V_{i}CA_{8})$, $cm3(V_{i})$, $cm3(V_{i}CD_{8})$, $cm3(V_{i}CA_{8})$, $E_{wave}(V_{i}CD_{8})$
Ofd-ij |  $var(V_{i}CD_{8})$, $var(V_{j}CD_{8})$, $skn(V_{i})$, $skn(V_{j})$, $skn(V_{i}CA_{8})$, $skn(V_{j}CA_{8}$, $cm3(V_{i})$, $cm3(V_{j})$, $cm3(V_{i}CD_{8})$, $cm3(V_{j}CD_{8})$, $cm3(V_{i}CA_{8})$, $cm3(V_{j}CA_{8})$, $mode(V_{i}CD_{8})$, $mode(V_{j}CD_{8})$, $E_{wave}(V_{i}CD_{8})$, $E_{wave}(V_{j}CD_{8})$
Ofd-ABC |  $var(V_{A}CD_{8})$, $var(V_{B}CD_{8})$, $var(V_{C}CD_{8})$, $var(V_{1}CD_{8})$, $var(V_{2}CD_{8})$, $var(V_{1}CA_{8})$, $var(V_{2}CA_{8})$, $skn(V_{A})$, $skn(V_{B})$, $skn(V_{C})$, $skn(V_{1})$, $skn(V_{A}CA_{8})$, $skn(V_{B}CA_{8})$, $skn(V_{C}CA_{8})$, $skn(V_{1}CA_{8})$, $skn(V_{2}CA_{8})$, $cm3(V_{A})$, $cm3(V_{B})$, $cm3(V_{C})$, $cm3(V_{1})$, $cm3(V_{A}CD_{8})$, $cm3(V_{B}CD_{8})$, $cm3(V_{C}CD_{8})$, $cm3(V_{A}CA_{8})$, $cm3(V_{B}CA_{8})$, $cm3(V_{C}CA_{8})$, $cm3(V_{1}CD_{8})$, $cm3(V_{2}CD_{8})$, $mode(V_{A}CD_{8})$, $mode(V_{B}CD_{8})$, $mode(V_{C}CD_{8})$, $mode(V_{A}CA_{8})$, $mode(V_{B}CA_{8})$, $mode(V_{C}CA_{8})$, $mode(V_{2}CA_{8})$, $E_{wave}(V_{A}CD_{8})$, $E_{wave}(V_{B}CD_{8})$, $E_{wave}(V_{C}CD_{8})$, $E_{wave}(V_{A}CA_{8})$, $E_{wave}(V_{B}CA_{8})$, $E_{wave}(V_{C}CA_{8})$, $E_{wave}(V_{1}CD_{8})$, $E_{wave}(V_{2}CD_{8})$, $E_{wave}(V_{1}CA_{8})$, $E_{wave}(V_{2}CA_{8})$
Ofp-i$H_{1}$ & Ofp-i$H_{2}$ |  $var(V_{i}CD_{8})$, $var(V_{i}CA_{8})$, $skn(V_{i})$, $skn(V_{i}CD_{8})$, $cm3(V_{i})$, $cm3(V_{i}CD_{8})$, $mode(V_{i})$, $mode(V_{i}CA_{8})$
Ofp-ij$H_{1}$ & Ofp-ij$H_{2}$ |  $var(V_{i})$, $var(V_{j})$, $var(V_{i}CD_{8})$, $var(V_{j}CD_{8})$, $var(V_{i}CA_{8})$, $var(V_{j}CA_{8})$, $skn(V_{i})$, $skn(V_{j})$, $skn(V_{i}CD_{8})$, $skn(V_{j}CD_{8})$, $cm3(V_{i})$, $cm3(V_{j})$, $cm3(V_{i}CD_{8})$, $cm3(V_{j}CD_{8})$, $cm3(V_{i}CA_{8})$, $cm3(V_{j}CA_{8})$, $mode(V_{i}CD_{8})$, $mode(V_{j}CD_{8})$, $E_{wave}(V_{i}CD_{8})$, $E_{wave}(V_{j}CD_{8})$, $E_{wave}(V_{i}CA_{8})$, $E_{wave}(V_{j}CA_{8})$
Ofp-ABC$H_{1}$ & Ofp-ABC$H_{2}$ |  $var(V_{2}CD_{8})$, $skn(V_{A})$, $skn(V_{B})$, $skn(V_{C})$, $skn(V_{0})$, $skn(V_{1})$, $skn(V_{A}CD_{8})$, $skn(V_{B}CD_{8})$, $skn(V_{C}CD_{8})$, $skn(V_{1}CD_{8})$, $skn(V_{2}CD_{8})$, $cm3(V_{A})$, $cm3(V_{B})$, $cm3(V_{C})$, $cm3(V_{0})$, $cm3(V_{1})$, $cm3(V_{A}CD_{8})$, $cm3(V_{B}CD_{8})$, $cm3(V_{C}CD_{8})$, $cm3(V_{1}CD_{8})$, $cm3(V_{2}CD_{8})$, $mode(V_{A}CD_{8})$, $mode(V_{B}CD_{8})$, $mode(V_{C}CD_{8})$, $mode(V_{2}CD_{8})$, $E_{wave}(V_{2}CD_{8})$


# Origin of the Ising Ferrimagnetism and Spin-Charge Coupling in LuFe2O4

H. J. Xiang  National Renewable Energy Laboratory, Golden, Colorado 80401, USA
E. J. Kan  Department of Chemistry, North Carolina State University, Raleigh,
North Carolina 27695-8204, USA  Su-Huai Wei  National Renewable Energy
Laboratory, Golden, Colorado 80401, USA  M.-H. Whangbo  Department of
Chemistry, North Carolina State University, Raleigh, North Carolina
27695-8204, USA  Jinlong Yang  Hefei National Laboratory for Physical Sciences
at Microscale, University of Science and Technology of China, Hefei, Anhui
230026, P. R. China

###### Abstract

The spin ordering and spin-charge coupling in LuFe2O4 were investigated on the
basis of density functional calculations and Monte Carlo simulations. The 2:1
ferrimagnetism arises from the strong antiferromagnetic intra-sheet Fe3+-Fe3+
and Fe3+-Fe2+ as well as some substantial antiferromagnetic Fe2+-Fe3+ inter-
sheet spin exchange interactions. The giant magnetocapacitance at room
temperature and the enhanced electric polarization at 240 K of LuFe2O4 are
explained by the strong spin-charge coupling.

###### pacs:

75.80.+q,71.20.-b,77.80.-e,64.60.De

Recently, multiferroics Kimura2003 ; Hur2004 ; Ikeda2005 ; Xu2008 ;
Subramanian2006 ; Zhang2008 ; Xiang2007A ; Xiang2007B ; Xiang2008 ; Angst2008
have attracted much attention because of their potential applications in novel
magnetoelectric and magneto-optical devices. Among the newly discovered
multiferroics, LuFe2O4 is particularly interesting due to its large
ferroelectric (FE) polarization Ikeda2005  and giant magnetocapacitance at
room temperature Subramanian2006 . In the high-temperature crystal structure
of LuFe2O4 with space group R3̄m, layers of composition Fe2O4 alternate with
layers of Lu3+ ions, such that there are three Fe2O4 layers per unit cell.
Each Fe2O4 layer is made up of two triangular sheets (hereafter, T-sheets) of
corner-sharing FeO5 trigonal bipyramids (Fig. 1). Below 320 K ($T_{CO}$)
LuFe2O4 undergoes a three-dimensional (3D) charge ordering (CO) (2Fe2.5+
$\Rightarrow$ Fe2+ \+ Fe3+) with the $\sqrt{3}\times\sqrt{3}$ superstructure
in each T-sheet; in each Fe2O4 layer, one T-sheet has the honeycomb network of
Fe2+ ions with a Fe3+ ion at the center of each Fe2+ hexagon (hereafter, the
type A T-sheet), while the other T-sheet has an opposite arrangement of the
Fe2+ and Fe3+ ions (hereafter the type B T-sheet).

LuFe2O4, with the novel CO-driven “electronic ferroelectricity”, Ikeda2005
presents several fundamental questions. First, LuFe2O4 shows strong Ising
behavior with the easy axis along $c$ Iida1993 ; Wu2008 . The spin anisotropy
of the non-CO state is understandable because the spin down electron of the
Fe2.5+ ion partially occupies the degenerate ($d_{x^{2}-y^{2}}$,$d_{xy}$)
orbitals Xiang2007A ; Dai2005 . However, the Ising behavior below $T_{CO}$ is
puzzling because the insulating $\sqrt{3}\times\sqrt{3}$ CO breaks the 3-fold
rotational symmetry hence lifting the degeneracy of the
($d_{x^{2}-y^{2}}$,$d_{xy}$) orbitals Xiang2007A . Second, LuFe2O4 undergoes a
ferrimagnetic spin ordering below 240 K ($T_{N}$) Iida1993 ; Tanaka1989 ;
Siratori1992 ; Christianson2008 . A number of experimental studies found this
spin ordering to be two-dimensional (2D) in nature Iida1993 ; Tanaka1989 ;
Funagashi1984 . In contrast, a recent neutron diffraction study observed a
finite spin correlation along $c$ and suggested a 3D spin structure without
considering CO Christianson2008 . The Mössbauer Tanaka1989  and neutron
diffraction Siratori1992  studies led to a detailed ferrimagnetic structure of
LuFe2O4, in which the majority spin lattice consists of all Fe2+ ions plus
one-third of the total Fe3+ ions while the minority spin sublattice consists
of the remaining Fe3+ ions. This 2:1 ferrimagnetic order was suggested to
originate from weak ferromagnetic (FM) interactions between the next-nearest
neighbor (NNN) Fe sites in the triangular antiferromagnetic (AFM) Ising
lattice Iida1993 . However, using the spin exchange parameters estimated from
the energy parameters of LaFeO3, Naka et al. Naka2008  predicted quite a
different spin structure that includes some Fe sites without unique spin
direction. Therefore, the detailed ferrimagnetic structure and its origin
remain unclear. Third, LuFe2O4 exhibits a giant magnetodielectric response at
room temperature Subramanian2006 , and a room-temperature dynamic
magnetoelectric coupling was also reported Park2007 . Furthermore, the FE
polarization of LuFe2O4 was found to increase around $T_{N}$ Ikeda2005 . These
observations suggest the occurrence of coupling between the CO and magnetism.
The understanding of the spin-charge coupling is crucial for future
magnetodielectric applications of LuFe2O4.

In this Letter, we explore these isuues on the basis of first principles
density functional calculations for the first time. A large spin anisotropy is
found along the $c$ direction due mainly to the Fe2+ ions of the B-sheet, the
spin ground state of the $\sqrt{3}\times\sqrt{3}$ CO state has the 2:1
ferrimagnetic spin arrangement proposed by Siratori et al. Siratori1992 , and
there occurs strong spin-charge coupling in LuFe2O4.

Our density functional theory calculations employed the frozen-core projector
augmented wave method PAW  encoded in the Vienna ab initio simulation package
VASP , and the generalized-gradient approximation (GGA) Perdew1996 . To
properly describe the strong electron correlation in the 3d transition-metal
oxide, the GGA plus on-site repulsion U method (GGA+U) Liechtenstein1995  was
employed with the effective $U$ value ($U_{eff}=U-J$ with $J=0$) of 4.61 eV
Xiang2007A . It is known experimentally Iida1993 ; Tanaka1989 ; Funagashi1984
that the interlayer magnetic interactions in LuFe2O4 are weak, which is
understandable due to its layered structure. In this work, therefore, we focus
on the 2D spin ordering within a single Fe2O4 layer. For the
$\sqrt{3}\times\sqrt{3}$ CO state of LuFe2O4, the FE ordering of the Fe2O4
layers will be assumed.

We first examine the magnetic anisotropy of the Fe ions by performing GGA+U
calculations, with spin-orbit coupling (SOC) included, for the FM state of
LuFe2O4 with the $\sqrt{3}\times\sqrt{3}$ CO. As shown in Fig. 1(a), there are
two kinds of Fe2+ ions and two kinds of Fe3+ ions in the
$\sqrt{3}\times\sqrt{3}$ CO state. We label the Fe2+ and Fe3+ ions of the type
A T-sheet as 2A and 3A, respectively, and those of the type B T-sheet as 2B
and 3B, respectively. In our GGA+U+SOC calculations with spins pointing along
several different directions, all Fe2+ and Fe3+ spins are kept in the same
direction. Our calculations show that the easy axis is along the $c$
direction, as experimentally observed Iida1993 ; Wu2008 ; the
$\parallel$c-spin orientation is more stable than the $\perp$c-spin
orientation by 1.5 meV per formula unit (FU). The orbital moments of 2A, 2B,
3A, and 3B for the $\parallel$c-spin orientation are 0.101, 0.156, 0.031 and
0.035, respectively, which are greater than those for the $\perp$c-spin
orientation by 0.019, 0.062, 0.015, and 0.018 $\mu_{B}$, respectively. As
expected, the Fe3+ ($d^{5}$) ions have a very small anisotropy, However, two
kinds of the Fe2+ ions also have different degree of spin anisotropy. The spin
down electron of the 2B Fe2+ ion occupies the ($d_{x^{2}-y^{2}}$,$d_{xy}$)
manifold Xiang2007A , therefore the 2B Fe2+ ion has the largest spin
anisotropy along $c$. Our calculations indicate a non-negligible orbital
contribution to the total magnetization, in agreement with the X-ray magnetic
circular dichroism result Wu2008 .

To determine the magnetic ground state of LuFe2O4 in the
$\sqrt{3}\times\sqrt{3}$ CO state, we extract its spin exchange parameters by
mapping the energy differences between ordered spin states obtained from GGA+U
calculations onto the corresponding energy differences obtained from the Ising
Hamiltonian whangbo2003 :

$H=\sum_{i,j}J_{ij}S_{iz}S_{jz},$ (1)

where the energy is expressed with respect to the spin disorder (paramagnetic)
state, $J_{ij}$ is the spin exchange parameter between the spin sites $i$ and
$j$, and $S_{iz}$ is the spin component along the $c$ direction ($|S_{z}|=2$
and $2.5$ for Fe2+ and Fe3+ ions, respectively). We consider all 15 possible
superexchange (SE) interactions and all 19 super-superexchange (SSE)
interactions with the O…O distance less than 3.2 Å. The intra- and inter-sheet
interactions within each Fe2O4 layer as well as the SSE interactions between
adjacent Fe2O4 layers are taken into account. To evaluate these 34 spin
exchange parameters reliably, we considered 111 different ordered spin states
leading to 110 energy differences. The 34 spin exchange parameters were
determined by performing a linear least-square fitting analysis. The SSE
interactions are generally much weaker than the SE interactions with the
magnitude of all SSE interactions less than 1.4 meV. The calculated SE
parameters are reported in Table 1. All intra-sheet SE interactions are AFM,
and the strongest interactions ($\sim 7.3$ meV) occurs between the 3B Fe3+
ions because of the large energy gain of the AFM configuration and almost zero
FM coupling. The inter-sheet SE interactions are weaker than the the intra-
sheet SE interactions, and are mostly AFM.

With the calculated spin exchange parameters, one can identify the spin ground
state of the CO state. The Metropolis Monte Carlo simulation of the Ising
model is performed to search for the ground state. Simulations with supercells
of several different sizes show that the spin ground state has the magnetic
structure shown in Fig. 2(a), which has the same cell as the
$\sqrt{3}\times\sqrt{3}$ CO structure. In this state, all Fe2+ ions contribute
to the majority spin, and the Fe3+ ions are antiferromagnetically coupled to
the Fe2+ ions in the type A T-sheet. In the honeycomb lattice of the type B
T-sheet, the Fe3+ spins are antiferromagnetically coupled. Thus, the spin
ground state is ferrimagnetic, as experimentally observed Iida1993 . This 2:1
ferrimagnetic structure is the same as the magnetic structure proposed by
Siratori et al. Siratori1992 , and differs from the structure proposed by Naka
et al. Naka2008 .

The observed ferrimagnetic ordering can be readily explained in terms of the
calculated exchange parameters. In the honeycomb network of the type B
T-sheet, the nearest-neighbor (NN) 3B ions are antiferromagnetically coupled
since their SE interaction is strongly AFM. In the type A T-sheet, the SE
interactions between the 2A ions are AFM, and so are those between the 2A and
3A ions, which leads to spin frustration. As a consequence, two possible spin
arrangements compete with each other in the type A T-sheet; the first is the
state in which the coupling between the NN 2A ions are AFM with the spin
direction of the 3A ion undetermined, and the second is the state in which all
2A ions are antiferromagnetically coupled to the 3A ions. The energies of
these two states (considering only the SE interaction) are
$E_{1}=-4(J_{2A1,2A2}+J_{2A1,2A4})$ per 3A ion, and
$E_{2}=-10(J_{3A1,2A1}+J_{3A1,2A2}+J_{3A1,2A3})+4(J_{2A1,2A2}+J_{2A1,2A4})$
per 3A ion, respectively. Due to the relatively strong AFM interactions
between the 3A and 2A ions (See Table 1) and the large spin of the 3A ions,
the second state has a lower energy, i.e., $E_{2}$ $<$ $E_{1}$. Without loss
of generality, we can assume the 2A (3A) ions constitute the majority
(minority) spin in the second state. Now, we examine the spin orientation of
the Fe2+ ions in the type B T-sheet. The intra-sheet interactions of the 2B
ion with 3B ions vanish due to the AFM ordering of the 3B ions. As for the
inter-sheet interactions involving the 2B ions, the dominant one is the AFM
interaction of the 2B ion with the 3A ion ($J_{3A1-2B1}$ in Table 1).
Consequently, we obtain the ferrimangetic ground state shown in Fig. 2(a), in
which the spin of the 2B ion contributes to the majority spin of the Fe2O4
layer. For the stability of the ferrimangetic ground state, the inter-sheet
interaction is essential. This was neglected in the model Hamiltonian study of
Naka et al. Naka2008 . The ferrimangetic state is not due to the FM
interactions between NNN Fe ions of the T-sheet because they must be
vanishingly weak and mostly AFM.

The electronic structure of the ferrimangetic state calculated for the
$\sqrt{3}\times\sqrt{3}$ CO structure of LuFe2O4 is shown in Fig. 3. Also
shown is the electronic structure calculated for the FM state. Both states are
semiconducting, and the highest occupied (HO) and the lowest unoccupied (LU)
levels of both states come from the spin-up Fe2+ and Fe3+ ions, respectively
Xiang2007A . In addition, the band dispersion from $\Gamma$ to A is rather
small, indicating a very weak interlayer interaction. However, there are some
important differences. First, the ferrimangetic state has a larger band gap
(1.68 eV) than does the FM state (0.77 eV). This is consistent with the
stability of the ferrimangetic state. Second, the FM state has an indirect
band gap with the HO and LU levels located at K and $\Gamma$, respectively. In
the ferrimangetic state, however, the LU level has the highest energy at
$\Gamma$ and the band dispersions of the HO and LU levels are almost flat from
M to K. This difference comes from the orbital interaction between the spin
down ($d_{x^{2}-y^{2}}$,$d_{xy}$) levels of the spin up Fe3+ and Fe2+ ions.

To probe the presence of spin-charge coupling in LuFe2O4, it is necessary to
consider the spin ordering in a CO state other than the
$\sqrt{3}\times\sqrt{3}$ CO state. The previous electrostatic calculations
Xiang2007A ; Naka2008  showed that the chain CO, in which one-dimensional (1D)
chains of Fe2+ ions alternate with 1D chains of Fe3+ ions in each T-sheet
[Fig. 2(b)], is only slightly less stable than the $\sqrt{3}\times\sqrt{3}$
CO, and has no FE polarization. We extract exchange parameters by mapping
analysis as described above. It is found that the intra-sheet SE between the
Fe3+ ions is the strongest ($J=6.7$ meV) as in the $\sqrt{3}\times\sqrt{3}$ CO
case. All intra-sheet SE’s are AFM with $J$(Fe3+-Fe3+) $>$ $J$(Fe2+-Fe3+) $>$
$J$(Fe2+-Fe2+). The inter-sheet SE between the Fe3+ ions is very weak
($|J|<0.3$ meV), and that between the Fe2+ and Fe3+ ions is FM with $J=-1.4$
meV. Interestingly, the inter-sheet SE between the Fe2+ ions is rather
strongly AFM ($J=6.3$ meV). Monte Carlo simulations using these spin exchange
parameters indicate that the spin state shown in Fig. 2(b) is the spin ground
state. In this spin ordering, the spins within each chain of Fe2+ ions or Fe3+
ions are antiferromagnetically coupled. The NN chains of Fe2+ ions belonging
to different T-sheets are coupled antiferromagnetically, whereas the
corresponding chains of Fe3+ are almost decoupled.

The above results show that the spin ordering of the chain CO state is
dramatically different from that of the $\sqrt{3}\times\sqrt{3}$ CO state. The
most important difference is that the total spin moments are 2.33 $\mu_{B}$/FU
for the $\sqrt{3}\times\sqrt{3}$ CO, but 0 $\mu_{B}$/FU for the chain CO. This
evidences a strong spin-charge coupling in LuFe2O4. The external magnetic
field will have different effects on the two CO states due to the the Zeeman
effect. It is expected that the magnetic field will further stabilize the
ferrimagnetic $\sqrt{3}\times\sqrt{3}$ CO state. Consequently, an external
magnetic field will reduce the extent of charge fluctuation and hence decrease
the dielectric constant. This supports our explanation for the giant
magnetocapacitance effect of LuFe2O4 at room temperature Xiang2007A  .

Without considering the inter-sheet interactions, Naka et al. Naka2008
suggested that the degeneracy of the spin ground state of the
$\sqrt{3}\times\sqrt{3}$ CO state is of the order O($2^{N/3}$)( N is the
number of the spin sites), which is much larger than the spin degeneracy
[O($2^{\sqrt{N}}$)] of the chain CO state. Thus, they proposed that spin
frustration induces reinforcement of the polar $\sqrt{3}\times\sqrt{3}$ CO by
a gain of spin entropy. However, our calculations show that there are
substantial inter-sheet spin exchange interactions between the 2B1 and 3A1
ions, which would remove the macroscopic degeneracy of the spin ground state
of the $\sqrt{3}\times\sqrt{3}$ CO state. The macroscopic degeneracy still
persists for the chain CO state. Thus, our work provides a picture opposite to
what Naka et al. proposed. Furthermore, we find that the
$\sqrt{3}\times\sqrt{3}$ CO state is more favorable for the spin ordering than
is the chain CO state; with respect to the paramagnetic state, the spin ground
state is lower in energy by $-78$ meV/FU for the $\sqrt{3}\times\sqrt{3}$ CO,
but by $-57$ meV/FU for the chain CO. The model of Naka et al. Naka2008
predicts that the polar $\sqrt{3}\times\sqrt{3}$ CO state is destabilized and
the electric polarization is reduced by the magnetic field, since it will lift
the macroscopic spin degeneracy. In contrast, our work predicts that the
magnetic field stabilizes the ferrimagnetic $\sqrt{3}\times\sqrt{3}$ CO state
due to the Zeeman effect, and provides an explanation for why the electric
polarization increases when the temperature is lowered below the Neel
temperature Ikeda2005 , because the charge fluctuation has an onset well below
$T_{CO}$ Xu2008 .

In summary, our first principles results explain the experimentally observed
Ising ferrimagnetism, and manifest the spin-charge coupling and
magnetoelectric effect in LuFe2O4.

Work at NREL was supported by the U.S. Department of Energy, under Contract
No. DE-AC36-08GO28308, and work at NCSU by the U. S. Department of Energy,
under Grant DE-FG02-86ER45259.

## References

  * (1) T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, and Y. Tokura, Nature (London) 426, 55 (2003).
  * (2) N. Hur, S. Park, P. A. Sharma, J. S. Ahn, S. Guha, and S-W. Cheong, Nature (London) 429, 392 (2004).
  * (3) N. Ikeda, H. Ohsumi, K. Ohwada, K. Ishii, T. Inami, K. Kakurai, Y. Murakami, K. Yoshii, S. Mori, Y. Horibe, and H. Kitô, Nature (London) 436, 1136 (2005).
  * (4) M. A. Subramanian, T. He, J. Chen, N. S. Rogado, T. G. Calvarese, and A. W. Sleight, Adv. Mater. 18, 1737 (2006).
  * (5) H. J. Xiang and M.-H. Whangbo, Phys. Rev. Lett. 98, 246403 (2007).
  * (6) Y. Zhang, H. X. Yang, C. Ma, H. F. Tian, and J. Q. Li, Phys. Rev. Lett. 98, 247602 (2007).
  * (7) M. Angst, R. P. Hermann, A. D. Christianson, M. D. Lumsden, C. Lee, M.-H. Whangbo, J.-W. Kim, P. J. Ryan, S. E. Nagler, W. Tian, R. Jin, B. C. Sales, and D. Mandrus, Phys. Rev. Lett. 101, 227601 (2008).
  * (8) X. S. Xu, M. Angst, T. V. Brinzari, R. P. Hermann, J. L. Musfeldt, A. D. Christianson, D. Mandrus, B. C. Sales, S. McGill, J.-W. Kim, and Z. Islam, Phys. Rev. Lett. 101, 227602 (2008).
  * (9) H. J. Xiang and M.-H. Whangbo, Phys. Rev. Lett. 99, 257203 (2007).
  * (10) H. J. Xiang, S.-H. Wei, M.-H. Whangbo, and J. L. F. Da Silva, Phys. Rev. Lett. 101, 037209 (2008).
  * (11) J. Iida, M. Tanaka, Y. Nakagawa, S. Funahashi, N. Kimizuka, and S. Takekawa, J. Phys. Soc. Jpn. 62, 1723 (1993).
  * (12) W. Wu, V. Kiryukhin, H.-J. Noh, K.-T. Ko, J.-H. Park, W. Ratcliff II, P. A. Sharma, N. Harrison, Y. J. Choi, Y. Horibe, S. Lee, S. Park, H. T. Yi, C. L. Zhang, and S.-W. Cheong, Phys. Rev. Lett. 101, 137203 (2008).
  * (13) D. Dai and M.-H. Whangbo, Inorg. Chem. 44, 4407 (2005).
  * (14) M. Tanaka, H. Iwasaki, K. Siratori, and I. Shindo, J. Phys. Soc. Jpn. 58, 1433 (1989).
  * (15) K. Siratori, S. Funahashi, J. Iida, and M. Tanaka, Proc. 6th Intern. Conf. Ferrites, Tokyo and Kyoto, Japan, 1992, p. 703.
  * (16) A. D. Christianson, M. D. Lumsden, M. Angst, Z. Yamani, W. Tian, R. Jin, E. A. Payzant, S. E. Nagler, B. C. Sales, and D. Mandrus, Phys. Rev. Lett. 100, 107601 (2008).
  * (17) S. Funahashi, J. Akimitsu, K. Siratori, N. Kimizuka, M. Tanaka, and H. Fujishita, J. Phys. Soc. Jpn. 53, 2688 (1984).
  * (18) M. Naka, A. Nagano, and S. Ishihara, Phys. Rev. B 77, 224441 (2008); A. Nagano, M. Naka, J. Nasu, and S. Ishihara, Phys. Rev. Lett. 99, 217202 (2007).
  * (19) J. Y. Park, J. H. Park, Y. K. Jeong, and H. M. Jang, Appl. Phys. Lett. 91, 152903 (2007).
  * (20) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994); G. Kresse and D. Joubert, ibid 59, 1758 (1999).
  * (21) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996); Phys. Rev. B 54, 11169 (1996).
  * (22) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * (23) A. I. Liechtenstein, V. I. Anisimov and J. Zaanen, Phys. Rev. B 52, R5467 (1995); S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
  * (24) M.-H. Whangbo, H.-J. Koo and D. Dai, J. Solid State Chem. 176, 417 (2003).

Table 1: Calculated superexchange parameters (in meV) in the $\sqrt{3}\times\sqrt{3}$ CO state of LuFe2O4 (For the spin sites of the 2A, 3A, 2B and 3B ions,see Fig. 1 ) A-A | $J_{3A1,2A1}$ | $J_{3A1,2A2}$ | $J_{3A1,2A3}$ | $J_{2A1,2A2}$ | $J_{2A1,2A4}$
---|---|---|---|---|---
| 3.2 | 4.0 | 4.7 | 1.9 | 3.6
B-B | $J_{3B1,3B2}$ | $J_{3B1,3B4}$ | $J_{2B1,3B1}$ | $J_{2B1,3B2}$ | $J_{2B1,3B3}$
| 7.0 | 7.6 | 1.5 | 2.8 | 1.3
A-B | $J_{3A1,3B1}$ | $J_{3A1,2B1}$ | $J_{2A1,2B1}$ | $J_{2A1,3B2}$ | $J_{2A1,3B3}$
| 2.0 | 1.9 | $\sim 0$ | $-0.6$ | 1.2

Figure 1: (Color online) Schematic representation of the
$\sqrt{3}\times\sqrt{3}$ CO structure. Large, medium, and small circles
represent the Fe2+, Fe3+, and O2- ions, respectively. The type A (type B)
T-sheet has the honeycomb network of Fe2+ (Fe3+) ions with a Fe3+ (Fe2+) ion
at the center of each hexagon. 2A and 3A (2B and 3B) refer to the Fe2+ and
Fe3+ ions of the type A (type B) T-sheet, respetively. The region enclosed by
dashed lines indicates the unit cell of the CO structure. There is a mirror
plane of symmetry, which is parallel to the $c$ axis and crosses the 3A1 and
2B1 sites. The inset shows an isolated FeO5 trigonal bipyramid. Figure 2:
(Color online) Schematic representations of (a) the spin ground state of the
$\sqrt{3}\times\sqrt{3}$ CO structure and (b) one of the macroscopic spin
ground states of the chain CO structure. The arrows denote the spin
directions. The region enclosed by the dashed lines on the bottom T-sheet
indicates the magnetic unit cell of the spin structure. Figure 3: (Color
online) Band structures calculated for (a) the FM state and (b) the
ferrimagnetic state of the $\sqrt{3}\times\sqrt{3}$ CO structure of LuFe2O4.
The solid and dashed lines represent the up-spin and down-spin bands,
respectively. The $\sqrt{3}\times\sqrt{3}\times 1$ hexagonal cell is used in
the calculations.