Load Margin Assessment of Power Systems Using Physics-Informed Neural Network with Optimized Parameters

Abstract

Challenges in the operation of power systems arise from several factors such as the interconnection of large power systems, integration of new energy sources and the increase in electrical energy demand. These challenges have required the development of fast and reliable tools for evaluating the operation of power systems. The load margin (LM) is an important index in evaluating the stability of power systems, but traditional methods for determining the LM consist of solving a set of differential-algebraic equations whose information may not always be available. Data-Driven techniques such as Artificial Neural Networks were developed to calculate and monitor LM, but may present unsatisfactory performance due to difficulty in generalization. Therefore, this article proposes a design method for Physics-Informed Neural Networks whose parameters will be tuned by bio-inspired algorithms in an optimization model. Physical knowledge regarding the operation of power systems is incorporated into the PINN training process. Case studies were carried out and discussed in the IEEE 68-bus system considering the N-1 criterion for disconnection of transmission lines. The PINN load margin results obtained by the proposed method showed lower error values for the Root Mean Square Error (RMSE), Mean Square Error (MSE) and Mean Absolute Percentage Error (MAPE) indices than the traditional training Levenberg-Marquard method.

1. Introduction

Currently, power systems face different and various operational challenges. The massive integration of renewable energy sources such as wind and photovoltaics with intermittent active power generation characteristics, the increase in demand for electrical energy and the most varied types of contingencies have hampered the operation of power systems and their purpose of supplying electrical energy. continuously, safely and reliably for everyone [1,2,3,4]. This scenario has required the development of different monitoring, control and protection tools to assist the power system operator in real-time operation.

Technological advances over the years have allowed the development of Synchronized Phasor Measurement Systems that provide voltage and current magnitude and phase measurements in real time with time synchronization and high sampling rates acquired by Phasor Measurement Units (PMUs) installed on system buses [5]. Data from PMUs have led to the development of different and promising tools for monitoring [6,7,8,9], protecting [10,11,12] and controlling [13,14,15,16,17,18] power systems.

Data from PMUs can also be very useful for power system stability studies and help the system operator in making decisions when necessary. Traditionally, the Voltage Stability Margin (VSM) is a fundamental index for evaluating how high the system load can be without causing a voltage collapse in the system. Different methods have been proposed in recent decades to calculate this VSM index such as Continuation Power Flow method [19], loop-analysis-based continuation power flow algorithm [20], methods based on maximum power transfer [21,22], coupled single-port Thevenin equivalent model [23], PQV curves [24], Jacobian matrix [25]. However, this method requires complete information from the entire system and calculation time is increased, making this method unfeasible for real-time applications.

Traditionally, VSM has always been the focus of voltage stability studies in static analysis of power systems [26]. However, the increase in load at extreme limits can also cause the emergence of low-frequency oscillation modes of electromechanical origin that can affect the angular stability of the system and, consequently, cause a blackout [27]. The system load level that causes the emergence of oscillation modes with low damping rates may occur before the load level that causes voltage collapse. These low oscillation modes can only be evaluated through linearization of the system’s dynamic equations, a procedure not performed in the VSM calculation [27].

Therefore, it is common to work with the Load Margin (LM) index instead of VSM. The LM consists of the difference in the system load level between the nominal case and the operating case where instability occurs, which could be voltage collapse in voltage stability studies or the oscillation modes of low damping rates in angular stability studies [27]. From the differential-algebraic equations of the system and the bifurcation theories of dynamic systems, it is possible to identify the LM of power systems. In dynamic systems, the Saddle-Node Bifurcation (SNB) is associated with the sudden loss of the equilibrium point and can be identified by the presence of a purely null eigenvalue in the linearized model of the system’s differential-algebraic equations [26]. In power systems, it is common to associate voltage collapse with SNB [26]. In dynamic systems, the Hopf Bifurcation (HB) is associated with the emergence of periodic orbits from the equilibrium point and can be identified by the presence of a pair of purely imaginary eigenvalues in the linearized model of the differential-algebraic equations of the system [26]. In power systems, it is common to associate oscillation modes with zero damping rate to HB [26].

In [28], the authors proposed a direct method to determine the LM associated with SNB and HB through the use of a determined set of equations and variables that represent the power system under the conditions of occurrence of SNB and HB. The preliminary results presented by the authors confirmed the feasibility of the method in correctly identifying the load margin at the voltage stability threshold or small signal stability associated with SNB or HB. After this proposal, different improvements were made over the years by the authors [29,30,31,32,33,34,35,36]. However, the challenges of this direct method still persist for real-time applications, such as the difficulty of converging initial conditions of variables far from the conditions of occurrence of SNB and HB, a complete and reliable model of power systems and many iterations may be necessary for the method to converge with a desirable error rate.

Machine learning is a data analysis method that allows the construction of analytical models in an automated way [37]. Machine learning techniques have been developed over the last few years for the most varied applications in power systems. Tools to calculate LM using machine learning techniques have been proposed. In [38], the authors use Artificial Neural Networks (ANN) to calculate VSM. In [39], the authors use Support Vector Machine (SVM) to determine the VSM. In [40], the authors use Decision Trees (DT) to determine the VSM. In [41], the authors use ANN and bio-inspired algorithms to calculate VSM. In [42], the authors combined Artificial Neural Networks and a network reduction technique to determine the VSM. In [43], the authors design and apply Radial Basis Function (RBF) neural network to determine VSM of power systems. In [44], the authors apply ANN and aim to determine a reduced set of measurements to be inputs to the ANN in order to correctly determine the VSM. In [45], the authors apply Convolutional Neural Network (CNN) to monitor the VSM of power systems in real-time. In [46], the authors use Graph Neural Networks (GNN) to calculate VSM with topology variations. In [47], the authors aim to apply energy functions in conjunction with ANN to calculate the VSM. In [48], the authors proposed a method based on Extreme Learning Machine (ELM) to determine the VSM of power systems. In [49], the authors propose a Genetic Algorithm based Support Vector Machine method to monitor the VSM of smart grids in real-time. In [50], the authors apply Recurrent Neural Networks (RNN) to determine the VSM for different typical scenarios of a power system. In [51], the authors use Random Forest to compute the VSM of smart grids. In [52], the authors apply decision tree to compute the VSM considering multiple operating conditions. The techniques are promising but the performance of these techniques may be unsatisfactory in large power systems due to the low generalization capacity of these techniques in critical situations. Most of these methods already proposed in the literature depend on a database of possible operating scenarios of the power system under study and scenarios not considered in the training stage can negatively affect the methods in decision making.

