**Advanced Approaches, Business Models, and Novel Techniques for Management and Control of Smart Grids**

Printed Edition of the Special Issue Published in *Energies* Pierluigi Siano and Miadreza Shafie-khah Edited by

www.mdpi.com/journal/energies

## **Advanced Approaches, Business Models, and Novel Techniques for Management and Control of Smart Grids**

## **Advanced Approaches, Business Models, and Novel Techniques for Management and Control of Smart Grids**

Editors

**Pierluigi Siano Miadreza Shafie-khah**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Pierluigi Siano University of Salerno Italy

Miadreza Shafie-khah University of Vaasa Finland

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Energies* (ISSN 1996-1073) (available at: https://www.mdpi.com/journal/energies/special issues/ Management Control Smart Grids).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03936-758-0 (Hbk) ISBN 978-3-03936-759-7 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


#### **Shi Chen, Hong Zhou, Jingang Lai, Yiwei Zhou and Chang Yu**


## **About the Editors**

**Pierluigi Siano** received an MSc degree in electronic engineering and a PhD in information and electrical engineering from the University of Salerno, Salerno, Italy, in 2001 and 2006, respectively. He is a Professor and Scientific Director of the Smart Grids and Smart Cities Laboratory with the Department of Management & Innovation Systems, University of Salerno. His research activities are centered on demand response, the integration of distributed energy resources in smart grids and the planning and management of power systems. He has co-authored more than 450 papers, including more than 250 international journal papers that received more than 8400 citations, with an h-index equal to 46. He received the Highly Cited Researcher award from the ISI Web of Science Group in 2019. He has been the Chair of the IES TC on Smart Grids. He is an Editor of the Power & Energy Society Section of *IEEE Access, IEEE Transactions on Industrial Informatics, IEEE Transactions on Industrial Electronics,* the *Open Journal of the IEEE IES and of IET Renewable Power Generation*.

**Miadreza Shafie-khah** received his first postdoc from the University of Beira Interior (UBI), Covilha, Portugal in 2015. He received his second postdoc from the University of Salerno, Salerno, Italy in 2016. Currently, he is a tenure-track Professor at the University of Vaasa, Vaasa, Finland. He is an Associate Editor for IET-RPG and an Editor of the IEEE OAJPE. He was considered one of the Outstanding Reviewers of the *IEEE Transactions on Sustainable Energy* in 2014 and 2017, one of the Best Reviewers of the IEEE Transactions on Smart Grid in 2016 and 2017, and one of the Outstanding Reviewers of the *IEEE Transactions on Power Systems* in 2017 and 2018. He has co-authored more than 310 papers that received more than 4500 citations, with an h-index equal to 39. He is also the editor of the book Blockchain-Based Smart Grids (Elsevier, 2020). His research interests include power market simulation, market power monitoring, power system optimization, demand response, electric vehicles, price and renewable forecasting and smart grids.

## **Special Issue on Advanced Approaches, Business Models, and Novel Techniques for Management and Control of Smart Grids**

#### **Pierluigi Siano 1,\* and Miadreza Shafie-khah 2,\***


Received: 27 April 2020; Accepted: 9 May 2020; Published: 26 May 2020

#### **1. Introduction**

The current power system should be renovated to fulfill social and industrial requests and economic advances. Hence, providing economic, green, and sustainable energy are key goals of advanced societies. In order to meet these goals, the latest features of smart grid technologies need to have the potential to improve reliability, flexibility, efficiency, and resiliency. This *Special Issue* aims to encourage researchers to address the mentioned challenges.

#### **2. Advanced Smart Grids**

Considering the aim of the special issue, ten papers have been published in the area of advanced approaches and novel techniques for the management and control of smart grids, from different standpoints. In reference [1], the classification of energy management systems has been carried out from a new viewpoint. Energy management systems have been classified into four categories based on the type of storage systems. Furthermore, energy management systems have been reviewed in terms of uncertainty modeling techniques, objective functions and constraints, optimization techniques, as well as simulation and experimental results. In order to perform dynamic optimizations in real-time, new solutions are required to consider data privacy and cybersecurity. In reference [2], a paradigm has been designed to meet these requirements based on a bottom-up method to lead to misinterpretations or wrong conclusions. Fractal analysis is also employed to identify unique and independent elements of smart grids required for the design of an architectural paradigm. The processes of demand response and conservation voltage reduction are also presented. In reference [3], the greedy smallest-cost-rate path first algorithm is presented to route power from sources to loads in a digital microgrid. In such a microgrid, one distributed energy resource may supply power to one or many loads, and one or many distributed energy resources may supply the power requested by a load. Reference [3] has addressed the high complexity by using heuristics to match sources and loads and to select the smallest-cost-rate paths in the digital microgrid. Reference [4] has studied the building energy and comfort management systems. These systems have a high potential to decrease energy consumption costs as well as to enhance the efficiency of resource exploitation, by implementing strategies and policies for resource control and demand side management. In order to prevent users' disaffection and the gradual abandonment of the system, reference [4] has developed a sensor-based system for the analysis of users' habits and the early detection and prediction of user activities.

Reference [5] has studied the phase synchronization of an islanded network with large-scale distributed generations in a non-homogeneous condition. The critical coupling strength formula is reduced for different weighting cases via the synchronization condition. On this basis, three cases of Gaussian distribution, power-law distribution, and frequency-weighted distribution have been analyzed. In reference [6], an islanding detection algorithm for distributed generations has been proposed. This work presents a passive islanding detection algorithm for all types of distributed generations by combining the measured frequency, voltage, current, phase angle and the rate of changes at the point of common coupling. The proposed algorithm can detect the islanding situation, even with the exact zero power mismatches. In reference [7], a novel approach has been presented for earth fault location by using a negative sequence current. The developed method is proposed for locating single-phase earth faults on non-effectively earthed medium voltage distribution networks, and it requires only current measurements. The method is based on the analysis of negative sequence components of the currents measured at secondary substations along with medium voltage distribution feeders.

In reference [8], an augmented Prony method has been proposed for power system oscillation analysis using synchrophasor data obtained from a wide-area measurement system. Intrinsic mode functions provide an intuitive representation of the oscillatory modes, and they are mainly calculated by Hilbert–Huang transform methods. However, those methods suffer from the end effects, mode-mixing and Gibbs phenomena. To overcome the drawback, the augmented Prony method has been employed in [8] and tested on an online oscillation-monitoring framework. Another important issue for the operation of islanded microgrids is frequency stability. An inherent delay in the communication system or other parts such as sensor sampling rates may lead microgrids to have unstable operation states. Therefore, it is essential to use an optimization algorithm to determine the optimum parameters. In reference [9], the communication system delay and ultra-capacitor size along with the parameters of the secondary controller have been obtained by using a Non-dominated Sorting Genetic Algorithm II (NSGA-II). Reference [10] has proposed a black-box modeling method to identify the voltage and frequency response model of microgrids. The microgrid has been considered a two-input two-output black-box system and the simplicity of the modeling can be met by input and output ports, while the model parameters can be adjusted by the change of microgrid operating structure, which makes the model more adaptable.

#### **3. Future Approaches, Models, and Techniques**

Although the *Special Issue* has been closed, there is a need for new studies in advanced approaches, business models, and novel techniques for management and control of smart grids. The studies should address the complexity of customers' privacy, climate change, and the role of smart grids for building a greener planet.

**Author Contributions:** Both guest editors have equally carried out the editorial work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** We would like to take this opportunity to record our gratitude to all the authors, hardworking and professional reviewers, and the dedicated editorial team of *Energies*.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Review* **A Survey on Microgrid Energy Management Considering Flexible Energy Sources**

#### **Hossein Shayeghi 1,\*, Elnaz Shahryari 1, Mohammad Moradzadeh <sup>2</sup> and Pierluigi Siano 3,\***


Received: 4 May 2019; Accepted: 31 May 2019; Published: 5 June 2019

**Abstract:** Aggregation of distributed generations (DGs) along with energy storage systems (ESSs) and controllable loads near power consumers has led to the concept of microgrids. However, the uncertain nature of renewable energy sources such as wind and photovoltaic generations, market prices and loads has led to difficulties in ensuring power quality and in balancing generation and consumption. To tackle these problems, microgrids should be managed by an energy management system (EMS) that facilitates the minimization of operational costs, emissions and peak loads while satisfying the microgrid technical constraints. Over the past years, microgrids' EMS have been studied from different perspectives and have recently attracted considerable attention of researchers. To this end, in this paper a classification and a survey of EMSs has been carried out from a new point of view. EMSs have been classified into four categories based on the kind of the reserve system being used, including non-renewable, ESS, demand-side management (DSM) and hybrid systems. Moreover, using recent literature, EMSs have been reviewed in terms of uncertainty modeling techniques, objective functions (OFs) and constraints, optimization techniques, and simulation and experimental results presented in the literature.

**Keywords:** microgrid; energy management system; demand-side management; uncertainty; energy storage; distributed generation

#### **1. Introduction**

There are strong incentives to utilize distributed generations (DGs) for reducing greenhouse gases, improving power system efficiency as well as its reliability, competitive energy policies and postponement of transmission and distribution system upgrading [1]. In fact, DGs are composed of renewable units such as wind turbines (WTs), photovoltaic (PV), fuel cells (FCs), biomass along with non-renewable ones such as micro-turbines (MTs), gas engines (GEs), diesel generators (DiGs), etc. [2]. DGs eliminate the need for the transmission system by being installed near the customers [3]. Integration and control of DGs along with storage devices and flexible loads can constitute a low voltage distribution network, called a microgrid, which can be operated in isolated or grid-connected mode [4]. The generic concept of a microgrid is shown in Figure 1.

**Figure 1.** Microgrid components.

Microgrids often face difficulties in supplying demand due to the lack of sufficient energy generation sources. This obstacle is caused by intermittent nature of loads and renewable energy sources [5]. As a result an energy management system (EMS) is necessary to tackle this problem. EMS for a microgrid represent relatively new and popular topics that attracted lots of attention, recently. Existing review papers in literature have investigated microgrid EMSs from a different aspects, as summarized in Table 1 [6–13].

**Table 1.** A survey on the existing review papers of microgrid energy management system.


The contributions of this review paper are:


The aim of an EMS is to determine the optimal use of DGs in order to feed the electrical loads [14]. An EMS can be operated in two modes, namely centralized and decentralize. In the centralized mode, the central controller aims to optimize the microgrid power exchanged based on the market prices and security constraints. In the decentralized mode, DGs and controllable loads have more degree of freedom [15]. As a result, the microgrid components are considered to be intelligent and try to

maximize the revenue of the microgrid by communicating with each other [16]. The initial duty of EMS in both centralized and decentralized mode is to ensure the microgrid of providing load-generation balance [17]. The EMS fail in matching the generation and load, whenever the total load is higher than the maximum capacity of DGs [18] and no other additional actions are taken. A solution for this drawback is to trade power with the utility or other micro-sources, however, this solution leads to an increment of pollution, costs and the need to solve a more complex unit commitment problem as a result of the additional units [19]. Various supporting systems such as DiGs, ESSs and DSM are employed to overcome the supply-demand mismatch of a microgrids. This paper provides the literature review of microgrid EMSs by classifying the existing articles into four categories as follows:


The remainder of the paper is organized as follows: in Section 2, the aforementioned categories are briefly introduced. A survey of uncertainty modeling techniques in EMSs is presented in Section 3. The mathematical formulation of objective functions along with constraints are presented in Section 4. The appropriate optimization techniques used by an EMS, microgrid test systems and obtained simulation results are reviewed in Sections 5–7, respectively. Finally, the conclusions are presented in Section 8.

#### **2. Classification of EMS Literature**

In this section, the necessity and utilization of the aforementioned categories are explained.

#### *2.1. Category 1: Non-Renewable Based EMS*

In case of failure or inaccessibility of ESSs, it is recommended to use non-renewable energy sources including diesel generators (DiGs), micro turbines (MTs), gas engines (GEs) or combustion turbines as a backup energy source in the microgrid. A DiG is made up of a combination of a diesel engine and an electric generator. The efficient selection of DiG depends on various factors such as load type, fuel cost, transportation cost, etc. [7].

#### *2.2. Category 2: ESS-Based EMS*

Microgrid EMSs face difficulties in the management of renewable energy sources such as wind and solar energy. This problem is due to the uncertain nature of the available energy which is caused by the difference between real-time and forecasted power production [20]. One of the solutions to tackle this problem is to utilize ESSs [21]. In most cases, ESSs maintain the power balance between generation and consumption by storing power during inexpensive or off-peak hours and discharging it during high-price or peak hours [22]. Diverse studies have been focused on the utilization of ESSs in a microgrid [23]. Utilizing ANN (Artificial neural network) for prediction of wind power generation, multi-objective energy management of microgrid considering uncertainties of wind generation in presence of ESS is studied in [24]. Karavas et al. [25] have studied the microgrid EMS considering a battery as ESS and solved the optimization problem based on distributed intelligence and MAS. Alavi et al. [26] solved the microgrid energy management problem considering a battery as a reserve energy source. In order to cover wind and solar power uncertainties, PEM has been used. Chen and Duan [27] have solved the optimal management problem of DGs with ESS in a microgrid by a MRCGA. Simultaneous capacity optimization of DGs and storage devices by considering the effects of weather condition and non-dispatchable power sources have been introduced in [28]. The EMS problem for multi-microgrid is proposed by Logenthiran et al. in [29] while the ESS being used as a reserve energy source.

#### *2.3. Category 3: DSM-Based EMS*

Another way to cope with unbalances of microgrid production and consumption is to utilize DSM programs [1]. All activities aiming to match the supply and demand by modifying time and/or shape of customers' demand profile is called DSM [30]. After the liberalization of electricity market, DSM is divided into following two categories [31]:


DG uncertainties and electricity price fluctuations cannot be managed and coped with for energy efficiency. As a result, special attention is given to DR programs in a microgrid and it is therefore crucial to study the impacts of various DR programs on the EMS of a microgrid.

DR programs are categorized into two groups namely: price-based and incentive-based DR programs [33]. A more detailed explanation of these categories can be found in [34].

Microgrid EMS in a deterministic case and with a multi-objective problem is studied in [35] wherein a price-quantity based DR package is utilized to manage the random nature of DGs. Falsafi et al. [36] have utilized both price-based and incentive-based DR programs to provide reserve, and to cover the wind power uncertainties in a multi-objective formulation. Nejhad et al. [37] have studied the EMS problem of microgrid as a stochastic problem in presence of DR programs. In addition, microgrid EMS taking into account effects of DR program is formulated as a MILP in [38]. A security-constrained EMS problem is considered in [39], in which frequency management is studied along with energy management and DR is utilized as a reserve energy source. The optimal operation of microgrid in the presence of electrical vehicles and responsive loads are discussed in [40] by considering wind and PV uncertainties.

#### *2.4. Category 4: Hybrid Systems Based EMS*

In this category, a combination of the abovementioned categories, two by two or all together, are surveyed to overcome microgrid EMS problems. Similar to the previous categories, this topic has also attracted lots of attention in the literature. In [41], the probabilistic coordination of DGs, ESS and DR are presented in a microgrid by performing a load reduction in the presence of security risks. Simultaneous implementation of ESS and DR is studied as a multi-objective problem by Marzband et al. [18]. Microgrid EMS in a system containing PV and wind turbines as DGs, DiG and small hydro generator as reserve power sources, DR and battery storage system is modeled by Zhao et al. [42] by using multi-agent systems (MASs). Talari et al. [43] have studied stochastic scheduling of microgrid components in presence of ESSs and DR programs. Pourmousavi et al. [44] have solved the management problem of an islanded microgrid which is composed of various DG units, storage and DR as a multi-timescale problem. The authors of [45] have utilized a pumped-storage unit and DR program in a new stochastic optimization framework to cover existed uncertainties of microgrid EMS problem. An overview of recent technologies in application of storage and DR operation of microgrids is presented in [33].

Solving the EMS problem of a microgrid is composed of various steps, including predicting uncertain parameters, modeling uncertainties, mathematical formulation of objective functions and constraints, choosing the optimization technique to solve the problem and selecting the understudying microgrid as the case study. Classification of the abovementioned steps along with sub-category of each step is shown in Figure 2. For instance, the existed uncertain parameters in microgrid EMS problem can be predicted in different time horizons such as short-term, mid-term and long-term time period. Prediction in short-term period can be performed using classical and intelligent techniques. Classical methods include ARIMA (auto-regressive integrated moving average), GHARCH (Generalized

auto-regressive conditional heteroskedasticity), DR and TF (Transfer function) while the intelligent ones are ANN, FNN (fuzzy neural network) and SVR (support vector regression). Theses explanations can be expanded for other steps, too.

#### **3. Prediction of Uncertain Parameters**

Lack of information which leads to a probability of a difference between real and forecasted values is defined as uncertainty [26]. The uncertain parameters in a microgrid EMS can be generally classified into the following two categories:


The prediction of uncertain parameters can be performed over various time horizons ranging from minutes to a couple of days as a short-term prediction, from several weeks to months as mid- term prediction and from multiple months to several years as long-term prediction [47]. Since the microgrid EMS problem is accomplished at hourly intervals, short-term forecasting methods are suitable to be used for this purpose. Classification of these short-term forecasting methods is represented in Figure 2 which is broken down into two categories namely classical and intelligent methods [37]. Classical techniques include well-known methods such as auto-regressive integrated moving average (ARIMA) [48], generalized auto-regressive conditional heteroskedasticity (GARCH) [49], dynamic regression (DR) and transfer function [50].

Intelligent methods are based on training historical data by mapping them on input-output of data-driven structures [37]. Artificial neural network (ANN) [51], fuzzy neural network (FNN) [52] and support vector regression (SVR) based on ARIMA can be classified in this category [53].

**Figure 2.** Overall EMS structure and underlying methods of each step.

#### **4. Uncertainty Modeling**

Uncertainty managing is one of the main difficulties decision makers have to deal with [46]. As a result, various methods have been employed in order to manage the uncertainty of the aforementioned parameters in Section 3 in microgrid energy management [46,54–57]. In this section, a review of all existing uncertainty modeling techniques used by an EMS are presented which are classified broadly into four categories as shown in Figure 2.

#### *4.1. Stochastic Methods*

Stochastic methods are used to approximate the PDF of random variables [58]. The first efforts in modeling uncertainties using stochastic methods was made by Dantzig [59]. These methods can be used, when the PDF of input parameters is known [60]. A brief introduction of three main stochastic technique are given below:

#### 4.1.1. Scenario Generation

Modeling uncertainties through scenario generation methods needs to know the PDF. In these methods, by dividing the PDF in multiple parts, each part is assumed as a scenario with an occurrence probability proportional to the PDF value of the preferred selected section [61]. The scenario generation method is used in [40] in order to cover uncertainties related to upper, lower and expected values of wind and PV systems. As reported by Alharbi and Raahemifar [41], uncertainties of wind, solar power and load are represented by discrete probability distribution sets as scenarios. In [39], load, wind and PV power random scenarios are firstly generated using MCS and roulette wheel mechanism. Then, to improve the computational speed, a scenario reduction algorithm is employed. Various scenarios for modeling wind and PV output powers are generated by Monte Carlo simulation by Talari et al. [43].

#### 4.1.2. MCS

In MCS-based methods, for each input parameter, a sample is generated using its PDF [62]. As MCS is a repetitive procedure, sample generation process is repeated for some iterations. Finally, histograms, statistic criteria or other approaches are used to analyze the outcome [63]. Caralis et al. [64] have used a Monte Carlo method to manage wind power uncertainties with the goal of earning profits from wind energy investment. Uncertainties related to solar radiation and load are investigated using MCS in [63] via optimization of the hybrid system.

#### 4.1.3. Point Estimate Method (PEM)

PEM is one of the methods used to estimate the value of uncertain parameters with high computational accuracy and is based on the concept of uncertain input parameters moments [65]. This method covers uncertainties by establishing a connection between input and output variables whose main steps are explained by Hong [66]. In [26], PEM is utilized to model uncertainties of wind and solar power. Peik-Herfeh et al. [67] employed PEM to model the uncertainties of market price and generation sources via a probabilistic price-based unit commitment problem.

#### *4.2. Fuzzy Method*

Based on fuzzy theory, a degree of membership can be attributed to each uncertain parameter using membership functions. Then, an α-cut method [68] can be used to determine membership functions of output variables according to membership functions of input parameters. An effort has been made in [69] to review the application of fuzzy methods in renewable energy sources. A fuzzy method has been used by Sourodi et al. [70] to study the effects of uncertain power output of DGs on power losses in distribution networks.

#### *4.3. Robust Optimization Method*

The robust optimization method [71] was introduced to explain parameter uncertainties using uncertain boundaries. Robust optimization is suitable when there is a lack of information about PDF of parameters. This technique is used to cover load uncertainties by Alavi et al. in [26]. Robust optimization is employed for wind and load uncertainty in Ref. [72].

#### *4.4. Information Gap Decision Theory (IGDT)*

As a decision-making approach, IGDT makes minor assumptions on the structure of uncertainty [73]. This technique makes robust decisions against uncertainty of the input parameters. The optimal bidding strategy in day-ahead electricity market considering uncertainties of market price is performed by implementing IGDT in [74]. Nojavan et al. [75] proposed a method based on IGDT to evaluate procurement strategy of large consumers.

#### **5. Mathematical Formulation of EMS**

Microgrid energy management is an optimization problem which aims to properly schedule short-term operation of DGs, ESSs and controllable loads with respect to various objective functions and constraints [76]. In this section, a literature review of existing objective functions as well as constraints considered by an EMS have been elaborated. Furthermore, Figure 2 represents a classification of objective functions utilized by the EMSs along with their constraints.

#### *5.1. Objective Functions*

The EMS can manage a microgrid by solving various objective functions. Objective functions may include capital or operational costs of the microgrid. Costs related to fuel, maintenance, start-up and shut-down, degradation as well as procurement from the utility in case of power deficiency, are considered as operational costs [77]. Table 2 provides a collection of utilized EMS objective functions in literature. In this table, objective functions are reviewed from being single objective and multi-objective perspectives, too.


**Table 2.** Survey through collection of EMS objective functions.


#### **Table 2.** *Cont*.

#### *5.2. Constraints*

Different constraints can have an effect on energy management of a microgrid. For instance, maximum and minimum limits of power generation units must be satisfied to secure their safe and

economic performance [82]. The balance between generation and consumption is another necessity of the system. The charge and discharge rates of ESS are also constrained. Violation of these constraints can lead to damaging effects on lifetime and efficiency of the ESS. Technical constraints of a microgrid include voltage at buses, feeder currents, frequency security aspects, start-up and shut-down reserve constraints, as well as ramping limits. In some of the studies which additionally consider responsive loads, constraints related to DR program must be satisfied. Table 3 provides a summary of the considered constraints used in the formulation of microgrid EMS.


**Table 3.** Survey through collection of EMS constraints.

#### **6. Optimization Techniques of a Microgrid Performed by an EMS**

The literature review reveals that various techniques have been utilized by researchers in order to solve the aforementioned optimization problems [17]. In this paper, these methods have been classified as follows:

#### *6.1. Heuristic Approaches*

Heuristic optimization techniques use exploratory approaches to solve the optimization problems while are unable to assure optimality of the obtained results [6]. A deterministic energy management problem is solved via multi-period gravitational search algorithm (MGSA) by Marzband and Ghadimi in [18]. In [35], the EMS optimization problem has been considered as a continuous multi-objective problem and is solved using multi-objective particle swarm optimization (MPSO) algorithm while it is solved by single-objective PSO in [26,40]. Motevasel and Seifi [24] have employed multi-objective Multi-objective Bacterial Foraging Optimization (MBFO) to solve complex and large-scale EMS problem. To deal with non-smooth cost functions, authors in [27,79] employed MRCGA. Implementation of multi-period ABC optimization algorithm to obtain an economic schedule of ESS and controllable loads of an islanded microgrid is studied in [83].

#### *6.2. MAS*

MAS is a smart system which is composed of multiple intelligent agents that are connected to each other within an environment [15]. The goal of MAS is to solve the optimization problems which are too complicated for a single agent. By assuming each member of microgrid as an agent, MAS has been used in [42] to find the optimal solution by the EMS problem. Optimal results of decentralized EMS problem have been studied using MAS approach by Karavas et al. [25]. A three-stage algorithm based on MAS is presented by Logenthiran et al. [29] to model the EMS problem in a multi-microgrid environment in which first stage schedules each microgrid to satisfy its load. The second and third stages determine microgrid bids and export power bids, respectively. Authors of Ref. [84] have utilized trust and reputation models to interconnect MAS effectively. In [85], a novel agent-based micro-storage management technique has been presented that allows all (individually-owned) storage devices in the system to converge a profitable and efficient behavior. A class of MAS technique is Distributed Constraint Optimization Problems (DCOP) in which a group of agents try to minimize the total cost related to a set of constraints. In a DCOP formulation the decision problem is translated to a set of constraints [86]. Long-term scheduling of microgrid is solved as an optimization problem using distributed constraint optimization (DCOP) by authors of [87]. Application of DCOP to solve large problems has difficulty as each agent can handle a single variable. So, a new multi-variable agent (MVA) DCOP decomposition method has been presented by athors of [88].

#### *6.3. Mathematical Methods*

Different kinds of software exist in order to solve microgrid EMS problems. Some of these solvers are mentioned below:

#### 6.3.1. CPLEX Solver

CPLEX is a solver of GAMS package which can solve integer and linear problems [6]. In [39], the EMS problem is modeled as a MILP considering uncertainties and scenario generation process, and has been solved by the CPLEX solver of the GAMS environment. The multi-scenario MILP model has been utilized to mathematically formulate the EMS problem in [41], and is solved by CPLEX solver. CPLEX solver has been utilized to solve MIP scheduling problem of microgrid by Talari et al. [43]. A large-scale MILP problem in [36] is solved also by CPLEX solver. Another EMS problem is modeled as a MILP in GAMS software and solved via CPLEX solver [79].

#### 6.3.2. SNOPT Solver

SNOPT is another software package of GAMS which is capable of solving nonlinear optimization problems [45]. The EMS mathematical formulation in [28] is a nonlinear programming problem, and has been solved by GAMS SNOPT solver.

#### 6.3.3. Gurobi Optimizer

Gurobi optimizer is used to solve MILP problems [89]. In [78], the EMS problem is formulated as a MILP problem and has been solved using Gurobi optimizer. Tenfen and Finardi [38] have employed the same solver for computational purposes.

#### **7. Microgrid Test Systems**

For the sake of evaluating the performance of the aforementioned EMS algorithms, they have been applied on the majority of microgrids. There are some efforts to sum up the studied microgrid test systems in the area of EMS in the literature [20]. Reference [90] has proposed the classification of microgrids based on their topologies. A review of existing real-world microgrid test systems worldwide is presented in [20,91]. Evaluation of microgrids which are applied in real-life as well as laboratory test systems are surveyed in [92]. Nonetheless, in this section we will solely focus on the test systems which are exploited in EMS studies. Table 4 has classified these test systems from various panoramas including single-microgrid, multi-microgrid, islanded and grid-connected systems.

**Table**

**4.**

Review

of

test

**Table4.***Cont*.

**Table**

**4.**

*Cont*.

18

*Energies* **2019**, *12*, 2156

20

*Energies* **2019**, *12*, 2156

**4.***Cont*. *Energies* **2019**, *12*, 2156

22

*Energies* **2019**, *12*, 2156

**Table 4.** *Cont*.

#### **8. EMS Simulation Scenarios**

The aforementioned EMS approaches are implemented on typical microgrid test systems in order to optimize the corresponding objective functions while considering various constraints in a 24-h time interval. Table 5 briefly summarizes the obtained simulation scenarios of reviewed papers to help researchers in directing their research efficiently. According to the first row of Table 5, three scenarios are considered in the results and simulation section of [24]. In the first scenario, all components of microgrid are operating within their limits and limited trade is available with the main grid. However, during the second scenario, WT operates at its maximum capacity while the third scenario studies effect of unlimited power exchange between microgrid and the upper utility. Same as the first row, remain rows of Table 5 explain the number and explanation corresponding to the simulation scenario of each reference.


**Table 5.** Scenarios of EMS optimization problem results.

#### **9. Conclusions**

This review paper has summarized energy management approaches for microgrids from a new perspective and classified it into four categories, namely non-renewable, ESS, DSM and hybrid- based EMS. This paper also provides a compendium on the modeling uncertainties associated with microgrid EMS, objective functions and constraints of microgrid formulation as well as many tools and techniques for solving this optimization problem. It is worth mentioning that considering various objective functions along with different technical and economic constraints has a large effect on the obtained EMS results. A brief review of test systems including single-microgrid and multi-microgrids is also described in this paper. A number of EMS simulation scenarios has been included while each one is an implemented simulation state. Furthermore, this work will help in identifying poorly researched areas in microgrid EMS for further investigation.

**Author Contributions:** H.S. supervised writing the paper, E.S wrote the paper, M.M. and P.S. checked and proofread the manuscript.

**Funding:** "This research received no external funding"

**Conflicts of Interest:** "The authors declare no conflict of interest."

#### **Abbreviations**


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **An Augmented Prony Method for Power System Oscillation Analysis Using Synchrophasor Data**

#### **Moossa Khodadadi Arpanahi 1, Meysam Kordi 2, Roozbeh Torkzadeh 3, Hassan Haes Alhelou 1,4 and Pierluigi Siano 5,\***


Received: 8 March 2019; Accepted: 29 March 2019; Published: 2 April 2019

**Abstract:** Intrinsic mode functions (IMFs) provide an intuitive representation of the oscillatory modes and are mainly calculated using Hilbert–Huang transform (HHT) methods. Those methods, however, suffer from the end effects, mode-mixing and Gibbs phenomena since they use an iterative procedure. This paper proposes an augmented Prony method for power system oscillation analysis using synchrophasor data obtained from a wide-area measurement system (WAMS). In the proposed method, in addition to the estimation of the modal information, IMFs are extracted using a new explicit mathematical formulation. Further, an indicator based on an energy and phase relationship of IMFs is proposed, which allows system operators to recognize the most effective generators/actuators on specific modes. The method is employed as an online oscillation-monitoring framework providing inputs for the so-called wide-area damping control (WADC) module. The efficacy of the proposed method is validated using three test cases, in which the IMFs calculation is simpler and more accurate if compared with other methods.

**Keywords:** power system oscillations; stability; Prony method; intrinsic mode functions (IMFs); wide-area measurement system (WAMS); phasor measurement units (PMUs)

#### **1. Introduction**

Power system oscillation analysis has received considerable attention over the past decades. Features such as frequency and damping of the power system oscillations provide critical information for system operators [1]. Analysis of such oscillations are conducted using: (i) model-based and (ii) measurement-based approaches. The former approach is based on the linearized differential equation of the system around a specific operating point [2,3], as accurate modeling of power systems for various operating conditions is complicated and these approaches may not be useful for practical applications. On the other hand, with developing phasor measurement units (PMUs) and wide-area measurement system (WAMS), the latter approach became more popular. Recently, several studies have attempted to extract oscillation parameters using synchrophasor data. Some of the well-known methods used for oscillation detection using measurement-based approaches include: Kalman Filter (KF) [4] estimation of signal parameters by rotational invariance technique

(ESPRIT) [5], Hilbert–Huang transform (HHT) [6–9], Yule–Walker algorithm [10], wavelet transform (WT) [11], recursive adaptive stochastic subspace identification (RASSI) [12], maximum likelihood method [13], and Prony analysis method [14–16]. Among them, the Prony method is one of the most commonly used measurement-based methods. The Prony method can directly estimate modal components of a measured signal. Recently, extensions of the Prony method have been developed that can calculate oscillatory modes of multiple signals in a central [14] or distributed [15] manner. Furthermore, a forward and backward extended Prony (FBEP) method is proposed in order [16] to reduce the sensitivity of the Prony method to noise.

A practical online oscillation monitoring framework should be able to determine the power system oscillations fast and accurately. Further, it should be capable of identifying contributions of different system buses in each oscillatory mode. Among the above methods, only the HHT method can determine the contribution of each state variable in various oscillatory modes. In general, each state variable of the system can be expressed as the sum of exponential sinusoids. According to the HHT method, each of the exponential sinusoids is called intrinsic mode functions (IMF). Using IMFs, participation of each system bus in an oscillatory mode can be determined by means of an indicator called node contribution factor (NCF). In the HHT method, IMFs can be extracted from the synchrophasor data using empirical mode decomposition (EMD), which is an iterative procedure. However, EMD process has some inherent shortcomings such as the end effects (EEs), mode-mixing and Gibbs phenomena. Some of those drawbacks have been addressed in improved versions of the HHT method in [6–9].

In this paper, we propose an augmented Prony based framework for real-time power system oscillation monitoring and analysis using synchrophasor data. The proposed framework estimates modal features (i.e., its damping and frequency) of a measured signal, while it is capable of determining contributions of different system buses in oscillatory modes by calculating IMFs and NCFs. Note that in the proposed framework, IMFs and thus NCFs are extracted using an explicit mathematical formulation based on an augmented Prony method. To the best of the authors' knowledge, the Prony method has never been used before for the calculation of IMFs and thus NCFs. Unlike the HHT method that relies on heuristic approaches such as EMD for the calculation of IMFs, in the proposed approach they are calculated using an explicit mathematical formulation, which is one of the novelties of the proposed method. The calculation of IMFs in the proposed framework is more accurate, while it requires a lower computational burden if compared to the existing HHT methodologies. The proposed online oscillation-monitoring framework can be used as an input for the so-called wide-area damping control (WADC) module, which is responsible for controlling the system oscillations. Note that obtaining NCFs associated with oscillatory modes in the system can assist power system operators to take better remedial actions enhancing system stability and security. The effectiveness of the proposed framework is validated using the Kundur's two-area test system and a practical power system. The main innovative contributions in this paper are:


The rest of the paper is organized as follows: the motivations behind the work and the formulations IMFs and NCFs calculation in the proposed framework are presented in Section 2. In Section 3, performance of the proposed framework is analyzed using the Kundur's two-area power system and a practical power system. Finally, Section 4 concludes the paper.

#### **2. Proposed Framework for WAMS-Based Oscillation Detection**

#### *2.1. Research Motivation*

Developing a real-time power system oscillation monitoring framework is an important part of small signal stability studies using WAMS technologies. An overview of Smart Grid components and technologies is shown in Figure 1. Observe that the smart grid consists of different state-of-the-art technologies enabling effective network management schemes in the presence of renewable energy systems, flexible loads, and energy storage systems. This new concept ultimately offers the efficient and reliable operation of future power systems via a two-way exchange of data and energy among different sectors in the electric power supply chainmn including generation, transmission, distribution, and consumption. Implementation of PMUs and WAMS to achieve a fully observable grid was defined as the first phase of a Smart Grid project [17,18].

**Figure 1.** General outline of a Smart Grid.

As shown in Figure 2, WAMS architecture includes five layers. The component layer consists of physical devices. It is responsible for providing data, functions, and communication services for the upper layers. The communication layer applies several methods and protocols to obtain data among various elements in the network. The information layer includes the data model and the communication rules that are used to exchange information. The function layer contains the logical functions or various WAMS applications, independent of the physical architecture; however, it is reliant on the business models and the organization's policies. Finally, the business layer defines the business models and the policies [17,18]

**Figure 2.** Wide-area measurement system (WAMS) architecture.

One of the main online applications associated with the proposed WAMS function layer is power system oscillation detection using synchrophasor data, which is the main scope of this paper. The oscillation detection algorithm is detailed in Sections 2.2–2.5.

#### *2.2. Motivation behind Methodology Selection*

In many Smart Grid projects, some well-known methods used for power system oscillation detection were analyzed and compared. This has been carried out to identify the best method in terms of: (1) accuracy in the estimation of the oscillatory mode features, (2) simplicity in the implementation, (3) ability to provide visual intuitiveness, (4) skill to determine the most associated state variable to each oscillatory mode, (5) robustness against the noise, and (6) self-verification ability. A survey of the literature showed that the extended Prony methods can be a reasonable choice for the proposed framework. Note that the basic Prony method has some weakness such as its sensitivity to noise and selection of the proper order of the system [14–16]. An overestimated order may lead to the generation of extraneous modes by the Prony method. These drawbacks were addressed in some extensions of the Prony method. For example, the study in [15] proposed a simple criteria for selecting the order of the system. The study in [16] suggested that the Prony method can be performed in forward and backward manners, where the final modes are selected from the collection of those extracted in such a process. This process enabled the Prony method to identify extraneous modes, while allowed elimination of the noise effects. Prony is inherently a self-verifiable method because it estimates the modes of a given signal and then reconstructs the main signal from the estimated mode. Comparing the main and the reconstructed signal is a convenient way to verify its performance.

The above features make the Prony method an effective one for practical applications; however, it does not provide visual intuitiveness on contribution various system modes in each oscillatory modes. This is important because power system operators are interested in visual information instead of numerical data. In this regard, the HHT method has advantages over the Prony method, in which the so-called IMFs plots represent the behavior of each mode intuitively, and therefore provide situational awareness on oscillations of the power system. Moreover, the HHT method can

determine the contribution of a state variable in a mode based on energy and phase relationship of IMFs. Nevertheless, as stated in the introduction, the HHT method suffers from the EEs, mode-mixing and Gibbs phenomena [6,9]. In this paper, an improved Prony method is proposed that is able to calculate IMFs and thus NCFs for a given signal.

Note that HHT has a much wider application than the Prony method and the authors don't mean the proposed method is generally better than HHT, especially in non-linear situations. We compare our method with HHT because it is the only method that can calculate IMFs and NCFs. To achieve the IMFs and NCFs, HHT should overcome some challenges such as end-effects, mode mixing, and Gibbs phenomena. These challenges are solved in different literature. For example, there are several ways to overcome end-effect issue [17]. In this paper, a direct method is proposed to obtain the IMFs and NCFs, in which the HHT challenges are avoided. Indeed, we are facing the alternative to obtain IMFs and NCFs as depicted in the following figure:


Therefore, both methods (HHT and proposed method) can handle the IMFs and NCFs with different approaches and can be used based on the situations. Figure 3, briefly shows the above discussions.

**Figure 3.** General outline of (**a**) Hilbert–Huang transform (HHT) method; (**b**) proposed method.

#### *2.3. Oscillation Analysis Framework*

Assuming a power system of order *n*, the zero-input time response is as follows [14]:

$$\underline{y}(t) = \sum\_{i=1}^{n} \gamma\_i e^{\lambda\_i t} \tag{1}$$

where *λ<sup>i</sup>* are called eigenvalues or modes of the system. The Prony method represents a given signal as a sum of exponential functions and then estimates the unknown parameters, i.e., *γ<sup>i</sup>* and *λ<sup>i</sup>* so that the real and the estimated output have a minimum difference. The estimated signal in the discrete time domain, and can be expressed as follows [14]:

$$\underline{\hat{y}}(k) = \sum\_{i=1}^{n} \beta\_i z\_i^{\;k} \quad k = 1, 2, \dots, N \tag{2}$$

where, *zi* are the eigenvalues of the system in the discrete time domain or *z*-domain, *β<sup>i</sup>* is the residue corresponding to *zi* and *N* is the number of the samples. Equation (2) can be described in a matrix form as [14]:

$$
\begin{bmatrix} z\_1^{1\_s} & z\_2^{1\_s} & \cdots & z\_n^{1\_s} \\ z\_1^{2\_s} & z\_2^{2\_s} & \cdots & z\_n^{2\_s} \\ \vdots\_1^{N\_s} & \vdots\_2^{N\_s} & \vdots & \vdots\_N^{N\_s} \\ z\_1^{N\_1} & z\_2^{N\_1} & \cdots & z\_n^{N\_s} \end{bmatrix} \begin{bmatrix} \beta\_1 \\ \beta\_2 \\ \vdots \\ \beta\_n \end{bmatrix} = \begin{bmatrix} y(1) \\ y(2) \\ \vdots \\ y(N) \end{bmatrix} \tag{3}
$$

According to the Prony method, for a discrete time system with an output as (2), the so called characteristic equation is as follows [14]:

$$z^n - a\_1 z^{n-1} - a\_2 z^{n-2} - \dots - a\_{n-1} z - a\_n = 0 \tag{4}$$

The coefficients of the above equation are obtained by solving the predictive equation as follows [14]:

$$
\begin{bmatrix} y(n) & \dots & y(1) \\ y(n+1) & \dots & y(2) \\ \vdots & \vdots & \vdots \\ y(N-1) & \dots & y(N-n) \end{bmatrix} \begin{bmatrix} a\_1 \\ a\_2 \\ \vdots \\ a\_N \end{bmatrix} = \begin{bmatrix} y(n+1) \\ y(n+2) \\ \vdots \\ y(N) \end{bmatrix} \tag{5}
$$

which can be given in a matrix form as:

$$D\underline{a} = \underline{d} \tag{6}$$

where, *D* is the Prony Matrix. For WAMS-monitored power systems, both *D* and *d* contain synchrophasor data (e.g., voltages, currents, or powers). Based on the above formulations, the algorithm of the Prony method can be summarized as follows [14–16]:


Based on the above algorithm, in Step 3, the *z*-domain eigenvalues are calculated. The eigenvalues of the main signal can be obtained as follows [14]:

$$z\_i = \epsilon^{\lambda\_i T} \Rightarrow \lambda\_i = \frac{\log(z\_i)}{T} \quad , i = 1, 2, \dots, n \tag{7}$$

where, *T* is the sampling rate of the main signal. In many cases, instead of complex values of eigenvalues, two well-known parameters of the modes, i.e., damping ratio *ξ* and damping frequency *fd* are represented, which are related to eigenvalues as follows [1]:

$$\lambda\_i = \sigma\_i \pm j\omega\_{d\_i} \text{ , } \xi\_i = \frac{\sigma\_i}{\left|\lambda\_i\right|} \text{ , } f\_{d\_i} = \frac{\omega\_{d\_i}}{2\pi} \quad , i = 1, 2, \dots, n \tag{8}$$

#### *2.4. Calculation of IMFs and NCFs Based on Augmented Prony Method*

As previously stated, an online oscillation monitoring framework should be able to determine the contribution of various system buses in each oscillatory mode. Here, as one of the novel contributions of this work, it is demonstrated that the Prony method can be used for such a purpose. Initially, the matrix *U* is defined as follows:

$$MI = \begin{bmatrix} z\_1^{1\_s} & z\_2^{1\_s} & \cdots & z\_n^{1\_s} \\ z\_1^{2\_s} & z\_2^{2\_s} & \cdots & z\_n^{2\_s} \\ \vdots\_1^{N\_s} & \vdots\_2^{N\_s} & \vdots & \vdots\_N^{N\_s} \\ z\_1^{N\_1} & z\_2^{N\_1} & \cdots & z\_n^{N\_s} \end{bmatrix} \tag{9}$$

where, *ui* is defined as column *i* of matrix *U*, thus, for a complex eigenvalue, there is another eigenvalue which is the conjugate of *zi*. Without the loss of generality, it is assumed that *zi*+<sup>1</sup> and *βi*+<sup>1</sup> are the conjugates of *zi* and *βi*, respectively (*zi*+<sup>1</sup> = *zi* <sup>∗</sup>; *βi*+<sup>1</sup> = *β<sup>i</sup>* ∗), thus:

$$
\beta\_i \ u\_i + \beta\_{i+1} \ u\_{i+1} = \beta\_i \ u\_i + \beta\_i \ ^\*\mu\_i \ ^\* \quad , i = 1, 2, \dots, n \tag{10}
$$

Substituting *ui* from (9) into (10) yields:

$$\beta\_i \, u\_i + \beta\_i^+ \, u\_i^+ = \beta\_i \begin{bmatrix} z\_i^1 \\ z\_i^2 \\ \vdots \\ z\_i^k \\ \vdots \\ z\_i^k \end{bmatrix} + \beta\_i^+ \begin{bmatrix} (z\_i^\*)^1 \\ (z\_i^\*)^2 \\ \vdots \\ (z\_i^\*)^k \\ \vdots \\ (z\_i^\*)^k \end{bmatrix} \quad , i = 1, 2, \dots, n \tag{11}$$

where, *β<sup>i</sup>* and *zi* are defined as follows:

$$z\_i = e^{\lambda\_i T} = e^{(\sigma\_i + j\omega\_{di})T}, \beta\_i = \frac{A\_i}{2} e^{j\varphi\_i} \qquad , i = 1, 2, \ldots, n \tag{12}$$

Thus, (11) can be rewritten as:

$$\begin{aligned} \left[\beta\_i \,\mu\_i + \beta\_i^\*\,\mu\_i^\* = \begin{bmatrix} A\_i e^{\delta\_i T} \cos[\omega\_{d\_i} T + q\_i] \\ A\_i e^{\delta\_i (2T)} \cos[\omega\_{d\_i} (2T) + q\_i] \\ \vdots \\ A\_i e^{\delta\_i (kT)} \cos[\omega\_{d\_i} (kT) + q\_i] \\ \vdots \\ A\_i e^{\delta\_i (NT)} \cos[\omega\_{d\_i} (NT) + q\_i] \end{bmatrix} \right] , i = 1, 2, \dots, n \end{aligned} \tag{13}$$

According to the definition of IMF (each of the exponential sinusoids contained in the main signal), the *k*th row of the vector in (13) shows the value of IMF of mode *i* at the time *t* = *kT*. Thus, IMFs can be calculated as in (14):

$$IMF\_i(t = kT) = \beta\_i \ u\_i(k) + \beta\_i {}^\* u\_i{}^\*(k) \quad , i = 1, 2, \dots, n \tag{14}$$

The above expression is valid for complex modes. For any real mode, its corresponding column in matrix *U* results in values of *IMFi*. Therefore, a generalized form for the calculation of IMFs can be expressed as:

$$IMF\_i(t = kT) = \beta\_i u\_i(k) + c\_i \beta\_i^\* u\_i^\*(k), \quad i = 1, 2, \dots, n \tag{15}$$

where, *ci* is a binary variable that shows either mode *i* is complex (*ci* = 1) or not (*ci* = 0). Note that it is assumed in (15) that columns of the matrix *U* are sorted based on their real part.

Having intrinsic function of each mode, it is possible to determine the most associated state variable corresponding to each oscillatory mode, using energy and phase relationship of IMFs. In other words, the generator or actuator that is the most associated ones with a specific mode can be identified. As mentioned in the Introduction, the participation of each system bus in an oscillatory mode can be determined using node contribution factor (NCF). This index is defined as follows:

$$NCF\_{s,i} = E\_{s,i} \cos(\varphi\_{s,i}) \tag{16}$$

where, *NCFs*.*<sup>i</sup>* represents the contribution of signal *s* in mode *i*. *Es*,*<sup>i</sup>* is defined as the energy of IMFi (related to mode *i*) obtained by signal *s*, which is defined as follows:

$$E\_{s,l} = \sum\_{k=1}^{N} |IMF\_{s,l}|^2 (t = kT) \tag{17}$$

In (16), *ϕs*,*<sup>i</sup>* is the average relative phase of *IMFs*,*<sup>i</sup>* and the first IMF (*IMFs*,1 that is associated with the dominant mode) in different times:

$$\varphi\_{s,i} = \left(\frac{1}{N}\right) . \sum\_{t=1}^{N} \theta\_{s,i}(t) - \theta\_{s,1}(t) \tag{18}$$

In order to calculate *ϕs*,*<sup>i</sup>* a phase angle is assigned to each data point of *IMFs*,*i*. For all IMFs extracted from a measured signal of the network, the initial phase angle of the first IMF is defined as the reference phase (*θs*,1(*t*) = 0). For a specific IMF, the phase angle of maximum, minimum, positive zero-crossing, and the negative zero-crossing points are set as π/2, 3π/2, 0 (or 2π) and π, respectively. According to the proposed formulation in (15)–(18), it is possible to directly calculate IMFs of the various signals of the system and therefore obtain their contributions in the different oscillatory modes. It is important to emphasize that so far, the Prony method has not been augmented for calculation of IMFs and thus NCFs. Our proposed formulation in this subsection bypasses the problems of nonlinear and the iterative method of HHT method for IMF calculation, such as the mode mixing of EMD, end effect and Gibbs phenomena [6–9].

Note that the understudy signal may contain noise that can be easily eliminated before performing the proposed method. It is worth mentioning that the mode mixing phenomena is not related to the noise. In the HHT method the IMFs are obtained using an iterative procedure with two loops, the outer loop repeats with the number of modes and for each mode an iterative procedure called sifting is performed to achieve the IMF associated with that mode. When the frequency of two modes of the system is close to each other, mode mixing may have occurred. In fact it is related to frequency heterodyne. Therefore, noise that is a high frequency signal may not mix with power system modes, which are low-frequency modes. In our proposed method, all sifting or decomposition procedures are not performed; therefore it is not prone to mode mixing. In the proposed method, each mode has a unique IMF. Thus, with an explicit mathematical formulation, each mode and its IMF are specifically extracted.

It is worth mentioning that power system oscillation detection and analysis is one of the topics of small signal stability assessment. In small signal stability assessment and analysis, it is generally accepted that the variation of the state variables of the system around the operating point are small, thus the differential equations of the system can reasonably be linearized at the operating point. The zero-input time response of such equations is seen in Equation (1). Therefore, despite the HHT method can handle signals without assumption of structure in Equation (1), in small signal stability analysis, which is the scope of this paper, we deal with the exponential sinusoids signals called intrinsic mode functions (IMFs). It is worth mentioning that the method is performed on many different signals

with different parameters (amplitude, damping, frequency, and phase angle), and in all cases, the performance of the method in estimation of the mode parameters and extraction of the IMFs, are approved. However, an example of such signal is presented in Section 3.1.

#### *2.5. Approximation System Order*

In practical situations, the number of modes incorporated in the signal is unknown. Several methods are proposed in the literature to approximate the effective order of the system [14–16]. In this paper, a new IMF-based index is introduced to determine the effective system order. The smallest index m that satisfies the inequality in (18) is the effective order of the system.

$$\frac{\sum\_{i=1}^{m} \, \, \, IMF\_{s,i} \, ^2}{\sum\_{i=1}^{n} \, \, \, IMF\_{s,i} ^2} \gtrsim \, \tau$$

where, the threshold value *τ* is close to 1 (e.g., 0.999) and *n* is the initial value of the system order.

#### **3. Simulation Results and Discussion**

#### *3.1. Test Case 1: Hypothetical Signal*

In this test case, a hypothetical signal give in (20) with four pairs of complex conjugate eigenvalues is considered. Parameters of this signal are shown in Table 1. The main signal (i.e., *y*(*t*)) and its components (i.e., *y*1(*t*) to *y*4(*t*)) are shown in Figure 4. It is assumed that signal *y*(*t*) is measured by PMU and resampled at 10 Hz used as an input of our proposed method. The main advantage of this test case is that the real modes and their corresponding IMFs are predefined, thus it is suitable for the validation of the proposed method. The jth mode associated with *y*(*t*) can be calculated according to (8). Also, *yj*(*t*) is by definition the jth IMF. The main signal is compared with the reconstructed signal by the Prony method in Figure 4. As shown in Figure 5, the reconstructed signal completely matches the real signal *y*(*t*), as expected. As stated in Section 2.3, Figure 5 illustrates the self-verification property of the Prony method.

$$y(t) = \sum\_{j=1}^{4} y\_j(t) = \sum\_{j=1}^{4} A\_j e^{\sigma\_j T} \cos \left[ 2\pi f\_{dj} t + \varphi\_j \right] \tag{20}$$


**Table 1.** Parameters of the Hypothetical Signal of Test Case 1.

**Figure 4.** Hypothetical signal of Test Case 1 and its components.

**Figure 5.** Comparison of real and reconstructed signals in Test Case 1.

The modes of the signal *y*(*t*) obtained by the Prony method are compared with real modes in Table 2. Observe that the Prony method is capable to extract modes accurately. The IMFs of the hypothetical signal are shown in Figure 6. Observe that extracted IMFs from our proposed mathematical method in Section 2.4 are exactly the same as real IMFs (*yj*(*t*); *j* = 1, ... , 4). Thus, our mathematical formulation computes IMFs of a given signal with high accuracy, without a need for heuristic algorithms such as those employed in the HHT method. Also, as illustrated in Figure 6, the proposed method represents IMFs of a given signal, in order of increasing damping for them. In other words, IMF of a more dominant mode or equivalently less damped ones will be represented first in our methodology. Therefore, in the Test Case 1, IMF1 to IMF4 are corresponding to *y*1, *y*4, *y*3, and *y*<sup>2</sup> in the hypothetical signal of (20), respectively, shown in Figure 4. 

**Figure 6.** Comparison of real and estimated intrinsic mode functions (IMFs) of signal of Test Case 1. (**a**–**d**) IMF1 to IMF4, respectively.


**Table 2.** Calculated Modes of Test Case 1.

#### *3.2. Test Case 2: Kundur's Two-Area Power System*

The second test case is carried out using Kundur's two-area test system [1], which is shown in Figure 7. These two areas are connected through weak transmission lines. The output active power of the generators denoted after the occurrence of a three-phase short-circuit fault is shown in Figure 8. The fault occurred at t = 1 s on Bus 8 and was cleared at t = 1.05 s. Assume that PG1 to PG4 are the output active power of the generators G1 to G4, respectively. Each of the four signals is analyzing using the proposed method in Sections 2.3 and 2.4 to extract modes of the system, their IMFs and consequently NCFs. The extracted modes based on each signal are depicted in Table 3. It is evident from the table that the extracted modes from the different signal (i.e., PG1 to PG4) result in the same modes. This is one way to validate the accuracy of the modes obtained by the Prony method. The main and reconstructed values of PG1 by the proposed method are compared in Figure 9, where again the accuracy of our framework is demonstrated. The extracted IMFs from generator G1 output power are demonstrated in Figure 10. Observe that IMF1 is associated with an unstable mode, while the other three modes are under-damped. While analyzing numerical results of Table 3 is important, our framework provides visual information using IMF plots (using the proposed formulation in Section 2.4) such as the one shown in Figure 10, which provides situational awareness for oscillations during the fault.

**Figure 7.** Kundur's two-area test system [1].

**Figure 8.** Power output of all generators related to Test Case 2.

**Table 3.** Comparison of Extracted Modes based on Active Power of Different Generators in Test Case 2; *ξ*(−), *fd*(*Hz*).


**Figure 9.** Comparison of main and reconstructed signals of PG1 in Test Case 2.

**Figure 10.** Extracted IMFs of PG1 in Test Case 2. (**a**–**d**) IMF1 to IMF4, respectively.

ሺ

Finally, based on energy and phase angle of IMFs, contribution factor of each generator output power in each of the modes can be determined, as demonstrated in Table 4.


**Table 4.** Node contribution factor (NCF) matrix in Case 2.

The matrix in Table 4 implies that the most associated generation in oscillatory modes 1, 2, and 4 is generator G4, while mode 3 is mainly affected by generator G1.

#### *3.3. Test Case 3: Interconnected Practical Power System*

In this section, the proposed method of Sections 2.3–2.5 has been used to analyze measured signals in a practical transmission power system. The output powers of the five generators resampled at 10 Hz are used to demonstrate performance of the proposed method. All presented results in this section are generated by the interconnected power system WAMS-based oscillation detection application utilizing our proposed method of Sections 2.3–2.5. A general schematic of this application is demonstrated in Figure 11. The signals related to five generators (here they are named G1 to G5) associated with the area that are affected by the fault are depicted in Figure 11. As shown in Figure 11, the information of oscillatory modes are provided to the system operator in a format of a table. Columns 1 to 4 in the table represent real and imaginary part of each oscillator mode, and their damping ratio and frequency, respectively. Moreover, the fifth column demonstrates stability status of each mode. For example, as shown in Figure 11, the first mode is an unstable mode shown with red color. The next columns show type of modes based on their frequency, i.e., local or inter-area. According to Figure 11, mode 2 is an inter-area mode with frequency of 0.4 Hz, while Modes 3 and 4 are local modes with frequency of 0.94 and 4.4 Hz due to oscillation of Generator G1 output power.

**Figure 11.** A schematic of the power system oscillation detection application.

As an example, the comparison of the original and reconstructed generation G5 signal based on the Prony-based method is shown in Figure 12, which demonstrates an acceptable matching between two signals. The accuracy of the modes and consequently IMFs and NCFs can be concluded from this matching. Having IMFs associated with all five generation signals, a NCF matrix with five rows and four columns is calculated by our framework (see Section 2.4), as shown in Table 5. This table shows that, the most associated signal to modes 1, 2, and 4 is G5 and mode 3 is mainly affected by G1. Thus, for controlling and stabilizing the system against the unstable mode (mode 1), also increasing the damping of the poorly damped mode 4, G5 should be controlled. Also, observe that generation G3 is the most suitable actuator for controlling mode 3. As seen from Table 5, the energy of the mode 2 is less than mode 1, while the energy of modes 3 and 4 relative to mode 2 is insignificant. Note that as the energy of the other modes is negligible, our oscillation detection framework discards them from the representation of their parameters as well as IMFs and NCFs. For better representing the contribution of each generator in each mode, the values of the NCF matrix are normalized and plotted in Figure 13. Other important plots in the power system oscillation detection application are related to IMFs of the measured signals, which are not shown here for brevity. The plot is similar to those shown in Figures 6 and 10 in Test Case 1 and 2, respectively. These provide an increase intuitiveness in the proposed oscillation detection framework, beside NCF plots. It was observed in IMF plots that the main or dominant mode (mode 1) has a growing IMF as it is unstable, while the other modes (mode 2 to 5) have exponential sinusoidal IMFs which are damped after 40, 20, and 7 s, respectively.

**Figure 12.** Comparison of real and reconstructed values of PG5 in Test Case 3.

**Figure 13.** Contribution factors of the generators G1 to G5 of the studied third system in modes 1–4.

44


**Table 5.** NCF Matrix in Test Case 3.

#### *3.4. Comparison with Related Methods*

In this section, we show the superiority of the proposed method for oscillation analysis. Due to space limitations, the proposed augmented Prony-based oscillation analysis is compared with the well-known and adopted techniques such as HHT and conventional Prony methods. It is clear from Table 6 that the proposed method is superior to others in terms of indices such as IMF, NCF, computational complexity, EE, Gibbs, and mode-mixing. For instance, the proposed method can eliminate the worse effects of EE, Gibbs and mode-mixing where other methods cannot. Considering the computation burden as an index, the proposed method is superior to HHT in which the complexity order of the proposed augmented method and HHT are O (1) and O (n × m), respectively. Where n is the number of the desired modes and m is the number of iterations in EMD process.

**Table 6.** Comparison of the Proposed and Other Related Works.


In order to further verify the proposed method, the results of the Test case 2 is compared with the results of PSAT software which is a MATLAB-based software, which can be used for small-signal stability analysis. Note that the analysis of PSAT is based on the detailed modeling of the generators, automatic voltage regulators (AVRs), power system stabilizers, PSSs, etc., whereas the proposed method is based on the synchrophasor measurements and the modeling of power system components is not needed. As seen from Table 7, the damping ratio and frequency of the dominant modes (Eig #21, Eig #22, Eig #13, and Eig #14) of the studied power system obtained by the detailed modeling in PSAT are comparable with the results of proposed method presented in Table 3. Both results, demonstrates that the dominant mode is an unstable mode. One of the capabilities of the PSAT toolbox is that it can determine the most-associated states related to each mode. The PSAT determines this using participation factors (PP) which be calculated using the elements of the state matrix A (if we consider the state equations as . x = Ax + Bu). Alternatively, in this paper we determine the most associated state variable corresponding to each oscillatory mode using the NCF index without the need to detail modeling of the system. As can be seen from Table 7, based on the PSAT results, the most associated generation in oscillatory modes 1, 2 is generator G4 which complies with the results obtained based on the NCF index in the proposed method.


**Table 7.** The Results of Test Case 2, Obtained by Detailed Modeling of Power System in PSAT Software.

#### **4. Conclusions**

In this paper, an online WAMS-based power system oscillation analysis framework was proposed. Features of the Prony method, such as simplicity in the implementation, self-verification, and robustness against the noise to estimate the modal information of a measured signal (i.e., its damping and frequency), nominated it for practical applications. The Prony method is then augmented here to be able to calculate IMFs and thus NCFs. IMFs provide an intuitive representation of the system modes and are mainly calculated using the iterative HHT-based methods. Our proposed framework calculates IMFs using an explicit mathematical formulation based on the augmented Prony method directly from PMU measurements. Therefore, it is simpler, faster and more accurate compared to those methodologies. The proposed Prony-based NCF index of this paper allowsthe identification of the most associated state variable to each mode accurately.

The simulation results on three test cases revealed that the proposed framework is able to determine IMFs and NCFs more accurately.

**Author Contributions:** All authors have contributed to the work. Conceptualization, Data curation, and Methodology: M.K.A., M.K., R.T.; Writing—review & editing: M.K.A., M.K., H.H.A., and P.S.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Improving Microgrid Frequency Regulation Based on the Virtual Inertia Concept while Considering Communication System Delay**

#### **Gholam Ali Alizadeh 1, Tohid Rahimi 2,\*, Mohsen Hasan Babayi Nozadian 2, Sanjeevikumar Padmanaban 3,\* and Zbigniew Leonowicz <sup>4</sup>**


Received: 6 April 2019; Accepted: 22 May 2019; Published: 26 May 2019

**Abstract:** Frequency stability is an important issue for the operation of islanded microgrids. Since the upstream grid does not support the islanded microgrids, the power control and frequency regulation encounter serious problems. By increasing the penetration of the renewable energy sources in microgrids, optimizing the parameters of the load frequency controller plays a great role in frequency stability, which is currently being investigated by researchers. The status of loads and generation sources are received by the control center of a microgrid via a communication system and the control center can regulate the output power of renewable energy sources and/or power storage devices. An inherent delay in the communication system or other parts like sensors sampling rates may lead microgrids to have unstable operation states. Reducing the delay in the communication system, as one of the main delay origins, can play an important role in improving fluctuation mitigation, which on the other hand increases the cost of communication system operation. In addition, application of ultra-capacitor banks, as a virtual inertial tool, can be considered as an effective solution to damp frequency oscillations. However, when the ultra-capacitor size is increased, the virtual inertia also increases, which in turn increases the costs. Therefore, it is essential to use a suitable optimization algorithm to determine the optimum parameters. In this paper, the communication system delay and ultra-capacitor size along with the parameters of the secondary controller are obtained by using a Non-dominated Sorting Genetic Algorithm II (NSGA-II) algorithm as well as by considering the costs. To cover frequency oscillations and the cost of microgrid operation, two fitness functions are defined. The frequency oscillations of the case study are investigated considering the stochastic behavior of the load and the output of the renewable energy sources.

**Keywords:** frequency stability; load frequency control; communication system delay; virtual inertia

#### **1. Introduction**

Microgrids are small power networks that are composed of loads, distributed generations, energy storage resources, telecommunication facilities and other infrastructure [1,2]. These networks are protected from damages resulting from the faults of the upstream network. The stochastic variations of generation units and load profiles [3,4] are common issues in these networks. The availability of different forms of energy storage in microgrids as well as their contribution to electricity markets at

different time periods justify microgrids more economically. These networks are suitable for critical applications such as hospitals and military centers, which require high levels of reliability. Natural disasters and system faults events can rarely affect the performance of these systems [5].

With the increasing penetration of renewable energy technologies in power grids, these grids experience new challenges. In a microgrid, power sources are connected to the network by power electronic converters. These interfaces have very fast transient responses, which reduces the overall inertia of the system. Consequently, low-inertia microgrids are less stable against any electrical disturbances. The imbalance between supply and load leads to frequency fluctuation. This problem, in turn, is aggravated by the presence of renewable energy sources. Renewable energy sources have unpredictable and fluctuating generation modes [6–9]. Fortunately, the energy storage units are capable of functioning as a power source or as load, and therefore can help attenuate the frequency fluctuations by absorbing or supplying power [10–14].

For the stable operation of microgrids, balancing between demand and supply in real time is an essential step, so that the frequency can be constant in its acceptable range. This task can be fulfilled using the load-frequency control mechanism [15]. The microgrid central controller divides the demand power between the power supply sources, thus maintaining the frequency of the system in the acceptable range [16]. Kalantar and Mousavi have studied the dynamical behavior of a hybrid system, which consisted of a wind turbine, microturbine, solar array, and battery according to the genetic algorithm and fuzzy controller [17]. Gao et al. have introduced active topology in its experimental and laboratory sense to control another hybrid system, including a battery and fuel cell [18]. Nayyeri Pour et al. have implemented the PI frequency control method in a local hybrid network composed of photovoltaic, a wind turbine, fuel cell, and super-capacitor. In addition, one PI controller is considered for each source, which determines the reference power of each unit for frequency stabilization by receiving frequency deviation [19]. In reference [20], a fuzzy set and Harmonic Search Algorithm is proposed to adjust the PI controller coefficients in the microgrid in the presence of an electric drive. In reference [21], an optimized fuzzy control system is used to overcome the frequency fluctuations of the microgrid against load and generation uncertainties. One of the modern methods to smooth frequency oscillation and improve the stability of microgrids is virtual synchronous generators [22,23]. In this method, pseudo inertial forces are realized by controlling DGs inverters to emulate the behavior of synchronous generators virtually into microgrids [24]. A Battery/Ultra capacitor hybrid energy storage system is introduced in reference [25] to cover power management and virtual synchronous generator emulation goals. In reference [26], an adaptive control method based on a virtual inertia system with MPC has been presented for frequency oscillation damping by emulating virtual inertia into the microgrid. The virtual inertia control is molded by a first-order transfer function in reference [26]. In reference [27], ultra-capacitor systems are controlled in a way that increases the virtual inertia of the studied microgrid. However, it is important to optimize the ultra-capacitor bank size to avoid extra costs.

On the other hand, with the development of computing and telecommunication systems, communication media are being used to exchange power between the central processing unit and the loads. However, with the increasing penetration of telecommunication networks in the control tasks, the delay between secondary controller and power resources cannot be ignored. Researchers have investigated the stability of the frequency control system in the presence of a time delay in several systems [27–29]. With decreasing communication delay, the cost of the communication infrastructure increases, however the frequency stability remains in a suitable state.

As discussed above, improving the inertia of microgrid and reducing the time delay of the communication system are effective methods to overcome frequency oscillations. However, this approach increases the operation cost of a microgrid. It is essential to select optimum values of time delay and Ultra capacitor size to have acceptable frequency oscillation of a microgrid and operation costs. To sum up, the challenges that discussed in this paper are given as follows:

• To reduce the installation costs, the size of the installed equipment should be minimal, so that the system will not face high costs during its long operation time.

• To enhance the frequency behavior of the microgrid, communicational infrastructure with minimum delays is required to update information quickly.

With considering the challenges that discussed above, this paper aim is determining optimal virtual inertia, comunication delay time and frequency control parameters simultaneously that are not investigated in privious papers. The proposed strategy guarantees frequency stability and decreases the cost of the microgrid's infrastructure. A Minimum ultra-capacitor size and communication delay value are selected in the proposed approach to decrease the cost of the microgrid's infrastructure. Also, the stable and unstable regions for different time delay and virtual inertia values are derived.

Multi-objective optimization methods are used to determine the parameters used in the load-frequency strategies to optimize the several objectives discussed previously [30]. The NSGA-II algorithm is a popular and justified optimizing algorithm [31,32]. Therefore, The NSGA-II algorithm is used in this paper to address the ultra-capacitor size as virtual inertia and the delay of the telecommunications system.

#### **2. Case Study**

Here, a microgrid system, which is isolated from the main grid, is considered for this case study. The islanded system consists of power generation, storage units, and telecommunication configurations. Diesel Generator Unit, Fuel Cell Unit, Solar Unit, and Wind Unit are responsible for power generation. The electrolyze unit is also use to absorb extra power. Ultra-capacitors are additionally used to generate and/or absorb power quickly in order to increase the virtual inertia of the network.

The central controller needs frequency information to appropriately order the microgrid power generation and absorption units. As discussed in reference [33], frequency and power information are to be sent or received with a specific time delay by a telecommunications platform, which is very effective for the stability of the microgrid. The communication delay with other delay sources are defined as a measurement delay in reference [34] and applied for primary and secondary controllers. The microgrid mentioned before is shown in Figure 1. Δ*f* is the frequency deviation and input signal of the controller, and the changes in the power generation and storage units as the output of this system. M represents the total inertia of a microgrid. When a difference between the generation and consumption occurs because of the load along with power generation increasing/decreasing, the frequency deviation of the microgrid depends on the inertia value and the load demand. It is assumed that the solar unit has variable and uncontrollable power. The central controller block sets a PI controller output for each power unit. The central controller unit block receives the network frequency information with a delay. The power and frequency information are passed through with a telecommunication system delay. The Diesel Generator Unit, Fuel Cell Unit, and electrolyze unit are equipped with primary and secondary frequency controllers. The Ultra-Capacitor unit participates in primary frequency control using a frequency signal. The parameters of this microgrid and power units are given in Table 1.


*Tdg* 2 (s) *fbase* 50 Hz

**Table 1.** Microgrid Constant Parameters.

**Figure 1.** Illustration of the Control Units of the Studied Microgrid.

#### *2.1. Modeling the Power Source of Diesel Generator*

The diesel generators are used to compensate the power required by the microgrid by using a governor and a speed drop system. The governor controls the fuel input of the fuel engine by adjusting the fuel throttle mechanism. The gasoline or gas engine acts as a turbine and moves a synchronous generator. The diesel generator governor can be modeled by a first-order conversion function that is described in Equation (1).

$$G\_{d\text{gs}}(\mathbf{s}) = \frac{1}{1 + sT\_{d\text{gs}}} \tag{1}$$

Similarly, the turbine of this generator is also modeled by the first order conversion function.

$$G\_{dt}(\mathbf{s}) = \frac{1}{1 + sT\_{dt}} \tag{2}$$

Due to the series performance of the mechanical and electrical unit of the diesel generator, the general conversion function is as follows.

$$G\_{d\lg}(s) = \frac{1}{1 + sT\_{d\lg}} \times \frac{1}{1 + sT\_{dt}} \tag{3}$$

The control commend to the diesel generator cannot be transferred to its output immediately because of the inherent delay of mechanical systems and generation rate constraint (GRCs). Therefore, these limitations are also applied to the Diesel Dynamic Model. The complete dynamic model of the diesel generator is shown in Figure 2.

**Figure 2.** Dynamic Model of the Diesel Generator [33].

#### *2.2. Fuel Cell*

The fuel cell is one of the static power generation equipment that converts the chemical energy stored in hydrogen into DC electrical energy. When the generated power of the renewable energy source is less than the power required for electric loads or the consumption of loads at their peak, the fuel cell begins to generate power by injecting hydrogen stored in its hydrogen tank. Since the ultimate performance of this system is similar to the voltage source inverter, the transfer function of this power source will be in the form of the first-order model as follows:

$$\mathcal{G}\_{fc}(s) = \frac{1}{1 + sT\_{fc}} \tag{4}$$

The generated power and its increase rate limit are also applied for the fuel cell. The dynamic model of a fuel cell is illustrated in Figure 3.

**Figure 3.** Dynamic Model of the Fuel Cell [33].

#### *2.3. Electrolyze Unit*

An electrolyze unit is a form of equipment that stores the extra energy of renewable energy sources. Water is decomposed into hydrogen and oxygen by passing the current between two separated electrodes. The hydrogen produced in the chemical reaction is stored in a tank and used as fuel for a fuel cell to meet load demand.

The water electrolyzes unit acts as a DC voltage source and requires a VSC inverter to be connected to a microgrid. Consequently, the transfer function will be as follows:

$$\mathbf{G}\_{AE}(\mathbf{s}) = \frac{\mathbf{K}\_{AE}}{1 + sT\_{AE}} \tag{5}$$

The power generation rate limits are also considered in the water electrolyze model. Figure 4 shows the dynamic model of the electrolyze unit.

**Figure 4.** Dynamic Model of the Electrolyzer [33].

#### *2.4. Model of Renewable Energy Sources*

In case of renewable energies, the MPPT strategy is implemented, so controllability will not be possible for wind/PV power generation [6,7]. These systems are not participating in frequency control.

#### *2.5. Ultra-Capacitor*

The Ultra-capacitors plays an important role in increasing the virtual inertia of the network, either by power injection or by absorption. The network exchange power of the Ultra-capacitor is directly related to the derivation of frequency fluctuations. These units have very little delay in response to frequency oscillations. Thus, the dynamic model of ultra-capacitors can be illustrated in Figure 5 [27].

**Figure 5.** Dynamic Model of Ultra-capacitor.

#### *2.6. Telecommunication System Delay*

The telecommunication system is used to transfer microgrid information to the central controller. In reference [34], a first order transfer function is considered for measurement of a delay that covers the communication delay. Therefore, for the simplification of analyses, all delays in communication units modeled as the first-order function in this paper, which is expressed in Equation (6). It should be noted that any other time delay model can be used and the optimizing algorithm can determine the optimum parameters of the selected model.

$$\mathcal{G}\_{\text{Comm}}(s) = \frac{1}{1 + sT\_{\text{comm}}} \tag{6}$$

#### *2.7. Optimization Algorithm*

The NSGA-II algorithm is used in this paper. The first goal is reducing the frequency fluctuations of the microgrid. Equation (7) is considered for this aim. The algorithm attempts to minimize (7) so the frequency oscillation will be in a minimum state. In addition, the second goal is determining the minimum ultra-capacitor size and the delay value of telecommunication systems to reduce the microgrid infrastructure cost. For a low ultra-capacitor size and the high delay value of telecommunication systems, the operation cost will be low. These two goals are expressed in two separate objective functions as below:

$$IAE = \sum\_{i=1}^{N} \left| \Delta f \right|\_i + \left| f \right|\_i. \tag{7}$$

$$\text{CostF} = \left(H\_{Vir\\_\text{average}} - H\right) + \frac{T\_{d\\_\text{average}}}{T\_d} = H\_{Vir} + \frac{T\_{d\\_\text{average}}}{T\_d} \tag{8}$$

In the above Equations, i.e., (7) and (8), Δ*f* stands for frequency fluctuations, *HVir* is the virtual inertia, *Td* is the delay of telecommunication systems, *Td*\_avarge is an average delay of telecommunication systems, and N is the number of samples in the simulation process.

The accuracy of resources and microgrid frequency response models has been investigated in the previous literature. These models have been used for the optimization process. The upper and lower limits of the decision variables of the optimization problem are listed in Table 2.


**Table 2.** The parameters limitation ranges.

The NSGA-II algorithm, like the conventional genetic algorithm, uses mutation and crossover operators at each stage. The algorithm uses repetitive loops to obtain the optimum population. The set of solutions for the right choice must first be ranked and each solution with fewer ratings is a better solution. This algorithm also begins its operation with the initial population, based on the values of the fitness functions. The solutions are ranked with the use of the mechanism of crowding distance and dominance, and based on this, a new population is generated with crossover and mutation functions [35]. As long as the two objective functions reach convergence, the algorithm should be continued. Figure 6 shows the convergence of these two conflicting functions in the the optimization algorithm steps.The multi-objective optimization is implemented in MATLAB/SIMULINK software

**Figure 6.** Pareto-based optimum solutions selection in the multi-objective optimization process.

#### **3. Simulation Results**

The frequency control system implemented in the MATLAB / SIMULINK software. The system parameters are obtained using the NSGA-II algorithm. The NSGA-II algorithm is used to optimize two conflict objectives, i.e. operation cost reduction and frequency oscillation reduction. The optimizing processes are showed as Pareto curves. The final Pareto curve is shown in Figure 7. The selected solution presented in the Pareto curve is considered an equally good solution. Low frequency oscillation with low operation costs can be derived in the selected point. The parameters values for the selected solution listed in Table 3. The optimum and non-optimum values are listed in Table 3. The non-optimum parameters were obtained by the trial and error method, reducing the infrastructure costs. The system was simulated in three scenarios:


**Figure 7.** Frontiers of Pareto of the two objective functions drawn from NSGA-II optimization.


**Table 3.** Optimum and non-optimum parameters of the control system, ultra-capacitor size, and telecommunications system delay values.

These scenarios show the effectiveness of the selected solution in overcoming frequency oscillation. In all the scenarios, the load and renewable energy source step changes are considered in the same profile. The magnitude of the load demand is equal to 1/1 per unit. When the simulation time reaches 1500 (s), it will reach 0.9 per unit. In contrast, the generation of renewable energy units is equal to one per unit. The frequency behavior of the microgrid in this section is shown in Figure 8. The microgrid frequency fluctuation has decreased with the application of optimum parameters and can be clearly seen in Figure 8. Meanwhile, the telecommunication system delay increases the undershoot and settling time of the frequency fluctuation. Up to t = 1500 s, diesel generator and fuel cell are expected to generate power, and after this instance the electrolysis unit is expected to absorb unbalanced power. The dynamic responses of the power sources and the electrolyze unit, while applying the optimum parameters, can be investigated using simulation results. Output powers of the source for the diesel generator, the fuel cell, and the electrolyze unit are shown in Figures 9–11, respectively. As can be seen from these figures, the electrolysis unit and diesel generator have the fastest and the slowest dynamic behavior responses, respectively. In the initial simulation time, the fuel cell injects power to grid with a high rate to compensate frequency drooping. Up to t = 1500 s, in steady state condition, the diesel generator and fuel cell contribute in power balancing and frequency stability.

Due to the high penetration of renewable energy sources, the actual inertia of the microgrid systems is small. This phenomenon shows its effect on steady frequency oscillation and high overshoot and undershoots in microgrid frequency responses. By increasing the communication delay, the stability margin of the microgrid can be improved. By increasing the microgrid inertia by increasing the virtual inertia value, the microgrid frequency is made less sensitive to the increasing commendation time delay. The studied microgrid behavior, for a different time delay value, is stable under the special condition that is shown in Figure 12. Two areas are distinguished in this figure. For a given Hvir value, if the time delay value is in the blue area, the system will be stable. The maximum virtual inertia margin that guarantees the stability is shown in Figure 13. To derive Figures 10 and 11, the conventional repetition process is used. The virtual inertia and time delay step changes are considered as being 0.5 Hz and 0.025 s, respectively.

With respect to general load and power generation profile, shown in Figure 14, the frequency behaviour of the grid while applying the variable profiles of the load power and power generated by renewable energy sources is shown in Figure 14. The frequency behavior of the grid with optimum

parameters (Figure 15) shows that the proposed strategy is more effective than the non-optimum based controller in damping frequency oscillation.

**Figure 8.** Comparison of the microgrid responses in terms of frequency in three scenarios.

**Figure 9.** The diesel Generator output power.

**Figure 11.** The electrolyzer output power.

**Figure 12.** The stable and unstable regains for the studied microgrid.

**Figure 13.** The maximum acceptable communication delay for given virtual values.

**Figure 14.** The renewable energy output and load demand powers.

**Figure 15.** Frequency deviation with optimum and non-optimum controllers.

In real grids, the output power of the renewable energy sources and load power consumption are stochastic. In this paper, the Gaussian distributions are assumed for the random output of the renewable energy sources and load profile, which expressed mathematically in Equations (9) and (10), respectively.

$$g\_1(\mathbf{x}) = \frac{1}{\sqrt{0.0005}} \exp\left(\frac{-1}{2} \left(\frac{\mathbf{x} - \mathbf{1}}{\sqrt{0.0005}}\right)^2\right) \tag{9}$$

$$g\_2(\mathbf{x}) = \frac{1}{\sqrt{0.001}} \exp(\frac{-1}{2} \left(\frac{\mathbf{x} - \mathbf{1}}{\sqrt{0.001}}\right)^2) \tag{10}$$

The frequency behavior of a microgrid while applying the specific uncertainties, optimum and non-optimum parameters is simulated and the results of these simulations are shown in Figure 16. The optimum-based controller can overcome power-unbalancing events and stabilize the frequency.

**Figure 16.** Frequency deviation under uncertainties.

#### **4. Conclusions**

A microgrid with controllable and uncontrollable sources, including solar, wind, a diesel generator, fuel cell, and an electrolyze unit, was investigated in this paper. The frequency stability of the grid was investigated while considering the ultra-capacitor unit and time delay of communication systems. The communication delay decreases the system function's ability to control the frequency. In this paper, tuning the virtual inertia constant and the time delay value in the communication

together with the parameters of the load-frequency controllers of the power sources was considered to improve the frequency stability of low inertia microgrids with the presence of the time delay in the communication network. In fact, the best way to tune these unknown parameters has been determined as a multi-objective optimization problem (NSGA-II algorithm), which considers both frequency stability and economic concerns. Reducing delay in telecommunication systems and increasing virtual inertia values, while applying the ultra-capacitor, can improve the damping of frequency fluctuations; however, it also increases the cost of grid operation. With the tuning of the control system parameters, time delay value and the virtual inertia value, the aforementioned challenges are solved. In the studied microgrid, the simulation was performed in MATLAB/Simulink numerical software. Investigation results prove the effectiveness of this network strategy to overcome frequency fluctuations.

**Author Contributions:** All authors involved equally and contributed for the final.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Black-Box Behavioral Modeling of Voltage and Frequency Response Characteristic for Islanded Microgrid**

#### **Yong Shi, Dong Xu \*, Jianhui Su, Ning Liu, Hongru Yu and Huadian Xu**

Anhui New Energy Utilization and Energy Saving Laboratory, School of Electrical Engineering and Automation, Hefei University of Technology, Hefei 230009, China; shiyong@hfut.edu.cn (Y.S.); su\_chen@hfut.edu.cn (J.S.); Ning.Liu@unb.ca (N.L.); 2017010029@mail.hfut.edu.cn (H.Y.); xuhuadian@mail.hfut.edu.cn (H.X.) **\*** Correspondence: 2017110310@mail.hfut.edu.cn; Tel.: +86-0551-6290-4042

Received: 8 May 2019; Accepted: 27 May 2019; Published: 29 May 2019

**Abstract:** The voltage and frequency response model of microgrid is significant for its application in the design of secondary voltage frequency controller and system stability analysis. However, most models developed for this aspect are complex in structure due to the difficult mechanism modeling process and are only suitable for offline identification. To solve these problems, this paper proposes a black-box modeling method to identify the voltage and frequency response model of microgrid online. Firstly, the microgrid system is set as a two-input, two-output black-box system and can be modeled only by data sampled at the input and output ports. Therefore, the simplicity of modeling steps can be guaranteed. Meanwhile, the recursive damped least squares method is used to realize the online model identification of the microgrid system, so that the model parameters can be adjusted with the change of the microgrid operating structure, which makes the model more adaptable. The paper analyzes the black-box modeling process of the microgrid system in detail, and the microgrid platform, including 100 kW rated power inverters, is employed to validate the analysis and experimental results.

**Keywords:** microgrid; recursive damped least squares; black-box modeling; online identification

#### **1. Introduction**

Since the microgrid technology can realize flexible and efficient applications of distributed power generations, more and more attention has been paid to its research [1,2]. However, as the scale of microgrid increases, its large-signal behavior becomes complex at the system level analysis. To deal with these problems, microgrid dynamic modeling is one of the most effective methods. At present, microgrid modeling methods mainly include two types: mechanism modeling and identification modeling [3]. The mechanism modeling is carried out on the premise of mastering microgrid control method and structure, and therefore the model has high precision accordingly. Using the constructed model, the small/large signal stability analysis can be performed and the high frequency response characteristics of microgrid system can be studied as well [4,5]. However, with the widespread use of commercial converters, the detailed information relating to converters' topology, parameters, and control methods is hardly to be obtained due to the commercial confidentiality and other reasons, which makes the mechanism modeling method extremely difficult to carry out. Compared with the mechanism modeling, the identification modeling method only needs to sample the input and output data to complete the identification modeling process, instead of mastering the detailed information of microgrid structure and control method. The identified model neglects to some extent the high frequency response characteristics, yet it is simple and convenient to use, thus increasing attention has been paid to it.

Identification modeling method was initially applied in the grid-connected photovoltaic (PV) distributed generation (DG) system modeling [6], and gradually extended to distributed power supplies and distributed power controllers [7–9]. For example, References [10,11] studied identification modeling method of DC/DC converter, and nonlinear models were selected to describe converters' output behaviors. In Reference [12], the black-box modeling process for a single three-phase inverter was proposed, and the dynamic parameters of the system were obtained through the step response. Moreover, the model also took the effects of cross-coupling into account, which further improved the accuracy of the model. In Reference [13], the parameters of the switched reluctance generator were identified, and the influence of nonlinear factors was considered as well, such as PI regulator and clamping functions. Currently, most of the research on microgrid identification modeling is aimed at the certain single converter integrated into the microgrid, such as DC/DC converter, photovoltaic inverter, and so on. By contrast, there are few studies on black-box identification modeling of microgrid at system level. Even for the identification modeling of the whole microgrid, the first step is to build non-linear models of the microgrid components, and then combine them into the microgrid model which is complex in model structure. For instance, Reference [14] described the relationships between the input and output voltage current of each DG by multi-segment nonlinear functions and then combined them together to build up the microgrid model. However, as the scale of the microgrid increases, the adaptability of this method is greatly limited. Reference [15] adopted the non-linear autoregressive method to construct the DGs integrated into the microgrid model by sampling the voltage data at the input port and the fault current at the output port, with up to 11 parameters to be identified.

On the other hand, most identification modeling methods are suitable for fixed microgrid structure, and off-line identification strategies can be employed to model microgrid system [16]. Although under some circumstances (such as similar power supply characteristics, operation mode is basically fixed, etc.), off-line identification has certain applicability, but for the common microgrid structure, there is a random change in power output and load switching, which makes the applicability of this method greatly being challenged. If the online identification method is adopted in identifying microgrid model, the real time update of the model parameters can make up for this deficiency. Nowadays, on-line identification methods are usually adopted in asynchronous motor parameters' identification [17,18], while application of these methods being adopted in microgrid system modeling has not yet been found in literatures.

Therefore, this paper proposes a novel black-box modeling method of microgrid system from a brand-new perspective. The model can be identified only by the voltage frequency data sampled at the point of common coupling (PCC) and the total active and reactive power references data obtained at the input port. The main features are as follows:


#### **2. Topology and Modeling Requirements of Microgrid**

#### *2.1. Typical Topology of Microgrid System*

The typical microgrid topology is shown in Figure 1. The system commonly comprises DGs such as photovoltaic and wind power generations, as well as large-capacity energy storage devices that maintain the stability of the microgrid voltage and frequency. These devices are connected to the AC bus via the inverters. Some active and reactive loads are also included in the microgrid system. Equipment with communication functions such as specific DGs and smart switches can be connected to the microgrid central controller (MGCC) which can also send control instructions to DGs and smart switches in turn. When the smart switch at PCC is disconnected, the whole microgrid system operates in islanded mode. At this time, the microgrid can adopt a peer-to-peer control structure, and then use the control strategy such as droop control [19] and virtual synchronous machine control [20] to form the voltage and frequency of the microgrid.

**Figure 1.** Topological Structure of Microgrid System.

#### *2.2. Modeling Requirements*

When studying the coordinated control of the microgrid system, usually the first step is to model and analyze the DGs included in the microgrid system, which is an important part before the overall modeling of microgrid. A rule of thumb is to use the mechanism modeling method to equivalent multiple energy storage inverters of the same type into one inverter so that the microgrid model can be greatly simplified. However, in most microgrid systems, besides the energy storage devices, there are various types of distributed power devices, which affect the system characteristics as well. These factors cannot be ignored, which lead to the model being too complex. Without the information of the internal structure of the microgrid system, the microgrid model can be equivalent to a black-box model and identification modeling method can be adopted to model the microgrid system.

According to the black-box model of microgrid as shown in Figure 2, the parameters optimization of secondary voltage and frequency controller can be carried out with less complexity. The black-box model of microgrid is considered as part of the open-loop transfer function of the overall control structure, and the corresponding root locus can be obtained according to the open-loop transfer function. Hence, the performance of microgrid system can be analyzed and the parameters of PI regulators can be optimized as well.

**Figure 2.** Microgrid Black-Box Model for Secondary Frequency and Voltage Regulation Control Structure.

#### **3. Identification of Black-Box Model**

#### *3.1. The Structure of Identification Model*

According to the modeling requirements discussed above, the equivalent structure of microgrid model can be obtained as shown in Figure 3. Instead of analyzing the internal structure of the microgrid, the black-box model can be directly constructed by sampling the active and reactive power references data at the input ports and the frequency and voltage amplitude data at the output ports.

**Figure 3.** Black-box model structure of microgrid.

In theory, if the line impedance and the output impedance of the inverters included in the microgrid system are purely inductive, the relationship between the frequency and the active power, as well as the relationship between the voltage and the reactive power are completely decoupled. However, in practical microgrid systems, the line impedance and the output impedance are not purely inductive, so the coupling term should be considered in the black-box model. Therefore, there are 4 parameters to be identified in the black-box model as shown in Equation (1):

$$
\begin{pmatrix} f\_{\text{PCC}} \\ E\_{\text{PCC}} \end{pmatrix} = \begin{pmatrix} G\_{fP}(s) & G\_{fQ}(s) \\ G\_{EP}(s) & G\_{EQ}(s) \end{pmatrix} \begin{pmatrix} P\_{0\text{total}} \\ Q\_{0\text{total}} \end{pmatrix} \tag{1}
$$

where *P*0total and *Q*0total refer to the total active and reactive power references in the microgrid respectively, and according to the power droop coefficients of each DG, MGCC assigns the total power reference values *P*0total and *Q*0total to each DG according to Equation (2):

$$\begin{cases} m\_1 P\_1 = m\_2 P\_2 = m\_3 P\_3 = \dots = m\_i P\_i \\\ n\_1 Q\_1 = n\_2 Q\_2 = n\_3 Q\_3 = \dots = n\_i Q\_i \\\ P\_{\text{Otatal}} = P\_1 + P\_2 + P\_3 + \dots \cdot P\_i \\\ Q\_{\text{Otatal}} = Q\_1 + Q\_2 + Q\_3 + \dots \cdot Q\_i \end{cases} \tag{2}$$

where *mi* and *ni* refer to the active and reactive power droop coefficients of *i*th inverter, *Pi* and *Qi* refer to the active and reactive power references of *i*th inverter, and *f* PCC and *E*PCC are frequency and voltage amplitudes of PCC. *GfP*(s) and *GEP*(s) are transfer functions of the active power reference with respect to the frequency and the voltage amplitude of PCC respectively, *GfQ*(s) and *GEQ*(s) are the transfer functions of the reactive power reference with respect to the frequency and the voltage amplitude of PCC respectively. Among them, *GEP*(s) and *GfQ*(s) are added to the model specifically in order to analyze the cross-coupling effect.

According to Equation (1), the microgrid system is regarded as a two-input and two-output black-box model. Without involving the characteristics of a single DG, the black-box identification model of the whole microgrid system can be obtained as shown in Figure 4.

**Figure 4.** Equivalent black-box model of voltage frequency response of microgrid.

By changing the total active and reactive power references issued by MGCC and sampling the voltage and frequency response data from the smart switch, the voltage and frequency response model of microgrid can be constructed by on-line identification method.

#### *3.2. Identification Experiments Design*

Since the step experiment process is easy to realize and has good identification performance, the step response experiment is adopted to identify the transfer function of the model. In order to simplify the identification process of the transfer functions mentioned above, this work divides modeled transfer functions into two groups. As shown in Figure 5, one group carries out the step experiment by performing a step on the active power reference value, the other group carries out the step experiment by performing a step on the reactive power reference value.

#### 3.2.1. Active Power Reference Step

The active power reference step experiment is applied to identify *GfP*(s) and *GEP*(s). The *Q*0total is set to 0, while the total active power reference is regarded as the input of the step experiment by performing a step on *P*0total. At the same time, the *f* PCC1 and *E*PCC1 are sampled at PCC as the output of the step experiment.

$$\begin{array}{l} G\_{fP}(s) = \frac{f\_{\text{PCC}}}{P\_{\text{Total}}} \Big|\_{Q\_{\text{total}} = 0} \\ G\_{EP}(s) = \frac{E\_{\text{PCC}}}{P\_{\text{Total}}} \Big|\_{Q\_{\text{total}} = 0} \end{array} \tag{3}$$

**Figure 5.** Input power references step experiments (**a**) Active power references step; (**b**) Reactive power references step.

#### 3.2.2. Reactive Power Reference Step

The reactive power reference step experiment is applied to identify *GfQ*(*s*) and *GEQ*(*s*). The *P*0total is set to 0, while the total reactive power reference is regarded as the input of the step experiment by performing a step on *Q*0total. Meanwhile, the *f* PCC2 and *E*PCC2 are sampled at PCC as the output of the step experiment.

$$\begin{array}{l} G\_{f\!Q}(s) = \frac{f\_{\text{PCC}}}{\text{Qtetaal}}|\_{P\_{\text{0total}} = 0} \\ G\_{EQ}(s) = \frac{F\_{\text{PCC}}}{\text{Qtetaal}}|\_{P\_{\text{0total}} = 0} \end{array} \tag{4}$$

#### *3.3. Identification Strategy and Process*

By sampling the input and output data of the microgrid system, the transfer function expressions of the microgrid model shown in Figure 4 can be obtained through an identification algorithm. The existing identification algorithms usually use the recursive least squares method to identify the system on-line and keep updating the parameters in real time. However, as the covariance matrix decreases in the recursive process, the parameters are prone to explode [21]. In order to enhance the stability of the recursive process, a damping term can be added to suppress the parameter explosion. Therefore, the recursive damped least squares method is used in this paper as the identification algorithm.

Moreover, in the classical control method for discrete systems, the bilinear transformation can ensure that the continuous system and the discrete system have the same stability and the same steady-state gain [22]. Therefore, the bilinear transformation is chosen to discretize the transfer functions shown in Equations (2) and (3), and the recursive damped least squares method is used for the parameter identification.

This paper takes *GfP*(*s*) as an example and analyzes the identification process in detail, other transfer functions *GEP*(*s*), *GfQ*(*s*) and *GEQ*(*s*) can adopt the same process.

#### 3.3.1. Data Preprocessing

The first task is to preprocess the sampled data *P*0total and *E*PCC1. It mainly includes two steps:


#### 3.3.2. Model Order and Initial Value Selection

Secondly, the number of coefficients of each polynomial (transfer function order) and the initial values of parameters should be determined. The first *n* data to be identified (in order to reduce the computation load, the selected data should not be too much) can be selected, and the least squares algorithm can be used to test a certain transfer function order iteratively. By comparing the fitness between the model output and the measured output in different orders, the order of the transfer function *GfP*(*s*) can be determined and remains constant in the process of recursion. The fitting performance can be calculated as follows:

$$f\_{\text{Rest fit}} = 100\% \times \left\{ 1 - \frac{\sqrt{\sum\_{k=1}^{N} \left( y(k) - \overline{y} \right)^2}}{\sqrt{\sum\_{k=1}^{N} \left( y(k) - \hat{y}(k) \right)^2}} \right\} = \frac{1}{N} \sum\_{K=1}^{N} y(k) \tag{5}$$

where *y*(*k*) is the measured output of microgrid, *y*ˆ(k) is the model output, *y* is the average value of the measured output. The higher the *f* Best fit, the better the model can reproduce the output characteristics of the microgrid port. Considering that the high model order affects the computing speed significantly, so as long as the fitness meets the requirements, choosing a lower transfer function order is recommended. Then, the order of the numerator of *GfP*(*s*) is set to 1, while that of the denominator is 2, the corresponding expression of *GfP*(*s*) can be obtained as

$$G\_{f\text{P}}(s) = \frac{f\_{\text{PCC}}(s)}{P\_{0\text{total}}(s)} = \frac{b\_1s + b\_2}{s^2 + a\_1s + a\_2} \tag{6}$$

The transfer function can be discretized by the bilinear transformation, as shown in Equation (7):

$$G\_{fP} = G\_{fP}(\mathbf{s})\Big|\_{\left(s = \frac{T\_2 - 1}{2z + 1}\right)} = \frac{B\_1 + B\_2 z^{-1} + B\_3 z^{-2}}{1 + A\_1 z^{-1} + A\_2 z^{-2}}\tag{7}$$

where

$$\begin{cases} A\_1 = \left(-T^2/2 + 2a\_2 - a\_1T\right)/\Delta \\ A\_2 = \left(T^2/4 + a\_1T/2 + a\_2\right)/\Delta \\ B\_1 = \left(b\_2 + b\_1T/2\right)/\Delta \\ B\_2 = \left(2b\_2 - b\_1T\right)/\Delta \\ B\_3 = \left(b\_2 + b\_1T/2\right)/\Delta \\ \Delta = T^2/4 + a\_1T/2 \end{cases} \tag{8}$$

#### 3.3.3. Online Recursive Algorithm

Then the corresponding difference equation of Equation (7) is:

$$f\_{\rm PCC}(k) = -A\_1 f\_{\rm PCC}(k-1) - A\_2 f\_{\rm PCC}(k-2) + B\_1 P\_{\rm Ototal}(k) + B\_2 P\_{\rm Ototal}(k-1) + B\_3 P\_{\rm Ototal}(k-2) \tag{9}$$

where *f* PCC (*k*) and *P*0total (*k*) are discrete values of *f* PCC and *P*0total at time *k*.

Let *k* = 3, 4,..., *m*, *m* is the total number of samples during the time of recursive computation. Equation (9) can be rewritten in the form of a matrix:

$$F\_m = H\_m \theta$$

where *F<sup>m</sup>* = [ *fPCC*(3), *fPCC*(4), ··· , *fPCC*(*m*)] T, *<sup>H</sup><sup>m</sup>* <sup>=</sup> [ϕ3,ϕ4, ··· ,ϕ*m*] T, <sup>θ</sup> = [*A*1, *<sup>A</sup>*2, *<sup>B</sup>*1, *<sup>B</sup>*2, *<sup>B</sup>*3] T, ϕ*<sup>m</sup>* = [ *f*PCC(*m* − 1), *f*PCC(*m* − 2), *P*total(*m*), *P*total(*m* − 1), *Ptotal*(*m* − 2)].

The idea of the recursive damped least squares method is to revise the estimated value θ recursively so that the minimum sum of squared errors between *<sup>F</sup><sup>m</sup>* and *F<sup>m</sup>* = *H<sup>m</sup>* θ*<sup>m</sup>* can be obtained:

$$\widehat{J(\theta)} = \left(\mathbf{F}\_{\text{ff}} - \widehat{\boldsymbol{H}\_{\text{m}}} \widehat{\boldsymbol{\theta}}\_{\text{m}}\right)^{\text{T}} \left(\mathbf{F}\_{\text{m}} - \widehat{\boldsymbol{H}\_{\text{m}}} \widehat{\boldsymbol{\theta}}\_{\text{m}}\right) = \text{min} \tag{11}$$

where *J*( θ) is the objective function. As discussed above, adding damping to the algorithm can suppress parameter explosion, so a damping term is added on the basis of equation (11):

$$f(\widehat{\boldsymbol{\theta}}) = (\boldsymbol{\mathcal{F}}\_m - \boldsymbol{H}\_m \widehat{\boldsymbol{\theta}}\_m)^T (\boldsymbol{\mathcal{F}}\_m - \boldsymbol{H}\_m \widehat{\boldsymbol{\theta}}\_m) + \mu \left[ \widehat{\boldsymbol{\theta}}\_m - \widehat{\boldsymbol{\theta}}\_{m-1} \right]^T \left[ \widehat{\boldsymbol{\theta}}\_m - \widehat{\boldsymbol{\theta}}\_{m-1} \right]^T \tag{12}$$

where μ is the damping factor, which is related to the linearity of the identification model.

This paper adopts the adaptive damping factor method to adjust the value of μ. Firstly, the matrix *HnH<sup>T</sup> <sup>n</sup>* is employed to determine the initial value of μ:

$$
\mu\_0 = 0.01 \times \text{mean}(H\_n H\_n^T) \tag{13}
$$

where the function mean is to compute the average value of the matrix. The linearity of the microgrid black-box model can be quantitatively evaluated by the ratio of the reduction and the change of the objective function *J* in the iterative process [22]:

$$\lambda = \frac{\overline{f(\widehat{\theta}\_{\mathfrak{M}})} - \overline{f(\widehat{\theta}\_{\mathfrak{M}-1})}}{\overline{f(\widehat{\theta}\_{\mathfrak{M}})} - \overline{f(\widehat{\theta}\_{\mathfrak{M}-1})}} \tag{14}$$

where *J*( θ*m*) is the calculated value of the objective function *J* at time *m*, *J*( θ*m*) is the second-order Taylor series expansion of the objective function *J* at time *m*:

$$\widehat{J}(\widehat{\boldsymbol{\theta}}\_{\rm m}) = \boldsymbol{J}(\widehat{\boldsymbol{\theta}}\_{\rm n}) + \boldsymbol{J}'(\widehat{\boldsymbol{\theta}}\_{\rm n})(\widehat{\boldsymbol{\theta}}\_{\rm m} - \widehat{\boldsymbol{\theta}}\_{\rm n}) + \frac{\boldsymbol{J}''(\widehat{\boldsymbol{\theta}}\_{\rm n})}{2}(\widehat{\boldsymbol{\theta}}\_{\rm m} - \widehat{\boldsymbol{\theta}}\_{\rm n})^2 \tag{15}$$

where *J* is the first derivative of *J*, *J*-is the second derivative of *J*.

Therefore, the damping factor can be corrected using the value calculated in Equation (14):


In order to achieve the adaptive correction of damping factor during the identification process, Proportional coefficient η is employed, whose value is generally between 2 and 10, as shown below:

$$\eta = \frac{f(\widehat{\boldsymbol{\theta}}\_{\rm m}) - f(\widehat{\boldsymbol{\theta}}\_{\rm m-1})}{(f\_{\rm PCC}(\boldsymbol{m}) - q(\boldsymbol{m}) \cdot \boldsymbol{\chi}\_{k}) \cdot \boldsymbol{\chi}\_{k}} + 2 \tag{16}$$

when μ is needed to be increased, μ can be taken as

$$
\mu = \mu\_0 \cdot \eta \tag{17}
$$

on the contrary, when μ is needed to be decreased, μ can be taken as

$$
\mu = \mu\_0 / \eta \tag{18}
$$

According to the extremum principle, the minimum value of *J* can be obtained by taking the differentiation to the right part of Equation (12). Through combing Equations (10)–(12), the recursive equation based on recursive damped least square method is derived as below:

$$\begin{cases} \boldsymbol{P}\_{m+1} = \left[ \boldsymbol{\mu} \boldsymbol{I} + \boldsymbol{H}\_{m-1}^{\mathrm{T}} \boldsymbol{H}\_{m-1} + \boldsymbol{q}\_{m}^{\mathrm{T}} \boldsymbol{q}\_{m} \right]^{-1} \\\ \boldsymbol{\theta}\_{m+1} = \boldsymbol{\theta}\_{m} + \boldsymbol{\mu} \boldsymbol{P}\_{m+1} \left[ \boldsymbol{\theta}\_{m} - \boldsymbol{\theta}\_{m-1} \right] + \boldsymbol{P}\_{m+1} \boldsymbol{h}\_{m+1}^{\mathrm{T}} \left[ \boldsymbol{f}\_{\mathrm{PCC}}(m) - \boldsymbol{q}\_{m+1} \boldsymbol{\theta}\_{m} \right] \end{cases} \tag{19}$$

Ever since a set of observation data ϕ*<sup>m</sup>* is added, *P*<sup>m</sup> and θˆ*<sup>m</sup>* are revised once, therefore new estimation values of parameters are derived and on-line identification of parameters are realized. The identification procedure is shown in Figure 6. Where ε is the allowable range of precision error, *N* is the total number of samples at the current moment. When the recursive algorithm compute to the sampling moment or the parameters estimation does not change, the recursive process is terminated until the new sampling data arrives. After obtaining the parameters of the transfer function *GfP*(*z*), the bilinear inverse transformation can be carried out to convert the transfer function from discrete domain to continuous domain. θ θ θθ θ

**Figure 6.** Flowchart of identification procedure.

#### **4. Experiments and Model Validation**

In this paper, a series of model identification experiments are carried out on a microgrid test platform to verify the accuracy of the proposed modeling method. The controllers of the two inverters integrated into the microgrid use MCU (TMS320F28335, Texas Instruments, Inc, Dallas, TX, USA), the pulse width modulation (PWM) carrier frequency of each controller is set to 6 kHz, and the integrated three-phase resistive load is configured to 30 kW. Both inverters operate in *P*-*f* droop control mode. Active and reactive power droop coefficients are chosen as 4 <sup>×</sup> 10−<sup>6</sup> and 4 <sup>×</sup> 10−<sup>4</sup> respectively. Under the normal operation of the experimental platform, the input and output data are sampled by the Yokogawa recorder. All devices are connected to the AC Bus and communication lines, the structure of the microgrid test platform is shown in Figure 7.

**Figure 7.** The structure of microgrid test platform.

#### *4.1. Identification Experiments of GfP(s) and GEP(s)*

To carry out the identification experiments of *GfP*(*s*) and *GEP*(*s*), parallel inverters and the three-phase active load are integrated into the microgrid. After the system reaches a steady state, the active and reactive power references are set to 0 by MGCC at first. Then the active power reference is changed to 80 kW at 9 s, and set to 0 again at 15.7 s. During this experiment, the reactive power reference remains at 0. Figure 8 shows the identification data *P*0total, *E*PCC1, and *f* PCC1 sampled by the waveform recorder at PCC.

According to the identification process shown in Figure 6, the sampled data *P*0total and *f* PCC1 can be used to identify the model transfer function *GfP*(*s*).

$$G\_{fP}(s) = \frac{4.229 \cdot 10^{-8}s^2 - 0.001s + 5.0184}{s^2 + 945.97s + 2.5172 \cdot 10^5} \tag{20}$$

By applying the same input on both the model and the microgrid system, the fitting performance can be obtained as shown in Figure 9. Since only the dynamic component is analyzed, the steady-state component contained in each measurement has to be removed. As it can be seen, the model output matches the microgrid output properly, and *f* Best fit reaches 95.36.

**Figure 8.** Active power step sample data for identification (**a**) Active power references; (**b**) Reactive power references; (**c**) PCC frequency; (**d**) PCC voltage amplitude.

**Figure 9.** Fitting curves of the transfer function model obtained from the active power step experiment.

Similarly, using the sampled data *P*0total and *E*PCC1, the transfer function G*EP*(*s*) can be obtained:

$$G\_{EP}(s) = \frac{0.06107s^2 + 64.24s + 152.9}{s^3 + 1249s^2 + 3.95 \cdot 10^5 s + 4.933 \cdot 10^8} \tag{21}$$

#### *4.2. Identification Experiments of GfQ(s) and GEQ(s)*

To carry out the identification experiments of *GfQ*(*s*) and *GEQ*(*s*), parallel inverters and the three-phase reactive load are integrated into the microgrid. After the system reaches a steady state, the active and reactive power references are set to 0 by MGCC at first. Then the reactive power reference is changed to 80 kVar at 9 s, and set to 0 again after 5.6 s. During this experiment, the active power reference remains at 0. Figure 10 shows the identification data *Q*0total, *E*PCC2, and *f* PCC2 sampled by the waveform recorder at PCC.

**Figure 10.** Reactive power step sampling data for identification (**a**) Active power references; (**b**) Reactive power references; (**c**) PCC frequency; (**d**) PCC voltage amplitude.

According to the identification process shown in Figure 6, the sampled data *Q*0total and *E*PCC2 can be used to identify the model transfer function *GEQ*(*s*).

$$G\_{EQ}(s) = \frac{1.0129 \cdot 10^{-5}s^2 - 1.288s + 7826.65}{s^2 + 2760.77s + 2.545 \cdot 10^6} \tag{22}$$

By applying the same excitation on both the model and the microgrid system, the fitting performance can be obtained as shown in Figure 11. Clearly, the model output also matches the microgrid output, and *f* Best fit reaches 90.72.

Similarly, using the sampled data *Q*0total and *f* PCC2, the transfer function *GfQ*(*s*) can be obtained as follows:

$$G\_{fQ}(\mathbf{s}) = \frac{8.085 \cdot 10^{-6} \text{s}^2 + 2.219 \cdot 10^{-4} \text{s} + 6.091 \cdot 10^{-4}}{\text{s}^3 + 8.47 \text{s}^2 + 710.6 \text{s} + 5187} \tag{23}$$

**Figure 11.** Fitting curves of the transfer function model obtained from the reactive power step experiment.

#### *4.3. Recursive Model Validation*

According to the power references step experiments, the voltage frequency response model is obtained under the certain microgrid structure. Furthermore, in order to verify the adaptability of the recursive damped least squares method applied in the model identification, that is, the online identification effect during the period of microgrid structure changing, several adjustments have been made: active power droop coefficients of two inverters are changed to 4 <sup>×</sup> 10−<sup>6</sup> and 8 <sup>×</sup> 10−<sup>6</sup> respectively; reactive power droop coefficients are changed to 2 <sup>×</sup> <sup>10</sup>−<sup>4</sup> and 4 <sup>×</sup> <sup>10</sup>−<sup>4</sup> respectively; the PV inverter which is supplied by Chroma solar PV array simulators is added to the microgrid structure, and the maximum output power of the PV inverter is set to 4 kW. During the power references step experiment, the number of inverters connected to the microgrid are changed, the active and reactive loads are set to different parameters, as listed in Table 1.


**Table 1.** Experimental parameters.

The input and output data during the experiment are recorded by the waveform recorder, as shown in Figure 12.

Meanwhile, the same experimental processes as shown in Table 1 are carried out on a corresponding simulation model under the MATLAB/Simulink environment (2014a, The MathWorks, Inc, Natick, MA, USA), which is constructed according to the black-box equivalent model shown in Figure 4. Initial values of *P* and θˆ can be figured out according to the transfer function parameters obtained in Sections 3.2 and 4.1. The frequency and voltage output of the simulation model can be obtained by using the recursive damped least squares method for online identification.

The experimental and simulation results are compared with each other, as shown in Figure 13. It can be seen that the simulation model is able to reproduce the system responses of the microgrid system during the experiment, especially when the frequency and voltage responses of the microgrid have step changes, the output of the simulation model tracks the changes accurately. Hence, the identification algorithm implemented in the simulation model is validated.

**Figure 12.** Data waveforms sampled at input and output ports of the microgrid (**a**) Active power references; (**b**) Reactive power references; (**c**) PCC frequency; (**d**) PCC voltage amplitude.

**Figure 13.** Identification model and measured output curves (**a**) PCC frequency; (**b**) PCC voltage amplitude.

#### **5. Conclusions**

In this paper, the black-box modeling technique of voltage frequency response model of microgrid system has been studied. In terms of the parameter identification algorithm of the recursive damped least squares, the model can be identified and is suitable for performing simulations of microgrid system on system level, especially when the internal information of a microgrid cannot be accessed. The procedure of constructing the black-box model of microgrid is analyzed emphatically. Compared with the traditional method of mechanism modeling, this method greatly simplifies the modeling steps and provides a general idea for microgrid system modeling.

The output data obtained from the simulation model has been compared with the measured output under the conditions of changing the microgrid structure. Depending on the recursive algorithm, the identification model can reproduce the frequency and voltage characteristics of the microgrid port accurately in all cases, which indicates that the black-box identification modeling method has a wide range of adaptability to realize on-line identification of the microgrid system.

In future, in order to improve the identification accuracy, the influence of different communication delays on the proposed identification strategy needs to be further studied. In addition to this, how to apply the proposed identification model to a secondary frequency and voltage regulation control system, so as to optimize parameters of secondary voltage and frequency controllers, are also set as the future research directions.

**Author Contributions:** This paper was a collaborative effort between the authors. Y.S. provided the original idea and reviewed the paper, D.X. and H.X. organized the manuscript and realized MATLAB-based simulation, H.Y. carried out experiment validation, J.S. and N.L. provided professional advice and technical comments.

**Funding:** This research was funded by Double First Class Project for Independent Innovation and Social Service Capabilities (45000-411104/012), National Key Research and Development Program of China (2017YFB0903503) and the Basic Operating Expenses of Central Scientific Research (PA2018GDQT0021).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Sensor-Based Early Activity Recognition Inside Buildings to Support Energy and Comfort Management Systems**

#### **Francesca Marcello 1,†, Virginia Pilloni 1,2,\*,† and Daniele Giusto 1,2**


Received: 20 June 2019; Accepted: 5 July 2019; Published: 9 July 2019

**Abstract:** Building Energy and Comfort Management (BECM) systems have the potential to considerably reduce costs related to energy consumption and improve the efficiency of resource exploitation, by implementing strategies for resource management and control and policies for Demand-Side Management (DSM). One of the main requirements for such systems is to be able to adapt their management decisions to the users' specific habits and preferences, even when they change over time. This feature is fundamental to prevent users' disaffection and the gradual abandonment of the system. In this paper, a sensor-based system for analysis of user habits and early detection and prediction of user activities is presented. To improve the resulting accuracy, the system incorporates statistics related to other relevant external conditions that have been observed to be correlated (e.g., time of the day). Performance evaluation on a real use case proves that the proposed system enables early recognition of activities after only 10 sensor events with an accuracy of 81%. Furthermore, the correlation between activities can be used to predict the next activity with an accuracy of about 60%.

**Keywords:** activity recognition; activity detection; activity prediction; smart building; energy and comfort management

#### **1. Introduction**

Smart buildings are characterized by the presence of sensors, actuators, and smart devices that give the opportunity to monitor and remotely control key equipment within buildings [1]. This is the concept behind Smart Building Energy and Comfort Management (BECM) systems [2,3]. In such an intelligent scenario, one of the major goals is to provide decision-support tools that support users in making cost-effective decisions in terms of energy consumption [4]. As a matter of fact, domestic electricity usage accounts for about 40% of the global energy consumption and contributes over 30% of total greenhouse gas emissions [5]. Nevertheless, user comfort is crucial when policies of Demand-Side Management (DSM) are put in place [6]. Indeed, a system that optimizes energy consumption without considering user preferences and habits about appliance usage quickly leads to user disaffection from the system and abandonment of it. Currently, most of the literature considers user comfort as a set of hard constraints on appliance usage, which are a priori set considering general statistics [7,8]. This approach neglects the fact that users are likely not only to have different subjective requirements with respect to the others, but they also dynamically change over time.

In this paper, user preferences and habits about appliance usage are inferred by monitoring them using sensors deployed inside the reference buildings. The first phase consists of analyzing the correlation between the information gathered from sensors and users' activities. The approach proposed in this paper is based on the classifier presented by Krishnan and Cook in [9], which recognizes activities based on sequences of sensor events. Nevertheless, as it will be better explained in the following, in order to improve accuracy, the classifier included in the proposed framework is enhanced with other significant information provided by statistics related to the monitored events. Accordingly, a profile specific to the monitored user is created, which enables early recognition of activities after only 10 sensor events with an accuracy of 81%. Furthermore, the correlation between activities can be used to predict the next activity with an accuracy of about 60%.

The main contributions provided by this paper can be summarized as follows:


The remainder of the paper is organized as follows. Section 2 presents past works and the required background. In Section 3 an overview of the considered system model is provided. Section 4 presents in detail the activity recognition algorithm that enables the early detection and prediction of users' activities. Section 5 describes the reference use case considered to test the performance of the system. Finally, in Section 6 a performance analysis is provided, and conclusions and final remarks are drawn in Section 7.

#### **2. Background**

#### *2.1. Smart Building Energy and Comfort Management Systems*

Smart technologies can be used in all kinds of different buildings (i.e., residential, office, and retail sectors) to improve the comfort and the safety of people in their home, concerning various topic, from healthcare and providing living assistance, to environmental monitoring and ensuring energy saving. Accordingly, BECM systems have the objective of combining power consumption minimization while preserving user comfort [2,10]. This issue has been addressed by researchers from many different perspectives. In [11] a system for intelligent energy management in buildings is proposed. It presents semantic modeling that integrates all the entities that constitute the environment of a Smart Building exploiting Internet of Thing (IoT) paradigm. The IoT-based system integrates different systems and makes use of various types of real-time data from different sources to achieve the common objective of the intelligent management of the building. In [12] a tool that provides effective automation and control of heating/cooling, ventilation/air conditioning and lighting and that uses optimization techniques to minimize energy consumption is proposed. The interactive system presented achieves a significant decrease in the operating cost of A/C system in a tertiary sector building, while maintaining desirable comfort taking into account two different time periods (peak hours and non-peak hours), a number of different zones and the end-user's preferences. The authors in [13] propose a multi-agent control system for integrated buildings and microgrids, which exploits Renewable Energy Sources (RES) efficiently among the agents. In [7], an algorithm for Distributed Energy Resources (DER) management is proposed to shave peak demand and increase energy efficiency in smart home environments. Customer comfort is considered to be a constraint on time preference ranges to run their appliances. A similar approach was presented in [8], where user preferences are expressed also in terms of indoor temperature and lightning. Furthermore, dynamic pricing is considered

to optimize DSM. The authors in [2], after a review of control systems for energy management and comfort in buildings, also present the architecture of a multi-agent control system that manages the user's preferences for thermal and lighting comfort, indoor air quality, and energy conservation. The system is based on a master-slave coordination mechanism to perform different tasks. In [14], an algorithm for thermostatically controlled household loads based on price and consumption forecasts of grid energy is presented. The optimization by means of the algorithm proposed in [14] takes into account the trade-off between customer comfort and cost of energy, by setting minimum and maximum boundaries for the thermal comfort. These boundaries are taken as hard constraints for a priori setting but a background on customer comfort profiling lacks. Collotta and Pau propose a BECM system in [15], where the management is based on a Fuzzy Logic Controller (FLC) which adaptively adjusts appliance starting times based on users' feedback. The issue of scheduling appliances according to user preferences was also addressed by [4], where Quality of Experience (QoE) is measured as a function of the interval between the preferred and proposed appliance starting time for switching controlled loads (e.g., washing machines and clothes dryers), and as a function of the interval between the preferred and proposed temperature for thermostatically controlled loads (e.g., Heating, Ventilation and Air Conditioning (HVAC) and water heaters). The system takes into account both dynamic pricing and RES production.

It is evident that user preferences and habits severely affect results of BECM systems. Indeed, BECM systems that only consider energy cost minimization may switch appliances on/off too early or too late. This is the case, for example, of a dishwasher that does not finish its cycle before dinnertime, or an HVAC that is switched off earlier than what the user considers thermal comfort. For this reason, in recent years researchers have started to observe users' behavior, in order to infer their habits and preferences.

#### *2.2. Activity Recognition in Smart Buildings*

As discussed in the previous subsection, BECM systems should be able to discover and predict users' habits and course of actions. The monitoring of activities of people in their home can be done by analyzing data that can be gathered with different technologies. Given the rapid development of sensor technology, wireless transmission technology, network communication technology, cloud computing, and smart mobile devices, large amounts of digitized information have been accumulated and the volume of data is growing rapidly with increasingly complex structures and forms [16]. The energy big data provides a new way to analyze and understand individuals' energy consumption behavior, improve energy efficiency, and promote energy conservation. Integration of big data technologies will help make the grid more efficient and it will fundamentally change the way in which regulators, utilities, grid operators, and end-users would interact [17,18].

In many cases, cameras and wearable sensors are used to collect all the information of interest and to understand what someone is doing [19]. These solutions present some problems because people are often not inclined to accept those devices [20]. Some studies are based on the data that are provided by phone accelerometer and gyroscope to understand repetitive body motions (walking, running, sitting) [21]. This solution is not very practical in home scenarios, where residents do not always take their phone with them. To monitor what activities people are performing in their house, non-intrusive sensors are often preferred: typical devices that are installed in the environment are motion sensors, door sensors or temperature and pressure sensors [22,23].

The authors in [24] propose the use of 3D depth sensors to detect user occupancy and profile their habits. The aim of the paper is to use this information to control HVAC and lighting, according to occupants' usual behavior. The relation between occupancy, energy consumption and users' comfort is also investigated in [25], where the Multi-Agent Comfort and Energy Management System is proposed. Along with building devices, MACES also considers occupants as active participants in the building energy reduction strategy and attempts to implement more energy conscious occupant

planning. This occupant planning is carried out using multi-objective Markov-Decision Problems (MDPs) to model the uncertainty of agent decisions and interactions.

An interesting approach for energy-consuming activity recognition is proposed in [26], where social media posts are analyzed to automatically extract information and describe energy-consuming activities.

In [27] the authors focus on the activity discovering problem, proposing an approach to build a model under the form of Hidden Markov Model (HMM), from a training database of observed events emitted by binary sensors, without the knowledge of actions really performed during the learning period. The discovering problem is presented also in [28,29], where a Discontinuous Varied-Order Sequential Miner (DVSM) algorithm is used to discover frequent activities that are continuously recorded in a smart environment and combined with a clustering algorithm to find frequent occurrences of activities and cluster familiar patterns together.

In [9], 4 algorithms that can identify the activities while they are being performed are proposed. In this work, activities are recognized even if they are done in an interleaved and concurrent manner. The models used in the proposed algorithms are a Naïve Bayes Classifier, an HMM with a time window, a frequency-based HMM with a sliding window and a frequency-based HMM with a shifting window. These algorithms have been later used in [30] to predict activities with the aim of controlling buildings to reduce energy consumption.

The problem of multi-resident activity recognition based on the use of non-intrusive sensors, along with smartphone-based sensed data, is also addressed in [31]. The approach used in this paper aims to exploit body-worn smartphone sensors to infer person-specific context that can be correlated with activity detected from ambient sensors. Other models for multi-resident activity recognition based on the use of non-intrusive sensors are proposed in [32]. In this paper the authors adopted three different directed graphical models including Poisson HMMs, coupled HMMs, and dynamic Bayesian networks, extended from coupled HMM by adding some vertices, to identify both individual and cooperative activities. A solution for multi-resident environments is given in [33] too, where they perform the tracking of people and recognition of activities by using different binary sensors and RFID to know the identity of the occupants as they enter or leave the environment.

#### *2.3. Background on Activity Recognition*

The data collected from sensors inside resident houses are analyzed using data mining and machine learning techniques to build activity models that are used as the basis of behavioral activity recognition. Feature extraction from the sequence of sensor events is a key step to better modeling and then recognizing human activities. In [34] four methods used to extract features for online recognition on streamed data are presented. With their approach, the authors can recognize activities while new sensor events are recorded.

A comparison of classification approaches for activity recognition is provided in [35,36]. As introduced in the previous subsection, with reference to modeling and classification methods researchers have investigated the recognition of resident activities using a variety of mechanisms, such as naïve Bayes classifiers, Markov models, and dynamic Bayes networks.

The various models and algorithms can give different performance results depending on the input data or on the implemented system. Usually, the results are adequately comparable, with one model giving better results in the recognition of some specific activities or using one specific dataset, and others giving better performance recognizing other activities or analyzing other datasets [37]. In multiple cases, in spite of its simple design and simplified assumptions, naïve Bayes classifiers often work much better than expected, especially when a specific group of sensors can easily be identified as characteristic of a certain activity [9].

Starting from the approach presented in [9] and from the analysis of sequences of sensor events for activity recognition solution, in this paper a system that adds statistics information about the context in which these activities are occurring, in order to improve their own recognition, is proposed. The ultimate goal of the system is to obtain a profile according to users' habits, which allows gaining personalized management of the whole system. Many systems in the literature propose management strategies focused on energy consumption and consider users' comfort, but they tend not to evaluate user activities in every daily aspect in order to learn specific habits and behavior and guarantee ad-hoc solutions.

#### **3. System Model**

The reference scenario considered in this paper is that of a BECM system that leverages distributed smart home sensor networks to elaborate user profiles that are later used to manage and control buildings. More specifically, sensors are used to make observations, report events that are detected in the building and that can be associated with users' actions, and learn which sets of detected events can identify specific activities. Therefore, the BECM system can predict users' activities based on their previous monitored actions, and make appropriate management decisions accordingly.

A fundamental step to achieve energy cost savings is the identification of possible causes of energy waste. The BECM system controls all energy consumption components (electric, gas, water) inside the buildings and plans specific actions aimed at reducing waste. A large number of home devices increase power consumption in two aspects: standby power and normal operation power. When high energy demand is detected, the less important loads, such as standby power mode devices, could be disconnected from the electricity grid. Around 10% of the total household power is consumed during the standby power mode [38]. The reduction of standby power is greatly necessary to reduce the electricity cost at home. Thanks to IoT networks, appliances such as washing machines and dishwashers can gather information about different energy prices related to different times of the day, thereby the system can automatically program them for the most convenient times. Air conditioning and heating systems should be monitored to guarantee desirable temperatures in all environments and to avoid unnecessary use. The BECM system records the inhabitants' preferences adapting itself to their way of life.

A system of this type must also be able to know the users' habits, in order to make coherent scheduling in the management of equipment and appliances, and because their activities and behavior have a considerable impact on energy consumption.

The whole system is therefore divided into several modules, each of them with specific features depending on the different tasks pertaining to them. Figure 1 shows the overall architecture of the system and how the modules interact with each other.

The "Acquisition Module" is the one that gets raw data from the sensors. The data are stored in a structured database with information about the sensors and their value, as well as the date and time of the relative information. These data are then processed and manipulated by the "Activity Recognition Module". This module can be organized into two separate tasks. At first, the data are used to understand users' behavior and to make models of the different activities usually performed; at a later time, the module has to recognize occurring activities based on previously created models. All the information about resident habits and preferences is provided as the output of this module. This information is then used by the "Energy and Comfort Management Module" to monitor and act on appliances and devices according to a precise scheduling algorithm based on user profiles and on activities performed. Indeed, each activity can be associated with a specific set of appliances that are turned on or adjusted accordingly. Table 1 includes the details of the most typical home appliances, along with their probability to be available at home in Italy [4]. Based on the activity recognition and/or prediction results, the appliances can be scheduled to improve energy savings and/or users' comfort [4]. The decisions taken by this last module can operate on the system automatically, via commands sent to actuators, or they can be translated in useful advice sent to the user through some interface. All the described modules are integrated into an intelligent device that oversees the data storage and the control of the building, either locally or in the cloud [39,40].

**Figure 1.** Integration of the activity recognition system inside a BECM system.


**Table 1.** Characteristic parameters of the most typical home appliances [4].

As an explanatory example, suppose that according to their profile, a user usually wakes up, then has a coffee watching TV, and later takes a shower with the bathroom heater on. Furthermore, suppose that the *wake up* activity is detected when the *turn on the bedroom light* event is identified (i.e., the light sensor inside the user bedroom detects some light) after the *sleep* activity occurred. Since the BECM system knows that the following activities are *have a coffee*, *watch TV* and *take a shower*, it can increase the user's comfort by: turning on the coffee machine as soon as the user wakes up, turning on the TV right after the coffee is made, and at the same time turning on the water and bathroom heater so that the water and room are warm when the user goes to the bathroom. Alternatively, the BECM system can choose to switch water and bathroom heater on earlier if the energy cost is lower, based on predictions about the usual morning routine of the user.

Figure 2 expresses the relationship between sensors, events, and activities. An event corresponds to a change in the state of a sensor. Each sensor has several possible states, depending on the significant values it measures; this is the case, for example, of a contact on a door, which has two different states (e.g., *OPEN*, *CLOSE*): an event is registered whenever the contact changes its state from *OPEN* to *CLOSE* and vice versa. On the other hand, a smart meter monitoring a washing machine can have multiple states: according to the measured power consumption, its related state can correspond to *OFF* or to any wash cycle. Accordingly, an event is registered whenever the smart meter changes its state.

**Figure 2.** Example of activities and relevant subdivision into actions and observations, in connection with related sensors.

Let E = {*ei*} be the set of events that can be detected by sensors inside a house. An event is then defined as the transition between possible states of the corresponding sensor monitoring that particular situation. According to this vision, activities are composed of several events. Even though the events that characterize an activity remain basically the same, their order may change when considering two different observations of the same activity. For example, when preparing a meal one can open the fridge to take the ingredients and then put a pan on the stove, or the same activity can be done in the opposite order. For this reason, a generic activity A*<sup>j</sup>* is modeled as the set of events that are usually observed when the reference user performs it. Accordingly, an instance *Ijk*(A*j*, *t s k*, *t e <sup>k</sup>*) of activity A*<sup>j</sup>* is defined as the array of events that are observed from the time *t s <sup>k</sup>* the activity started to the time *t e <sup>k</sup>* the activity finished. In the following sections, the process used to recognize the activities based on the observed events will be described in detail.

#### **4. Activity Recognition Algorithm**

The set of events E is the basic set of information needed to understand what users are doing, and as seen previously they are strongly connected to sensors activation. The proposed system is then based on recognizing activities performed in smart environments from sequences of collected sensor readings. The choice of which sensor types to use leads to different models and algorithms for solving activity recognition problems because of different kinds of data, produced according to the type of sensor that has to be manipulated in order to extract important information.

As seen in Section 2, the solutions for monitoring people in their home include: 1. cameras, 2. wearable sensors and 3. different kinds of ambient sensors. The first solution gives the type of information with the highest accuracy, but it requires the heaviest process to manipulate images and videos. Moreover, some people find problematic the idea of having cameras in their homes. Wearable sensors give information about the physical state of a person, so it is possible to understand simple actions, but it is harder to recognize more complex and complicated activities. Besides, there could be problems if the users forgot to wear them. Ambient sensors produce data with a low level of semantic information, but choosing only non-intrusive binary sensors is a better option for experiments in real life so that there is no need for people to remember to always wear wearable sensors or to be monitored with cameras.

The flowchart in Figure 3 expresses and clarifies all the steps of the proposed activity recognition algorithm that are described in detail in the following lines. The algorithm is constituted by two main phases: the training phase, during which the characteristics and statistics that describe how the observed user performs the activities are created; the running phase, where the events are observed and processed so that activities can be recognized by the classification module.

**Figure 3.** Flowchart of the steps for the activity recognition algorithm.

#### *4.1. Training Phase*

During the training phase, each activity instance is observed within a time window <sup>O</sup>*<sup>A</sup>* included in {*t s k*, *t e <sup>k</sup>*} and is defined by the sequence of sensor events, i.e., sensors that have changed their state within the considered window. For each observed activity, a feature vector *F jk*(*Ijk*)=[ *f*1*jk*, *f*2*jk*,..., *fijk*,... ] is computed with the rates of event occurrences for each sensor that is the number of events related to one specific sensor with respect to the total number of events observed considering all the sensors within the time window <sup>O</sup>*A*. Then, for each type of activity A*j*, a model vector *mj* = mean*k*(*F jk*)=[ *f* <sup>1</sup>*jk*, *f* <sup>2</sup>*jk*, ... , *f ijk*] is defined such that the rates of event occurrences of its sensors is the average rate for all the observed instances associated with the same activity. This model vector *m<sup>j</sup>* is representative about the probability that a sensor is connected to every activity. When an event associated with a sensor is counted for a certain number of times during the observed sequence, it is possible to understand which activity is statistically more probable.

From raw data acquired by sensors, feature vectors that can be analyzed to perform classification of activities are constructed. Every activity is strongly connected to a specific group of sensors that change their states during a definite time span, as described in the previous section.

Furthermore, statistics Γ*<sup>j</sup>* about relevant conditions that can be associated with the performed activity A*<sup>j</sup>* are evaluated and stored by the system. Indeed, activities are often performed within the same time window (e.g., sleeping, preparing meals). Therefore, Γ*<sup>j</sup>* includes the following statistics about A*j*: the average starting and ending times *t s* <sup>A</sup>*<sup>j</sup>* and *t e* A*j* , and their standard deviations *σ<sup>s</sup>* A*j* and *σ<sup>s</sup>* A*j* . If there is more than one time window, statistics are generated for each observed time window.

#### *4.2. Running Phase*

After the probabilistic model is obtained, the system must recognize the activities performed by evaluating which the most likely to be happening is. To this aim, a sensor-based windowing implementation is considered [9]. This approach consists of dividing a sequence of incoming events into subsequences using an observation window O<sup>W</sup> (*t*) starting at time *t*, which contains a certain number of events equal to the size W of the aforesaid window. Another possible approach would have been the time-based windowing, which consists of segmenting the incoming sequence of events using a window of fixed temporal length [30]. This second approach is mostly used when analyzing data coming from sensors such as gyroscope or accelerometer, because there is a constant amount of data over time, while with binary sensors there could be moments without any sensor readings and some stall phases. Since the proposed system uses this type of binary sensors, a sensor-based window has been chosen. The result is that every sensor is treated as a feature and is associated with a particular activity based on its distribution probability to be in a sequence that is labeled with the name of that activity. This is done by implementing a Naïve Bayes Classifier (NBC) [41]. This type of classifier is based over the Bayes' theorem of independence between features, so that the model is constructed to find, for each activity A*j*, the probability *p* that given a certain number of features, i.e., the set of events, the class being observed is A*j*:

$$p(\mathcal{A}\_{\bar{j}}|\mathcal{E})\tag{1}$$

Accordingly, for each sequence of events observed in the observation window <sup>O</sup>*W*, a feature vector *<sup>F</sup><sup>W</sup> <sup>z</sup>* is computed similarly to the previous vector *F jk* of the training phase.

Finally, in the classification phase, the sequences of observed events are classified based on their probability to belong to a given activity. To this, the possible activities to be associated with the observed sequence of events are first filtered based on statistics Γ*j*: only the activities that are usually observed within a time window that includes the current observation window's starting time *t* are considered in the following step. Accordingly, the following condition needs to be fulfilled:

$$
\widetilde{\mathfrak{f}}\_{\mathcal{A}\_{\flat}}^{s} - \mathfrak{a} \cdot \sigma\_{\mathcal{A}\_{\flat}}^{s} \le t \le \widetilde{\mathfrak{f}}\_{\mathcal{A}\_{\flat}}^{s} + \mathfrak{a} \cdot \sigma\_{\mathcal{A}\_{\flat}}^{s} \tag{2}
$$

where *α* is a weighting factor that adjusts the time window to be considered. Finally, the cosine similarity *Sjz* between model vector *mj* and feature vector *f<sup>W</sup> <sup>z</sup>* is calculated with the equation below over the remaining activities, to evaluate which the more likely to be observed is:

$$\mathcal{S}\_{\bar{j}z} = \frac{\sum\_{i} (m\_{ji} \cdot f\_{zi}^{W})}{\sqrt{\sum\_{i} (m\_{ji})^2} \cdot \sqrt{\sum\_{i} (f\_{zi}^{W})^2}} \tag{3}$$

The observed activity is then labeled as the activity that corresponds to the highest cosine similarity.

#### **5. Reference Use Case**

The algorithm for modeling the activities and then discovering what the resident is doing is implemented and tested using the Aruba real-word dataset from the CASAS smart environment project of the Washington State University [37]. The data were collected from one smart apartment provided with motion sensors, contact sensors in the doors or cabinets and temperature sensors. Figure 4 shows the house plan of the apartment and the exact position of every sensor in the rooms. The description provided by the CASAS project did not give information about the specifics of the used sensor. Table 2 explains the number of sensors per type placed in each room. The values provided by motion and contact sensors are Boolean, whereas the ones provided by temperature sensors are numbers. There are two more sensors not listed in this table, one motion sensor and one contact sensor because they are not located in a specific room, but they are linked to the entrance of the house.

To correctly evaluate the correlation between the sets of events and the observed user's activities, without interference from other people, a dataset with only one resident living in the home was considered. The events decoded by these sensors are significant for recording elementary actions that people are performing, for example, door sensors are easily associated with opening and closing medical cabinet, food storage, or the entrance door, while with motion sensors it is possible to monitor the presence of the resident in one room and the proximity with a specific object or piece of furniture. The aggregation of these elementary actions defines one activity of interest.

**Figure 4.** House plant of the apartment for the Aruba dataset [37].


**Table 2.** Number of sensors per type in every room of the apartment.

The gathered data are presented with information about the date and time of every sensor event registered, the id of the activated sensor with its value and the beginning or end of each activity that is monitored. The dataset has the structure presented in Table 3. In the dataset, 10 different activities performed by the resident are noted. Table 4 shows the details of the number of times each activity appears in the dataset, as indicated by the user. The "Relax" activity is the one that the user has denoted as the activity occurring while staying in the living room, and it involves the set of sensors arranged in that room, as shown in the house map (Figure 4). The "Work" activity is the one performed in the office room and involving the specific group of sensors placed in that area. Lastly, the "Housekeeping" activity involves a great number of all the sensors of the house, due to the intrinsic dynamism of this type of activity. Sensors detect even the activities that are not registered, that correspond to "Other activity" with no label in the dataset. Since they cannot be classified accurately, this has been ignored in the proposed framework.


**Table 3.** Data extracted from the Aruba dataset [37].

**Table 4.** Activities and Statistics of the Aruba dataset.


Part of the dataset has been first used to train the system, whereas another part has been used to test it, by simulating the running phase. Most parts of the monitored activities have a long representation in terms of sequences of events. Only three of them generally are concluded in less than 15 sensor events, and they are the activities of "Enter Home", "Leave Home" and "Bed to Toilet Transition". The others are more variable and can last longer. Because of this variability, it is difficult to find a common size W to consider for the observation window OW. The size of this window is essential because for the activity recognition problem the algorithm evaluates the sensors' events that occur within this window. Based on that, the most probable activity being performed is found by searching for the minimum distance between the modeled feature vector and the new instances that are happening at the moment and that has to be classified, as defined in the previous Section.

#### **6. Performance Evaluation**

To evaluate the algorithm, an assessment of the classification accuracy that shows the percentage of correctly classified sequences of events for each class is used. The accuracy is obtained by observing the number of times an activity is correctly labeled, compared to all the occurrences of that activity in the dataset. Accordingly, the accuracy for activity *l* is expressed by

$$Accuracy\_I = \frac{T\_l}{T\_l + F\_l} \tag{4}$$

where *Tl* and *Fl* are respectively the number of times the activity was correctly and erroneously labeled.

The test is performed considering a training time of 2 months, one week of test data from the dataset and taking a window of W = 10 sensor events as the sequence that must be classified. The choice of the window is done by taking into account that there are three activities among those indicated in Table 4 that are shorter than the others. These activities often last for several events included between 5 and 15. Furthermore, it is preferable that activities are recognized as early as possible so that relevant management actions can be promptly started. Nevertheless, choosing a window that is too small could lead to errors when longer activities, which are the most numerous, are being performed, because the algorithm frequently confuses possible activities that are similar.

#### *6.1. Accuracy of the Activity Recognition Algorithm*

To first evaluate the accuracy of the recognition algorithm by itself, the system was first run excluding the filtering with respect to the Γ*<sup>j</sup>* statistics, as it is proposed in [9]. With the current choices, simulation results achieve an average accuracy of 70.43%. The results are presented in Figure 5a where the confusion matrix with the accuracy percentages of the classification for all the activities is presented.

(**a**) Classification without filter (**b**) Classification with filter

**Figure 5.** Confusion matrices of the activity recognition algorithm proposed in [9] (**a**), which corresponds to a classification without the proposed filter, and of the proposed algorithm that includes the statistics-based filter (**b**).

It is evident that almost half of the times there are errors discerning the "Meal Preparation" activity from the "Wash Dishes" activity. This is due to the fact that both activities are performed in the kitchen and involve the same sensors. The reason for the frequent errors between the "Relax" and "Housekeeping" activities is explainable remembering how this second activity was described in the previous Section 5. Due to the fact that in several cases the sensors involved are the ones placed in the living room, the algorithm often confuses the two activities.

Results improve consistently when adding the filter before the cosine similarity calculation, as it is shown in Figure 5b. Indeed, in this second case, the algorithm gave an average accuracy of 81.14% and it is possible to observe how all the activities are better recognized. This is even more evident in Figure 6, where the comparison of the percentage accuracy for every class between the two examined cases is presented.

**Figure 6.** Comparison between activity recognition accuracy in percentage with and without the filter. The classification without the filter corresponds to the algorithm in [9].

#### *6.2. Accuracy Results for Different Sizes of the Observation Window* O<sup>W</sup>

Table 5 shows the accuracy results for the different sizes W of the observation window with the indication of the 95% confidence interval for the three examined cases. Changing the size of the observation window to 15 and 20 consecutive sensor events, the performance of the algorithm are quite better, giving results of average accuracy equal to 87.02% and 84.09% respectively, making the choice of the size window of 15 sensor events the preferable one. This is due to the reason explicated before, for which longer activities are more easily observed with bigger windows while shorter activities are usually less frequent and affect the performance less significantly. Results also confirm that performance shrinks when considering observation windows that are longer than the shortest activities, i.e., those that have less than 15 events as mentioned at the beginning of this Section. The confidence interval is narrower for smaller windows than for bigger windows, because in these cases there is a greater number of samples, due to the fact that the same instance is divided into more parts compared to the use of a bigger window.


**Table 5.** Accuracy and confidential interval using different size of observation window *W*.

#### *6.3. Accuracy Results for Different Training Periods*

In Figure 7 the different overall accuracy values obtained with increasing training days, from 2 weeks to 3 months, are presented. The horizontal bars represent the 95% confidence interval with values that go from ±1.95% for 2 weeks of training to ±1.68% for 3 months of training.

**Figure 7.** Percentage overall accuracy for different training days.

The increase in the performance from 2 months to 3 months is not as prominent as the increasing noticed during fewer days of training. Analyzing this result, the conclusion is that it is probably better to have a training phase that does not last too long so that is not necessary to wait a long time period to get results. Instead, it is more important to evaluate the need for a new training phase after a while, in order to check if there have been changes in the habits of the user after some time. Furthermore, the various activities under consideration are not carried out with the same frequency every day or every week. This explains why, with a training set of a few days, the accuracy obtained is so stable and low, due to the fact that during that time some activities have only a few samples that could be used for making the models, making it harder to recognize those same activities later in the running phase.

#### *6.4. Prediction Results for Subsequent Activities*

The current activity can be used to predict the activities that are going to be performed in the next future. For this purpose, an evaluation of the probability of transition from an activity in a column to an activity in a row is reported in Table 6, whereas Table 7 shows the probabilities that the activity "C" is happening given the fact that activities "A" and "B" occurred. In the last table, only some of the results are presented, i.e., the ones with the greatest probability value for each combination of the first two activities. When two consecutive activities indicated in the tables turn out to be the same repeated activity, it is due to the fact that between these two instances of the same activities there was a gap in which the user performed some actions that were unknown and that were part of the "Other activity" group that was not considered during these tests as it was specified in Section 5. The last column reports the average time duration expressed in minutes of the sequence of the 3 activities indicated as activity "A", "B" and "C". This information could be useful for future works for making prediction and consideration of possible sequences of activities over a long time.


**Table 6.** Transition Probability of Activities.

**Table 7.** Conditional probability for a sequence of three activities.


#### **7. Conclusions and Future Works**

This paper focuses on the problem of users' activity recognition inside Smart Buildings to support BECM systems, by exploiting the acquisitions made by sensors deployed inside the building. To this aim, a system that analyses the user behavior to profile it and later perform early recognition and prediction of users' activities is presented. The proposed system is a sensor-based activity

recognition system that analyzing raw data from binary sensors arranged in an apartment can model users' activities and understand known behaviors. In addition to the data collected by sensors, other significant information, provided by statistics Γ*<sup>j</sup>* related to relevant external conditions that have been observed to be correlated (e.g., time of the day), are incorporated in the system, allowing the obtaining of better results than the simple analysis case of sensors data. The system is proved to enable activity recognition with an accuracy of more than 80% after only 10 sensor events are registered. Furthermore, the prediction of the following activity is achieved with an accuracy of about 60%.

It should be mentioned, however, that there are some limitations to the proposed method. As already highlighted explaining the obtained results, there are some problems related to the fixed size of the observation window, which hold the performance of the algorithm depending on the activity examined. A wider window is better to recognize long activities but, on the other hand, a longer time is required to start the assignment of an activity. Another important aspect concerns the period of the training phase. A longer time span gives better accuracy results, with more available samples to better model each activity, but it leads to exaggerated waiting times to first get any results. Lastly, this approach was thought and tested only for consecutive activities performed by a single user. Some changes must be evaluated for an extension to more than one resident and for recognize concurrent and interleaved activities.

Future works will be focused on improving the system accuracy by including statistics of other relevant external conditions that can be correlated, such as the weather, or the fact that it is a working day or not. Furthermore, other devices, such as smartphones, will be included to enlarge the number of sensors available for the analysis. The introduction of measurements coming from personal devices is expected to provide a more thorough insight into users' habits. Moreover, the system will be expanded to consider cases with more residents and to recognize different contemporary actions. Finally, the proposed system will be first tested using commercial software, to be later included on a real BECM system scenario, so that the actual convenience on energy cost savings and the quality of experience perceived by users can be assessed.

**Author Contributions:** Conceptualization, F.M., V.P.; validation: F.M.; writing–draft preparation, F.M., V.P., D.G.; writing–review and editing, F.M., V.P., D.G.

**Funding:** This work was partially supported by MIUR, within the Smart Cities framework (Project CagliariPort2020, ID: SCN 00281), by the Italian Ministry of Economic Development (MiSE, Project INSIEME, HORIZON 2020, PON 2014/2020 POS. 395), by "Fondazione di Sardegna" within the research project "SUM2GRIDS—Solutions by mUltidisciplinary approach for intelligent Monitoring and Management of power distribution GRIDS"—Convenzione triennale tra la Fondazione di Sardegna e gli Atenei Sardi, Regione Sardegna—L.R. 7/2007 annuity 2017—DGR 28/21 of 17 May 2015, within the project Design and Implementation of a Novel Hybrid Energy Storage System for Microgrids, which is funded by the Sardinian Regional Government (Regional Law no. 7, 7 August 2007) under the Grant Agreement no. 68 (Annuity 2015), and within the project LEAPH – anaLytics and data Enrichment plAtform for Pharma and pHarmacy owner, funded by the Sardinian Regional Government under the Sardinian POR FESR 2014-2020 (ID: RICERCA\_1C-124).

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Greedy Algorithm for Minimizing the Cost of Routing Power on a Digital Microgrid**

#### **Zhengqi Jiang 1, Vinit Sahasrabudhe 1, Ahmed Mohamed 2, Haim Grebel <sup>1</sup> and Roberto Rojas-Cessa 1,\***


Received: 22 June 2019; Accepted: 30 July 2019; Published: 9 August 2019

**Abstract:** In this paper, we propose the greedy smallest-cost-rate path first (GRASP) algorithm to route power from sources to loads in a digital microgrid (DMG). Routing of power from distributed energy resources (DERs) to loads of a DMG comprises matching loads to DERs and the selection of the smallest-cost-rate path from a load to its supplying DERs. In such a microgrid, one DER may supply power to one or many loads, and one or many DERs may supply the power requested by a load. Because the optimal method is NP-hard, GRASP addresses this high complexity by using heuristics to match sources and loads and to select the smallest-cost-rate paths in the DMG. We compare the cost achieved by GRASP and an optimal method based on integer linear programming on different IEEE test feeders and other test networks. The comparison shows the trade-offs between lowering complexity and achieving optimal-cost paths. The results show that the cost incurred by GRASP approaches that of the optimal solution by small margins. In the adopted networks, GRASP trades its lower complexity for up to 18% higher costs than those achieved by the optimal solution.

**Keywords:** digital microgrid; power grid; integer linear programming; routing energy; distributed energy resources; Dijkstra algorithm; integer linear programming

#### **1. Introduction**

The power grid infrastructure has been dramatically evolving from the traditional electrical grid to the smart grid to make it more manageable and efficient [1–4]. Moreover, the interest in integrating renewable resources into the grid to decrease greenhouse gas emissions into the atmosphere continues to increase. It is no surprise that the next-generation smart grid will use optimized electricity generation and transmission infrastructure, provide elastic electricity distribution, reduce peaks in power usage, and sense and react to power outages by incorporating the Internet and information technologies along with intelligent control algorithms [5,6].

Microgrids are becoming a functional and building block of the next-generation power grids [1]. They enable local control for the distribution of power and configurations tailored to satisfy local needs. These features are critical for improving the reliability and resilience of power distribution in the occurrence of failures or natural disasters [7]. The flexibility of microgrids to embrace one or multiple energy resources, configure one or multiple paths for redundancy, and the incorporation of communications to monitor and control the configuration of distribution legs enable such features [8]. Moreover, microgrids have access to the main power line in addition to the local energy resources. Therefore, microgrids can operate in grid-connected and islanded modes [9]. In a grid-connected mode, DERs may increase supply capacity, and in islanded mode, DERs may energize critical loads as needed during extenuating circumstances, such as in emergencies.

Several works propose the concept of the digital grid (DG) as means to fulfill the requirements of the next-generation power grid [10–19]. There are several approaches to the DG. They all share the use of internet protocol (IP) addresses to identify lines and grid equipment, such as power routers, to make them be able to communicate and participate in a the management and the distribution of power. The DG not only enables a robust operation and a seamless integration of renewable energy sources into a grid or microgrid, but it also leverages the incorporation of energy storage. The use of energy storage habilitates the aggregation of power generated by different sources and the supply of specific amounts of power to specific loads. The transmission of data in the Internet have inspired the development of the DG [20–23].

The DG supplies energy in a specific amount and delivers it to a unique destination, such as a piece of grid equipment, link interconnecting a substation, energy storage unit, or the final consuming load. In the DG, the energy carries both the IP address of the recipient and information of the amount of power delivered [14,24,25]. The transmission of the recipient address and the amount of granted power may be performed by either power lines that also transmit data or by a parallel data network. On a higher operating layer, a request-grant protocol is performed between loads and sources to schedule the delivery of power [14,24,25]. In such a protocol, loads issue a request for the needed amount of power, and one or several sources grant such requests and supply the requested amount of power. The DG concept is particularly suited for managing a microgrid, which becomes a digital microgrid (DMG). A service provider, acting as a DMG controller, may exercise centralized management for a DMG with numerous sources.

Digitization of power (i.e., the supply of finite and exact amounts of power to a specific destination or interconnection leg) makes it possible to route power from one or multiple energy sources to a load of a DMG. This routing of power, however, is difficult to realize on a traditional microgrid. In a DMG, the distribution network, formed by the interconnected transmission lines and cables, may yield more than one path for a source to supply power to a load. Shinabo [26] considers a similar feature, although the focus is the optimization of the power generation by multiple generators. Here, each generator is independent and interconnected by a power router [27]. However, this approach suffers from both limited scalability and uncontrollable power delivery because the power router only sets the connection or disconnection of lines rather than controlling the amount of power delivery.

However, the efficiency of a DMG may be affected by the costs of power distribution. In this paper, we are interested in minimizing the total cost, which we can model as power loss. We define routing here as the selection of the smallest-cost-rate links and as an optimization problem that calls for integer linear programming (ILP) as the solving approach [28]. The cost-rate is the rate that each unit of power costs to distribute (i.e., cost= cost rate \* power).

Because microgrids may count with multiple distributed energy resources (DERs), the routing problem includes a matching between a load and one or multiple DERs to satisfy the load's power demand. We show the optimal selected paths on different power networks, including different IEEE test buses, where each of them includes a different number of energy sources and loads. While this optimal method works well in the networks tested in this paper, we envision a power network with a vast number of loads for which the high complexity of ILP would make the selection of paths infeasible. Because the optimal ILP method is NP-hard, we propose the greedy smallest-cost-rate path first (GRASP), a heuristic approach but with lower complexity than ILP. GRASP matches the loads to sources and finds the smallest-cost-rate paths to route power in an iterative and greedy approach, and it uses the Dijkstra algorithm as the path searching mechanism [29]. We analyze the cost of GRASP and compare it with that of the optimal method. The discrepancy between the two approaches is the cost associated with the reduction of algorithm complexity. We show that GRASP approaches the optimal solution by 18% additional cost, but with a reduced computation complexity.

We organize this paper as follows. Section 2 introduces an optimal routing and source assignment method and GRASP. Section 3 presents the evaluation results of the optimal method and GRASP for different distribution networks, including IEEE test buses. This section also includes a practical example of the application of GRASP. Section 4 presents our conclusions.

#### **2. Methods: Optimal and GRASP**

In a conventional microgrid, the delivery of power takes place over a single shared bus, which directly connects the power source to each load. However, a large number of DERs in a DMG ensemble a distributed power generation. Examples of such DERs are rooftop solar panels, wind turbines, and conventional fossil-fuels-based generators. We model the set of interconnections formed by generators, grid nodes (or power routers), and loads as a graph where edges represent distribution links and nodes represent sources, power routers, or loads. Power routers may be multi-leg devices; they enable building a distribution network with multiple distribution paths. Therefore, a DMG may distribute power to loads through disjoint lines. Here, several loads may share a bus and the supply from many DERs but they may use diverse paths to interconnect to DERs. The DMG uses a power router to transfer energy from one link (or smaller grid) to another [30]. Such a power router has multiple inputs and outputs, or legs. Such routers can interconnect multiple segments and form an extended microgrid. The end topology of the microgrid would be a meshed network [11]. Figure 1 shows the architecture of a power router. This power router comprises a power path, which delivers energy to loads, and a data path, which transmits the information to monitor the demand and supply of energy. The inputs and outputs of the power router handle data and power. Here, energy is aggregated or partitioned by the router to supply amounts that match those requested.

**Figure 1.** Architecture of a power router. Reproduced from [14]. Publisher: IEEE.

In this paper, we focus on operating expenses (OPEX) of the use of power routers in a DMG. We consider losses on the transmission of power to be the associated costs. We consider two sources of losses: transmission lines and power routers. The lines may have a constant loss or one proportional to the energy transmitted. We adopt the latter one. The routers have power losses inherent to the transfer of energy from one input to an output, associated with the power router operation. Also, because the microgrid may be extended, there may be a cost associated with these extensions, as worst-case scenarios. For simplicity, and without losing generality, we associate these cost-rates to each link. Because some links and paths may incur costs as power flows through them, it is of interest to minimize the cost (or overhead) of power distribution in the DMG. We also address the practical issue of how links can be interconnected to enable multiple supply paths between sources and loads, and are also useful for fault-tolerant purposes. A router with many legs may enabled us to set up a complex network that includes redundant paths in the microgrid.

Moreover, the presence of multiple DERs allows the supply to a load by a single or multiple energy sources and an energy source to supply single or multiple loads. Thus, matching loads to sources was also a process included in the routing of energy in a DMG. The minimum cost routes include the assignment of sources to loads through minimum cost paths.

#### *2.1. Problem Definition*

The following example depicts the path and source assignment problem in the supply of power requests in a DMG. Figure 2a shows an example of a simple DMG. This network has two sources (nodes 1 and 2), two loads (nodes 3 and 4), and two power routers (nodes 5 and 6). The two power routers provide a path with links that share the power transmitted to both loads. Each source has a supply capacity indicated by the labels (0.5R units for node 1 and R units for node 2), and the number on each link indicates the cost per unit of transmitted power and the capacity (i.e., the amount of power that the link can carry) as (*p*, *q*), respectively. Each load requests an amount of power, as indicated by the labels beside the node.

In this example, node 3 requests R units of power, and node 4 requests 0.5R units. Figure 2b shows a path for transmitting power from node 2 (source) to node 3. This path has a low cost but leaves a load request unsatisfied as there is no other path to connect to node 1 and because the shared link (link from node 5 to node 6) has reached its capacity with the flow from node 2 to node 3. The figure shows the flow that supplies power to Node 3 as a solid blue line. Figures 2c,d shows the flow that supplies node 4 as a dashed red line. Figure 2c shows that the selected paths satisfy all requests, but the total cost is high (7R). Figure 2d shows a solution where the two loads are also satisfied, and the cost is the lowest. This solution also shows an example of a source supplying two loads (node 2) and a load receiving energy from two sources (node 3).

**Figure 2.** Example of problem definition of source and path selection on a small digital microgrid (DMG) model: (**a**) network model with paths, two sources and two loads; (**b**) Low cost path but unsatisfied requests. Total cost is 4R; (**c**) satisfied requests but high routing cost. Total cost is 7R; and (**d**) Lowest cost path and satisfied requests. Total cost is 6.5R.

As the example shows, the optimal routing of power in a DMG comprised of (a) the selection of the lowest-cost path to supply each load, (b) the selection of sources that supply energy according to their generation capacities, and (c) the consideration of the link transmission capacities.

#### *2.2. Optimal Solution: Integer Linear Programming Formulation*

We model the selection of paths with the smallest cost as an optimization problem to minimize the energy cost on the distribution of energy. The cost and losses in the distribution network of the DMG were linear and expressed with integer numbers, therefore, ILP is suitable for solving this optimization problem. We first modeled the DMG as a distribution network. Consider *G* = (*V*, *E*) be a graph with two special vertices: a source *s* and a sink *t*. Here, *u*, where *u* ∈ *V*, is a node in the set of nodes *V* and edge (*i*, *j*), such that (*i*, *j*) ∈ *E*, in the set of edges *E*. Every edge has a capacity *cap*(*i*, *j*) ≥ 0 and a cost or weight associated with that edge. The capacity and cost are two main features of the edges in the graph that indicate the maximum capacity of the link and cost per power unit (kW) transmitted, respectively. One or multiple flows may use one edge. A flow is a real-valued function on edges that occupies a portion of the power capacity of a link, and it is indicated as *f*(*i*, *j*). In practice, one can set the transmitted capacity as the amount of power (or current) that the link can carry and that the source can provide. The properties of a flow are as follows:


The set of ILP equations describe the network model, the selection of paths and the matching between loads and sources. This set of equations consists of an objective function and a set of constraints. The objective function is defined as follows:

$$\min \sum\_{(i,j)\in E} d\_{ij} \mathbf{x}\_{ij}. \tag{1}$$

Subject to:

$$0 \le s\_{i\prime} \forall i \in V\_s \tag{2}$$

$$0 \le t\_{\rangle}, \forall j \in \mathcal{V}\_t \tag{3}$$

$$\sum\_{j:(i,j)\in E} \mathbf{x}\_{ij} - \sum\_{j:(j,i)\in E} \mathbf{x}\_{ji} = s\_{i\prime} \,\forall i \in V\_s \tag{4}$$

$$\sum\_{i:\,(i,j)\in E} x\_{ij} - \sum\_{i:\,(j,i)\in E} x\_{ji} = t\_{j'} \,\forall j \in V\_t \tag{5}$$

$$\sum\_{j:(i,j)\in E} \mathbf{x}\_{ij} - \sum\_{j:(j,i)\in E} \mathbf{x}\_{ji} = \mathbf{0}\_{\prime} \forall i \in V - V\_{\prime} - V\_{t} \tag{6}$$

$$0 \le \mathbf{x}\_{ij} \le \text{cap}(i, j), \forall (i, j) \in E,\tag{7}$$

where *xij* is the amount of power the generator transmits to the load, *dij* is the cost per unit of power edge (*i*, *j*) carries, *si* and *tj* are the amount of power source *i* generates and sink *j* requests, respectively. *Vs* is a subset of *V* that contains all the source nodes. Similarly, *Vt* contains all the sink nodes. Here, (1) defines the cost minimization function, (2) and (3) are constraints corresponding to the amount of power generated by the sources and required by the sinks, (4) and (5) describe the conservation of energy flow of the power sources and sinks, respectively, (6) describes the conservation of energy flow of the transport nodes, and (7) shows the capacity constraint for each link. The complexity of the optimal method for the selection of the smallest cost path is *<sup>O</sup>*(*n*(2|*V*|+2)(|*V*|*a*)(|*V*|+1)(|*V*|+2)) where *<sup>n</sup>* is the largest number of legs of a network node, *a* = *max*{|(*i*, *j*)|, |*Vs*|, |*Vt*|}; the largest link cost, source capacity, and load demand, and |*V*| is the number of nodes in the network.

#### *2.3. Greedy Smallest-Cost-Rate Path First (GRASP) Algorithm*

Although the optimal method converges to an optimal solution for energy allocation in the presence of multiple energy sources and loads, the method is NP-hard, so that its complexity is prohibitively high for a network with a considerable number of nodes [31]. GRASP aims to match sources and loads and, in combination, find the smallest cost-rate paths in between, restricted to the source and link capacities, and the losses associated to links. By using the cost rates, GRASP aims to obtain a small routing cost. It should be noted that this smallest cost rate may not be the smallest overall cost but the use of this parameter reduces the complexity of finding the paths in comparison with the optimal method. The algorithm operates in an iterative and greedy approach, and its objective is also described by the objective function and constraints in Section 2.2. GRASP finds the *k* shortest cost-rate paths, one per each load-source pair and matches loads (demands) to as many sources as needed with an aggregate capacity generation that is able to satisfy the power demands of the loads.

There are other constraints in this problem: (a) the amount of power that the selected path can carry was set by either (1) the link with the smallest of power capacity; (2) the capacity of the power source and; (3) the amount of power the customer requests. (b) The selection of a path cannot select a congested link (i.e., a link carrying an amount of power equal to its capacity) for accommodating another request. Because different sources may partially fulfil the power requested by a load, the number of paths needed to supply a demanding load may be many. A routing algorithm must select the sources needed to satisfy a load demand and whose path to the source(s) incurs in the minimum transmission costs. This property is a significant difference between routing in a DMG and that in computer networks.

GRASP works as follows:


The pseudocode of this algorithm, presented in Algorithm 1, shows a greater level of detail.

Therefore, it is easy to see that the complexity of GRASP is the complexity of the smallest-cost-rate link search or *O*(|*E*|*log*|*V*|) [32] times the number of loads times the number of sources, or *O*(*Ns NL* |*E*|*log*|*V*|), where *Ns* is the number DERs, *NL* is the number of loads whose power requests are satisfied, and |*E*| is the number of links in the network.


#### **3. Results and Analysis**

Here, we evaluate the total cost incurred by the optimal method and GRASP. We show the application of both methods on two different networks. One was a simple network and the other was National Science Foundation (NSFNET). While NSFNET was not a power network, was is an example of a network with a different number of paths between two nodes (as in computer networks). It also shows how a power distribution bus may look after the adoption of power routers. We also consider different IEEE test buses, which found their use in various power tests. Some of these networks had a small number of power sources and loads, and others had a large number of both, thus presenting different scenarios. These networks offered different examples of the application of multi-leg nodes and a variety of paths.

#### *3.1. Example DMG Networks*

#### 3.1.1. A Simple Network

First, we consider a simple network, as Figures 3 and 4 show. These figures show the cost and capacity link assignment and source capacity for case 1. In the remainder of this paper, we show diagrams for case 1 of the considered networks. In this network, there are three power sources, nodes 1–3, depicted in green color. Each source node had a generation capacity of 35 kW, 45 kW, and 20 kW, respectively, resulting in a total generation of 100 kW available to the two loads; nodes 10 and 11, depicted in yellow color. Each load requested 50 kW. We indicate both the cost (e.g., power loss) per unit of power transmitted through it and *power capacity* of each link on the network.

**Figure 3.** Paths selected by the optimal method on a simple network.

We consider three test cases per network. We assigned an arbitrary cost and capacity per link in case 1. We modified the costs and link capacity of the links in case 2 to test the selection of a different path. We changed the amount of power capacity of sources and demand (requests) for driving the loads to match with different sources, leading to a varied selection of paths, in case 3.

**Figure 4.** Paths selected by greedy smallest-cost-rate path first (GRASP) on a simple network. (**a**) Paths selected by GRASP on a simple network indicating each power flow. (**b**) Paths selected by GRASP on a simple network indicating the aggregated power.

In general, the GRASP algorithm greedily selects the smallest-cost-rate path of all the source-load pairs in each iteration. As shown in line 4 of the pseudocode, the smallest-cost-rate paths of all source-load pairs are pre-calculated and cached in a list path. For each iteration of the loop, in line 13, path *p*, which is the smallest one amount the paths in path, will be selected and the source of this path is assigned to the load of this path. To show how the proposed GRASP algorithm works, we took the simple network as an example. As shown in Figure 4a, the blue arrows show the selected paths by using GRASP. First, GRASP calculates the smallest-cost-rate path of all the source-load pairs. It selects the path 1 → 5 → 9 → 11, where *a* → *b* means a link selected from node *a* to node *b*, that has a cost rate equal to 0.08. Then source 1 transmits 35 kW of power through as limited by its capacity. Then, it selects 2 → 5 → 9 → 11 as the second smallest-cost-rate path, to transmit 5 kW because this is the maximum power that could be transmitted form node 5 to 9 as 35 kW were already reserved on this link. GRASP keeps going in this way until all the power demands are satisfied.

A path from a source to a load can be represented by the path cost, which is the sum of link costs, and the path capacity, which is the capacity of the link with the smallest power capacity of all links comprising the path. The total cost is the sum of incurred flow costs on each link:

$$\mathbf{C}\_{t} = \sum\_{(i,j)\in E} d\_{ij} \mathbf{x}\_{ij\prime} \tag{8}$$

where *Ct* is the total cost, *dij* is the cost per unit of power transmitted on edge (*i*, *j*), and *xij* is the amount of power transmitted on edge (*i*, *j*).

Figures 3 and 4 also show the selected paths and the amounts of power granted over each selected path, as blue arrows, selected by the optimal method and GRASP, respectively. As Figure 4a shows, some power flows share some of the links in the selected paths as they are associated with small costs. Figure 4b shows the total power flowing on each link. The power on the links is the sum of the individual power per flow shown in Figure 4. We use this representation of results in the remainder of this paper. In the network, the total cost of routing power from all the sources to all the loads was 12.7 kW. However, the total cost achieved by GRASP was 13.5 kW or an increase of 6.3% over the optimal solution. We also performed a test for cases 2 and 3 (we do not show the figures for these two cases because of the limited space). In case 3, the power capacity of nodes 1–3 were 40 kW, 40 kW, and 20 kW, respectively, and the amount of power requested by nodes 10 and 11 was 70 kW and 30 kW, respectively. Figure 5 shows the comparisons of the total cost incurred by the optimal method and GRASP for each the three cases. As the figure shows, GRASP achieved a cost of 8.33% above the optimal costs as the worst case of the tested scenarios.

**Figure 5.** Cost comparison between the optimal solution and GRASP on the simple network.

We also performed route calculations on more complex networks, such as NSFNET [33], IEEE14 test [34] and IEEE39 test buses [35,36]. In this section, we present the total cost recurred by GRASP and the optimal method. Similarly, as in the simple network discussed above, we consider three different cases for each network.

#### 3.1.2. NSFNET

This network is the result of a program of coordinated and evolving projects sponsored by the National Science Foundation (NSF) that was initiated in 1985 to support and promote advanced networking among U.S. research and education institutions. Figures 6 and 7 show this network. This network had three sources and two loads, as in the simple network, but it also had a larger

number of links and a higher node degree. The increased number of links in this network enabled a wider selection of independent paths (i.e., paths with fewer or no shared link) of source-load pairs. Particularly to this network, the cost of each link is directly proportional to the geographical distance between two nodes (each node in NSFNET is a U.S. city) and the capacity of each link was set arbitrarily.

**Figure 6.** Optimal path selection and source-load matching on the NSFNET network.

**Figure 7.** Path selection and source-load matching by GRASP on the NSFNET network.

Figures 6 and 7 show the results of case 1 in NSFNET as resolved by the optimal method and GRASP, respectively. Figure 8 shows the overall cost results obtained by these two methods. These results indicate that GRASP fell 3.9% short of achieving the optimal solution by an increase of about 12.6 kW in case 2. Similarly to the simple network, case 2 of the NSFNET network used different link costs and link power capacities from case 1, while case 3 used different power source capacities and demands of power. Case 2 then is used to observe different route selection from case 1 as link costs are modified, and case 3 may lead to the use of different paths but triggered by a different distribution on power generation and demand. From these two cases, case 2 draws the larger discrepancy between these two methods. However, the differences between the optimal solution and GRASP in the NSFNET network are smaller than those achieved in the simple network. A possible cause of this improvement is the larger number of routing options provided by the larger number of links and node order in NSFNET.

Considering the complexity of NSFNET, the matching of sources and loads and the selected distribution paths in this network is also represented by a power routing table, as Tables 1 and 2. These tables correspond to the path selection presented in Figures 6 and 7, as obtained by the optimal method and GRASP, respectively. These tables show the amounts of power assigned to the different links. In each table, *xij* means the amount of power transmitted from node *i* to node *j*. A positive *xij* means that the power flows from node *i* to node *j* and a negative *xij* means that the power flows from node *j* to node *i* (opposite direction). Also, if *xij* = 0, there is no power flow between nodes *i* and *j*.

**Figure 8.** Cost comparison between the optimal method and GRASP on the NSFNET network. **Table 1.** Power assigned to the links of the NSFNET network by using the optimal method.


**Table 2.** Power assigned to the links of the NSFNET network by GRASP.


#### 3.1.3. IEEE14 bus

The IEEE14 bus system [34] represents a simple approximation of the American Electric Power system as of February 1962. This test bus has 14 buses, five generators (i.e., sources), and 11 loads.

However, to stress the performance of GRASP in this network, we assume only three loads from those in the bus. In this way, we increase the search space to match the loads to the sources as the number of sources is much larger than the number of loads. Figures 9 and 10 show the matched pairs and resulting paths by the two considered schemes. The data in the figures indicate the amount of power provided by the sources and requested by the loads and the cost and capacity of each link of case 1.

**Figure 9.** Optimal path selection and source-load matching on the IEEE14 bus.

**Figure 10.** Path selection and source-load matching performed by GRASP on the IEEE14 bus.

The results in this figures show that the total cost of the paths selected by GRASP is comparable to that achieved by the optimal method, as Figure 11 shows. The largest discrepancy in this scenario was 18.0%, which was the additional cost incurred by the selections GRASP make.

**Figure 11.** Comparison of power costs achieved by the optimal method and GRASP on the IEEE14 bus system.

#### 3.1.4. IEEE39 bus

This was a large bus, it had 39 buses, 46 transmission lines, 10 power sources, and 19 loads. This network had two nodes, nodes 31 and 39, represented in blue color, that performed a combined functional mode; each node generated and spent energy. These two nodes played both roles; a power source and a load, at the same time.

Figures 12 and 13 show the results obtained by the optimal method and GRASP on this test bus. With the larger number of sources and customers, the IEEE39 test bus showed a much more complex scenario than the three previous network topologies. Here, loads had a larger number of choices of sources to select, but the number of links and node order is equal to or smaller than those considered in NSFNET. It was expected that there were fewer choices for routing power in this test bus than in NSFNET and at the same time, more flows needed to be routed by the larger number of loads in this test bus.

**Figure 12.** Optimal path selection and source-load matching on the IEEE39 bus.

Figure 14 shows a comparison of the different cases tested on this network. This graph shows that GRASP achieves 9% as the most substantial additional cost in comparison to the optimal solution in case 3. Figure 15 shows the largest additional cost incurred by GRASP in comparison to those obtained by the optimal method for the three considered networks and cases.

**Figure 13.** Path selection and source-load matching performed by GRASP on the IEEE39 bus.

**Figure 14.** Comparison of power costs for the optimal solution and GRASP on the IEEE39 bus system.

**Figure 15.** Summary of the largest additional costs required by GRASP for the different considered distribution networks.

#### *3.2. An Example of a Practical Application of GRASP*

We used low-frequency power data collected from six houses from the Reference Energy Disaggregation Data Set (REDD) [37]. The data of each house consisted of power consumption information measured from different appliances, such as a refrigerator, lighting, microwave, etc. The data contained power consumption from real homes, for the whole house as well as for each circuit in the house (labeled by the main type of appliance on that circuit).

We used the power consumption profile of six different houses as loads and five different energy sources. Figure 16 shows the power profiles for 24 h of the six loads. The energy resources were three photovoltaic (solar) generators [38] and two wind turbine generators [39]. We used the same power profile for the solar and wind sources. Figure 17 shows the profile of the solar and wind generators.

**Figure 16.** Practical case with the consumption of six loads in IEEE14 bus.

**Figure 17.** Generation of power in practical sources for solar and wind turbine.

We used the IEEE 14 testbus to interconnect the sources and loads. In this practical example, we used one generator of constant power (power line) as node 1, two identical Solar generators as nodes 2 and 3, and two identical wind-turbine generators, as nodes 6 and 8. The power capacity of the constant power generator is 10 kW. The loads in this testbus were Nodes 12, 13, 14, 11, 5, and 4 as houses 1–6, respectively.

We performed power routing for this example, and the summary of per-hour distribution provides the results Figure 18 shows. These results show that the largest cost discrepancy is 0.1 kW at 6 p.m. Figure 19 shows the routing of power at 6 p.m. for the optimal method and GRASP.

**Figure 18.** Summary of the costs difference required by GRASP and optimal method in a practical example.

**Figure 19.** *Cont.*

**Figure 19.** Routing of power at 6 pm on IEEE14 testbus by (a) the optimal method and (b) GRASP on the practical example.

#### **4. Conclusions**

In a source-dense digital microgrid, the distribution of power consists of matching loads to the power sources that can satisfy their power request, and of the selection of the lowest-cost paths interconnecting them. Such paths transfer the requested amounts of power but with minimum total costs. We considered the costs incurred by a function of the power transferred through the link. To realize power routing, we first modelled the problem as an integer linear programming approach. While this method provides an optimal solution, its complexity is NP-hard. Therefore, the application of this optimal method on a digital microgrid with a large number of sources is prohibitive. As a solution, we have proposed the greedy smallest-cost-rate path first algorithm, or GRASP in short, to assign sources to loads and obtain the smallest cost routes between sources and loads. We tested GRASP on four different network topologies, including IEEE test buses. We compared the total cost obtained by the optimal method and GRASP. The results show that GRASP achieves comparable costs as does the optimal solution, with a discrepancy of up to 18 % on the considered scenarios. We also present a practical application of GRASP where we use the power consumption profile of six houses and the actual power generation profile of five sources on IEEE14 testbus and show the differences between the optimal method and GRASP. We showed that the largest discrepancy cost in this example is 0.1 kW. Our analysis also shows that GRASP achieves these additional costs in exchange for a computation complexity smaller than that of the optimal method.

**Author Contributions:** conceptualization, Z.J., V.S., and R.R.-C.; methodology, V.S. and Z.J.; investigation, Z.J., H.G. and R. R.-C.; writing, review and editing, Z.J., H.G., A.M., and R.R.-C.; project administration, R.R.-C.

**Funding:** This research was supported in part by the NSF Award CNS 1641033.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **An E**ff**ective Passive Islanding Detection Algorithm for Distributed Generations**

#### **Arash Abyaz 1, Habib Panahi 1, Reza Zamani 2, Hassan Haes Alhelou 3, Pierluigi Siano 4,\*, Miadreza Shafie-khah <sup>5</sup> and Mimmo Parente <sup>4</sup>**


Received: 1 July 2019; Accepted: 12 August 2019; Published: 16 August 2019

**Abstract:** Different issues will be raised and highlighted by emerging distributed generations (DGs) into modern power systems in which the islanding detection is the most important. In the islanding situation, a part of the system which consists of at least one DG, passive grid, and local load, becomes fully separated from the main grid. Several detection methods of islanding have been proposed in recent researches based on measured electrical parameters of the system. However, islanding detection based on local measurements suffers from the non-detection zone (NDZ) and undesirable detection during grid-connected events. This paper proposes a passive islanding detection algorithm for all types of DGs by appropriate combining the measured frequency, voltage, current, and phase angle and their rate of changes at the point of common coupling (PCC). The proposed algorithm detects the islanding situation, even with the exact zero power mismatches. Proposed algorithm discriminates between the islanding situation and non-islanding disturbances, such as short circuit faults, capacitor faults, and load switching in a proper time and without mal-operation. In addition, the performance of the proposed algorithm has been evaluated under different scenarios by performing the algorithm on the IEEE 13-bus distribution system.

**Keywords:** Islanding detection; distributed generation; modern power system; anti-islanding; protection

#### **1. Introduction**

Modern power systems are facing high complexity and sensitivity. One of the aspects of these grids and also future power system is the high penetration of distribution generations (DGs) into the network, which raises several issues. Islanding detection is a serious problem that emerges in this environment. In the islanding situation, one or more of the DGs and a portion of the network become separated from the main grid and continue to energize the islanded network. Due to several concerns, such as safety hazard, personal safety, and out-of-phase reclosing, islanding detection has vital importance. The islanding should be detected less than two seconds, according to the IEEE 1547 standard [1] and also with the respect of the re-closer operation. The re-closer constant time usually is less than 500 ms, hence the time of detection is a critical issue to discriminate islanding before the operation of the re-closer. Furthermore, some grid codes [2] imply that the DG should have the capability to continue to operate in the autonomous islanded mode. Therefore, the islanding is detected to the purpose of disconnecting the DG or preparing the DG to operate in the new control mode after the occurrence of the islanding.

In general, islanding detection methods can be classified into remote and local techniques. Remote techniques detect the islanding situation using communication links between the DG and the point of common coupling (PCC), which connects the network to the upstream grid [3–6]. In Reference [5], islanding can be detected using PMU measurements. Power line signaling utilized in Reference [6] to detect the islanding by transmitting and an insignificant signal from the DG bus and receiving it at the circuit breaker point. Remote methods detect the islanding with a negligible or even zero none detection zone (NDZ). In addition, these techniques are fast in the detection of the islanding situation and discriminating islanding from grid-connected disturbances. However, the main drawbacks of the remote methods are their high cost, risk of loss of communication links, and requiring to the backup protection problem.

Local methods can be divided into active and passive ones. A disturbance is added to the system deliberately to detect the islanding situation in the active methods. Consequently, active methods can be designed based on the type of DG. Hence, the active methods are used for specified types of DGs. In addition, in some active methods, the disturbance signal is not added to the system by means of the DG. For instance, inserted capacitor or condenser to intentionally disturbing the equality power mismatches to detect the islanding situation [7,8]. Some of the active methods for inverter-based generators are Sandia frequency shift [9], current harmonic injection [10], active frequency drift [11], voltage drifting technique [12], and positive feedback method [13]. Capability issue of current injection method is investigated in Reference [14] and a solution to solve this problem is also introduced based on considering only certain high frequency components as injected current. Likewise, in References [15–17], islanding detection methods for synchronous based generators are presented. In Reference [15], two positive feedback for active and reactive power control loops are presented to make the system unstable in the islanding situation. Afterwards, in Reference [17], new active and passive power control loops proposed to simultaneously improve the ride-through capability of the synchronous generator, as well as detecting islanding-situation. Integral controllers are added to the governor, and excitation system of the synchronous generator are added in Reference [16] to make the system marginally unstable in the islanding situation. In Reference [18], two modified active and reactive power control loops of synchronous DG is presented to make the system marginally unstable in islanding situation. Moreover, a signal processing technique is used to detect a feature which indicates the islanding conditions. A method to control the active and reactive power outputs of inverter-based DGs using two probabilistic phasing neural network controllers has been presented in Reference [19]. Active techniques can interfere with other devices which inject the disturbance to the system and also non-linear loads. Also, power quality degradation is another drawback of the active islanding methods.

Passive methods are based on monitoring the electrical quantities of the system. In this manner, the islanding is detected by violating the presets thresholds. Over/under voltage or frequency relays are the most common passive loss of mains detection [20]. Vector surge relays [21] are another type of passive islanding detection technique. Rate of change of frequency (ROCOF) is another common islanding detection method [22] which has large NDZ. In addition, setting the threshold of the ROCOF is a direct compromise between reducing the NDZ and maloperation in non-islanding events. Also, the voltage measurement-based methods, such as the rate of change of voltage (ROCOV) [8] have an unsuccessful operation to distinguish islanding condition, especially in the case of reactive power imbalance. Afterwards, the rate of change of phase angle difference (ROCPAD) is presented in Reference [23]. The maloperation in short circuit fault is the main drawback of these methods. Signal processing techniques are other passive methods that have been introduced recently [24–26]. The wavelet transform is used in References [25,26] to detect the islanding and minimize the NDZ. A passive method for inverter-based DGs introduced in Reference [27] by combining the help of close loop frequency control and a high frequency impedance detection method. In Reference [28], learning methods have been developed to extract the features which illustrate the difference between islanding and grid connected conditions. At first, different features are analyzed using signal processing methods

and then these features applied to a deep learning-based algorithm to classified islanding and grid connected disturbances. A voltage index has been used in Reference [29] to detect the islanding condition for large power imbalances, also for small power imbalances, the line current measured to disconnect some loads, due to transfer small imbalances to large power mismatches.

This paper proposes a passive islanding detection technique to eliminate the NDZ and improve the detection time. The measured values of the frequency, voltage, and current are used in the proposed technique to discriminate between islanding and non-islanding conditions. Also, a method to accurately calculate the phase angle difference between voltage and current is proposed in this paper for real-time applications. In addition, the proposed islanding detection technique can be used for all types of the DGs and without depending on the load type of the system. Performance of the proposed method can be retained when multiple DGs are connected to the islanded network. The maximum islanding detection time using the proposed islanding technique is 6 ms even though the power mismatch is zero. The detection time is reduced by increasing the power imbalances based on the proposed technique. The main contributions of this paper can be summarized as follows:


#### **2. Proposed Algorithm**

Many methods have been proposed for islanding detection. The main concept of most islanding detection methods is that some system parameters (like the voltage, frequency, etc.) change greatly in islanding mode, while they do not change significantly when the distribution system is grid connected. In some cases, existing methods do not recognize islanding mode and cause mal-operation. In addition, these methods may have maloperation in power system disturbances like short circuit faults. This paper, by combining different methods that are described below, provides a new algorithm that has a better operation for islanding detection.

One method for generator's islanding detection is the use of under frequency (UF) and over frequency (OF) relays (Figure 1). When islanding mode happens due to the power mismatch between load and generation, generator's speed changes according to Equation (1) [30].

$$2 \times \text{H} \times \frac{d\omega}{dt} = P\_m - P\_\varepsilon = P\_{mismatch} \tag{1}$$

where:


Frequency depends on speed; therefore, when speed changes, the frequency will also change. If the generation is greater than demand, the frequency will increase, and it will decrease if the generation is less than load. The islanding mode will be detected if the frequency drops below or increases above a defined setting for a time greater than a defined delay.

In this paper, frequency is calculated using the PCC voltage. First, the frequency range is estimated using the zero-crossing method. Then, using the least squares error (LSE) method, it is calculated accurately.

Consider the main signal as follows:

$$V(t) = V\_p \text{Sim}(\omega t + \theta),\tag{2}$$

where:

*V*(*t*) Measured voltage (V)

ω Angular speed (rad/s)

*VP* Peak of voltage (V)

θ Angle (rad)

*t* Time (s)

By sinusoidal expansion, Equation (2) becomes Equation (3).

$$V(t) = V\_p \text{Cost}\theta \text{Sim}(\omega t) + V\_p \text{Sim}\theta \text{Cas}(\omega t). \tag{3}$$

Equation (3) is expended by using Taylor's sinus and cosine expansions to obtain Equation (4).

$$V(t) = \text{Sim}(\omega t)V\_p\text{Co}\Theta + (\Delta\omega)t\text{Co}\text{os}(\omega t)V\_p\text{Co}\Theta + (\Delta\omega)tV\_p\text{Si}\text{in}(\omega t)V\_p\text{Si}\text{in}\Theta - \frac{1}{2}\sin(\omega t)(\Delta\omega)^2\text{Co}\Theta - \frac{t^2}{2}\cos(\omega t)(\Delta\omega)^2\text{Sin}\Theta$$

where Δω is equal to change in angular speed. There are six unknowns in Equation (4), so at least seven equations are needed to calculate the unknowns by LSE method. In other words, at least seven samples are required to estimate the frequency. To calculate angular speed variations, the known variable and unknown variable of Equation (4) are separated and written as Equation (5). Finally, by adding frequency changes to the base frequency obtained from zero crossings, the frequency is calculated according to Equation (7).

$$V(t) = S\_{11}X\_1 + S\_{12}X\_2 + S\_{13}X\_3 + S\_{14}X\_4 + S\_{15}X\_5 + S\_{16}X\_{6\prime} \tag{5}$$

Where,

$$\begin{array}{cccc} S\_{11} = V\_p \text{Cos} \theta & S\_{12} = tV\_p \text{Cos} \theta & S\_{13} = V\_p \text{Sin} \theta\\ S\_{14} = -tV\_p \text{Sin} \theta & S\_{15} = -\frac{t^2}{2} \text{Cos} \theta & S\_{16} = -\frac{t^2}{2} \text{Sin} \theta \end{array}$$

$$\begin{array}{cccc} X\_{11} = \text{Sin}(\omega t) & X\_{12} = (\Delta \omega) \text{Cos}(\omega t) & X\_{13} = \text{Cos}(\omega t) \\ X\_{14} = (\Delta \omega) \text{Sin}(\omega t) & X\_{15} = \text{Sin}(\omega t) \left(\Delta \omega\right)^2 & X\_{16} = \text{Cos}(\omega t) \left(\Delta \omega\right)^2 \end{array}$$

$$\Delta \omega = \frac{X\_2 + X\_4}{X\_1 + X\_3} \tag{6}$$

$$\Delta \omega$$

$$f = f\_{zcro} + \frac{\Delta \omega}{2\pi f\_0} \tag{7}$$

**Figure 1.** Islanding detection using under frequency (UF) and over frequency (OF) relays.

Rate of change of frequency (ROCOF) is an important technique for islanding detection. In low power mismatches, the frequency decreases slowly, and it takes a long time to exceed the relay set point. Therefore, islanding mode is detected faster using the ROCOF method. The ROCOF method is based on the theory that there is a mismatch between generation and load when islanding mode

happens. Immediately after islanding, because of the power mismatch, the frequency will change dynamically, which neglecting the governor action can be approximated by Equation (8).

$$\frac{df}{dt} = -\left(\frac{P\_L - P\_G}{2H \times S\_{GN}} \times f\_r\right) \tag{8}$$

where:

$$\text{( \newline \text{ \newline \text{ \textquotedblleft} Prequencary (Hz)})}$$


*fr* Rated frequency (Hz)

ROCOF is determined by deriving of the frequency which is obtained by the LSE method. The ROCOF is measured to compare it with the set threshold limit. If its value exceeds the pre-specified threshold limit, islanding is detected. The block diagram of the ROCOF method is shown in Figure 2.

**Figure 2.** Islanding detection using rate of change of frequency (ROCOF) method.

Another method of islanding detection is ROCPAD. The fundamental of this method is based on the rate of change of phase angle difference at DG end. The process starts with measuring voltage and current signal at the DG end, and then from the phase angle information of the voltage and current signals, the ROCPAD is computed for islanding detection. The five-sample LSE method is used to find the current and voltage angle.

The voltage signal is sampled at t = t0:

$$V(t\_0) = V\_p \text{Sim}(a\_0 t\_0 + \theta\_V),\tag{9}$$

where θ*<sup>v</sup>* is the phase angle of the voltage signal. Equation (9) can be written as Equation (10):

$$V(t\_0) = V\_p \text{Cost} \theta \nu \text{Sim}(a \rtimes t\_0) + V\_p \text{Sim} \theta \nu \text{Cost} (a \rtimes t\_0). \tag{10}$$

Consider Equation (10) as Equation (11):

$$V(t\_0) = a\_{01}X\_1 + a\_{02}X\_2. \tag{11}$$

where,

$$a\_{01} = \text{Sim}\left(a\_0 t\_0\right) a\_{02} = \text{Cov}\left(a\_0 t\_0\right).$$

$$X\_1 = V\_p \text{Cov}\theta\_V\\X\_2 = V\_p \text{Sim}\theta\_V.$$

It is done for four other voltage samples, and five equations are obtained similar to Equation (11). The matrix form of these equations is given by Equation (12),

$$[A]\_{\mathfrak{F}\times 2}.[X]\_{\mathfrak{Z}\times 1} = [V]\_{\mathfrak{F}\times 1}.\tag{12}$$

The unknown values of matrix X can be calculated from Equation (13),

$$[X]\_{2\times 1} = [A]\_{2\times 5}^{+}. \ [V]\_{5\times 1}.\tag{13}$$

$$\left[\left[A\right]\right]^{+} = \left[\left[A\right]^{T}.\left[A\right]\right]^{-1}.\left[A\right]^{T}.\tag{14}$$

Finally, voltage angle is calculated from Equation (15). It is also done for the current signal to get the current angle.

$$\theta = \tan^{-1} \left( \frac{X\_2}{X\_1} \right) \tag{15}$$

The ROCPAD can be calculated by deriving of the difference between the voltage and current angle. Islanding happened when its value exceeds a defined threshold. The block diagram of the ROCPAD method is shown in Figure 3.

**Figure 3.** Islanding detection using rate of change of phase angle difference (ROCPAD) method.

The rate of change of voltage, which is known as ROCOV index is one of the islanding detection methods. Most of DGs are required to operate at unit power factor. Thus, the lack of active power will be possible. Change of voltage is a function of reactive power, and as a result of the change of reactive power, the voltage will be changed. Capacitor banks may be the only available source of reactive power in the islanded network with DG operating in unit power factor. The amount of reactive power generated by capacitor banks is a function of the voltage. When the voltage changes, as a result of islanding, the reactive power which capacitor banks generate will also change, and it will change the voltage further. Thus, when reactive power mismatch has a large value, this causes a change in the voltage, so ROCOV value exceeds a preset threshold and islanding mode is detected.

Different methods of Islanding detection are described above. Each of these methods has problems that cause mal-operation or maloperation. The operation of OF/UF relays depends on the frequency. If the mismatch between the generation and load in the islanded network is low, then the frequency will have small changes. Therefore, it may not exceed the threshold limit and leads to false detection. On the other hand, when short circuit faults happen in power systems, the frequency decreases, due to a sudden decrease in load. In this situation, this method may operate incorrectly.

ROCOF method measures the rate of change of frequency, and once the rate of change of frequency exceeds the pre-determined setting, a trip signal is initiated. Similar to OF/UF relays, when active power imbalance in the islanding system is low, the frequency changes slightly. Therefore, the ROCOF method may become ineffective. In addition, network transient events may also cause changes in system frequency, resulting in the incorrect operation of the ROCOF method.

The reactive power imbalance is the detection criterion of the ROCOV method. If a load with a low power factor is disconnected from the network in a normal situation, the reactive power imbalance will have a large value. Therefore, the voltage changes and the ROCOV method may operate incorrectly. On the other hand, when reactive power imbalance in the islanding mode has small variation, the ROCOV value does not sense the threshold and cannot detect the islanding mode.

The ROCPAD method was proposed to solve the maloperation of the ROCOF, and OF/UF relays in the low power mismatch. This method even detects islanding mode with 0% power mismatch. However, the ROCPAD method may have mal-operation during a short circuit fault.

This paper, by combining ROCOF, ROCPAD, ROCOV, and OF/UF relays, presents an algorithm that solves the problems of previous methods. All previous methods can detect islanding mode. However, a more secure operation is achieved using the proposed algorithm. The proposed islanding detection algorithm is shown in Figure 4.

**Figure 4.** Proposed islanding detection algorithm.

In the proposed algorithm, the combination of ROCOF, ROCPAD, ROCOV and OF/UF relays is used to detect islanding mode fast. The operation time is an important parameter in protection algorithms. There is a direct relationship between the time of detection and the threshold value in classic methods. In order to fast operation, the ROCOV, ROCOF, OF/UF and ROCPAD methods with suitable setting are combined in the proposed method. As the operation speed of the classic algorithm increases, the probability of mal-operation will also increase. The proposed method can provide security and fast operation simultaneously.

The proposed algorithm can tackle the problems of classic algorithms by a logical combination of their output. The ROCPAD relay may have mal-operation in power system short circuit faults. While the mal-operation problem is solved by raising the threshold of this relay, the operation time increases. In the proposed method, two ROCPAD relays with different threshold are used for fast and secure operation. The first ROCPAD relay (ROCPAD-threshold1) has a lower threshold and is faster.

As mentioned above, the ROCOF, ROCPAD and OF/UF relays may have mal operation in power system short circuits. In short circuit faults, the voltage suddenly drops and reaches to a new value after a short time, but the voltage drops slowly in islanding mode. Using this rule, short circuit faults are separated from islanding modes. The ROCOV relay which measures the rate of change of voltage can separate short circuit faults from islanding mode. Therefore, the output of the ROCOV, ROCOF, ROCPAD and OF/UF relays are connected together by a logical AND gate to limit their operation in short circuit faults.

The second ROCPAD relay (ROCPAD-threshold2) has a higher threshold. Hence, it does not recognize short circuit faults as islanding mode. Its trip signal is connected to the output with a logical OR gate. If the first ROCPAD relay is limited, islanding will be detected through the second ROCPAD relay.

There is a direct relationship between the time of detection and the threshold of ROCOF relay. By reducing the threshold, as the islanding detection speed is increased, the probability of wrong detection increases. Therefore, Similar to ROCPAD relay, two ROCOF relays with different threshold are used in the proposed method. The first ROCOF relay (ROCOF-threshold1) has a lower threshold and operate faster. Noises and the transient faults change the frequency quickly. In this situation, the ROCOF value may exceed its threshold, while the frequency does not change significantly. In order to block the operation of the first ROCOF relay under transient faults and noises, the trip signal of it and OF/UF relay combined with a logical AND gate. The second ROCOF relay (ROCOF-threshold2) with higher threshold operates in parallel with second ROCPAD relay to enhance security.

The operation of the proposed algorithm in three general network modes is as follows:


#### **3. Simulation and Results**

In this paper, the voltage and current measured at the PCC by using capacitor voltage transformer (CVT) and current transformer (CT). Moreover, the least square error algorithm with five samples and seven samples, respectively, has been used to find the voltage and current phase angle and frequency of power system. The sampling rate is 1000 Hz. IEEE 13-bus distribution test feeder, shown in Figure 5, is used to simulate the proposed algorithm. The power system frequency is 60 Hz. Two gas generators are connected to the bus 12. The rated power of each generator is 2 megawatts. At bus 10, there is a 1 megawatt solar power plant connected to the network by an inverter. More details about generators are described in the Appendix A. There is also a 1.2 MVAr compensation capacitor on the network connected to the bus 13.

There is a fuse in the transmission line between the bus 6 and bus 7. The islanding mode happens due to burning the fuse. In all the cases in this paper, the islanding mode has happened, the transmission line between the bus 6 and bus 7 is disconnected. Simulations are performed using MATLAB and PSCAD software. In this paper, IEEE 13-bus distribution system is simulated by PSCAD software. The voltage and current obtained from the PSCAD at PCC location are used for islanding detection. The frequency estimation and the voltage and current phase angle extraction algorithms are developed in MATLAB software. Moreover, the proposed islanding detection method is simulated in MATLAB software. Specifically, the PSCAD is used to simulate the studied system and to generating the data considering different simulation scenarios as explained in the paper, while the MATLAB is used to implement the proposed islanding detection algorithm. In this case, the real power system behavior and the implementation of islanding detection are modeled the same to the reality. As a result, by the parallel execution of the PSCAD and MATLAB, the detection time and other implementation technical issues can be well-considered and determined.

**Figure 5.** IEEE 13-bus distribution test feeder.

#### *3.1. Islanding with High Power Mismatch*

In the first case study, the power mismatch is very high. At t = 3, islanding mode happens with 50% power mismatch. The output of the ROCPAD-1, ROCPAD-2, ROCOF-1, ROCOF-2, ROCOV, and OF/UF relays is shown in Figures 6–9.

**Figure 6.** The output of OF/UF relay in islanding mode with high power mismatch.

**Figure 7.** The output of ROCOF relay in islanding mode with high power mismatch.

**Figure 8.** The output of ROCPAD relay in islanding mode with high power mismatch.

**Figure 9.** The output of ROCOV relay in islanding mode with high power mismatch.

As shown in Figure 6, after the islanding happens, the frequency begins to increase immediately and after a short time exceeds the threshold. Figure 7 shows that the ROCOF value increases rapidly and exceeds both thresholds, due to high power mismatch. The ROCPAD value is shown in Figure 8. The phase angle difference between the voltage and current begin to increase rapidly, and ROCPAD value exceeds both thresholds in a short time. Figure 9 shows that the ROCOV value increases rapidly and exceeds its threshold of less than 10 ms.

The trip signal of relays and the proposed algorithm is shown in Figure 10. It can be seen that islanding has been quickly detected by all relays and the proposed algorithm.

**Figure 10.** Trip signal of relays and the proposed algorithm in islanding mode with high power mismatch.

#### *3.2. Short Circuit Fault*

In the second case study, a short circuit fault occurs in the middle of the transmission line between bus 11 and bus 12 at t = 3. The short circuit resistance is equal to 1.5 Ω. The output of the ROCPAD-1, ROCPAD-2, ROCOF-1, ROCOF-2, ROCOV and OF/UF relays is shown in Figures 11–14.

**Figure 11.** The output of OF/UF relay in short circuit fault.

**Figure 12.** The output of ROCOF relay in short circuit fault.

**Figure 13.** The output of ROCPAD relay in short circuit fault.

**Figure 14.** The output of ROCOV relay in short circuit fault.

As shown in Figure 11 in the short circuit fault, the frequency remains in an acceptable range. Therefore, the OF/UF relay does not work. Figures 12 and 13 show the ROCPAD and ROCOF value, respectively. It can be seen that the ROCPAD-threshold1 and ROCOF-threshold1 relays exceed their threshold. Thus, they operate incorrectly in short circuit faults. As mentioned above, after short circuit faults, the voltage drops immediately and reaches to a new value. On the other hand, the voltage drops with a lower slop in islanding modes. Thus, the ROCOV relay can separate islanding modes from non-islanding modes. As shown in Figure 14, the ROCOV value does not exceed the threshold and works correctly. Therefore, the output of the combination of these 4 relays with AND gate is zero. In this situation, if the ROCPAD and ROCOF are used, the generator is wrongly separated from the power grid.

Raising the ROCOF and ROCPAD relays threshold makes them secure in short circuit faults. As seen in Figures 12 and 13 the ROCOF and ROCPAD value in ROCOF-threshold2 and ROCPAD-threshold2 relays does not exceed the threshold. The trip signal of these relays is equal to zero. So, the output of these two relays with a logical OR gate is zero. Finally, the output signal of the proposed algorithm is equal to zero, and it does not work in short circuit faults.

The trip signal of relays and the proposed algorithm is shown in Figure 15. It can be seen that the trip signal of the proposed algorithm is zero. So, the proposed algorithm does not work in short circuit faults.

**Figure 15.** Trip signal of relays and the proposed algorithm in short circuit fault.

#### *3.3. Islanding with Zero Power Mismatch*

In the third case study, islanding mode happens with 0% power mismatch at t = 3. The trip signal of relays and the proposed algorithm is shown in Figure 16.

**Figure 16.** Trip signal of relays and the proposed algorithm in islanding mode with 0% power mismatch.

The frequency and voltage are not significantly low regarding power mismatch. The ROCOV, ROCOF, and OF/UF relays do not exceed their threshold. So, these relays do not detect islanding mode with low power mismatch. The ROCPAD relays can detect islanding mode even in 0% power mismatch. As seen in Figure 16, the trip signal of the ROCOV, ROCOF, and OF/UF is zero. The ROCPAD-threshold1 relay finds the islanding mode sooner than the ROCPAD-threshold2. The operation of ROCPAD-threshold1 is limited by ROCOV, ROCOF, and OF/UF relays. So, when the ROCPAD-thereshold2 exceeds its threshold, the proposed algorithm operates, and the generator will be separated from the power grid.

#### *3.4. Load Shedding*

In the fourth case study, at bus 7, a 1 MVA load with a power factor of 0.8 is disconnected from the power network at t = 3 s. The trip signal of relays and the proposed algorithm is shown in Figure 17.

**Figure 17.** Trip signal of relays and the proposed algorithm in load shedding.

As seen in Figure 17, load shedding does not affect the operation of the relays. Thus, the output of the proposed algorithm is zero.

#### *3.5. Capacitor Bank Switching*

In the fifth case study, at bus 13, a 1.2 MVAr capacitor is connected to the power network at t = 3 s. The trip signal of relays and the proposed algorithm is shown in Figure 18.

As seen in Figure 18, capacitor switching does not affect the operation of the relays. Thus, the output of the proposed algorithm is zero.

#### *3.6. Motor Switching*

In the sixth case study, at bus 11, a 2.2 MW motor is connected to the power network at t = 3 s. The trip signal of relays and the proposed algorithm is shown in Figure 19. As seen in Figure 19, due to starting, the inrush current the voltage of PCC decreases. Therefore, the ROCOV value exceeds the pre-defined threshold. In this situation the ROCOV relay has mal-operation, but the proposed method operates correctly.

The proposed algorithm has been investigated in different islanding mode, different short circuit faults with different fault resistance, different load connecting, load shedding and capacitor bank disconnecting. The results show that the proposed algorithm separates clearly islanding mode from the non-islanding mode. Table 1 shows the comparison of detection time between the proposed algorithm and classic algorithms at different power mismatches. It can be seen that the proposed algorithm detects the island condition in less than 6 ms time delay. Moreover, in all cases, the speed of the proposed algorithm is faster than the classic algorithm.

**Figure 18.** Trip signal of relays and the proposed algorithm in capacitor switching.


**Table 1.** Result of islanding detection using the proposed algorithm and classic algorithms for different active and reactive power mismatches.

**Figure 19.** Trip signal of relays and the proposed algorithm in motor switching.

Table 2 shows the mal-operation of the proposed algorithm and classic algorithms in non-islanding events. It can be seen that the proposed algorithm, unlike classic algorithms, has not maloperation on non-islanding events.


**Table 2.** Result of operation of the proposed method and classic methods in the non-islanding event.

#### **4. Conclusions**

Islanding detection becomes very important by employing distributed generation in the power network. If the islanding is not detected correctly, it will damage to the network equipment. This paper presents a new algorithm for islanding detection in distributed generation. The proposed algorithm by combining different islanding detection methods presents a new method that detects islanding mode very fast. It does not suffer from drawbacks, which may be seen in other methods. It works effectively under a wide range of power mismatch. The proposed algorithm does not require much data. It detects islanding modes only by using the voltage and current at DG's locations. Mal-operation during short circuit faults is one of the major problems of most of the islanding detection methods. The proposed algorithm can distinguish between short circuit faults and islanding mode. The operation of the proposed algorithm was investigated in various conditions, including different short circuit faults with different fault resistance, different load connection and rejection and disconnection of compensator capacitor. In all cases, the proposed algorithm works correctly.

**Author Contributions:** All authors worked on this manuscript together. All authors read and approved the final manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**


**Table A1.** The parameters of synchronous DG


#### **Table A2.** The parameters of inverter-based DG

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Earth Fault Location Using Negative Sequence Currents**

#### **Amir Farughian \*, Lauri Kumpulainen and Kimmo Kauhaniemi**

School of Technology and Innovations, University of Vaasa, 6500 Vaasa, Finland; lauri.kumpulainen@uva.fi (L.K.); kimmo.kauhaniemi@uva.fi (K.K.)

**\*** Correspondence: amir.farughian@uva.fi

Received: 27 August 2019; Accepted: 27 September 2019; Published: 30 September 2019

**Abstract:** In this paper, a new method for locating single-phase earth faults on non-effectively earthed medium voltage distribution networks is proposed. The method requires only current measurements and is based on the analysis of the negative sequence components of the currents measured at secondary substations along medium voltage (MV) distribution feeders. The theory behind the proposed method is discussed in depth. The proposed method is examined by simulations, which are carried out for different types of networks. The results validate the effectiveness of the method in locating single-phase earth faults. In addition, some aspects of practical implementation are discussed. A brief comparative analysis is conducted between the behaviors of negative and zero sequence currents along a faulty feeder. The results reveal a considerably higher stability level of the negative sequence current over that of the zero sequence current.

**Keywords:** sequence components; earth fault location; negative sequence current

#### **1. Introduction**

Feeder automation is one of the distinguishing features of smart grids. It aims at developing self-healing systems, able to locate faults and perform isolation and supply restoration automatically. Reliable fault location and indication is the key to this functionality. For short circuit faults, there are established methods for locating them, whereas, for locating earth faults, there is not a universally accepted, reliable and cost-effective method in the market for isolated neutral or compensated networks. However, a number of methods have been put forward to address this matter. A comprehensive review of the state-of-the-art methods for locating single-phase earth faults in medium voltage (MV) distribution networks is provided in Reference [1].

There are different types of fault location methods. Impedance-based methods [2–7] work based on calculating the apparent impedance seen when looking into the line from an end (measuring point) during the fault condition [8]. Since the line length between the measuring point and the fault location is proportional to the calculated impedance, by knowing the line impedance per unit length, the fault distance can be estimated. In Reference [2], the fault location and resistance are estimated using an iterative method. In Reference [3], a model developed for high-impedance faults consisting of two antiparallel diodes is proposed. A parameter estimation procedure is performed which uses the proposed model along with voltage and current signals to estimate the fault distance and other parameters. In Reference [4], the whole load of the feeder is modelled as one equivalent load tap located where the voltage drop due to the fault is at a maximum. Two assumptions are made i.e., once the load tap is assumed to be located before the fault point and once it is assumed to be after the fault. By solving the equations corresponding to those assumptions, the fault distance is estimated in unearthed networks. The method presented in Reference [5] scans the data obtained from power quality monitors and relays to estimate the fault distance in terms of the line impedance. The concept introduced in Reference [4] can be adapted for

compensated distribution networks as presented in Reference [6]. In Reference [7], a method for fault distance calculation for compensated networks is put forward which requires a fewer number of settings compared to the method presented in Reference [4]. Despite all the efforts, these impedance-based methods are subject to errors as the fault impedance is unknown. Moreover, to locate single-phase to ground faults, the phase to ground voltages must be available [4], which is cost prohibitive. Generally, these methods are more suitable in transmission lines.

Travelling-wave based methods [9–13] work based on the transient voltages and currents (impulses) which are generated at the fault location. These waves propagate away from the fault point in both directions towards the ends. As these types of waves travel at the speed of light, by measuring the wave and their reflection arrival times at either terminal of the line, the fault distance can be determined [8]. A concept to create travelling waves is proposed in Reference [9], which is based on earthing the neutral via a controlled thyristor that provides a short period of high fault current. The method in Reference [10] works based on frequency spectrum components of waves generated as a result of fault occurrence. Reference [11] is based on the wave velocity of zero and aerial mode components and their arriving timestamps. Some results from field tests are reported in Reference [12] indicating the limitations of applying travelling wave fault location on distribution networks. In Reference [13], the frequency information of travelling voltage waves is integrated with their time domain information using continuous wavelet transform to locate faults. However, the main problem when using travelling-wave based methods is that due to the large number of branches and laterals in distribution networks, the reflections captured at the beginning of the feeder, may come from different sources and not necessarily from the fault location only. This makes the use of traveling wave-based methods problematic in distribution networks.

Signal injection-based methods are another type of fault location method [14–16]. The idea is to inject a pulsating current into the neutral of the system following a fault occurrence. The injected signal is detectable only in the faulted feeder and from the injection point up to the fault point. Although, this method has been field tested, the results show that it can be used for fault resistances only up to 300 Ω (without voltage measurement). In addition, the use of this method is limited to only compensated networks where the injection device can be coupled with the Petersen coil.

In References [17–20], earth fault location using zero sequence components is proposed. Despite the simplicity this approach offers, it suffers from a shortcoming i.e., false indication in networks where the length of the faulty feeder is long and the fault occurs close to the beginning of the feeder. In addition, this method in the manner proposed in Reference [20], when applied to compensated networks, requires connection and disconnection of the resistor in parallel with the Petersen coil in the right time, which could cause complexity.

Despite all these efforts, earth fault location in isolated and compensated neutral distribution networks remains a problem for which there does not seem to be a widely adopted solution and thus it is worth studying. Fault indication and location is a more challenging task in isolated and compensated neutral networks because their fault currents are so low that the simple threshold detection will often lead to poor results. One solution is to utilize neutral voltage measurement, in addition to current measurement. However, this neutral voltage measurement is not usually preferred at the MV side due to additional costs [18]. Therefore, there is a need for an earth fault indication method which is based on current measurement only. On the other hand, recent developments in smart grids have enabled communication between different parts of the network. By taking advantage of this feature, better solutions can be achieved as measurements from different parts of the network can be utilized to locate the fault.

In this paper, a new method for identifying the faulted segment (the segment between two secondary substations on which fault has occurred) in case of a single-phase to ground fault is presented. It is purely based on current measurements and no voltage measurement is required, which can be considered as an advantage. The proposed method employs the negative sequence current components to locate the fault. First, In Section 2, the negative sequence current along a faulted feeder is analysed in depth to gain insights into its behaviour. Then, based on that analysis, the fault location procedure is established in Section 3. In Section 4, the effectiveness of the method is evaluated by means of simulation using various types of

medium voltage distribution networks. Lastly, some aspects of implementing the proposed method are discussed in Section 5.

#### **2. Theoretical Background of the Proposed Method**

In this section, to gain some insight into the proposed fault location method, the negative sequence currents on different locations on a feeder following a fault occurrence are discussed. The focus is on the difference between the negative sequence currents before and after the fault point.

To have a better understanding of the fault current, an MV distribution network with one healthy feeder and one with a single-phase to ground fault on it is shown in Figure 1. Earth capacitances and the flow of currents following a fault occurrence are shown in the figure. Fault current, especially in isolated neutral networks, is determined by network capacitances and mostly by the undamaged parts of the network.

**Figure 1.** A single phase to ground fault in an isolated neutral system [21].

For studying the fault location problem with more details, the MV distribution network shown in Figure 2 is used. In this paper, the style of single line diagrams is derived from Reference [22]. A typical MV network consists of several feeders but from a theoretical viewpoint, only the faulty feeder needs to be studied with more details and others can be represented by a simple electrical equivalent circuit This aspect is considered in the figure using the single block "Background network" representing all the healthy feeders. A single-phase to ground fault occurs on the feeder under study between two secondary substations H and J at point F with the fault resistance *RF*. The dashed lines indicate that the feeder can be of any length and can have any number of secondary substations. In Figure 2a, two arbitrary secondary substations are shown which indicates that the following analysis is valid for any point on the feeder. The goal here is to analyse the negative sequence current at point B, which is before the fault point and the negative sequence current after the fault location at point A. Note that points B and A are not electrically the same location. This is shown, for clarification, in Figure 2b.

**Figure 2.** *Cont*.

**Figure 2.** Single line diagram of a distribution network. (**a**) two arbitrary secondary substations are shown which indicates that the following analysis is valid for any point on the feeder. (**b**) B and A are not electrically the same location.

The sequence equivalent networks of the MV distribution system shown in Figure 2 can be obtained, as shown in Figures 3–5 (see also the notations list presented at the beginning of the paper). Note that line (to earth) capacitances only appear in the zero sequence network. In case of a single phase to ground fault, the sequence networks will be connected in series. The complete sequence network of the faulted network under study is shown in Figure 6. It can be simplified, as shown in Figure 7 where *ZG* and *ZJ* are equivalent impedances of the parts marked on Figure 6. The currents *IB* and *IA* are the fundamental frequency components of the negative sequence currents at points B (before the fault) and A (after the fault) in Figure 2, respectively.

**Figure 5.** Zero sequence circuit.

The values of load impedances are normally rather large compared to the system impedances, so that they have a negligible effect on the faulted phase current. Therefore, it is common practice

in fault analysis to neglect load impedances for shunt faults [22]. However, in the following, both scenarios are discussed, i.e., one where the effect of load is neglected and another where its effect is considered.

#### *2.1. With No Load*

Consider Figure 6. Shunt branches consist of transformers' reactance in series with load impedances. With no load assumption, the values of all the load impedances in the shunt branches are infinite and hence no current flows through the shunt branches in the positive and negative sequence networks. As a result, the negative sequence current measured at any point from beginning of the feeder to the fault point is constant and it is zero after the fault point. Therefore:

$$I\_B > I\_A \tag{1}$$

where, *IB* and *IA* are the magnitudes of the phasors *IB* and *IA*, respectively.

#### *2.2. With Load*

In Figure 6, when taking load impedances into account, the negative sequence current measured along the feeder from S to F is not constant anymore as some currents will flow through the shunt branches as well. In addition, similarly, due to the flow of current through shunt branches after the fault point, the negative sequence current after the fault point is not zero anymore.

In practice, the series impedances in Figure 6 are negligible when compared to the shunt impedances, with the exception of *Z*2*<sup>B</sup>* and *X*2*s*. Indeed, *Z*2*<sup>B</sup>* is comparatively low as it is the equivalent negative sequence impedance of several feeders i.e., the background network.

Now, consider Figure 7 which shows the simplified network where *ZG* and *ZJ* are the total negative sequence impedances before and after the fault point, respectively. In Figure 7, the following is valid.

$$
\tilde{I}\_B = \tilde{I}\_A + \tilde{I}\_F \tag{2}
$$

$$
\tilde{I}\_B = \frac{V\_f}{Z\_G + Z\_{2HF}} \tag{3}
$$

$$
\tilde{I}\_A = -\frac{V\_f}{Z\_{2\text{F}\text{I}} + Z\_{\text{I}}} \tag{4}
$$

As *ZG* represents the part of the network which consists of negligible series impedances and the shunt branch *Z*2*B*, it is a comparatively low impedance. In contrast, *ZJ* is a comparatively high impedance so that:

$$\text{abs}(Z\_{\text{G}} + Z\_{2\text{HF}}) < \text{abs}(Z\_{2\text{F}} + Z\_{\text{I}}) \tag{5}$$

Equations (3)–(5) yield:

$$I\_B > I\_A \tag{6}$$

Therefore, the negative sequence current before the fault point is higher than the negative sequence current after the fault point in both scenarios i.e., with and without considering the load.

**Figure 6.** Sequence networks and interconnections for a phase-to-ground fault in the system of Figure 2.

**Figure 7.** Simplified form of the network shown in Figure 6.

#### **3. Outline of the Implementation of the Method**

Considering the theory discussed in the previous section, when a single-phase-to-ground fault occurs, the negative sequence current exists and varies only a little in the section between the primary substation and the fault location whereas after the fault point, it is very low. This is simply illustrated in Figure 8.

**Figure 8.** Negative sequence current magnitude along the feeder following a single-phase to ground fault occurrence.

Therefore, the proposed new fault location procedure is structured as follows:


For example, in Figure 8, the fault information is sent to the control room from those measuring points at which the magnitude of the negative sequence current has exceeded the threshold i.e., measuring points 1–4. The last three secondary substations (points 5–7) send no signal to the control room and therefore the faulted segment is identified as the one linking measuring points 4 and 5.

#### **4. Simulation Results**

To evaluate the validity of the proposed method, a set of simulations was carried out with PSCAD™/EMTDC™. The effectiveness of the proposed method was studied with different fault resistances in each of the following network types:


The network shown in Figure 9 is an urban compensated medium voltage distribution network with compensation degree 0.95. The voltage level in the medium voltage side is 20 kV which is fed from a 110-kV supply network using a 40 MVA transformer. The length of the feeder under study is 5.4 km consisting of AHXAMK underground cables. The feeder under study has five secondary substations equipped with current measurements. The current measurement points 2–6 are the secondary substations along the feeder, points 1 is at the beginning of the feeder and H refers to the measuring point at the beginning of an adjacent healthy feeder. A single-phase-to-earth fault, with the resistance of 30 Ω, occurs at 2.5 km from the beginning of the feeder.

**Figure 9.** Cabled urban compensated network.

The magnitudes of the fundamental frequency phasors of the negative sequence currents with respect to time for faulty and healthy feeders are shown in Figure 10. The stabilized values obtained at *t* = 0.2 *s* are shown in a separate figure (Figure 11). Using the proposed method, one can readily determine the faulted section i.e., the one between secondary substations with measuring points 3 and 4. There is a considerable difference in the amplitude of the negative sequence currents before and

after the fault points so that the negative sequence current is practically negligible after the fault point. It is worth noticing that the negative sequence current of the healthy feeder behaves in a similar way as that of the points after the fault point.

In addition, the variations of the negative sequence currents along the feeder from the beginning of it up to the fault point are so insignificant that the first three graphs (corresponding to locations before the fault point) are almost matched. This highlights the strength of using the negative sequence current in locating faults.

**Figure 10.** Negative sequence current magnitude.

**Figure 11.** Negative sequence current magnitude.

To enable a comparison between negative and zero sequence currents, similar types of graphs are obtained for the zero sequence currents and shown in Figures 12 and 13. Contrary to the negative sequence current, the zero-sequence current varies considerably along the feeder and is not negligible after the fault point. This reveals the shortcoming of using the zero-sequence current over the negative sequence current.

**Figure 12.** Zero sequence current magnitude.

**Figure 13.** Zero sequence current magnitude.

To further investigate the validity and performance of the proposed method, five fault resistances are studied ranging from 0.01 Ω up to 1 kΩ. The changes in the negative sequence currents, caused by the fault, at the beginning of the healthy and faulty feeders and at each secondary substation of the faulty feeder are provided in Table 1. The choice of the threshold might vary from network to network. Developing a procedure to decide the proper threshold for a given network is not the focus of this paper. For simulations presented in this paper, the threshold is chosen to be 3 A for urban networks and 2 A for rural networks. For fault resistances between 0.01 Ω and 1000 Ω, it is clear from the table that for points 1 to 3, negative sequence currents exceed the threshold (3 A) whereas these values are below the threshold for points 4 to 6. Therefore, the faulted segment is the one between the second and third secondary substations i.e., between points 3 and 4.

**Table 1.** Negative sequence currents in urban compensated network.


The same procedure can be applied to isolated networks, as shown in Figure 14. The changes in the negative sequence currents, following the fault occurrence, are provided in Table 2. Using the table and the same 3 A threshold, the faulted segment can be determined successfully for all the 5 fault resistances.

**Figure 14.** Urban isolated network.


**Table 2.** Negative sequence currents in urban isolated network.

Now consider the rural compensated network shown in Figure 15. The first 5 km of the feeder under study consists of cables and the rest 45 km is overhead lines. The cable and overhead lines types modelled in the simulations are shown on the figure. A single-phase to ground occurs at 28 km from the beginning of the feeder. Similarly, using a threshold of 2 A and Table 3, the faulted segment is successfully identified, as the one between the second and third secondary substations (points 3 and 4).

**Figure 15.** Rural compensated network.

**Table 3.** Negative sequence currents in rural compensated network.


Lastly, in Figure 16, a rural isolated neutral network is shown. The change in the magnitude of the negative sequence currents for the points marked on the figure are shown in Table 4. In a similar fashion, the faulted segment is determined as the one between points 3 and 4 using the same 2 A threshold.

**Figure 16.** Rural isolated network.


**Table 4.** Negative sequence currents in rural isolated network.

#### **5. Discussion**

In practice, current measurements at each secondary substation can be made using Rogowski coils, which are also suitable for retrofit installations. In the tables presented in Section 4, there are current values which are very low. It should be noted that the accuracy of the negative sequence current depends on the accuracy of CTs or current sensors (Rogowski coils) used and that decides the lowest possible value of the negative sequence current that can be reliably measured.

When a fault occurs on a distribution feeder, the information on the fault conditions must be delivered to the control room where SCADA system or DMS processes the information to visualize the fault location to the operator or to generate an automatic FLIR (fault location, isolation and restoration) switching sequence. The proposed method requires some level of communication between secondary substations and the primary substation. However, as communication will be widespread in future distribution networks, the proposed method is worth considering.

In an ideal symmetrical network, the negative zero sequence measured at any point is zero under normal condition. However, in practice, there is some level of negative sequence current even under no fault condition, for example due to asymmetry of loads. Therefore, it is important to note that the proposed method is based on the change in the negative sequence current and not the negative sequence current alone.

Compared to the current injection method presented in Reference [14], the proposed method in this paper has a clear advantage as there is no need for any extra injecting device. In addition, the pulse injection method works only with compensated neutral networks where the injection device can be connected in parallel with the Petersen coil. For isolated neutral networks where there is no possibility for connecting the injection device, the pulse injection method cannot be used.

In case of a single-phase-to-ground fault, unlike the zero sequence current, the negative sequence current remains almost constant with negligible variations along the feeder as shown in Section 4. Moreover, in compensated networks, the difference in the zero sequence current level for points before and after the fault could be too low for finding suitable and reliable criteria for locating the fault. Therefore, there is a need to increase this deference level somehow. Reference [20] suggests to connect and control an additional resistor in parallel with the Petersen coil. However, that action may lead to more complexity when it comes to practical implementation. In addition, there is a risk that in case of a long feeder, the line to ground capacitances of that part of the feeder, which is located after the fault point, provide a considerable amount of zero sequence current so that the zero-sequence based method fails to operate. This highlights the advantage of using the negative sequence current over the zero sequence current in locating faults in distribution networks.

In summary, the main advantages the negative sequence current based method offers are:


#### **6. Conclusions**

Earth fault location on MV distribution networks using the negative sequence currents in the manner proposed in this paper appears to be promising. The proposed method is applicable to both compensated and isolated neutral networks. When using this method, only current measurements are required and not voltage measurements, which is an advantage. The theoretical reasoning to justify the proposed method was discussed in detail and simulation results obtained from different types of practical MV distribution networks confirmed the validity and performance of the method. Some further studies are still needed to address the aspects relating to practical implementation and for example the effects of decentralized compensation.

**Author Contributions:** Conceptualization, L.K.; Formal analysis, A.F.; Funding acquisition, K.K.; Investigation, A.F.; Project administration, K.K.; Software, A.F.; Supervision, K.K.; Writing—original draft, A.F.; Writing—review & editing, L.K. and K.K.

**Funding:** This work was supported by the European Regional Development Fund (ERDF) (Business Finland grant No. 4332/31/2014 and Council of Tampere Region grant No. A73094).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **List of Notations**

*Y*0*<sup>B</sup>*

*X*0*<sup>S</sup>* Zero sequence source reactance

*X*1*<sup>S</sup>* Positive sequence source reactance *X*2*<sup>S</sup>* Negative sequence source reactance *X*0*<sup>T</sup>* Zero sequence reactance of the main transformer *X*1*<sup>T</sup>* Positive sequence reactance of the main transformer *X*2*<sup>T</sup>* Negative sequence reactance of the main transformer *Z*0*GH* Zero sequence impedance of the feeder between nodes G and H *Z*1*GH* Positive sequence impedance of the feeder between nodes G and H *Z*2*GH* Negative sequence impedance of the feeder between nodes G and H *X*0*TH* Zero sequence reactance of transformer H *X*1*TH* Positive sequence reactance of transformer H *X*2*TH* Negative sequence reactance of transformer H *Z*0*LH* Zero sequence impedance of load H *Z*1*LH* Positive sequence impedance of load H *Z*2*LH* Negative sequence impedance of load H *Z*0*HF* Zero sequence impedance of the feeder between node H and the fault location *Z*1*HF* Positive sequence impedance of the feeder between node H and the fault location *Z*2*HF* Negative sequence impedance of the feeder between node H and the fault location *RF* Fault resistance *Z*0*FJ* Zero sequence impedance of the feeder between the fault location and node J *Z*1*FJ* Positive sequence impedance of the feeder between the fault location and node J *Z*2*FJ* Negative sequence impedance of the feeder between the fault location and node J *X*0*TJ* Zero sequence reactance of transformer J *X*1*TJ* Positive sequence reactance of transformer J *X*2*TJ* Negative sequence reactance of transformer J *Z*0*LJ* Zero sequence impedance of load J *Z*1*LJ* Positive sequence impedance of load J *Z*2*LJ* Negative sequence impedance of load J Equivalent zero sequence admittance representing the phase to earth capacitances of the background network *Z*1*<sup>B</sup>* Equivalent positive sequence impedance of the background network *Z*2*<sup>B</sup>* Equivalent negative sequence impedance of the background network *Y*0*GH* Zero sequence admittance representing phase to earth capacitances between nodes G and H *Y*0*HF* Zero sequence admittance representing phase to earth capacitances between nodes H and the fault location *Y*0*FJ* Zero sequence admittance representing phase to earth capacitances between the fault location and node J

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Design of the Smart Grid Architecture According to Fractal Principles and the Basics of Corresponding Market Structure**

#### **Albana ILO**

TU Wien—Institute of Energy Systems and Electrical Drives, 1040 Vienna, Austria; albana.ilo@tuwien.ac.at; Tel.: +43-1-58801-370114

Received: 7 October 2019; Accepted: 28 October 2019; Published: 31 October 2019

**Abstract:** Nowadays, there is a dramatic upsurge in the use of renewable energy resources, ICT and digitalization that requires more than the straightforward refinements of an established power system structure. New solutions are required to perform dynamic optimizations in real time, closed loops and so on, taking into account the high requirements on data privacy and cyber security. The *LINK*-paradigm was designed to meet these requirements. It was developed on the basis of the bottom-up method that can lead to misinterpretations or wrong conclusions. This work mainly deals with the verification of the authenticity and correctness of *LINK*. Fractal analysis is used to identify the unique and independent elements of smart grids required for the design of an architectural paradigm. The signature of the fractal structure, the so-called fractal pattern, is founded and referred to as electrical appliances (ElA). The latter has proven to be the key component of the architectural *LINK* paradigm. The definition of the *LINK* paradigm is finally validated: It consists of unique and independent elements that avoid misinterpretation or the need for any changes in its definition. Additionally, the fractal analysis indicates two fractal anomalies in the existing power system structure, while the fractal dimension calculation insinuates the highest complexity in the fractal level of electrical devices. The *LINK*-based holistic architecture is given a finishing touch. A compact presentation of the control chain strategy is provided that should facilitate its practical implementation. The basis for the harmonization of the market structure with the grid link arrangements is established. The processes of demand response and conservation voltage reduction are presented under the new findings.

**Keywords:** holistic power system architecture; smart grid; fractal design; fractal grid; *LINK*-paradigm; market design; local electricity market

#### **1. Introduction**

In the last 15 years many papers and books developed to investigate and design smart grids have been written. For the most part, concepts are introduced and various models are developed that lead to extremely ramified and complex schemas. Various projects are focused on certain parts of power systems without considering the integrity of smart grids. As a result, all efforts remain to the level of prototypes or isolated model regions.

So, what is the problem? "An extremely diverse and complex topic" is the answer. "... Each time we get into this logjam of too much trouble, too many problems, it is because the methods that we are using are just like the ones we have used before ... " teaches us Richard Feynman, the greatest [1]. This means the problem may be in is in the origin, in the definitions of the circulating smart grid concepts. Without elaborate concepts or paradigms characterized by unique and independent elements, serious design flaws and unclear operating procedures may result [2].

#### *1.1. Literature Review on the Popular Concepts of Smart Grids*

The traditional structure of power systems became questionable at the beginning of the liberalization; the market rules overcome the technical. Several blackouts and electricity crises were the consequence [3]. After the blackout that plagued the United States and Canada, on August 14th 2003, the term "smart grid" was introduced. It was mostly related to the increase of the transmission capacity and level of automation in the grid [4]. The electricity crisis in California (2000–2001) gave the rise of distributed generation (DG) a major boost [5]. The installation of the small photo voltaic (PV) plants on house roofs transformed consumers during the day into electricity producers. A new category of customers appeared: The prosumers. The meaning of the term smart grids have evolved and now stands for the modernization of power systems and meeting all the requirements of the time. According to [6], the scope of smart grids covers the entire power system (high, medium and low voltage level, HV, MV and LV, respectively) right down to the individual electrical appliances in customer plants (CP).

The need to design smart grids has brought onstage various smart grid concepts: Virtual power plants (VPP), microgrids, cellular approach (CA), web of cells (WoC).

The VPP concept as an aggregation of a number of distributed generators (DG) was first introduced in 2001 [7,8] to enable their participation in the electricity market. It soon turned out that an adaptation was necessary to take into account the voltage and frequency performance [9]. VPP optimization focuses on internal energy dispatch and external market participation [10]. The definition of VPP, in particular its technical aspect, continues to be in discussion [11–14]. The basic cell controller architecture [15] uses the VPP concepts. It has a layered control hierarchy that uses distributed agent technology. It defines three control modules: Local, regional and enterprise, requiring a tremendous amount of data to be exchanged between its modules [16]. Market researches focus on details of how they function, such as the development of different participation models in the balancing market [17] or the provision of secure market products [18], but not on the harmonization of the market structure with the power grid architecture.

Microgrids were developed with a focus on the technical issues of DG integration [19,20]. They lead to very complex architectures, setting the microgrids, nanogrids, etc., according to the Matryoshka-doll principle [21]. The detailed definition of microgrids is still under discussion in technical forums [22]. It focuses on distribution level and includes non-unique and undefined elements. "Load" is an important component of the microgrid definition that according to [23] has an ambiguous meaning. It can present the consumption of a device, a part of the grid or the total active and/or reactive power consumed by all devices of the power system. The scope of "the host power system" is undefined. All this complicates the definition of the size of microgrids and consequently the design of clear and standardized structures. The lack of a definitive definition gives to any project or initiative the opportunity to adapt the basic definition to their specific needs and develop individual processes. The latter cannot guarantee financial revenue streams, cannot be reliably audited, impedes pooling of multiple microgrid projects into a financial asset class and does not allow for wide-spread and attractive microgrid and distributed energy resource project deployment, [24].

CA can be seen as an attempt to further develop the microgrid concept taking into account the distribution and transmission grid. It is a self-controlled small microgrid, which is integrated with a modular smart grid ICT infrastructure [25]. The CA-architecture is based on five different energy cell types: Residential, commercial and industrial-plant, -area and –park. The distribution grid is considered through "connection corridors" while the transmission grid through "transmission corridors" [26,27].

WoC is the latest concept introduced into the landscape of the smart grid concepts. The cell definition has evolved very dynamically in a short time [28–32], is generic and does not appear unique in the literature [31,32]. In the cell-based decentralized control framework WoC, cells are defined as "non-overlapping topological subsets of a power system associated with a scale-independent operational responsibility" [31] without specifying the criteria for the practical definition of its size. The design of market structure is handled separately, not in combination with the technical architecture to enable their harmonization [32–34]. There are also attempts to design the architecture of smart grids without using a particular concept [35,36].

As a result all architectures based on the VPP, microgrids or their combinations with VPP, CA and WoC are very complex and hardly practicable. The harmonization of the technical smart grid architectures with the market structures does not seem to be in focus yet.

An ultra-large scale comprehensive control framework is needed to avoid the growing and unmanageable patchwork of grid implementations that are not sustainable on a large scale [37]. *LINK*-paradigm enables the holistic consideration of the power system (HV, MV and LV levels) and CP [38] and of the market [39]. However, it was developed on the basis of the bottom-up method that can lead to misinterpretations or wrong conclusions. The goal of this paper is the verification of the authenticity and correctness of *LINK*. The identification of unique and independent construction elements of smart grids is done using fractality principles.

#### *1.2. Literature Review on the Use of Fractal Analysis in Smart Grids*

In the past few years, fractality is increasingly being used to describe complex structures in nature [40]. It has recently been considered as a useful tool for analyzing, understanding and designing smart grids. According to [41] the scale-invariant properties of power grids are due to their design. Self-similarity is already identified in load flow characteristics [42,43], power demand behaviour [44], fault diagnostics [45,46], and protective relays [47,48]. Similarities have also been discovered in blackouts of interconnected power grids and their sub grids [49]. Fractal principles have been studied in grid structures and fractal dimensions have been calculated for different grids [50,51] to develop a cluster power system philosophy with cluster network structure [52,53]. Other authors have organized a smart grid based on a fractal model and have associated it with the holonic concept and holarchy [54]. Others consider fractals as an instrument that facilitates the realization of smart grids in a decentralized fashion [55]. However, the identification of a fractal pattern describing smart grids in its entirety has not yet been a research focus.

Section 2 gives an overview of the methodology used in this work. The in-depth investigation of the fractality features throughout the smart grid structure, including the calculation of fractal dimension, is presented in Section 3.1. The identified signature of the fractal structure, the so-called fractal pattern, in smart grids constitutes the foundation of the *LINK*-paradigm and of the holistic model, described in Section 3.2. Section 3.3 presents the holistic architecture and some applications of the *LINK*-solution, including a compact presentation of the control chain strategy and some applications as the harmonization of the technical architecture with the market structure, demand response, etc. Finally, some concluding remarks are summarized.

#### **2. Methodology**

Figure 1 shows the methodology used to design the smart grid architecture. In the first phase a superordinate connection between some well-known operation processes as "volt/var control in medium voltage grids" was established and the holistic technical model was predicted [38]. According to the prediction, the holistic technical model should include the entire power system, i.e., high, medium and low voltage levels and the customer plants. The bottom-up or induction method was initially used to derive the paradigm of smart grids "*LINK*" and the corresponding holistic architecture. The latter have shown the fractality feature of similar structural details.

The definition of the architectural paradigm of smart grids should be clear and contain unambiguous elements so that experts or users cannot make different interpretations [2]. Since misinterpretations or wrong decisions are preprogrammed in the bottom-up method, the top-down or the deduction method is used to verify the authenticity of the paradigm *LINK* and to exclude any suspicion or misinterpretation. After a detailed analysis, the fractal pattern of smart grids is revealed. The architectural *LINK*-paradigm is confirmed: Among others, it consists of the fractal pattern. Thereafter, the holistic technical and market model and the *LINK*-based holistic architecture have been given a finishing touch. The new architecture enables the description and realization of all processes required for the operation of smart grids, such as demand response, load/generation balance and so on.

**Figure 1.** The methodology used to design the smart grids architecture and the market structure.

#### **3. Results**

#### *3.1. Fractality in Smart Grids*

Fractals are constructs characterized by a cascade of similar structural details that appear on all levels [40]. The cascade is never-ending in theoretical mathematical cases. However, in our study, we are restricted to an ending cascade starting with the highest voltage level and ending with the customer device level, as shown in Table 1. This means that we consider the power system holistically, regardless of the voltage levels, equipment sizes and technologies used.


**Table 1.** The fractal structure of smart grids based on different fractal levels.

#### 3.1.1. Fractal Pattern

The fractal of the entire power system (including the customer devices sporadically connected to the mains) is conceived as five fractal levels shown in Table 1. To effectively meet the ever-increasing demand for electricity, power engineers have designed and developed power systems based on different voltage levels. The electrical appliances connected in each of the three voltage levels, high, medium and low (HV, MV and LV, respectively) grids, are contained in fractal levels 1, 2 and 3, respectively. The electrical appliances connected in the CP grid correspond to level 4. Level 5 takes into account all internal elements of the electrical devices that are connected by wires/conducting media. Power system appliances are categorized in two groups: The grid itself and active power appliances (APAs).

The grid is a construction that provides the connections amongst the points of power production and consumption. The relevant grid elements are wires, transformers and reactive power devices (RPDs) [56]. It is a fact that power, whether using alternating or direct current, needs wires to conduct it from the point of production to one of consumption. Thus, conductors are present in all fractal levels. In level 1, the high voltage level, they are called lines, while in the medium and low voltage levels, levels 2 and 3, they are more commonly referred to as feeders. High voltage lines are the longest, ranging from 100 to 2000 km. The world's longest power transmission line is the 600 kV direct current Rio Madeira transmission link in Brazil [57], with a length of 2385 km. Feeders in medium voltage grids are long, ranging from a few km to 70 km, while in low voltage grids they are short and range between several tens of meters and two kilometers. The CPs included in level 4 have short wires installed in building walls. They range between a few to several tens of meters.

Electrical devices connected to the customer plant grid have internal small/tiny wires that range from a few to several tens of centimeters. Thus, as the fractal level increases, the wire length decreases. The same trend also characterizes the wire diameters. The wire length and diameter repeatedly become smaller in the branches from level 1 upwards. Thus, wires are constructs that repeat themselves at progressively smaller scales.

The use of different voltage levels in the power system design dictated the introduction of transformers as connecting elements. Very large, large and medium sized transformers are used in HV, MV and LV (or fractal levels 1, 2 and 3), respectively. Transformers do not exist at the CP level, the gray elements in Table 1. Thus, the fractal structure presented in Table 1 shows a design anomaly in level 4. This anomaly may indicate further optimization potential of the power system structure. In the past, losses were one of the key optimization criteria in the power system design process. Only in the 1970s was it discovered that supplying a load using lower voltages reduces demand and energy consumption [58]. Since then, many studies based on conservation voltage reduction (CVR) have been conducted to assess the effectiveness of the implementation [59]. The current CVR implementation is based on controlling the voltage through distribution grids. Apart from the fact that its implementation in wide areas carries the risk of exposing some customers to unacceptable under-voltage conditions, the radial structures in distribution grids are subject to decreasing voltage profiles. The latter means that not all customers can be supplied with the lowest permitted voltage. This can be solved by installing 1:1 transformers with voltage control ability at the CP level, as indicated in the fractal structure given in Table 1. Their adjustment can guarantee the end customer's supply with the lowest possible voltage, thus ensuring maximum load reduction and energy savings. Additionally, this may lead to the liberalization of the voltage limits in super-ordinate grids, which in turn would make operation of the entire power system more efficient. However, basic economic studies are needed to substantiate this hypothesis. Tiny transformers exist in several customer devices (level 5). As in the wires case, the transformer size decreases with voltage level, meaning that it decreases with increasing fractal level. Transformers are thus constructs that repeat themselves at progressively smaller scales.

RPDs are mainly locally controlled (LC). They are widely used in the high and medium voltage grids to control the voltage. Very large and large sized RPDs are used in HV and MV (fractal levels 1 and 2), respectively. RPDs are not present in the LV level. The fractal structure presented in Table 1 thus shows a design anomaly in level 3, the gray element. In contrast, in level 4, the table is already filled, because with the integration of rooftop PVs, associated inverters are used to support the voltage control in low voltage grids by providing reactive power (almost purely inductive) [60]. Studies have shown that the use of customers' own devices to control the voltage on low voltage grids is technically not effective [61,62] and can lead to social problems [63]. The use of distribution system operators' (DSOs) reactive power sinks to control the voltage in LV grids is found to be more effective [61]. Thus, level 3 of the fractal structure table can be filled by DSOs and customer inverters should be used only

to compensate for their reactive power requirements [62]. Tiny RPDs exist in several customer devices (level 5). The RPD size follows the same trend as discussed above; it decreases with increasing fractal level. RPDs are constructs that repeat themselves at progressively smaller scales.

APAs are appliances that mainly produce, consume or store active power (rotating machines, photo voltaic, batteries, etc.). Very large sized APAs present in hydro, nuclear, etc., power plants as well as in pumped hydro power plants (storage) are connected to HV grids (level 1). With the emergence of distributed energy resources, large, medium and small sized APAs are already connected in MV and LV grids and in CPs (levels 2, 3 and 4), respectively. Moreover, several devices such as computers have tiny APAs in the form of batteries (level 5). The column of APAs in Table 1 is completely filled. The APA size follows the same trend as discussed above; it decreases with increasing fractal level. Thus, APAs are constructs that repeat themselves at progressively smaller scales.

As a result, although very heterogeneous, the construction of smart grids is self-similar. Figure 2 shows its fractal pattern. Both groups of self-similar constructs (grid and APA) are figuratively wrapped in an ellipse, see Figure 2a.

By definition, the fractal pattern of the power system consists of electrical appliances (ElA) designed for a pre-defined level.

Each ElA-ellipse denotes a separate chain link. It includes the grid of a specific level and all APAs connected to it. Figure 2b illustrates ElA in different levels: HV-ElA, MV-ElA, LV-ElA, CP-ElA and Dev-ElA for levels 1, 2, 3, 4 and 5, respectively. Figure 2c magnifies the CP- and Dev-ElA patterns.

**Figure 2.** Link chain fractal of power systems: (**a**) The fractal pattern electrical appliances (ElA); (**b**) ElA in different fractal levels; (**c**) magnification of levels 4 and 5.

#### 3.1.2. Fractal-Set of Smart Grids

From the bird's-eye view, all identified fractal levels are in the same territory, T. The overall pattern of level 1, HV-ElA, spreads over the entire T of a country or of a part of it, Figure 3. The patterns of the higher fractal levels propagate in smaller areas of the same territory T, repeating themselves many times. The repeating number of patterns in different fractal levels is derived from the relationship between different zones of power systems in real conditions and the number of appliances comprised in each of them. Table 2 shows a structural overview of the power system in real conditions and of its fractal.

Since the MV and LV grids have radial structures, the number of MV and LV zones are defined by the number of suppling substations (HV/MV) and distribution transformers (MV/LV), respectively.

A real European power system, operated by one transmission system operation (TSO), have about 300 supplying and 100,000 distribution substations to supply about 4 million customer plants.

**Figure 3.** The fractal of high voltage grids (HVGs).


**Table 2.** The fractal structure of smart grids based on different fractal levels.

It is assumed that each distribution substation supplies on average about 40 customer plants [64]. Additionally, it is assumed that in each customer plant exists on average for ten electrical devices. Figure 4a shows a graphical overview of the structure of smart grids. The graph increases exponentially as the number of customer plants, 4 million, increases drastically compared to the number of the HV zones, that is one. Thus, the fractal set cannot be assembled using only one scaling factor, because the number of zones in different levels is very different. In Table 2 are shown the different scaling factors used to define the number of ElAs for each fractal level and the derived numbers of ElAs used. The scaling factor for level 1 is set to one, while for levels 2 and 3 they are set to 50 and 5,555.56, respectively. The same scaling factor of 55,555.56 is used for levels 4 and 5 because the number of zones included in them differs only by a factor of ten. Figure 4b shows a graphical overview of the structure of the fractal set. The graph behaves similarly to the one of the power system in real conditions, Figure 4a. Using this approach, each of the six MV-ElAs touches the single HV-ElA at one point and tangentially to each other, Figure 5.

**Figure 4.** Overview of the structure of smart grids: (**a**) Real conditions; (**b**) fractal.

**Figure 5.** The fractal set of HV and medium voltage graphs (MVGs).

The latter symbolizes the tie switches (normally open points) connecting feeders of different MV-subsystems. For the 18 LV-ElAs is used the same principle as for MV-ElAs, Figure 6.

**Figure 6.** The fractal set of HV, MV and low voltage graphs (LVGs).

Four CP-ElAs are added into every LV-ElA, each touching only one point of the LV-ElA pattern (every customer plant has one electrical connection point with the LV grid and no connection with each other), Figure 7.

**Figure 7.** The fractal set of HV, MV and LVGs and customer plants (CP).

The same principle is used to add the ten Dev-ElA**s** into every CP-ElA, Figure 8. The fractal pattern ElA of different levels is assembled in a common fractal set; a link chain fractal limited to five interlinks.

**Figure 8.** The fractal set of smart grids.

The significance of the fractal set of smart grids is verified by the calculation of the fractal dimension. It contains information about the geometrical structure at multiple scales reflecting the complexity of the fractal set: The greater the complexity, the larger the fractal dimension.

#### 3.1.3. Fractal-Dimension

Fractals can be characterized by the fractal dimension *D*. For a fractal located in a bi-dimensional space, *D* is larger than 1 and less than 2,

$$1 < D < 2.\tag{1}$$

The larger *D*, the greater the fractal complexity.

*D* is calculated using the box-counting method, which is based on the counting of boxes (of side length (*s*) occupied by the fractal. The number of occupied boxes, *N*, increases by making the grid finer (*s* → 0). D is estimated as the exponent of a power low as follows:

$$D = \lim\_{s \to 0} \frac{\log N(s)}{\log \left(\frac{1}{s}\right)} \tag{2}$$

The software Fractalyse 2.4 [65] is used to calculate the fractal dimension of the entire power system, customer plants and electrical devices. The fractal dimension is obtained by using a logarithmic linear regression to fulfil the following objective function

$$\log N(s) = D \cdot \log \left(\frac{1}{s}\right) + \mathbb{C},\tag{3}$$

where *C* is the equation constant.

Figure 9 shows the separate fractal sets of smart grids in different levels. While the logarithm diagram for measuring D of each level is shown in Figure 10. The segments' slope increases steadily with increasing fractal level. Figure 11 shows the fractal dimension D for different separate fractal levels. D increases from 1184 for level 1 (HV) to 1272 and 1308 for the levels 2 and 3 (MV and LV), respectively. The D-increase from level 1 to level 2 is significant (ΔDL1→L2 = 0.088), which means that the complexity of smart grids at MV level increases dramatically compared to the HV level. Another significant D-increase, even bigger than the previous one, is observed between levels 3 (LV) and 4 (CP), (ΔDL3→L4 = 0.108). The complexity of smart grids at the CP level increases considerably compared to the LV level. The slope of the trend line of levels 4 and 5 (dashed line) is clearly bigger than in the case of the trend line for levels 2 and 3 (dotted line), 0.156 and 0.0347, respectively.

**Figure 9.** Smart grid's separate fractal sets in different levels: (**a**) Level 1; (**b**) level 2; (**c**) level 3; (**d**) level 4; (**e**) level 5.

The increase in complexity at the level of electrical device compared to the customer plant is stronger than that between the LV and MV levels.

Usually, the resources required to develop and purchase the different parts of a complex system are proportional to their complexity. The fractal dimension analysis indicates that, to realize smart grids, the highest global resources should be provided for the development and purchase of electrical devices, followed by a continuous reduction in resources for CPs, LV and MV up to the lowest resource allocation for the HV level.

**Figure 10.** Logarithmic diagram for measuring the fractal dimension of different separate fractal levels.

**Figure 11.** The fractal dimension D for different separate fractal levels.

#### *3.2. LINK-Paradigm and the Associated Holistic Model*

The realization of smart grids is related to a growing demand for more sensors, more communication, more computation and more control [66]; these represent an ever-growing cornucopia of new technologies [67,68]. Therefore, the ElA fractal pattern is combined with control schemes and interfaces to create the smart grid paradigm, Figure 12.

**Figure 12.** Overview of the smart grid paradigm.

By definition, the *LINK*-paradigm is a set of one or more ElAs, i.e., a grid part, storage device or a producer device, the controlling schema and the interface.

The *LINK*-paradigm is used as an instrument to design a *LINK*-based holistic architecture. It facilitates modelling of the entire power system from high to low voltage levels, including CPs and includes the description of all power system operation processes such as load-generation balance, voltage assessment, dynamic security, price and emergency driven demand response, etc. [69]. The *LINK*-paradigm is the fundament of the holistic, technical and market-related model of smart power systems with large distributed energy resource (DER) shares. Figure 13a shows the technical holistic model (the "energy supply chain net"), which extends to fractal level 4. It illustrates the links' compositions and their relative position in space, both horizontally and vertically. In the horizontal axis, the interconnected high voltage grids (HVGs) are arranged. They are actually owned and operated by TSOs. The medium and low voltage grids (MVGs and LVGs, respectively) and the customer plant grids (CPGs) are arranged on the vertical axis. MVGs and LVGs are actually owned and operated by the DSOs, while CPGs are operated by customers.

**Figure 13.** Overview of the holistic models: (**a**) Technical, the "energy supply chain net"; (**b**) market.

By definition, an "energy supply chain net" is a set of automated power grids intended for chain links (abbreviated as links), which fit into one another to establish a flexible and reliable electrical connection. Each individual link or link bundle operates autonomously and has contractual arrangements with other relevant boundary links or link bundles [70].

The holistic model associated with the energy market is derived from the technical holistic model, the "energy supply chain net", as shown in Figure 13b. The whole energy market consists of coupled market areas (balancing groups) at the horizontal and vertical axes. TSOs operate on the horizontal axis of the holistic market model, while DSOs operate on the vertical. Based on this model, not only TSOs but also DSOs will communicate directly with the whole market to ensure a congestion-free distribution grid operation and to take over the task of load-production balance. The owner of the

distributed energy resources as well as the prosumers (producers and consumers of electricity) may participate directly in the market or may do so via aggregators or local energy communities (LECs) [71]. The creation of the local retail markets (LRM) attracts the demand response bids and stimulates investment in the LEC areas (see Section 3.3.2).

#### *3.3. LINK-Based Holistic Architecture*

When the structure of electricity supply changes radically because of many DG units, each with the possibility of interfering with the system operation at all voltage levels, then a new architecture is needed to utilize this added flexibility so that power system operation can remain reliable. The new operational architecture of power systems should guarantee that it performs as expected by unifying their whole structure and by systematizing the execution of operational tasks. The basis for the design of the *LINK*-based holistic architecture is in its definition [39,72]:

Definition: A holistic power system architecture is an architecture in which all relevant components of the power system are merged into one single structure. These components could consist of the following: Electricity producer (regardless of technology or size, e.g., large power plants, DGs, etc.); electricity storage (regardless of technology or size, e.g., pumped power plants, batteries, etc.); electricity grid (regardless of voltage level, e.g., high, medium and low voltage grids); customer plants and electricity market.

Fleshing out the holistic model into a schematic holistic architecture requires specification of the main independent architecture elements, Figure 14. There exist three independent main power system components, i.e., producer, storage and grid, which are the three basic elements of the architecture and create the basis for the definition of the architecture elements. They are defined as follows:


To overcome data privacy and cyber security challenges, a distributed *LINK*-based holistic architecture is chosen [38]. Its key principle is to prohibit access to all resources by default, allowing access only to a minimum of data. Each link has its own operator. Based on the architecture elements, operators can be classified into three types:


(3) The grid link operator operates the grid regardless of voltage level (excluding the customer grid). Customers themselves are responsible for the operation of all their home elements.

**Figure 14.** Main elements of the *LINK*-based holistic architecture and the corresponding operators.

The technical/functional stage includes the conjunction of all three architecture elements in the four fractal levels (Figure 15a)) and a standardized architecture structure encompasses the four fractal levels. The different links communicate with each other via well-defined technical interfaces [38] "T". This is a detailed architecture level facilitating all technical/functional processes, which are needed to reliably and economically operate the decentralized power system. All processes can be described by using the unified modelling language (UML) diagrams (see Section 3.3.3). Additionally, it enables the step into the level of management systems, i.e., the application level, to develop concrete applications, (see details in the "Methods" section). The standardized structure of the technical/functional architecture allows the transition to the generalized architecture level (Figure 15b), where the four fractal levels are represented very compactly. The generalized architecture is the core of the *LINK*-based holistic architecture (Figure 15c). The market surrounds it and communicates with it through the market interfaces, "M", by exchanging aggregated meter readings, external schedules, etc. [39]. At the holistic level, the grid links of customer plants are removed from the generalized presentation because they are too small to participate directly in the whole market. They may participate in the common market through aggregators or through energy communities [71]. For the sake of privacy and cyber security, the technical interfaces "T" are designed apart from the market interfaces "M". The new designed holistic architecture facilitates two operating modes: (1) Autonomous—each individual link or link bundle operates independently by respecting the contractual arrangements with other relevant boundary links or link bundles. All links are connected together creating a large power system; (2) autarkic or self-sufficient—an optional operating mode, which may be applied in any link bundle that consists of at least one grid link and one producer link or storage link as long as it is self-sufficient and self-sustaining without any dependency on electricity imports. Restoration—an option of the autarkic operating mode, which may be applied after a blackout during the restoration process to supply electricity to at least the communication appliances. In order to successfully switch the operating mode from autonomous to autarkic, a familiar resynchronization process should be established, where each grid link has a secondary control on active and reactive power that supports the synchronization process. Depending on the properties of the links, the resynchronization with other links may be automatic or manual. However, resynchronization philosophies should be initially investigated to determine the most appropriate approach.

**Figure 15.** Different stages of the *LINK*-based architecture: (**a**) Technical/functional, (**b**) generalized and (**c**) holistic.

Operators are responsible for safe, reliable and efficient operation of the power grid at all times. The *LINK*-based holistic architecture supports them to achieve these goals, facilitating the execution of a set of functions and tasks that are encapsulated in various technical/functional operational processes, such as: Monitoring, generation-load balance [38], voltage assessment and var management [67], static and dynamic security [38], emergency [38], price driven demand response (see Section 3.3.3), etc. Some applications in the context of the holistic architecture are given in the following.

#### 3.3.1. Compact Representation of the Control-Chain Strategy

The standardized architecture of *LINK*-based architecture enables a compact and clear representation of the control strategy, which is designed as a chain net. According to Section 3.2, the latter is conceived on the X- and Y-axes. The different units of secondary control form a chain, which is expressed below in terms of volt/var control (VVC) chain.

The VVC chain in the X-axis extends to very HV and HV grids and is presented by the following generalized equation:

$$\begin{aligned} & \quad VV\mathbb{C}^{X-\text{axis}}\\ &= \left\{ VV\mathbb{S}\mathbb{C}^{HV} \Big( Volt\mathbb{P}\mathbb{C}^{HV}\_{\text{OLTC}^{V}}, var(\mathbb{P}\mathbb{C}or\mathbb{L}\mathbb{C})^{HV}\_{\text{RD}}, var\mathbb{P}\mathbb{C}^{HV}\_{\text{G}}, var\mathbb{N}\mathsf{g}\mathbb{b}Grid\mathbb{L}ink^{HV}\_{HV,MV,LV} \Big) \right\} \end{aligned} \tag{4}$$

where *VVSCHV* calculates:


While, the VVC chain in the Y-axis extends to HV, MV and LV grids and CP. It is presented by the following generalized equation:

$$\begin{array}{l} \text{VVC}^{\text{C}} \text{-axis} = \{\text{VVC}\mathbf{C}^{\text{HV}} \text{(}\text{VoltP}\mathbf{C}^{\text{HV}}\_{\text{OLTC}}, \text{var}(\text{PCorrL}\mathbf{C})^{\text{HV}}\_{\text{RD}}, \text{var}\mathbf{P}\mathbf{C}^{\text{HV}}\_{\text{G}}, \text{var}\mathbf{N}\mathbf{g}\text{bGrid}\text{L}\text{im}^{\text{HV}}\_{\text{HV}\text{MV}\mathbf{L}\mathbf{V}}\},\\ \text{VVC}\mathbf{C}^{\text{MV}} \text{(}\text{VoltP}\mathbf{C}^{\text{MV}}\_{\text{OLTC}}, \text{var}(\text{PCorrL}\mathbf{C})^{\text{MV}}\_{\text{RD}}, \text{var}\mathbf{P}\mathbf{C}^{\text{MV}}\_{\text{DG}}, \text{var}\mathbf{N}\mathbf{g}\text{bGrid}\text{L}\text{im}^{\text{MV}}\_{\text{MV}\mathbf{L}\mathbf{V}\mathbf{C}^{\text{LV}}}, \text{var}\mathbf{C}\mathbf{n}^{\text{MV}}\_{\text{HV}\mathbf{C}})\},\\ \text{VVC}\mathbf{C}^{\text{LV}} \text{(}\text{var}\mathbf{P}\mathbf{C}^{\text{LV}}\_{\text{RD}}, \text{var}(\text{PCorrL}\mathbf{C})^{\text{MV}}\_{\text{RD}}, \text{var}\mathbf{P}\mathbf{C}^{\text{LV}}\_{\text{DG}}, \text{var}\mathbf{N}\mathbf{g}\text{bGrid}\text{Link}^{\text{LV}}\_{\text{CP}}, \text{var}\mathbf{C}\mathbf{n}^{\text{LV}}\_{\text{MV}\mathbf{f}}),\\ \text{VVC}\mathbf{C}^{\text{CP}} \text{(}\text{var}\mathbf{P}\mathbf{C}^{\text{CP}}\_{\text{inv}}, \text{var}\mathbf{C}^{\text{mV}}\_{\text{$$

where *VVSCMV* calculates:


*VVSCLV* calculates:


*VVSCCP* calculates:

• the var set point *varPCCP inv* of the inverter connected on CP link grid; while respecting the var constraint on the border with LV *varCnsCP LV* Some practical applications of Equation (5) are shown in Section 3.3.4 and in reference [74].

#### 3.3.2. Smart Grid and Market Harmonized Structures

The large-scale DER integration increases the efficiency of the power grid and helps with its decarburization. In addition, their existence contributes to the LRM design, enabling their market-based control, which in turn stimulates investment in the LEC areas.

The market-based control must take into account the technical behaviour of the smart grids and support their safe reliable and resilient operation, while at the same time attracting the demand response bids. This is illustrated below by describing the harmonization and coordination of the market structure with the grid link arrangement, Figure 16. The implementation of LRMs requires the creation of balancing group areas (BGA) within the whole market, the day-a-head market that harmonizes with the grid links as follows:


**Figure 16.** Harmonization of the market structure with the grid link arrangement.

The electricity market is thus similar to the physical reality of power flows, as neighboring grids links are physically connected anyway and electricity always takes the shortest route from producer to consumer—across grid links and market area boundaries. The rules and methods for LRM operation should be in line with national policies and are therefore the subject of other studies.

The harmonized electricity market structure facilitates the introduction of the demand response process in large scale as follows.

#### 3.3.3. Demand Response Process

This is illustrated below by describing the process of price driven demand response. The activation of residential, commercial and small business sectors, which join the real-time pricing demand response through already concluded contracts, may be triggered at any time. Their degree of participation in the demand response process may be different depending on the time of the day, duration interval, price value, etc. In the case of a surplus of electricity in the market, the electricity price decreases. All market participants and market operators will be notified to allow them the possibility to act on time. Figure 17 shows the information flow during price driven demand response. This enables residential, commercial and small business sectors to perceive transparent energy prices and to contribute in the reliable and efficient operation of electric power systems.

**Figure 17.** Unified modelling language (UML) diagram of the process price driven demand response.

#### 3.3.4. Conservation Voltage Reduction in MV Level

The proof of the concept is performed in the frame of the industrial research project ZUQDE (central volt/var control in presence of DGs), Salzburg, Austria, based on one of the major entities of power systems, the voltage [75,76]. The scope was fractal level 2 (MV), where the secondary control link was realized by means of the volt/var control, Figure 18.

**Figure 18.** The implemented technical/functional *LINK*-based architecture in the ZUQDE project.

The applied algorithm calculated the set points by respecting the constraint, which was set to the HV/MV transformers by means of a constant cosΦ. The set points, reactive power Q and voltage, were sent to all four "run of river" distributed generators (DG) and to the feeder head bus bar, respectively, Figure 19. All relevant generators were upgraded along with the primary control, thus building up the producer component. All distributed transformers were modelled loads. As a result, the voltage in the Lungau region was automatically controlled and the grid was dynamically optimized in real time for more than one year.

$$VVC^{Y-\text{axis}} = \left\{VVSC^{MV}\left(Volt\mathbf{PC}^{MV}\_{OLTC'}\circ var\mathbf{PC}^{MV}\_{DG'}\circ \cos(\varphi)\mathbf{Cns}^{MV}\_{HV}\right)\right\}\tag{6}$$

**Figure 19.** The realized volt/var control scheme in management system level in ZUQDE project.

The voltages, power and current during the CVR switching process are shown in Figure 20. The snapshot is made on the SCADA system of the control center. The curves of voltage in one of the 30 kV bus bars and of the active, measured on the supplying transformer level, are highlighted.

Significant in the consideration are the curves shown in blue, green and yellow. The blue curve shows the voltage set point on the 30 kV side of the supplying transformer. At 10:18 a.m., the VVC set in the CVR mode calculated new set points for the 30 kV bus voltage calculated new set points for the 30 kV bus voltage and the reactive power of DGs. The voltage set point changed from 31.7 kV to 30.5 kV. This corresponds to an adjustment of the on-load tap changer by two steps. The green curve shows the measured voltage on the 30 kV bus bar. Due to the reaction time of the OLTC, the new voltage set point was reached after about one minute. The yellow curve shows the active power, which

decreases from 19.9 MW to 19.0 MW. This corresponds to a reduction of the active power by 4.7%. The light blue or turquoise curve shows that the voltage reduction results in a slight increase of the current.

**Figure 20.** Voltage and active power course during the execution of the conservation voltage reduction (CVR) process in the ZUQDE project [76].

#### **4. Conclusions**

The definition of the architectural paradigm *LINK* is finally validated by the fractal analysis: It consists of unique and independent elements that avoid misinterpretation or the need for any changes in its definition. Its use allows the convergence of the results of various scientific works and the implementation of smart grids on a large scale.

The fractal analysis is a suitable instrument for characterizing the heterogeneity of the smart grid structure, as it is capable of characterizing the spatial differences within the smart grid. Based on fractal analysis, two fractal anomalies are identified in the existing power system structure. Firstly, 1:1 transformers with voltage control capability are missing in fractal level 4. Their presence could contribute to increased energy savings, which is a key objective under actual climatic conditions. However, fundamental studies are needed to assess this change in the current power grid structure. Secondly, reactive devices are missing in fractal level 3. Studies have shown that eliminating this anomaly produces essential benefits. The fractal dimension analysis indicates that to realize smart grids, the highest global resources should be provided for the development and purchase of electrical devices, followed by a continuous reduction in resources for CPs, LV and MV up to the lowest resource allocation for the HV level.

The *LINK*-based holistic architecture reflects the principles of the fractal. Repeating the same structures in ever smaller dimensions greatly improves the practical relevance of the *LINK* solution. The compact and clear representation of the control strategy, which is designed as a chain net, is straightforward und understandable for every electrical engineer. The *LINK*-based holistic architecture considers the entire power system from high, medium and low voltage levels, including customer plants and the market. It facilitates the description of all power system operation processes such as load-generation balance, voltage assessment, dynamic security processes, etc. It enables the large-scale DER implementation, respecting data privacy and minimizing cyber security risks. The new market structure harmonizes with the grid link arrangement, facilitating the demand response process. The CVR process is also easy to implement.

The realization of significant processes required for the operation of smart grids, such as demand response and load/generation balance need further investigation to show the effectiveness using the proposed architectural *LINK*-paradigm.

**Funding:** This research received no external funding.

**Acknowledgments:** Open Access Funding by TU Wien.

**Conflicts of Interest:** The author declare no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Phase Synchronization Stability of Non-Homogeneous Low-Voltage Distribution Networks with Large-Scale Distributed Generations**

#### **Shi Chen 1, Hong Zhou 1, Jingang Lai 2, Yiwei Zhou 3,\*, and Chang Yu <sup>1</sup>**


Received: 9 February 2020; Accepted: 3 March 2020; Published: 9 March 2020

**Abstract:** The ideal distributed network composed of distributed generations (DGs) has unweighted and undirected interactions which omit the impact of the power grid structure and actual demand. Apparently, the coupling relationship between DGs, which is determined by line impedance, node voltage, and droop coefficient, is generally non-homogeneous. Motivated by this, this paper investigates the phase synchronization of an islanded network with large-scale DGs in a non-homogeneous condition. Furthermore, we explicitly deduce the critical coupling strength formula for different weighting cases via the synchronization condition. On this basis, three cases of Gaussian distribution, power-law distribution, and frequency-weighted distribution are analyzed. A synthetical analysis is also presented, which helps to identify the order parameter. Finally, this paper employs the numerical simulation methods to test the effectiveness of the critical coupling strength formula and the superiority over the power-law distribution.

**Keywords:** large-scale; distributed generations; low-voltage active distribution network; islanded mode; non-homogeneous model; synchronization; stability

#### **1. Introduction**

Renewable energy sources have a broad prospect in power systems these days—more and more distributed generations (DGs) are integrated into the modern electricity grid. The distributed generation has the advantages of being environmentally-friendly and sustainable. In addition, two weaknesses include low output power per unit and energy output in remote areas that can not be ignored. For one thing, the low output power per unit of solar power and wind power results in a low capacity for one DG unit. This is why the distribution network requires a large number of installed DGs to provide effective output power. For another thing, the remote distribution leads to power systems usually working in the islanded operation. Different from the normal generators, there is no rotor in the DGs such as batteries and photovoltaic units. That is to say, the frequency stability is easy to be broken in the islanded microgrid composed of DG units without the rotating mass [1–3]. Hence, we will address the issue of phase and frequency stability in the islanded system.

Previous studies have based their criteria for selection on the local control of a few DGs in the ideal model [4–6]. Reference [7] shows the change of various parameters of generators in a large scale photovoltaic system and presents two different operations to research the voltage stability in an IEEE 30-bus model with a Large-Scale photovoltaic module. Indeed, the large-scale photovoltaic module can be considered as one high-capacity DG in the case of IEEE 30-bus, and no method of improving stability was given in [7]. Conventionally, when the number of composed DG increases, the active

low-voltage islanded network will face the issue of high-order, multi-dimensional, and strong coupling. In brief, with a sharp increase in the number of DGs, the system will encounter the problem named "Curse of Dimensionality". Therefore, it is difficult for previous solutions that consider a few DGs to address relevant problems in current research effectively.

We note that recently there has been some advanced complex network theory that provides the stability analysis of large-scale network systems. It has wide applications in cyber science [8], medical science [9], social science [10], and so on. Nevertheless, there are a few related studies on the frequency stability of an active distribution network composed of DG units.

In addition, previous research mainly assumes that the networks are unweighted and undirected; this means that the coupling relationship between nodes has no direction, and the strength of all couplings is consistent. However, this assumption mismatches the characteristics of the real network of many cases. For instance, in an active distribution network, the coupling relationship between DGs which takes droop control is determined by the node voltage, line impedance, and droop coefficient, where the line impedance between DGs is generally determined by the distance between DGs. Apparently, the assumption that all DGs in the active distribution network is equidistant is undoubtedly harsh and impractical. On the contrary, the assumption that the coupling relationship between DGs is non-homogeneous is more realistic. Fortunately, complex network theory has made a difference in weighted networks. Reference [11] investigates the system of Kuramoto oscillators whose weight is related with a degree in the asymmetrical situation and comes to one conclusion that the critical coupling strength has a positive correlation with the degree exponent. Reference [12] indicates that the degree sequence contributes to the ability to synchronize in the mean field. Reference [13] shows that the degree distribution of the fixed system has a negative degree correlation between its nodes, which improves its synchronizability. Besides the weighted methods with the degree, some other papers pay attention to the weighted methods with frequency. Reference [14] introduces a weighting method which cares about the relationship between frequency and link to achieve the phase synchronization. Linear and nonlinear weighting procedures are applied in mean and heterogeneous systems to verify the method. Frequencies of symmetrical oscillators are studied in the bipolar model built by [15], while the relationship between frequency and coupling strength is analyzed. Reference [16] comes to the conclusion that asymmetric weights that depend on loads make contributions to the stability of the whole system. Overall, weighted methods that are based on the inherent qualities of the inverter nodes [17,18], i.e., node degrees [19], node loads [20], and network's global properties in the connection links have been discussed widely to enhance the network synchronization in recent studies.

However, the general impact of weighted coupling on network synchronization is not studied comprehensively, especially the weighted method. These weighted procedures above, unfortunately, are designed only concerned with the inherent attributes of the network, i.e., the size of the nodes degree or the edges between vertices, which limits their applications. Thus, these weighted procedures may have a significant difference from real systems. Therefore, it is important to describe such real weighted distribution methods in weighted networks.

Thereby, it is of great importance to study different weighted distribution methods of a distribution network, which has not been discussed yet. This paper employs some weighted distribution methods that are found in many realistic networks such as the link weight and investigates the phase synchronization performance in a large-scale islanded network to pave the way for the solution of the similar weighted distribution issues in the future. The innovations of this paper are as follows.


(3) Different from other studies focusing on the weighted distribution methods which only associate with the inherent properties of the network, this paper adopts several realistic distribution styles as weighted distribution methods to research the synchronization stability in the large-scale distribution network.

The rest of this paper is organized as follows. In Section 2, we introduce correlation complex network theory. Section 3 focuses on the modeling of weighted distribution and derivation of the critical coupling strength formula. The numerical results and simulation analyses in Section 4 are used to compare the network stability of different weighting methods. Section 5 concludes this paper.

#### **2. Correlation Theory**

The structure of an active distribution network can be achieved by mirroring algebraic graph theory. Furthermore, the synchronization condition of an active distribution network will be analyzed by order parameters that are introduced thoroughly in the following. Then, the critical coupling strength of network stability can be obtained with the help of synchronization condition.

#### *2.1. Graph Theory*

The microgrid can be modeled as a graph *G*, where *V*= {*θ*1,*θ*2, ··· *θN*} is the set of edges while *E* ∈ *V* × *V* is the set of edges. Generally, the connected graph can be called an undirected graph when (*θi*, *θj*) and (*θj*, *θi*) represent the same edge. The value of weight *aij* defines whether the graph is weighted (*aij* = 1) or not (*aij* = 1) and hence we have the associated adjacency matrix *A* = [*aij*]. Moreover, the Laplace matrix *L* is a symmetric zero-sum matrix of undirected graphs, where *lii* = <sup>∑</sup>*<sup>i</sup>* <sup>=</sup>*<sup>j</sup> aij*, *lij* = −*aij*. The matrix *<sup>L</sup>* satisfies:

$$\begin{cases} l\_{ij} \le 0, i \ne j \\ \sum\_{j=1}^{N} l\_{ij} = 0, i = 1, \dots, N \end{cases} \tag{1}$$

Whether the graph *G* is a directed graph or an undirected graph, its associated matrix is based on a directed label, which means edge *e* = {*i*, *j*} ∈ *E* is represented by an ordered pair (*i*, *j*). In the undirected graph, we assign an arbitrary direction to it. For the ordered pair (*i*, *j*), we define *bje* = 1 if *j* is the sink node of edge *e*, *bie* = −1 if *i* is the source node of edge *e*, and *bke* = 0 otherwise in the incidence matrix *<sup>B</sup>* = [*bie*] <sup>∈</sup> *<sup>R</sup>n*×*m*. In addition, the edge (*i*, *<sup>j</sup>*) satisfy (*BTX*)*<sup>e</sup>* <sup>=</sup> *xj* <sup>−</sup> *xi*.

#### *2.2. Synchronization Conditions*

Specifically, the phase and frequency of all DG units tend to stabilize at a fixed value when *t* → ∞. Therefore, we focus on the synchronization process of phase and frequency to research the system stability intuitively.

If the coupling strength *K* between the DG nodes in a complex network is larger than the critical value *Kcr*, the DG nodes will divide into two groups, the phase of one cluster is fixed at a stable specific frequency while the other oscillators have a drift.

The order parameter is a visual signal to show the change of phase, whose definition is:

$$r e^{i\psi} = \frac{1}{N} \sum\_{j=1}^{N} e^{i\theta\_j} \tag{2}$$

The size of *r* shows the percentage of stable nodes. The value is:

$$|r| = \frac{1}{N} \left| \sum\_{j=1}^{N} e^{i\theta\_j} \right| \tag{3}$$

where *r* = 1 means that a system enters a stable state absolutely while *r* = 0 means that a system becomes unstable.

**Criterion 1 [21]:** For a complex network graph *G* with an incidence matrix *B*, each pair of interconnected oscillators {*i*, *j*} ∈ *ε* could enter frequency synchronization and phase convergence *θ<sup>i</sup>* <sup>−</sup> *<sup>θ</sup><sup>j</sup>* <sup>≤</sup> *<sup>γ</sup>* <sup>&</sup>lt; *<sup>π</sup>*/2 when satisfying

$$\left\| L^{\dagger} \omega \right\|\_{\varepsilon,\infty} := \left\| B^{T} L^{\dagger} \omega \right\|\_{\infty} \leq \sin(\gamma) \tag{4}$$

where *<sup>L</sup>*† is the inverse of Laplace matrix *<sup>L</sup>*, *<sup>ω</sup>*= (*ω*1, ··· *<sup>ω</sup>n*) *<sup>T</sup>* is the vector of the natural frequency of each oscillator in the network and · *ε*,<sup>∞</sup> evaluates the worst-case dissimilarity of this weighted projection.

This paper will use this synchronization condition to obtain the critical coupling strength of network synchronization under weighted conditions.

#### *2.3. Critical Coupling Strength*

In the previous subsection, the synchronization conditions recently mentioned in [22] have been introduced. In different occasions, these conditions have been proven to be applicable to oscillator models. This paper takes the Kuramoto oscillator to express the nonlinear ordinary differential equations:

$$\theta\_i = \omega\_i + \sum\_{j \in N\_i} a\_{ij} \sin(\theta\_j - \theta\_i), i = 1, 2, 3...N \tag{5}$$

The model can be expressed as a matrix:

$$\dot{\theta}\_{\text{i}} = \omega - B\mathcal{W}\sin(B^T\theta\_{\text{i}}) , \text{i} = 1, 2, 3... \text{N} \tag{6}$$

*W* := *diag*(*w*) is the diagonal matrix of weights. In a complex network, which determines the stability of the system is not the absolute value of the phase but the phase difference. The network entry synchronizes only when two phases of coupled nodes *ϕ*<sup>1</sup> and *ϕ*<sup>2</sup> are locked at a certain ratio, which means <sup>|</sup>*nϕ*<sup>1</sup> <sup>−</sup> *<sup>m</sup>ϕ*2<sup>|</sup> <sup>&</sup>lt;constant. Therefore, the phase difference vector is *<sup>ϕ</sup>*(*t*) :<sup>=</sup> *<sup>B</sup>Tθ*, and the Equation (6) can be expressed by the phase difference:

$$\dot{\varphi} = B^T \omega - B^T B \mathcal{W} \sin(\varphi) \tag{7}$$

In the above equation, the synchronous solution can be expressed as *ϕ*˙=0 or *ϕ*˙=*ϕ*∗ for a fixed point *<sup>ϕ</sup>*<sup>∗</sup> <sup>∈</sup> *<sup>R</sup>n*.

In an analysis of [23], Fazlyab found that all fixed points of Equation (7) can be expressed as:

$$\sin(\varphi^\*) = \mathcal{W}^{r-1} B^T (B\mathcal{W}^r B^T)^\dagger \omega + \mathcal{W}^{-1} F y \tag{8}$$

where *<sup>F</sup>* <sup>∈</sup> *<sup>R</sup>m*×*m*−*n*+<sup>1</sup> is a matrix whose columns span the null space of *<sup>B</sup>*, i.e., *BF* <sup>=</sup> 0, *<sup>y</sup>* <sup>∈</sup> *<sup>R</sup>m*−*n*<sup>+</sup>1, and *r* is an arbitrary real number.

When *F* = 0 and the fixed point *ϕ*<sup>∗</sup> is unique, the particular solution *Wr*−1*BT*(*BWrBT*)†*ω* is independent of *r*. In the special case of *r* = 0, we have that

$$\sin(\varphi^\*) = B^T L^\dagger \omega \tag{9}$$

provided that *B<sup>T</sup> <sup>L</sup>*†*<sup>ω</sup>* <sup>∞</sup> ≤ 1.

This section has proposed the theory about synchronization. The next section will analyze the situation of different weighting methods in star topology through algebraic graph theory, and obtain the critical coupling strength by using the synchronization condition.

#### **3. Weighted Distribution**

The weighted distribution is not only helpful to describe the features of real networks but are also crucial for controlling network synchronization. Three different weighting methods will be discussed to research the phase synchronization of a large-scale islanded network in this section. This issue takes an ideal star topology as the research object. Furthermore, a critical coupling strength formula for the weighted star network would be proposed first in this paper.

#### *3.1. Star Topology*

The purpose of this investigation is to explore the synchronization stability of large-scale grid-tied DG in islanded operation. Considering that most of the distribution networks are star topology [24,25], this paper chooses star topology as the object, as shown in Figure 1. The star network is shaped by one primary node and many branch nodes connected with it. It is characterized by high reliability, easy operation, and easy maintenance when the branch lines fail. However, this structure will consume a lot of metal and switchgear. Therefore, it is suitable for important construction sites with large capacity equipment or special places with concentrated load and wet, corrosive environments [21,26,27].

**Figure 1.** Star topology.

Due to the practical application of the power grid, the traditional research object of multi-machine (three to five machines) can no longer meet the actual needs. This is why this paper delivers an explicit analysis of the network stability with large-scale grid-tied DG. The number of DG *N* in a Large-scale model generally refers to *N* ≥ 200, where this paper takes *N* as 201. There are no conventional generators with rotors used in this model, and all DGs are photovoltaic units. The photovoltaic community electrical diagram with star topology is depicted in Figure 2. We assume that the value of the line impedance between load nodes and the master node is consistent, and the inverter of each load node has the same impedance, so the impedance value between load node and the master node is equal. Therefore, *Z* = *Zline* + *Zinv* and the admittance value is *Y* = <sup>1</sup> *<sup>Z</sup>* <sup>=</sup> <sup>1</sup> *Zline*+*Zinv* .

**Figure 2.** Star topology electrical diagram.

In Figure 2, each load node is only connected to the master node *L*1, and the impedance value between nodes *Zij* = *Z* + *Zinv* are equal. Therefore, the admittance matrix with star topology is

$$\begin{aligned} \left| Y\_{\bar{\mathcal{Y}}} \right| &= \begin{bmatrix} 0 & \left| Y\_{12} \right| & \cdots & \left| Y\_{1N} \right| \\ \left| Y\_{21} \right| & 0 & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots \\ \left| Y\_{N1} \right| & \cdots & \cdots & 0 \end{bmatrix}\_{\left( N \times N \right)} \\ &= \begin{bmatrix} 0 & \left| \frac{1}{Z + Z\_{\text{inv}}} \right| & \cdots & \left| \frac{1}{Z + Z\_{\text{inv}}} \right| \\ \left| \frac{1}{Z + Z\_{\text{inv}}} \right| & 0 & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots \\ \left| \frac{1}{Z + Z\_{\text{inv}}} \right| & \cdots & \cdots & 0 \end{bmatrix}\_{\left( N \times N \right)} \end{aligned} \tag{10}$$

#### *3.2. The General Formula for Weighted Distribution*

The DG units with droop control can be described as a Kuramoto oscillator [28–30]:

$$\begin{split} \frac{d\theta\_i}{dt} &= \theta\_t = \omega\_i - \omega\_\theta = -m\_i \left( p\_i - p\_{i\theta} \right) = m\_i p\_{i\theta} - m\_i p\_i \\ &= m\_i p\_{i\theta} - m\_i E\_i \sum\_{j=1}^n \left| \begin{smallmatrix} \mathbb{Y}\_{ij} \end{smallmatrix} \right| E\_j \sin \left( \theta\_i - \theta\_j \right) \end{smallmatrix} \tag{11}$$

where *mi piθ*=*ωi*, *ω<sup>i</sup>* is the rated frequency of inverter in the *i*th DG node. In addition, the value of *ω<sup>i</sup>* is chosen from the probability density function *g*(*ω*). When the standard frequency is 50 Hz, we take *g*(*ω*) as the symmetric normal distribution of 50, therefore *g*(*ω*)∼ (50, 1).

In the above equation, the weight between *DGi* and *DGj* nodes is obtained by the droop coefficient *di*, the voltage of node *Ei*, *Ej*, and the admittance *Yij* . At this point, the coupling strength

$$K\_i = d\_i \left| \mathcal{Y}\_{i\bar{j}} \right| E\_i E\_{\bar{j}} = k m\_i \bar{z}$$

For the sake of simplicity, we divide the weight between the DG nodes into two parts. One is the coupling strength *k*, and the other is the weight index *mi* that conforms to some kind of weighted distribution. Therefore, the weight of the directed edge {*i*, *j*} is determined by the weight index *mi* while the weight of the edge {*j*, *i*} is determined by the weight index *mj*, so the weight distribution of one network can be decided by the distribution of weight index *mi*. In addition, we can synchronously adjust the size of all *mi* values to change the size of the weight. Thus, we have the weight diagonal matrix *W* = *k* × *diag*(*m*), where *m*<sup>1</sup> is the weight when a load node is the sink node and the master node is the source node, *m*<sup>2</sup> ··· *mn* is the weight when the master node is the sink node and a load node is the source node. Thus, the weight distribution of the network can be simulated by adjusting the value of the weight index:

$$\begin{aligned} A &= \begin{bmatrix} 0 & km\_1 & \dots & km\_1 \\ km\_2 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ km\_n & 0 & \dots & 0 \end{bmatrix}\_{(201\times 201)} \\ D &= \begin{bmatrix} (n-1)\times km\_1 & 0 & \dots & 0 \\ 0 & km\_2 & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & \dots & \dots & km\_n \end{bmatrix}\_{(201\times 201)} \\ L &= \begin{bmatrix} (n-1)\times km\_1 & km\_1 & \dots & km\_1 \\ km\_2 & -km\_2 & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ km\_n & 0 & \dots & -km\_n \end{bmatrix}\_{(201\times 201)} \end{aligned} \tag{12}$$

In a star topology whose *N* is 201, the incidence matrix is

$$B = \begin{bmatrix} 1 & 1 & \cdots & 1 & 1 \\ -1 & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & -1 & 0 \\ 0 & 0 & \cdots & 0 & -1 \end{bmatrix}\_{(201 \times 200)} \tag{13}$$

#### *3.3. Weighted Distribution Model*

Different weighting methods have different performances in different topology, but the topology of an actual power grid is fixed and difficult to adjust. Thus, the weighting methods, which are associated with the degree distribution and node strength, are not considered in this paper. Only the weighting method associated with the distribution of coupling strength is considered in this paper.

In practice, the error occurs in the installation process often obey Gaussian distribution, the weight distribution of complex network is mostly in line with the power-law distribution, and the frequency-weighted distribution can reflect some peculiarities of the actual network. Thus, this paper discusses the weighting methods of Gaussian distribution, power-law distribution, and frequency-weighted network, and the aforementioned methods can be simulated by adjusting the weight index.

#### 3.3.1. Gaussian Distribution

At first, this paper considers the case that the weight distribution method is the same as the rated frequency of an inverter distribution method so that the distribution function *g*(*m*) of *mi* is the same as the distribution function *g*(*ω*) of *ωi*; both are Gaussian distribution. Unlike the rated frequency of the inverter being a symmetric normal distribution of 50, the weight distribution is a symmetric normal distribution of the *y*-axis. Its Kuramoto oscillator model is:

$$\begin{aligned} \theta\_i &= \omega\_i + k \lg(m) \sum\_{j=1}^N a\_{ij} \sin(\theta\_i - \theta\_j), i = 1, \dots, N \\ \theta\_i &= \omega\_i + k \left| m\_i \right| \sum\_{j=1}^N a\_{ij} \sin(\theta\_i - \theta\_j), i = 1, \dots, N. \end{aligned} \tag{14}$$

*aij* = 1 when DG*<sup>i</sup>* has coupling relationship with DG*j*, otherwise *aij* = 0.

#### 3.3.2. Power-Law Distribution

Power-law distribution has extensive application in complex network theory. Barabasi and Albert pointed out that the degree distribution of complex networks is the power-law distribution in [31]. Then, Chunguang Li and Guanrong Chen verified that the weight distribution is a power-law distribution in a general complex network. Furthermore, according to the research of Barrat et al. in [32], the node strength distribution in the weighted network obeys the power-law distribution as well. In this paper, only the weighted distribution of weights is discussed. Therefore, it is worth mentioning that the assumed distribution function *g*(*m*) of *mi* satisfies the power-law distribution, which means *<sup>g</sup>*(*m*) <sup>∼</sup> *<sup>m</sup>*−*β*, and

$$\begin{aligned} \theta &= \omega\_i + k \lg(m) \sum\_{j=1}^{N} a\_{ij} \sin(\theta\_i - \theta\_j), i = 1, \dots, N \\ \theta &= \omega\_i + k m\_i \sum\_{j=1}^{N} a\_{ij} \sin(\theta\_i - \theta\_j), i = 1, \dots, N \end{aligned} \tag{15}$$

#### 3.3.3. Frequency-Weighted Network

The frequency-weighted network reflects the feature of several natural and social networks. In fact, a power grid could be expressed as a Kuramoto oscillator model, and the weight of the coupling strengths are obtained by their natural frequency. In this case, the weighted coupling coefficients are related to DGs' own inverter rated frequency rather than the topology of a network. We added the inverter rated frequency to the dynamic model as follows:

$$\dot{\theta}\_{i} = \omega\_{i} + \frac{k \, |\omega\_{i}|}{\sum\_{j=1}^{N} a\_{ij}} \sum\_{j=1}^{N} a\_{ij} \sin(\theta\_{i} - \theta\_{j}) , i = 1, \ldots, N \tag{16}$$

#### *3.4. Critical Coupling Strength Formula*

Reference [23] pointed out that any fixed point in Equation (7) has local exponential stability when *ϕ*∗ <sup>∞</sup> ≤ *π*/2. For the sake of simplicity, we decompose the matrix *A* = *k* × *A* , where *A* is the adjacency matrix entirely determined by weight indexes, and *k* is the coupling strength. The in-degree matrix *D* is a diagonal matrix, whose values in each row are the row sum of corresponding rows of *A* matrix. Certainly, the *D* matrix can be decomposed into *D* = *k* × *D* , where *D* is the in-degree matrix entirely determined by weight indexes. At the same time, the Laplacian matrix *L* = *D* − *A* can be decomposed into *L* = *k* × *L* , where *L* is the Laplacian matrix entirely determined by weight indexes. Furthermore, Equation (9) can be decomposed into:

$$\sin(\varphi^\*) = k^{-1} B^T (L')^\dagger \omega \tag{17}$$

There is sin(*BTθ*∗) <sup>∞</sup> = 1 when the network enters a critical stable state, and then the critical coupling values of star topology in different weighting methods can be transformed into the following formula:

$$k\_{\mathcal{E}} = \left\| B^T (L')^\dagger \omega \right\|\_{\infty} \tag{18}$$

The calculation of the critical coupling value in different weighting methods is described in detail in Section 4.

#### **4. Simulation**

In this section, three weighted methods from above are simulated, and the proposed formula of critical coupling (18) is verified.

#### *4.1. Case 1*

It is inevitable that deviations will generate in the process of installing the grid-tied DG inverters in an active distribution network, and the distribution of deviations often conforms to the Gaussian distribution.

In this case, the assumed actual parameters are reported in Table 1:

**Table 1.** Related parameters in Case 1.


#### 4.1.1. Weighted Undirected Network

Let us now turn to the case of the undirected situation, the weight distribution of edges in this situation conforms to the Gaussian distribution *X* ∼ *N*(0, 1), and source point and sink point of each edge are arbitrary. When *N* takes 201, and the number of edges takes 200 in this case; meanwhile, the distribution is depicted in Figure 3. Furthermore, Figure 4 shows the probability density function of this Gaussian distribution, and it is obvious that the distribution law conforms to the (0,1) distribution.

**Figure 3.** Nodes in Gaussian distribution.

**Figure 4.** Probability density function of Gaussian distribution.

Thus, the ordinary differential equation of oscillator *θ<sup>i</sup>* is

$$\dot{\theta} = \omega\_i + kX \sum\_{j=1}^{N} a\_{ij} \sin(\theta\_i - \theta\_j), i = 1, \dots, N \tag{19}$$

where *X* is the weight of the corresponding edge.

Turn now to verifying the critical coupling strength formula of the weighted situation proposed in Section 3. In this case,

$$L = \begin{bmatrix} \begin{bmatrix} -\sum \mathbf{X} & \mathbf{X}\_1 & \cdots & \mathbf{X}\_{200} \\ \mathbf{X}\_1 & -\mathbf{X}\_1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ \mathbf{X}\_{200} & 0 & \cdots & -\mathbf{X}\_{200} \end{bmatrix}\_{\left(201 \times 201\right)} \tag{20}$$

Furthermore, the related parameters in this case are depicted in Table 1, and the weight *k* = *<sup>a</sup> <sup>Z</sup>*+*Z*inv <sup>×</sup> *<sup>E</sup>*<sup>2</sup> <sup>=</sup> 0.1×<sup>2202</sup> <sup>100</sup>×0.153+10×1.15 <sup>=</sup> 180. The phase and frequency stability are shown in Figures 5 and 6, and it can be seen that the phases of DG nodes change at the same rate and the frequency of DG nodes fixed in one constant value.

**Figure 5.** Phase stabilization process of undirected network in case 1.

**Figure 6.** Frequency stabilization process of undirected network in case 1.

In this paper, order parameters are used to reflect the synchronization of DGs in the network. Figure 7 shows the relationship between the coupling strength *K* and the order parameter *r* when *N* has different values, which reflect the change of network synchronization with the coupling strength. It is obvious that the network stability is the best when the value of *N* is 51, followed by 101. When *N* takes 201, the network has the worst stability. It can be seen that the stability of the network decreases with the increase of *N*, which is not consistent with the situation in which the stability does not change with the value of *N* when the star topology is an undirected network.

**Figure 7.** The stability of the weighted undirected network when *N* is 51,101,201.

#### 4.1.2. Weighted Directed Network

Thus far, this paper has focused on an ideal situation that the weighted distribution of the undirected edge in the network is Gaussian distribution. Let us now turn to consider the weighted directed case. In a power grid, each DG can become not only the power transmitter but also the power receiver. Therefore, it is of great practical value to consider the different weights of edges when they take different directions between coupling DG nodes. That is, the edge weight is the weight index of the master node when the master node is the source point, while the edge weight is the weight index of the load node when the load node is the source point:

Since the inverter rated frequency of the DG unit also conforms to the Gaussian distribution, this paper assumes that the weight distribution of edges is the same type as the distribution of inverter rated frequency. Namely, the ordinary differential equation with the oscillator *θ<sup>i</sup>* is

$$\theta = \omega\_i + k \left| \omega\_i \right| \sum\_{j=1}^{N} a\_{ij} \sin(\theta\_i - \theta\_j), i = 1, \dots, N \tag{21}$$

This paper takes the weight index of the load node as the Gaussian distribution and the weight index of the master node as the mean value of the weight index of each load node to realize the power balance. Thus, there is

$$L = \begin{bmatrix} (n-1) \ast m\_1 & m\_1 & \cdots & m\_1 \\ & m\_2 & -m\_2 & \cdots & 0 \\ & \vdots & \vdots & \vdots & \vdots \\ & m\_n & 0 & \cdots & -m\_n \end{bmatrix}\_{(201\times 201)}$$

$$= \begin{bmatrix} -\sum\_2^{201} m\_i & \frac{\sum\_2^{201} m\_i}{200} & \cdots & \frac{\sum\_2^{201} m\_i}{200} \\ & m\_2 & -m\_2 & \cdots & \cdots \\ & \vdots & \vdots & \vdots & \vdots \\ & m\_n & 0 & \cdots & -m\_n \end{bmatrix}\_{(201\times 201)} \tag{22}$$

Frequency deviation will lead to power imbalance and even instability of the whole system. The phase and frequency stabilization processes in the directed situation are shown in Figures 8 and 9 when *k* = 180. Obviously, the coupling strength *k* in this case is much greater than the critical coupling strength required for network stability. Let us now turn to the critical coupling strength in the case of Gaussian distribution. When *N* = 201, the value of L matrix can be substituted for formula (18) to obtain the critical coupling strength *KC* = 6.7; when *N* is 51 and 101, the critical coupling strength *KC* = 3.7 and 4.46 can be obtained from the calculation formula, respectively. Figure 10 shows that, when *N* takes 201, the network enters synchronization rapidly before *K* < 6.7, while the network is nearly stable after *K* > 6.7. This result is basically consistent with the result obtained by formula (4-1). Similarly, when the value of *N* is 51 and 101, the network becomes stable after *KC* = 3.7 and 4.46, respectively. Thus far, two conclusions have been shown. The first is that the validity of the critical coupling strength formula (18) in the directed weighted network is verified, and the second is that the network stability decreases with the increasing number of nodes in the island star topology directed weighted network.

**Figure 8.** Phase stabilization process of directed network in case 1.

**Figure 9.** Frequency stabilization process of directed network in case 1.

**Figure 10.** The stability of the weighted directed network when *N* is 51,101,201.

#### *4.2. Case 2*

In particular, previous studies show that most of the nodes have a degree that is a less than average value while a few nodes have large connections in the real network. The study of Chunguang Li and Guanrong Chen [17] confirmed that the distribution of connection weights would follow the power-law distribution in general actual complex networks. Thus, it is important to research the stability of a weighting case that conforms to the power-law distribution.

In this case, the assumed actual parameters are reported in Table 2:

**Table 2.** Related parameters in Case 2.


#### 4.2.1. Weighted Undirected Network

In this case, we assume that the weight distribution of 200 edges conforms to the power-law distribution, as shown in Figure 11. Then, Figure 12 shows the probability density function diagram. It can be seen that, with the increase of abscissa, the probability of weight distribution is lower. Furthermore, the probability of randomly selecting a node with greater weight is negatively correlated with the number of the same weight. The power-law index in Figure 11 is 2. The distribution of nodes depicted as a line whose slope is negative in the double logarithmic diagram, as demonstrated in Figure 12. This distribution has an asymptote of the black dotted line, where the exponent determines the rate of decay.

**Figure 11.** Nodes in power-law distribution.

**Figure 12.** Probability density function of power-law distribution.

The weighted distribution in this case follows *<sup>g</sup>*(*X*) <sup>∼</sup> *<sup>X</sup>*−*β*, which means that the ordinary differential equation of oscillator *θ<sup>i</sup>* is

$$\dot{\theta} = \omega\_{\bar{i}} + kX \sum\_{j=1}^{N} a\_{\bar{i}j} \sin(\theta\_{\bar{i}} - \theta\_{\bar{j}}),\\ \dot{i} = 1, \dots, N \tag{23}$$

where *X* is the weight of the corresponding edge.

Like the weighted undirected network in case 1, the weighted distribution in this case is reflected by the value of the edge weight, so the Laplace matrix is

$$L = \begin{bmatrix} \begin{bmatrix} -\sum \mathbf{X} & \mathbf{X}\_1 & \cdots & \mathbf{X}\_{200} \\ \mathbf{X}\_1 & -\mathbf{X}\_1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ \mathbf{X}\_{200} & 0 & \cdots & -\mathbf{X}\_{200} \end{bmatrix}\_{\left(201 \times 201\right)} \tag{24}$$

The same as the coupling strength value of case 1, the phase and frequency stabilization process diagram are shown in Figures 13 and 14. In addition, the influence of different power exponents on network stability is shown in Figure 15. Figure 15 shows the network synchronization process when power-law index *β* grades –5, –3, –2, 2, 3, 5. It is obvious that the critical coupling strength of the system when the power-law index is positive does not differ greatly, while the positive power-law index values have significantly better stability than the negative values. This paper comes to a conclusion that the positive and negative value of power-law index rather than the size of the Power-law index has a key influence on the stability of the network when ignoring the deviations generated in the calculation process.

**Figure 13.** Phase stabilization process of undirected network in case 2.

**Figure 14.** Frequency stabilization process of undirected network in case 2.

**Figure 15.** The stability of the weighted undirected network when *β* is ±2, ±3, ±5.

The simulation analysis of the different number of nodes *N* is shown in Figure 16. When the index is 2, the value of *N* is negatively correlated with network stability. The conclusion also reflects the analysis of the relationship between the node number and the network stability in the previous case.

**Figure 16.** The stability of the weighted undirected network when *N* is 51,101,201(*β* = 2).

#### 4.2.2. Weighted Directed Network

Turn now to the directed weighted case. Similar to case 1, the method of adjusting the weight index is used to complete the weight distribution conforms to the power-law distribution *<sup>g</sup>*(*m*) <sup>∼</sup> *<sup>m</sup>*−*β*, and the ordinary differential equation is given in the formula (25):

$$\dot{\theta} = \omega\_i + k \lg(m) \sum\_{j=1}^{N} A\_{ij} \sin(\theta\_i - \theta\_j), i = 1, \dots, N \tag{25}$$

Similar to the weighted directed network in case 1, the weight index of the master node in this case is the mean value of the weight index of each load node considering the problem of power balance. Therefore, the Laplace matrix in the network with 201 nodes is:

$$L = \begin{bmatrix} (n-1) \ast m\_1 & m\_1 & \dots & m\_1 \\ m\_2 & -m\_2 & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ m\_n & 0 & \dots & -m\_n \end{bmatrix}\_{(201\times 201)}$$

$$L = \begin{bmatrix} -\sum\_2^{201} m\_i & \frac{\sum\_2^{201} m\_i}{200} & \dots & \frac{\sum\_2^{201} m\_i}{200} \\ m\_2 & -m\_2 & \dots & \dots \\ \vdots & \vdots & \vdots & \vdots \\ m\_n & 0 & \dots & -m\_n \end{bmatrix}\_{(201\times 201)}\tag{26}$$

The phase and frequency stabilization process in the directed situation are shown in Figures 17 and 18 when the value of power-law index beta is 2, and coupling strength *k* is 180. When the value of the power-law index respectively takes –5, –3, –2, 2, 3, 5, the critical coupling strength *KC* of network is respectively 1.4976, 1.4353, 1.9639, 1.4976, 1.4353, and 0.7771 by formula (18). These critical coupling strengths are nearly consistent with the critical value depicted in Figure 19. In addition, the weighted directed network has similar performance to the weighted undirected network, which means that the positive power-law index value has better network stability performance than the negative one, while the size of the power-law index value has almost no effect. Apparently, in the case of power-law distribution, the stability of the network is obviously better than that in the case of Gaussian distribution. When the Power-law index is 2, the network can be stable only on the condition that the coupling strength is greater than the critical value *KC* = 0.5475.

Furthermore, the relationship between the number of DG nodes and network stability is shown in Figure 20. The negative correlation between DG number and critical coupling strength is still evident.

**Figure 17.** Phase stabilization process of directed network in case 2.

**Figure 18.** Frequency stabilization process of directed network in case 2.

**Figure 19.** The stability of the weighted directed network when *β* is ±2, ±3, ±5).

**Figure 20.** The stability of the weighted directed network when *N* is 51,101,201 (*β* = 2).

#### *4.3. Case 3*

The following case will discuss the Frequency-weighted situation. The natural frequency is inverter rated frequency in the Kuramoto oscillator model. Thus, case 3 chooses the absolute ratio between the inverter rated frequency of one DG node and its in-degree as the weight; there is the ordinary differential equation:

$$\theta = \omega\_i + \frac{k \, |\omega\_i|}{\sum\_{j=1}^{N} A\_{ij}} \sum\_{j=1}^{N} A\_{ij} \sin(\theta\_i - \theta\_j), i = 1, \dots, N \tag{27}$$

The weighted distribution in this case is different from the previous two cases, whose weight index is not a random value which obeys a certain distribution. The Laplace matrix in this case is

$$L = \begin{bmatrix} -\left|\omega\_1\right| & \frac{\left|\omega\_1\right|}{200} & \cdots & \frac{\left|\omega\_1\right|}{200} \\ \left|\omega\_2\right| & -\left|\omega\_2\right| & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots \\ \left|\omega\_{n-1}\right| & 0 & \cdots & 0 \\ \left|\omega\_n\right| & 0 & \cdots & -\left|\omega\_n\right| \end{bmatrix} \tag{28}$$

In case 3, the phase and frequency stabilization process diagram are shown in Figures 21 and 22, which are similar to the above two cases. According to formula (18), the critical coupling strength of this case is *KC* = 12.89 when *N* = 201, and the network synchronization progress is shown in Figure 23. It can be seen that the weighted method in this case is not as stable as the Gaussian distribution and the power-law distribution.

**Figure 21.** Phase stabilization process in case 3.

**Figure 22.** Frequency stabilization process in case 3.

**Figure 23.** Synchronization progress when *N* is 51,101,201 in case 3.

#### *4.4. Stability Analysis*

Through comparing the critical coupling strength of three cases as illustrated in Figure 24, we find that the power-law distribution has the best stability in different weighted distribution methods. Moreover, the second largest eigenvalue *λ*<sup>2</sup> or the eigenvalue proportion *λN*/*λ*<sup>2</sup> in adjacency matrix A could describe the synchronization ability of the system. The network coherence is proportional to *λ*<sup>2</sup> and inversely proportional to *λN*/*λ*2. Table 3 shows that the *λ*<sup>2</sup> of power-law distribution has the smallest value, which means it has the best synchronization ability, followed by the Gaussian distribution and frequency-weighted distribution. Thus, the results of Table 3 and Figure 24 confirm each other. In addition, the contrastive analysis of weighted directed network and weighted undirected network show that the weighted directed network which conforms to the Gaussian distribution or power-law distribution is easier to get into synchronization than the weighted undirected network. Furthermore, by studying the network with different numbers of DG, this paper finds that the stability of the weighted network decreases with an increase of *N*, which is different from the average field. For the weighting method of power-law distribution, the size of the power-law index *β* is not significant to the stability of the network, but the sign of the power-law index *β* is not negligible to the stability of the network.

**Figure 24.** Synchronization progress of different weighting distributions.

**Table 3.** Spectrum analysis.


#### *4.5. Further Analysis of Weighted Network in the Actually Distributed Network*

Actually, the topology structure of an actual power grid network often isn't a standard star topology. The standard star topology is an ideal assumption. In order to verify the stability analysis from the above subsection, we take the dynamic IEEE test systems [33] and the Iceland 189-node distributed grid [34] as a research object. The synchronizing processes of Gaussian distribution, power-law distribution, and frequency-weighted network are shown in Figures 25 and 26.

#### 4.5.1. Dynamic IEEE Test Systems

Three different weighted distribution methods have been applied to the dynamic IEEE models, including 9-bus, 14-bus, 57-bus, and 118-bus modified test systems. The diagram of coupling strength and order parameter is shown in Figure 25, where the labels G, P, and F represent Gaussian distribution, Power-law distribution and frequency-weighted distribution, respectively, and the network with the Power-law distribution method always enter the stable state at first. The value of the critical coupling strength of three weighted distribution methods in the dynamic IEEE test systems is shown in Table 4, where each value is calculated by Equation (18). It shows that the network with the power-law distribution method always has the least critical coupling strengths, which corresponds to the stability performance in Figure 25.

**Table 4.** *Kc* of three weighted distribution methods in different IEEE systems.


**Figure 25.** Synchronization progress of different weighting methods in the dynamic IEEE test systems.

#### 4.5.2. Iceland 189-Nodes Distribution Grid

Figure 26 shows the results of the stability performance of three distribution methods in an Iceland 189-nodes distribution grid. It can be seen that the power-law distribution undoubtedly has the best stability performance. Second is the frequency-weighted network. Moreover, the Gaussian distribution needs a long process to enter synchronization. These results further support the conclusions of the above subsection.

**Figure 26.** Synchronization progress of different weighting methods in the Iceland 189-node distributed grid.

#### **5. Conclusions**

Above all, this paper analyzes the stability of large-scale grid-tied DG in a non-homogeneous active power network by using complex networks theory.

Then, three different weighting distribution methods, including Gaussian distribution, power-law distribution, and frequency-weighted distribution, are analyzed. In addition, the critical coupling strength formula for the weighted situation *kc* = *BT*(*L*) † *ω* <sup>∞</sup> has been proposed and verified. The critical coupling strength obtained from this formula conforms to the reality when ignoring tiny deviation.

Among the three different weighting methods, the power-law distribution has the best stability. Moreover, the difference of power-law index value has little effect on the network stability, but the sign of power-law index has quite an effect on the system stability. The second is the Gaussian distribution, which is obviously less stable than the power-law distribution. Finally, the frequency-weighted network has the worst synchronization performance, which is significantly inferior to the other two. In addition, further analysis of the IEEE dynamic test systems and the Iceland 189-node distributed grid has proved the better stability of power-law distribution in this paper.

**Author Contributions:** Conceptualization, S.C. and Y.Z.; methodology, H.Z.; software, J.L.; validation, C.Y., S.C. and J.L.; formal analysis, Y.Z.; investigation, C.Y.; resources, H.Z.; data curation, Y.Z.; writing–original draft preparation, S.C.; writing–review and editing, J.L.; visualization, S.C.; supervision, Y.Z.; project administration, J.L.; funding acquisition, H.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China under Grants Nos. 61873195 and 51707071.

**Acknowledgments:** The authors like to express sincere gratitude to School of Electrical Engineering and Automation, Wuhan University for making this research work for execution and technical implementation in real-time.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Energies* Editorial Office E-mail: energies@mdpi.com www.mdpi.com/journal/energies

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-03936-759-7