The generalization capacity of machine learning techniques can be improved by incorporating physical knowledge of the problem. When it comes to Neural Networks, there is a consensus that Physics-Informed Neural Networks (PINNs) can improve this generalization capacity through the incorporation of rules and/or physical laws in their formulation and design. In [53], the authors present a set of PINN applications in power systems. In [54], the authors use PINN to calculate the AC-Optimal Power Flow and the results showed good generalization capacity of the application. In [55], the authors apply PINN to determine the rotor angle and frequency of test systems. In [56], the authors use PINN to study the thermal dynamic behavior of power transformers. In [57], the authors apply PINN to identify power systems. In [58], the authors apply PINN for control-oriented modeling of existing buildings. In [59], the authors propose a framework based on PINDA to issue a Li-ion battery prediction. In [60], the authors use PINN to solve parametric magnetostatic problems. In [61], the authors propose a PINN-based framework for energy-efficient food production in existing buildings. In [62], the authors apply PINN to evaluate management and prognostics of Lithium-Ion Batteries. In [63], the authors use PINN to solve transient heat conduction problems. In [64], the authors apply PINN with other devices to solve the task of temperature field inversion of existing and new heat-source systems. In [65], the authors use PINN to model electric water heaters more efficiently. In [66], the authors apply PINN to solve heated tube equations. In [67], the authors apply PINN to reconstruct a good flow field. In [68], the authors apply PINN in river silting simulation. In [69], the authors use PINN for the heat equation. In [70], the authors propose a Physics-Guided Neural Network (PGNN) for Load Margin Assessment of Power Systems and the main difficulty encountered by the authors was the tuning of the PGNN parameters, which can compromise the results and performance of the proposed method. Thus, there is a recent set of PINN applications in the literature and preliminary results show superior performance to traditional Neural Networks obtained only with empirical knowledge and no physical knowledge. In the applications presented in the literature, different physical knowledge was applied in the Neural Networks training stage and aimed to adapt to the problem itself and which characteristics or physical laws must be obeyed.

Bio-inspired algorithms are optimization algorithms that draw on the principles and inspiration of nature’s biological evolution to develop new search tools for optimization problems. Different bio-inspired algorithms have been proposed such as Particle Swarm Optimization [71], Coati Optimization Algorithm [72], Pelican Optimization Algorithm [73], Marine Predators Algorithm [74], Electric Eel Foraging Optimization [75], Hippopotamus Optimization Algorithm [76], Several applications of bio-inspired algorithms in power systems have been proposed over the years such as wide-area damping control design [77,78], electricity theft detection [79], integrated energy system optimization [80], optimization of HVAC systems [81], power system stabilizer design [82], load dispatch for microgrid [83], energy management [84], load profile generation [85], power system state estimation [86], short-term hydrothermal scheduling [87], distributed power generation planning [88], reactive power optimization [89], maximum power point tracking [90], wind turbine placement [91], coordination of directional overcurrent relays [92], placement of electric vehicle charging station [93], optimal DG unit placement [94], power quality disturbances identification [95], optimal power flow [96]. The use of bio-inspired algorithms to tune traditional ANN is possible and some authors have already pointed out this benefit to improve the generalization capacity of the ANN.

Based on all of the above, this article proposes a design method for Physics-Informed Neural Networks (PINNs) whose parameters will be tuned by bio-inspired algorithms in an optimization model. Three bio-inspired algorithms existing in the scientific community were chosen to solve this optimization model: Particle Swarm Optimization (PSO) [71], Coati Optimization Algorithm (COA) [72] and Pelican Optimization Algorithm (POA) [73]. These three bio-inspired algorithms will be applied, compared and discussed in case studies that will be carried out on the IEEE 68-bus system composed of 16 generators and considering the N-1 criterion for disconnection of transmission lines. The main contributions of this article can be described as

• A new method for calculating the Load Margin (LM) of power systems using machine learning techniques;
• A new Neural Network architecture called Physics-Informed Neural Network (PINN) as a machine learning technique for the LM calculation problem;
• An optimization model was developed to tune the PINN weights in such a way as to minimize the training error given by the Root Mean Square Error (RMSE) index.
• Three different bio-inspired algorithms will be applied and evaluated in the proposed optimization model in both the training and testing stages of PINN.

This article has the following organization: Section 3 presents the PINN for evaluating the LM of power systems and which variables must be determined, Section 5 presents the four bio-inspired algorithms that will be applied separately to determine the PINN parameters, Section 4 presents the optimization model proposed to determine the PINN, Section 6 presents the case studies and discussions of the results achieved, Section 6 presents the conclusions of the article.

2. Direct Method to Obtain LM

The successful application of the LM monitoring method based on Physics-Informed Neural Networks requires the construction of a very consistent database that is representative of the set of typical scenarios of a power system. With this in mind, the direct method proposed in [28] will be applied in this research to build the database. The benefit of this method is that it provides an LM that meets voltage stability and small-signal stability thresholds. The direct method consists of constructing and solving a determined set of equations and variables under operating point conditions at the threshold of small-signal stability and voltage stability associated with Hopf Bifurcation and Saddle-Node Bifurcation, respectively. This determined system of equations and variables is formed by the following vector of variables

$$K={\left[\begin{array}{cccccccc}{x}^T& y^T& \mu & v_R^T& v_I^T& w_R^T& w_I^T& {\omega }_0\end{array}\right]}^T$$

(1)

where the vector that describes the state variables such as angle and speed of the generators is x, the vector that describes the algebraic variables such as voltage and current is y, the system load level is $\mu$ and is the information of interest in this system of equations because it defines the LM, the imaginary and real part of the eigenvector relative to the vector x are $v_I$ and $v_R$ respectively, the imaginary and real part of the eigenvector relative to the vector y are $w_I$ and $w_R$ respectively, and the frequency of the pair of purely imaginary eigenvalues of HB is described by ${\omega }_0$.

The system proposed in [28] is made up of several equations. The first two equations involve the existence and uniqueness of an equilibrium point of the power system and is related to the power flow [97]. These two equations are described in (2) and (3) by differential and algebraic functions. Thus, the convergence of the method provides the system’s equilibrium point at the stability thresholds with the voltage, current and power information on the system buses. The other equations require knowledge of linear algebra rules such as eigenvalues and eigenvectors and how this knowledge is associated with voltage stability and small-signal stability thresholds. The Hopf bifurcation is associated with the small-signal stability threshold as there is the emergence of periodic orbits around the equilibrium point of the dynamic equations. In linear algebra concepts, HB is linked to the presence of a pair of purely imaginary eigenvalues with zero real part $({\lambda }_{1,1}=\pm {\omega }_0)$ and, thus, Equations (6) and (7) can be used relative to this purely imaginary eigenvalue. The Saddle-Node bifurcation is associated with the voltage stability threshold as there is a sudden loss of the equilibrium point of the dynamic equations. In concepts of linear algebra, the SNB is linked to the presence of a single eigenvalue that is purely null $(\lambda =0)$ and, therefore, Equations (4) and (5) relative to this null eigenvalue can be used. Furthermore, two Equations (8) and (9) associated with the eigenvectors were used to obtain a system with the same number of variables and equations. The larger the size of the power system, the greater the number of variables to be determined. Therefore, systems with many generators, transmission lines and buses can make the process of searching for the values of the variables in this system difficult.

$$f(x,y,d,\mu )=0$$

(2)

$$g(x,y,d,\mu )=0$$

(3)

$$J·\left[\begin{array}{c}{v}_R\\ w_R\end{array}\right]=0$$

(4)

$$J·\left[\begin{array}{c}{v}_I\\ w_I\end{array}\right]=0$$

(5)

$$J·\left[\begin{array}{c}{v}_R\\ w_R\end{array}\right]+{\omega }_0·\left[\begin{array}{c}{v}_I\\ 0\end{array}\right]=0$$

(6)

$$J·\left[\begin{array}{c}{v}_I\\ w_I\end{array}\right]-{\omega }_0·\left[\begin{array}{c}{v}_R\\ 0\end{array}\right]=0$$

(7)

$$\left[\begin{array}{cc}{v}_R^T& w_R^T\end{array}\right]·\left[\begin{array}{c}{v}_I\\ w_I\end{array}\right]=0$$

(8)

$$\left[\begin{array}{cc}{v}_R^T& w_R^T\end{array}\right]·\left[\begin{array}{c}{v}_R\\ w_R\end{array}\right]+\left[\begin{array}{cc}{v}_I^T& w_I^T\end{array}\right]·\left[\begin{array}{c}{v}_I\\ w_I\end{array}\right]=1$$

(9)

This complete system can be solved by any method offline. The authors [28] employed Newton method to solve this set of equations for various scenarios. However, the same authors found obstacles to applying the method online: (i) difficulty in converging to initial conditions far from the correct solution, (ii) need for precise knowledge of all model parameters, (iii) uncertain number of iterations for final convergence. This method is ideal for building a database of possible power system scenarios, but it is not an ideal method for real-time applications. Therefore, a machine learning technique called Physics-Informed Neural Network was proposed and is presented in detail below.

3. Physics-Informed Neural Network (PINN)

Em [70], the author proposed a Physics-Informed Neural Network (PINN) or Physics-Guided Neural Network (PGNN) to calculate the load margin of power systems considering voltage stability and small-signal stability thresholds. Traditional Artificial Neural Networks can perform poorly due to overfitting and low generalization capacity for scenarios not considered in the training stage. With this in mind, the author [70] incorporated physical characteristics of the operation of power systems as a requirement for the existence of a convergent power flow from the output data of the second ANN during the training stages to obtain better Root Mean Square Error (RMSE) results.

The author proposed PINN whose architecture is described in [70]. PINN is composed of two ANNs, the output data of one ANN is the input data of the second ANN. The input vector ${\mathbf{X}}_{\mathbf{1}}$ of the first ANN is described in (10) and is composed of voltage magnitude $\left(V_{i,0}\right)$ and angle measurements $\left({\theta }_{i,0}\right)$ collected by PMUs installed on system buses $(i=1,\cdots ,N_b)$. PMUS are measuring equipment developed to measure voltage and current with high precision and high sampling rates, typically 60 samples/s. Initial attempts have shown that voltage measurements promote better PINN performance in correctly calculating the LM of smart grids. This first ANN provides as output the vector ${\mathbf{Y}}_{\mathbf{1}}$ described in (11) and composed of the Load Margin (LM) measurement $\mu$ at the stability threshold and the magnitude $\left(V_{i,B}\right)$ and voltage angle $\left({\theta }_{i,B}\right)$ measurements of the buses at the stability threshold. In the database, the output data $\left({\mathbf{Y}}_{\mathbf{1}}\right)$ from the first ANN is known and wil be compared with the estimated values $\left({\widehat{\mathbf{Y}}}_{\mathbf{1}}\right)$ and the objective is for these values to be as close as possible. Thus, the errors associated with the estimated value and the database value must be minimized by the first error function $l_1$ given by (1)

$${\mathbf{X}}_{\mathbf{1}}=\left[\begin{array}{ccccc}{V}_{1,0}& {\theta }_{1,0}& \cdots & V_{N_b,0}& {\theta }_{N_b,0}\end{array}\right]$$

(10)

$${\mathbf{Y}}_{\mathbf{1}}=\left[\begin{array}{cccccc}\mu & V_{1,B}& {\theta }_{1,B}& \cdots & V_{N_b,B}& {\theta }_{N_b,B}\end{array}\right]$$

(11)

$$l_1={\widehat{\mathbf{Y}}}_{\mathbf{1}}-{\mathbf{Y}}_{\mathbf{1}}$$

(12)

These output measurements at the stability threshold must present data consistent with the existence of a convergent power flow. With this in mind, the author used a second ANN to try to imitate the calculation of a power flow [97]. The inputs of this second ANN are the ${\mathbf{X}}_{\mathbf{2}}$ vector itself described in (13) and the output is the ${\mathbf{Y}}_{\mathbf{2}}$ vector described in (14) and composed of the estimates of active and reactive power at the stability threshold. From the system bus voltage data, it is possible to calculate the active $\left(P_{i,B}\right)$ and reactive $\left(Q_{i,B}\right)$ power of the buses using the power flow equations described by Equations (15) and (16). Again, the database has the actual active and reactive power data. Thus, the second ANN must be designed in the training stage by minimizing the error function $\left(l_2\right)$ described by Equation (17). Thus, in the process of minimizing the error function l${ }_2$, the programmer is physically guaranteeing the need for a convergent power flow. The application of this second function proves useful in the process of generalizing this machine learning technique in dealing with scenarios or situations not previously studied in the training stage with a database.

$${\mathbf{X}}_{\mathbf{2}}=\left[\begin{array}{ccccc}{V}_{1,B}& {\theta }_{1,B}& \cdots & V_{N_b,B}& {\theta }_{N_b,B}\end{array}\right]$$

(13)

$${\mathbf{Y}}_{\mathbf{2}}=\left[\begin{array}{ccccc}{P}_{1,B}& Q_{1,B}& \cdots & P_{N_b,B}& Q_{N_b,B}\end{array}\right]$$

(14)

$$P_{i,B}=\sum _{k=1}^{N_b}{V}_{i,B}{V}_{k,B}\left(G_{i,k}cos({\theta }_{i,B}-{\theta }_{k,B})+B_{i,k}sin({\theta }_{i,B}-{\theta }_{k,B})\right)$$

(15)

$$Q_{i,B}=\sum _{k=1}^{N_b}{V}_{i,B}{V}_{k,B}\left(G_{i,k}sin({\theta }_{i,B}-{\theta }_{k,B})-B_{i,k}cos({\theta }_{i,B}-{\theta }_{k,B})\right)$$

(16)

$$l_2={\widehat{\mathbf{Y}}}_{\mathbf{2}}-{\mathbf{Y}}_{\mathbf{2}}$$

(17)

PINN will be trained from a database. The training of the two ANNs will be done with the objective of minimizing two error functions $l_1$ and $l_2$ described in Equation (18). The first error function $l_1$ is typical of training ANNs and trying to find a map that relates inputs and output. The second error function $L_2$ attempts to incorporate physical knowledge into the training process to ensure the existence of a power flow. Weights ${\alpha }_1$ and ${\alpha }_2$ must be chosen carefully to ensure effective PINN training.

$$l={\alpha }_1·l_1+{\alpha }_2·l_2={\alpha }_1·({\widehat{\mathbf{Y}}}_{\mathbf{1}}-{\mathbf{Y}}_{\mathbf{1}})+{\alpha }_2·({\widehat{\mathbf{Y}}}_{\mathbf{2}}-{\mathbf{Y}}_{\mathbf{2}})$$

(18)

In [70], the author used Levenberg-Marquardt [98] as a training algorithm for PINN whose stopping criterion is the maximum number of epochs defined as 500. The preliminary results were satisfactory as they met previously established expectations. However, the PINN design can be improved with more effective training techniques capable of providing a more adequate training error in less time. In this research, three bio-inspired algorithms will be used as training algorithms for PINN and this method will be described in the next section.

4. Proposed Method

The objective of this research is the design of a PINN capable of monitoring the LM of power systems equipped with PMUs taking into account voltage stability and small-signal stability requirements with the lowest possible error between measured data and estimated data. Therefore, the training process of the PINN proposed in this research will be carried out through the iterative resolution of the optimization model by different bio-inspired algorithms. Thus, the weights of the two ANNs described by the vector $\mathbf{W}$ in (19) that make up the PINN will be tuned through an optimization model using bio-inspired algorithms. The motivation for applying bio-inspired algorithms is their ability to solve optimization problems of different complexity levels with good success rates. Each bio-inspired algorithm is composed of its own operators who try to imitate behaviors of animals in the wild such as hunting and escape processes. Thus, the vector $\mathbf{W}$ consists of the variable of the optimization problem.

$$\mathbf{W}=\left[w_{a,b}^{c,d}\right]$$

(19)

where $a=1,\cdots ,N_a$ represents the number of neurons in one layer and $b=1,\cdots ,N_b$ represents the number of neurons in the other layer, $c=1,2$ represents the number of connections between layers and $d=1,2$ represents one of the two ANNs. In this research, the first ANN will have an input layer composed of $2·N_{bPMUs}$ neurons where $N_{bPMUs}$ is the number of buses with PMUs, a hidden layer with $N_{bPMUs}$ neurons and an output layer with $2·N_{bPMUs}+1$ neurons. The second ANN will have an input layer composed of $2·N_{bPMUs}$ neurons, a hidden layer with $N_{bPMUs}$ neurons and an output layer with $2·N_{bPMUs}$ neurons.

The objective function $\left(F_{obj}\right)$ of this optimization model is described in (20) and consists of minimizing the training error l based on a fixed architecture of the two ANNs. The error functions associated with empirical knowledge and physical knowledge will then be minimized simultaneously in this process to guarantee satisfactory performance indices in the PINN testing stage.

$$F_{obj}\left(\mathbf{W}\right)=l={\alpha }_1·l_1+{\alpha }_2·l_2={\alpha }_1·({\widehat{\mathbf{Y}}}_{\mathbf{1}}-{\mathbf{Y}}_{\mathbf{1}})+{\alpha }_2·({\widehat{\mathbf{Y}}}_{\mathbf{2}}-{\mathbf{Y}}_{\mathbf{2}})$$

(20)

The optimization model can then be defined as a minimization problem as described in (21) where the weights of the two ANNs that make up the PINN must be found that minimizes the two error functions associated with empirical knowledge and physical knowledge.

$$\begin{array}{cc}\mathrm{Find}\hfill & \mathbf{W}=\left[w_{a,b}^{c,d}\right]\hfill \\ \mathrm{Minimize}\hfill & F_{obj}\left(\mathbf{W}\right)=l={\alpha }_1·l_1+{\alpha }_2·l_2={\alpha }_1·({\widehat{\mathbf{Y}}}_{\mathbf{1}}-{\mathbf{Y}}_{\mathbf{1}})+{\alpha }_2·({\widehat{\mathbf{Y}}}_{\mathbf{2}}-{\mathbf{Y}}_{\mathbf{2}})\hfill \end{array}$$

(21)

The stopping criterion for this optimization model will be the number of epochs dedicated to training the PINN. Furthermore, the weights ${\alpha }_1$ and ${\alpha }_2$ must be chosen carefully. If ${\alpha }_1=1$ and ${\alpha }_2=0$, it is a traditional ANN because only empirical knowledge will be used in the training process. If ${\alpha }_1=0.5$ and ${\alpha }_2=0.5$, both empirical knowledge and physical knowledge will be used in the same proportion for the PINN training process and will be beneficial in obtaining a PINN with greater ability to generalize across scenarios unknown and not considered in the training stage. It is also important to highlight that the PINN input signals will be from a select group of buses $\left(N_b\right)$ with PMUs because the greater the number of PINN inputs, the greater the number of problem variables would be and could hinder the convergence of the optimization model.

The optimization model can be solved by any metaheuristic and it was decided to use three different bio-inspired algorithms for its resolution: Particle Swarm Optimization (PSO), Coati Optimization Algorithm (COA) and Pelican Optimization Algorithm (POA). These three algorithms will be applied and evaluated separately in the minimization problem. The next section describes each of these three bio-inspired algorithms.

After the convergence of the proposed optimization model, the PINN weights were defined and the training stage was successful. The next step is to carry out the testing stage and evaluate the ability of the proposed method to deal with scenarios not considered in the training stage. Three indices were applied in the PINN test stage: Root Mean Square Error (RMSE), Mean Square Error (MSE) and Mean Absolute Percentage Error (MAPE) whose formulas are detailed below

$$RMSE=\sqrt{\frac{1}{n}\sum _i^n{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(22)

$$MSE=\frac{1}{n}\sum _i^n{\left(y_i-{\widehat{y}}_i\right)}^2$$

(23)

$$MAPE=\frac{1}{n}\sum _i^n\frac{|y_i-{\widehat{y}}_i|}{y_i}$$

(24)

where $y_i$ is the true value, ${\widehat{y}}_i$ is the predicted value and n is the number of samples present in the database for the PINN test stage.

5. Bio-Inspired Algorithms

Three different bio-inspired algorithms were applied and evaluated separately in the proposed method in the search for the best parameters in PINN in order to correctly determine the LM of power systems equipped with PMUs. These three bio-inspired algorithms are Particle Swarm Optimization (PSO), Coati Optimization Algorithm (COA) and Pelican Optimization Algorithm (POA). Each of them are succinctly described below.

5.1. Particle Swarm Optimization (PSO)

The Particle Swarm Optimization (PSO) metaheuristic attempts to imitate the behavior of particles in nature [71]. This is an old metaheuristic that is well used by the scientific community. It is a bio-inspired algorithm with few operations, simple to implement and effective in obtaining satisfactory solutions to optimization problems. There are several applications of PSO in the literature [99,100,101,102,103]. The positions of the particles are the variables of the search problem and are updated according to the velocity values of each particle, which can be mathematically formulated as

$$v_i^{k+1}=\omega v_i^k+c_1{r}_1(x_i^L-x_i^k)+c_2{r}_2(x^G-x_i^k)$$

(25)

$$x_i^{k+1}=x_i^k+v_i^{k+1}$$

(26)

where $x_i$ is the position of the i-th particle, $v_i$ is the velocity of the i-th particle, $\omega$ is a parameter that evaluates the percentage of speed that must be considered for the next iteration, $r_1$ and $r_2$ are random values of a uniform distribution between 0 and 1, $c_1$ and $c_2$ are fixed parameters chosen by the programmer, $x^G$ is the best position so far according to the objective function, $x_i^L$ the best position for each particle so far. Particle positions are updated according to iterations and respecting the maximum epoch limit $N_E$$(k=1,\cdots ,N_E)$. After reaching the maximum iteration limit, the final solution is the vector $x^G$. Over the years after the release of this bio-inspired algorithm, different improvements have been applied to traditional PSO through initialization processes, operator changes or hybrid approaches with other existing methods. Found several improvements in different optimization problems. In this research, the traditional PSO available in [71] was applied.

5.2. Coati Optimization Algorithm (COA)

The authors [72] created the metaheuristic called Coati Optimization Algorithm (COA) which mimics the behavior of coatis hunting and attacking iguanas in nature. Some recent applications of COA in optimization problems are described in [104,105,106,107,108]. Two typical coatis behaviors are transformed in two steps of the algorithm proposed by the authors: (i) processes of hunting and attacking iguanas by coatis and (ii) process of escape of coatis from predators.

The positions of the coatis in nature are the variables of the search problem of this metaheuristic. These positions assume random values according to a uniform distribution with minimum and maximum limits as follows

$$X_i:x_{i,j}=l{b}_j+r·(u{b}_j-l{b}_j)$$

(27)

where $X_i$ is the position of the i-th coati, $x_{i,j}$ is the value of the j-th decision coati, for $i=1,\cdots ,N$ and $j=1,\cdots ,m$, N is the number of coatis, m is the number of decision variables, r is a random value arising from a uniform distribution in the interval [0,1], the maximum and minimum values of the decision variables are called $u{b}_j$ and $l{b}_j$ respectively. Considering all coatis in the search problem, we have a population of coatis that can be defined by the following matrix X

$$X=\left[\begin{array}{c}{X}_1\\ ⋮\\ X_i\\ ⋮\\ X_N\end{array}\right]=\left[\begin{array}{ccccc}{x}_{1,1}& \cdots & x_{1,j}& \cdots & x_{1,m}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ x_{i,1}& ·& x_{i,j}& \cdots & x_{i,m}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ x_{N,1}& ·& x_{N,j}& \cdots & x_{N,m}\end{array}\right]$$

(28)

In this matrix of coatis, each line of X represents a coati and a possible candidate in the problem of searching for an optimal value of an objective function. Thus, throughout the iterations, each coati position has its objective function $F\left(\right)$ calculated as follows

$$F\left(X\right)=\left[\begin{array}{c}{F}_1\\ ⋮\\ F_i\\ ⋮\\ F_N\end{array}\right]=\left[\begin{array}{c}F\left(X_1\right)\\ ⋮\\ F\left(X_i\right)\\ ⋮\\ F\left(X_N\right)\end{array}\right]$$

(29)

Throughout the iterations, the best candidate will be the one that presents the best objective function value. The operators of this metaheuristic will always change the values of these positions in order to have better candidates according to the objective function. The interactions of coatis in the wild were interpreted by the authors through two steps that are described in detail in the following sections.

5.2.1. Phase 1: Hunting and Attacks on Iguanas

In this first phase, the positions of the coatis will be updated imitating the process of attacking the iguanas, their prey. In the process of updating positions, the best candidate will be in the iguana position representing the coati with the best attack and performance. Half of the coatis will attack the iguanas while the other half wait. Half of the coatis that will attack the iguanas can be represented by the following modeling

$$X_i^{P1}:x_{i,j}^{P1}=x_{i,j}+r·(Iguan{a}_j-I·x_{i,j})$$

(30)

for $i=1,\cdots ,N/2$ and $j=1,\cdots ,m$.

When the iguana, prey of the coatis, falls to the ground in a random position, the other half of the coatis update two positions according to the following formulas

$$Iguan{a}^G:Iguan{a}_j^G=l{b}_j+r·(u{b}_j-l{b}_j)$$

(31)

$$X_i^{P1}=x_{i,j}^{P1}=\left\{\begin{array}{cc}{x}_{i,j}+r·(Iguan{a}_j^G-I·x_{i,j}),\hfill & \mathrm{If}{F}_{Iguan{a}^G}<F_i\hfill \\ x_{i,j}+r·(·x_{i,j}-Iguan{a}_j^G),\hfill & \mathrm{else}\hfill \end{array}\right.$$

(32)

for $j=1,\cdots ,m$ and $i=N/2+1,\cdots ,N$.

These new position updates $\left(X_i^{P1}\right)$ will be accepted into the population matrix $\left(X_i\right)$ only if they actually improve the objective function. Otherwise, the position matrix remains at the current values. This analysis can be summarized as

$$X_i=\left\{\begin{array}{cc}{X}_i^{P1},\hfill & \mathrm{If}{F}_i^{P1}<F_i\hfill \\ X_i,\hfill & \mathrm{else}\hfill \end{array}\right.$$

(33)

5.2.2. Phase 2: Escape from Predators

The second stage of the process of updating the coatis’ positions attempts to imitate the behavior of coatis in the presence of predators. If a predator tries to attack a coati, the coati changes its position to escape the predator. In order for this imitation to be successful, a random position is generated for each coati close to its current position according to the following formulation

$$X_i^{P2}:x_{i,j}^{P2}=x_{i,j}+(2-2r)(l{b}_j^{local}+r(u{b}_j^{local}-l{b}_j^{local}))$$

(34)

where $i=1,\cdots ,N$, $j=1,\cdots ,m$ and

$$l{b}_j^{local}=\frac{l{b}_j}{t}$$

(35)

$$u{b}_j^{local}=\frac{u{b}_j}{t}$$

(36)

where $t=1,\cdots ,T$.

These new calculated positions will be effectively used only if they improve the objective function values. Otherwise, the positions of the coatis will be maintained in this iteration. This analysis can be summarized as

$$X_i=\left\{\begin{array}{cc}{X}_i^{P2},\hfill & \mathrm{If}{F}_i^{P2}<F_i\hfill \\ X_i,\hfill & \mathrm{else}\hfill \end{array}\right.$$

(37)

Thus, COA is a bio-inspired algorithm with recent development and composed of several operations that aim to find the best solution to a minimization or maximization problem. It should be noted that there are few COA parameters to be defined by the programmer and this can facilitate the implementation of the algorithm.

5.3. Pelican Optimization Algorithm (POA)

The authors [73] proposed Pelican Optimization Algorithm (POA) that tries to imitate the behavior of pelicans during hunting. Some recent applications of COA in optimization problems are described in [109,110,111,112,113,114]. The positions of the pelican will be the variables of this problem of searching for the optimal value of the objective function. The position of each pelican will be randomly generated according to a uniform distribution with minimum and maximum limits respecting the following formulation

$$x_{i,j}=l_j+rand·(u_j-l_j)$$

(38)

where $x_{i,j}$ is the position of the i-th pelican for $i=1,\cdots ,N$ and $j=1,\cdots ,m$, N is the number of pelicans, m is the number of decision of variables, $rand$ is a random value from a uniform distribution in the interval [0,1], the maximum and minimum limits of the decision variables are $u_j$ and $l_j$ respectively.

Each pelican represented by (38) can be grouped into a population and thus we have a population matrix X where each line represents a pelican and each column represents a decision variable. This matrix X is formulated as

$$X=\left[\begin{array}{c}{X}_1\\ ⋮\\ X_i\\ ⋮\\ X_N\end{array}\right]=\left[\begin{array}{ccccc}{x}_{1,1}& \cdots & x_{1,j}& \cdots & x_{1,m}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ x_{i,1}& \cdots & x_{i,j}& \cdots & x_{i,m}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ x_{N,1}& \cdots & x_{N,j}& \cdots & x_{N,m}\end{array}\right]$$

(39)

Each pelican in this population is a candidate for the search problem and, therefore, presents an objective function value for the optimization problem. The objective function $\left(F\right(\left)\right)$ calculated in the following way defines whether each pelican will have its position changed or not for the next iteration.

$$F=\left[\begin{array}{c}{F}_1\\ ⋮\\ F_i\\ ⋮\\ F_N\end{array}\right]=\left[\begin{array}{c}F\left(X_1\right)\\ ⋮\\ F\left(X_i\right)\\ ⋮\\ F\left(X_N\right)\end{array}\right]$$

(40)

This metaheuristic presents two distinct stages for updating the pelicans’ positions: (i) hunting for prey, and (ii) flying on the water surface. Each of these two stages are described in the following subsections.

5.3.1. Phase 1: Hunting for Prey

In this phase, pelicans find the position of the prey and move to this position. The mathematical formulation of this process is given by

$$x_{i,j}^{P1}=\left\{\begin{array}{cc}{x}_{i,j}+rand·(p_j-I·x_{i,j}),\hfill & \mathrm{If}{F}_p<F_i\hfill \\ x_{i,j}+rand·(x_{i,j}-p_j),\hfill & \mathrm{else}\hfill \end{array}\right.$$

(41)

where $p_j$ is the position of the prey, $x_{i,j}^{P1}$ is the new position of the pelican, I is a parameter that can take the value 1 or 2, $F_p$ is the objective function of the new position.

After calculating the new positions $\left(X_i^{P1}\right)$ and objective functions, the new positions will be effectively applied if they improve the values of the current objective functions, otherwise they maintain the current positions. This analysis then respects the following formulation

$$X_i=\left\{\begin{array}{cc}{X}_i^{P1},\hfill & \mathrm{If}{F}_i^{P1}<F_i\hfill \\ X_i,\hfill & \mathrm{else}\hfill \end{array}\right.$$

(42)

5.3.2. Phase 2: Flying on the Water Surface

In this second phase, the pelicans go to the surface of the water, they spread their wings on the surface of the water to attract the fish upwards and thus attack the prey. This behavior is described by the following mathematical formulation

$$x_{i,j}^{P2}=x_{i,j}+R·\left(1-\frac{t}{T}\right)·(2·rand-1)·x_{i,j}$$

(43)

After calculating the new positions $\left(X_i^{P2}\right)$ and objective functions, the new positions will be effectively applied if they improve the values of the current objective functions, otherwise they maintain the current positions. This analysis then respects the following formulation

$$X_i=\left\{\begin{array}{cc}{X}_i^{P2},\hfill & \mathrm{If}{F}_i^{P2}<F_i\hfill \\ X_i,\hfill & \mathrm{else}\hfill \end{array}\right.$$

(44)

Therefore, POA is also a bio-inspired algorithm with recent development and composed of several operations whose objective is to find the best solution to an optimization problem. Furthermore, it is also worth highlighting that there are few POA parameters to be defined by the future programmer and this can facilitate and speed up the implementation of this algorithm.

6. Results

The proposed method to determine LM using PINN and bio-inspired algorithms will be evaluated in this section through its application in the IEEE 68-bus system. The nominal operating point of this test system is described in [115]. This test system is evaluated for the most varied studies on power systems. It is composed of 16 synchronous generators and 68 buses and, therefore, static and dynamic analyzes are possible to be carried out. Loads are modeled as constant power type. Now that the test system has been defined, the proposed optimization model-based method for determining PINN weights for calculating LM in power systems can be evaluated from different aspects.

The first step in applying the proposed method is the creation of a reliable operating database for the test system. From the nominal operating point of the test system described in [115], new 10,000 operating points were obtained through random variation of active power and reactive power of the system load level between $-15\%$ and +15%. The load growth direction adopted in this research consisted of increasing the load level in the same proportion while keeping the power factor constant. The active power generation redispatch model consists of each generator in the system providing an active power value proportional to the inertia constant with a greater participation of generators with greater inertia. Thus, for each of these new 10,000 operating points, an LM is obtained by applying the direct method described in Section 2, totaling 10,000 LMs. Figure 1 presents a histogram with the distribution of these 10,000 load margins of the test system. Some important aspects can be evaluated from the results of this database histogram. There is a uniform set of cases with different load margins (LM) that are in the range between 0% and 15%. This facilitates the division of cases for the training and testing stages as there are very representative groups. It is also observed that there are cases that provide very low and very high LM. High LMs are desirable for safe and reliable system operation. Therefore, detecting low LM is essential to inform the system operator to take control actions to try to effectively increase LM.

From the database, 70% of the data was applied to the PINN training stage, and 30% of the data was applied to the PINN testing stage. PINN input data are measurements of voltage magnitude and angle coming from PMUs that may be subject to noise. Thus, an error of 1% standard deviation of a Gaussian distribution was incorporated into the measurements.

The three bio-inspired algorithms Particle Swarm Optimization (PSO), Coati Optimization Algorithm (COA) and Pelican Optimization Algorithm (POA) were applied to the optimization model with the parameters defined in the original articles in a population defined in 40 particles, coatis or pelicans, and a maximum number of epochs defined as 5000. In the objective function, the weights assumed the following values ${\alpha }_1=0.5$ and ${\alpha }_2=0.5$. Thus, there is an equal proportion between both error functions during the training process and the PINN to be tuned will have an equal weight of empirical knowledge and physical knowledge. The test system has 68 buses and it was assumed that 20 buses are equipped with a PMU and collect voltage measurements that will be used by the PINN to be designed by the proposed optimization problem. The PINN training step was executed 100 times to evaluate the uniformity of convergence.

Table 1. Results of 100 simulations of the optimization model for designing a PINN by different training algorithms.

Algorithm	Minimum RMSE	Average RMSE	Maximum RMSE	Standard Deviation
Levenberg-Marquard	$5.452\times {10}^{-3}$	$6.128\times {10}^{-3}$	$6.889\times {10}^{-3}$	$4.247\times {10}^{-4}$
Particle Swarm Optimization (PSO)	$4.621\times {10}^{-3}$	$5.226\times {10}^{-3}$	$5.778\times {10}^{-3}$	$3.539\times {10}^{-4}$
Coati Optimization Algorithm (COA)	$4.293\times {10}^{-3}$	$4.854\times {10}^{-3}$	$5.396\times {10}^{-3}$	$3.149\times {10}^{-4}$
Pelican Optimization Algorithm (POA)	$4.016\times {10}^{-3}$	$4.556\times {10}^{-3}$	$5.088\times {10}^{-3}$	$3.034\times {10}^{-4}$

presents the root mean square error (RMSE) results achieved in the PINN training stage for different bio-inspired algorithms and the traditional Levenberg-Marquard method for training ANNs.

Now, the PINNs designed by different algorithms will be evaluated in the testing stage for different scenarios not evaluated in the training stage.

Table 2. RMSE results of the test stage of the 100 PINNs obtained by different training algorithms.

Algorithm	Minimum RMSE	Average RMSE	Maximum RMSE	Standard Deviation
Levenberg-Marquard	$8.572\times {10}^{-3}$	$9.779\times {10}^{-3}$	$10.961\times {10}^{-3}$	$7.118\times {10}^{-4}$
Particle Swarm Optimization (PSO)	$5.632\times {10}^{-3}$	$6.368\times {10}^{-3}$	$7.091\times {10}^{-3}$	$4.333\times {10}^{-4}$
Coati Optimization Algorithm (COA)	$5.344\times {10}^{-3}$	$6.010\times {10}^{-3}$	$6.707\times {10}^{-3}$	$4.056\times {10}^{-4}$
Pelican Optimization Algorithm (POA)	$5.094\times {10}^{-3}$	$5.682\times {10}^{-3}$	$6.276\times {10}^{-3}$	$3.846\times {10}^{-4}$

,

Table 3. MSE results of the test stage of the 100 PINNs obtained by different training algorithms.

Algorithm	Minimum MSE	Average MSE	Maximum MSE	Standard Deviation
Levenberg-Marquard	$7.348\times {10}^{-5}$	$9.613\times {10}^{-5}$	$12.014\times {10}^{-5}$	$1.3954\times {10}^{-5}$
Particle Swarm Optimization (PSO)	$3.172\times {10}^{-5}$	$4.074\times {10}^{-5}$	$5.023\times {10}^{-5}$	$5.513\times {10}^{-6}$
Coati Optimization Algorithm (COA)	$2.856\times {10}^{-5}$	$3.629\times {10}^{-5}$	$4.499\times {10}^{-5}$	$4.899\times {10}^{-6}$
Pelican Optimization Algorithm (POA)	$2.594\times {10}^{-5}$	$3.243\times {10}^{-5}$	$3.939\times {10}^{-5}$	$4.379\times {10}^{-6}$

and

Table 4. MAPE results of the test stage of the 100 PINNs obtained by different training algorithms.

Algorithm	Minimum MAPE	Average MAPE	Maximum MAPE	Standard Deviation
Levenberg-Marquard	$3.704\times {10}^{-4}$	$4.842\times {10}^{-4}$	$5.969\times {10}^{-4}$	$6.875\times {10}^{-5}$
Particle Swarm Optimization (PSO)	$1.591\times {10}^{-4}$	$2.049\times {10}^{-4}$	$2.494\times {10}^{-4}$	$2.605\times {10}^{-5}$
Coati Optimization Algorithm (COA)	$1.443\times {10}^{-4}$	$1.828\times {10}^{-4}$	$2.244\times {10}^{-4}$	$2.310\times {10}^{-5}$
Pelican Optimization Algorithm (POA)	$1.183\times {10}^{-4}$	$1.566\times {10}^{-4}$	$1.961\times {10}^{-4}$	$2.177\times {10}^{-5}$

provide the RMSE, MSE and MAPE values achieved in the test stage by the different algorithms.

The results achieved in the PINN training and testing stages obtained by different training algorithms allow the following evaluations:

• The RMSE index results from PINN training by the optimization model show that the three bio-inspired algorithms provided the least values compared to the traditional Levenberg-Marquard method. The bio-inspired algorithm called Pelican Optimization Algorithm (POA) was developed recently and provided the PINN with lowest RMSE values.
• The RMSE, MSE and MAPE indices obtained from the PINN test stage again show the superior performance of the bio-inspired algorithms compared to the traditional Levenberg-Marquard method. Furthermore, it is observed that the index results are low, showing the good ability of PINN to correctly calculate the LM of power systems for data not considered in the training stage.
• Only voltage measurements (magnitude and angle) of some system buses equipped with PMUs were sufficient to provide PINN with low errors in the test stage for scenarios not evaluated in the training stage. Therefore, the installation or selection of some PMUs is sufficient for the correct calculation of the LM of power systems.
• The presence of two ANNs improved the performance of PINN in determining the LM of power systems with greater precision. Incorporating physical aspects of power systems such as ensuring a convergent power flow is beneficial for PINN to have greater generalization capacity and thus provide more accurate results for unseen scenarios.
• The optimization model for the training stage considered equal weights of the error functions, indicating that empirical knowledge and physical knowledge should have the same weight in the search process for PINN parameters. This consideration of equal weights provided satisfactory results of the training and testing errors of the PINN projected by the different methods.

7. Conclusions

This work proposed an optimization model for the design of a Physics-Informed Neural Network (PINN) for monitoring the Load Margin of power systems equipped with PMUs in some system buses. Any bio-inspired algorithm can be applied to the optimization model and three were evaluated in this research: Particle Swarm Optimization (PSO), Coati Optimization Algorithm (COA) and Pelican Optimization Algorithm (POA). Some conclusions can be drawn from the work developed in this article:

• The results achieved by the proposed method show low error values compared to the traditional training Levenberg-Marquard method. Bio-inspired algorithms proved to be able to find PINN parameters with better efficiency than the traditional method.
• The use of voltage data (magnitude and angle) provided by PMUs of some system buses is sufficient to monitor the Load Margin (LM) of power systems that meet the requirements of small-signal stability and voltage stability. Therefore, a massive number of PMUs are not needed in the system, just a few strategically located PMUs are sufficient for the satisfactory performance of the proposed method.
• Incorporating physical characteristics of power system operation such as power flow in the training process increases the generalization capacity of PINN and promotes greater precision in obtaining more precise LMs for scenarios not considered in the training stage.