Linearized-Trajectory Model-Predictive Controller for Improving Microgrid Short-Term Stability With Real-Time Implementation

A model-predictive controller (MPC) for synchronous generator-based microgrids is introduced to enhance voltage, frequency, and transient stability in fast timescales. The centralized controller uses mathematical models of power system dynamics to predict the evolution of system states over a finite horizon based on information from local state observers and relays. The MPC dynamically adjusts the setpoints of existing primary controls to maximize system stability. To address the computational challenges of a real-time MPC implementation due to nonlinear power system dynamics and short sampling intervals, a trajectory linearization technique is applied to the MPC formulation. The proposed controller is validated via a hardware-in-the-loop (HIL) testbed. The testbed includes an implementation of the controller in hardware, which interfaces with a high-fidelity microgrid simulation model running in a real-time simulator. State observers based on the Extended Kalman Filter (EKF), which provide the controller with estimates of the system state, are also implemented. The HIL testbed is used to verify the real-time operation of the controller and evaluate its performance under modeling uncertainties.


I. INTRODUCTION
E NSURING system stability is a crucial aspect of the safe and reliable operation of power systems.Herein, we focus on short-term stability, where the primary controls of Dakota Hamilton, Loraine Navarro, Dionysios Aliprantis, and Steven Pekarek are with the Elmore Family School of Electrical and Computer Engineering, Purdue University, West Lafayette, IN 47907 USA (e-mail: hamilt89@purdue.edu;navarr50@purdue.edu;dionysis@purdue.edu;spekarek@purdue.edu).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TIA.2023.3304621.
synchronous generator excitation and turbine systems, including speed governors, automatic voltage regulators (AVR), and power system stabilizers (PSS), react in timescales on the order of a second.These controls, which typically rely on proportionalintegral-derivative (PID) feedback loops, have become standard in the power system industry due to their effectiveness and simplicity.However, tuning such systems to maintain system-wide stability under a wide variety of disturbances and operating conditions is a significant engineering challenge.In microgrids, which is the focus of this work, poor PID tuning has been shown to lead to instability due to lower inertia and higher R/X ratios of conductors [1].
In this article, we propose a linearized-trajectory modelpredictive control (LTMPC) architecture that improves shortterm microgrid stability by dynamically adjusting the reference setpoints of existing primary controls.Furthermore, this reduces the engineering effort needed to tune PID controls, which can be time-consuming and tedious.Model-predictive control (MPC) uses a dynamic model to predict the evolution of system states over a finite horizon and determines an optimal control that minimizes a given objective.In certain microgrid applications (e.g., in industrial settings), details of generation, loads, and the network are relatively well-known, and an adequate communication infrastructure is typically in place (see Fig. 1).Thus, an MPC approach that employs information about the network topology (from relays) and system states (from observers in individual machines) becomes practicable since a reasonably accurate system model can be developed.However, due to the combination of short sampling intervals with the nonlinear power system dynamics, a real-time MPC implementation at such timescales is computationally challenging.The main benefit of the proposed LTMPC approach is its computational efficiency, which opens new possibilities for power system control at fast timescales.
A brief review of recent work on MPC-based power system control is provided next.MPC-based secondary control structures to ensure long-term voltage and frequency stability have been proposed [2], [3], [4], [5], [6], [7], [8].These instability mechanisms develop over timescales of several minutes, and are primarily due to power flow imbalance in a network.Thus, the control inputs determined by these MPC approaches include generation dispatch, load shedding, and switching of shunt capacitor banks.In contrast, we focus on short-term stability, often associated with electric machine rotor dynamics.The proposed LTMPC controls the inputs of generator excitation systems and turbines on timescales of fractions of a second.In addition, distributed or decentralized MPC approaches, which reduce the number of optimization variables and facilitate a nonlinear MPC formulation, have been proposed [9], [10], [11], [12].These methods can be advantageous in bulk power systems with high dimensionality and communication constraints; however, their performance may suffer from lack of a system-wide perspective.In contrast, a nonlinear MPC formulation can be intractable in a centralized approach due to the problem dimensionality.In [11], [13], piecewise-affine approximations to the system dynamics are computed offline and used to predict future system states.Alternatively, a state-space transformation, where nonlinear functions are used to map control signals to their physical counterparts, can be used to obtain a linear MPC [14].To decrease computational effort, reduced-order generator dynamic models are commonly used.
The key idea behind the LTMPC is the separation of nonlinear power system dynamics from the control optimization problem through trajectory linearization [15].Trajectory sensitivities have been used to handle nonlinear dynamics in a number of power systems applications from dynamic security assessment and parameter identification [16], [17], [18], [19] to operational decision making, such as load shedding and generation dispatch [20], [21].Trajectory sensitivities have also been used within MPC formulations to address stability concerns at slower timescales [2], [4], [5].In the proposed LTMPC, the use of trajectory linearization leads to a computationally efficient linearly-constrained quadratic program (QP).It also enables higher-fidelity models of power system components, incorporating nonlinearities and constraints such as those found in realistic excitation and turbine systems, to be used in a real-time LTMPC implementation.
In summary, the main features of the proposed LTMPC are: 1) Speed and voltage setpoints of conventional PID loops in turbine prime movers and field winding excitation systems are dynamically adjusted, which enhances short-term system stability.
2) A centralized control approach utilizes system state and network topology information, while accounting for realistic communication and timing delays.3) Trajectory linearization around the nonlinear power system dynamics reduces computational cost, and enables real-time implementation.4) Detailed models of nonlinear power system components (e.g., control loops including limiters) are incorporated.The theoretical underpinnings of the LTMPC were introduced in [22], where the LTMPC was embedded in a simulation environment, and real-time implementation was not a concern.Here, we extend our previous work by implementing the LTMPC in physical hardware, verifying its real-time operation, and testing its performance in a hardware-in-the-loop (HIL) testbed.The following practical implementation considerations are addressed.First, we implement a state observer based on the Extended Kalman Filter (EKF) to estimate values of dynamic states from measurements available in practice (e.g., voltages, currents).Second, realistic communication and timing delays are accounted for.Finally, we verify the LTMPC performance under a certain degree of uncertainty in knowing the exact parameters of the physical plant.
The LTMPC and EKF are introduced in Sections II and III, respectively.The LTMPC implementation for a notional industrial microgrid is discussed in Section IV.Sections V and VI describe the HIL testbed setup, and provide case study results, respectively.Section VII concludes.

II. GENERAL LTMPC FORMULATION
The power system dynamics can be modeled as a system of nonlinear differential-algebraic equations (DAE) of the form Here, x(t) ∈ R n is a (column) vector of system states, u(t) ∈ R m is a vector of control inputs, and y(t) ∈ R p is a vector of system outputs tracked in the MPC objective function.The vector I qd (t) ∈ R 2r consists of the qd stator currents (in individual rotor reference frames) injected into the network by r rotating machines (generators plus motors).The nonlinear differential equations associated with the system state dynamics are described by the vector-valued function f .The function g represents the algebraic network equations and reference frame transformations.The algebraic system output equations are represented by h, where we assume no direct feedthrough from input to output.The overarching goal is to develop an MPC that accepts measurements or estimates of system variables at regular sampling intervals, utilizes a mathematical model of the nonlinear power system dynamics to predict the development of future states and outputs over a finite horizon, and determines the optimal control actions that minimize an appropriate cost function.Nevertheless, the real-time implementation of a nonlinear MPC is impractical due to several factors.First, solving a nonlinear program in each Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
time step is prohibitively expensive with today's computational technology.Second, a relatively short sampling time is required to accurately capture fast dynamics, such as those of the flux linkages of synchronous generator rotor windings.
The proposed LTMPC overcomes these challenges by separating the nonlinear nature of the power system equations from the optimization as follows.In each sampling step, the LTMPC first collects information (with some communication delay) about system states from local state observers, and determines any changes in network topology from circuit breaker relays.Then, a nonlinear simulation of the nominal system dynamics is performed over the prediction horizon, given the latest pieces of information gathered about the present system state, and based on a reasonable control input trajectory forecast (i.e., the optimization result from the previous time step).Next, the system equations are linearized about this trajectory, so linear relationships between changes in inputs and outputs are obtained.The optimal control input trajectory is then determined by solving a QP.The first of these control actions is taken, and the remainder of the trajectory is saved to serve as the control input forecast for the subsequent sampling interval.

A. LTMPC Timing and Delays
Every T s seconds, the LTMPC samples the system state and initiates its calculations.First, the LTMPC predicts system states and outputs over a prediction horizon of N p steps, as shown in Fig. 2. Based on this prediction, the optimal control input is computed over a control horizon of N c steps.To account for delays associated with the LTMPC computation and communication with devices dispersed around the power system, we introduce a delay interval of N d ≥ 1 steps such that N d T s is longer than the worst-case sum of these delays.Hence, an LTMPC calculation that begins at time t 0 includes decision variables starting at time t 0 + N d T s , and maintains preceding control inputs at their initial value (obtained from the optimization of the previous time step).If the control horizon is selected such that N c < N p − N d , the final control input of the control horizon is repeated until the end of the prediction horizon.
Fig. 2(c) shows the delays associated with LTMPC and state observer computations, and delays due to communication between devices.At time t 0 , the state observers take sensor measurements and estimate the dynamics states of interest.The computation time required is shown as the interval (A).State observers then send the state estimates to the centralized LTMPC, with communication delay (B).The LTMPC uses this information to calculate the optimal control inputs.The time required to perform the LTMPC computation is shown as the interval (C).The optimal control setpoints are then sent from the LTMPC to the local PI controllers, with communication delay (D).Finally, the update to the PI controller setpoints is applied at time t 0 + N d T s , which is accounted for in the LTMPC formulation.We emphasize that the proposed control architecture is robust to failures of the communication network since it does not replace the existing PI-based primary controls.That is, if communication between the LTMPC and PI controls was severed, the PI controls could continue to operate normally using their default reference setpoints.
The LTMPC may also receive occasional updates of the network topology following actions of protection equipment.For example, following the clearing of a fault by opening a line, protective relays will send signals to the LTMPC indicating the opened line.The associated communication delays are accounted for in the simulations of Section VI.

B. Nonlinear Simulation and Trajectory Linearization
The LTMPC computation begins by predicting the trajectories of future system states and outputs via simulation.For a computation initiated at time t 0 , the LTMPC is given a measurement (or estimate) of the current system state x(t 0 ) and a nominal input trajectory over the prediction horizon u * (t), which is a staircase signal of the form (Throughout this article, an asterisk is used to denote variables that are predicted by simulation.)A continuous-time solution to the nonlinear DAE ( 1)-( 3) is obtained via an ODE solver.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Hence, the system state x * k , output y * k , and qd currents I * qd,k are obtained by evaluating the continuous-time solution at time t 0 + kT s for each discrete step k of the prediction horizon.
The system dynamics are linearized at each time step around the nominal input trajectory and the predicted trajectories found via simulation.Hence, Jacobian matrices are evaluated at (x * k , u * k , I * qd,k ).In particular, matrix elements corresponding to nonlinear terms are recalculated for each k.Thus, a linear state-space equation for the perturbed system dynamics for each interval k of the prediction horizon is obtained: where δ denotes a perturbation in the variable, and (7) Next, the continuous-time ( 5)-( 6) are discretized: where By recursively applying ( 8)-( 9), the perturbed system outputs for each step of the prediction horizon are written as a function of the changes in control input.This procedure is fairly standard for developing a discrete MPC [24].Thus, the perturbed control inputs, ) are the decision variables.Performing algebraic manipulations, we obtain a linear map between input and output perturbations, (11) where δY ∈ R p(N p −N d ) , and Z is a block lower-triangular matrix.Fig. 2 illustrates the trajectory linearization process.

C. Control Input Optimization
The LTMPC has two objectives: (i) to stay as close as possible to desired reference outputs, and (ii) to return the system to a steady-state operating point as quickly as possible.Thus, the following objective function is defined: where and Y, Y * , Y ref , U, U * , U min , and U max are defined similarly to δY and δU. 1 For example, Y is defined by replacing every term δy k in ( 11) by y k .The first term in ( 12) plays a key role in stabilizing the system during a transient by penalizing deviations between system outputs and reference values.The remaining two terms play an auxiliary role, as they help reduce oscillations in the decision variables once the system gets close to a steady-state.The second term penalizes changes in control between consecutive time steps, whereas the third term equals In steady-state, the combined effect of the two terms is that u k = u * N d for all k ≥ N d .These final two terms also ensure a unique solution of the optimization problem [22].The weight matrices Q, V 1 , and V 2 are diagonal and positive definite.In general, the elements of Q should be chosen such that the first term dominates at the beginning of a transient.The constraints (15) capture limits on the control input magnitude.
This optimization problem can be cast as a standard QP: where

III. STATE OBSERVER BACKGROUND
In practice, not all dynamic states used in the LTMPC model can be directly measured.Hence, state observers are required.In this section, we provide a general formulation of the EKF [25] for a subsystem/component of the microgrid (e.g., a synchronous machine, turbine, or induction machine).
We model the continuous-time subsystem dynamics by where x γ (t) is a vector of subsystem states, u γ (t) is a vector of actual subsystem inputs, and ũγ (t) is a vector of measured subsystem inputs.The vector w 1,γ (t) represents continuoustime noise in the input measurements.Similarly, the process noise w 2,γ (t) is introduced to capture modeling uncertainties.
The subscript γ denotes a variable or equation corresponding to a specific subsystem/component, rather than the entire system as in (1).Here, x γ (t) ⊂ x(t).Next, we discretize the continuoustime model ( 22)-( 23) using the forward Euler method, leading to a discrete-time model for the subsystem dynamics at each discrete time κ: where T f is the sampling time of the EKF.The discrete-time noise vectors w 1,γ,κ and w 2,γ,κ follow zero-mean Gaussian distributions with covariance matrices Q 1,γ,κ and Q 2,γ,κ , respectively. 2We model the subsystem outputs measured by the observer at each discrete time κ by where ỹγ,κ is a vector of the measured outputs, and v γ,κ is the output measurement noise.We assume v γ,κ follows a zero-mean Gaussian distribution with covariance matrix R γ,κ .

IV. LTMPC IMPLEMENTATION
The notional microgrid shown in Fig. 1 is representative of the power system of an industrial facility.This 13.8-kV system is capable of operating in either grid-connected or islanded mode.The loads on each bus consist of a constant impedance and two dissimilar induction motors, the first being a large induction motor and the second representing an aggregation of smaller motors.Each synchronous generator is equipped with a brushless excitation system, and is controlled by a PI-based AVR and speed governor.The proposed LTMPC will dynamically adjust the voltage and speed setpoints of these PI controllers to improve the system response following a disturbance.Key details of the combined controller and state observer implementation for the notional microgrid are discussed herein.System parameters are provided in [22].

A. Synchronous Machine Model
The LTMPC uses a 1.1 synchronous machine model neglecting stator transients [26].While a reduced-order model is used here, higher fidelity machine models (e.g., a 2.2 model) can be used in the LTMPC.The differential equations governing the 1.1 model are: All variables are in per unit (pu).The rotor angle δ (not to be confused with the perturbation operator) is defined with respect to a common synchronously rotating reference frame with frequency ω s .Fig. 3 shows a voltage-behind-transient-reactance (VBR) circuit which represents the algebraic stator-side equations.The field winding voltage, E fd , and the mechanical torque, T m , are inputs to the machine model.The rotor field winding current, I fd , is given by: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. Brushless Exciter Model
Fig. 4 shows a block diagram of the brushless exciter model [27].The field winding voltage of the exciter alternator, E fe , is controlled by a power amplifier circuit.The power amplifier is represented by a first-order transfer function and the input voltage signal, V a , is determined by the AVR.The voltage of the exciter 3-phase winding, which is supplied to the rotating rectifier, is represented by the state V e .The exciter dynamics depend on the main alternator field winding current, I fd , and are nonlinear due to saturation effects.The rectifier modes of operation are captured by the piecewise function When linearizing the exciter dynamics in the LTMPC, the piecewise nature of (39) is addressed by assuming that the rectifier remains in the same mode of operation as predicted during nonlinear simulation.

C. Turbine Model
The turbine model shown in Fig. 5 represents a simplified version of both the gas turbine model GGOV1 and the steam turbine model IEEEG1 [28].The speed governor provides a reference, fsr, for the position of an actuator valve that controls the fuel flow, W f .The model includes nonlinear slew rate limiters.A first-order transfer function is used to model the conversion of fuel into mechanical power, P m .

D. PI Control Model
The AVR and speed governor of each generator are modeled by PI controllers, as shown in Fig. 6.For the AVR, z PI is the q + V 2 d , and y PI is the exciter input V a .For the speed governor, z PI is the rotor speed ω, and y PI is the turbine input fsr.The input, u PI , is the voltage or speed reference setpoint determined by the LTMPC.
Thus, in the four-generator microgrid under consideration, the LTMPC controls the following eight PI setpoints (voltage and speed for each generator), which represent the inputs of the system model ( 1): (40)

E. Load Model
The induction motor loads are modeled in the LTMPC using a 1.1 model in the rotor reference frame [29] and an algebraic VBR circuit similar to Fig. 3.The mechanical load torque is represented as a quadratic function of rotor speed.All constant impedance loads are added into the network admittance matrix.

F. Network Model
The power system network model in the LTMPC neglects electromagnetic transients.A system bus admittance matrix reduced to "internal nodes," Ỹ, is formed by combining the admittances of transmission lines, circuit breakers, constant impedance loads, and VBR machine models [26].The stator currents for the r machines are given by the nonlinear, complex algebraic equation [30], where I q = [I q,1 . . .I q,r ], I d = [I d,1 . . .I d,r ], S ∈ R r×r is a constant diagonal matrix accounting for rotor saliency, and T(x) = diag{[e jδ 1 . . .e jδ r ]} is the rotor reference frame transformation.The VBR voltage sources are (42)

G. LTMPC Formulation for Microgrid Short-Term Stability
To define the system outputs of interest, we consider three types of microgrid short-term stability: frequency stability, transient (rotor angle) stability, and voltage stability.We define the center-of-inertia (COI) speed as where H g and ω g are the inertia constant and rotor speed of generator g, respectively.Frequency instability occurs when the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
COI frequency diverges from the desired system frequency (e.g., 60 Hz).Transient instability occurs when synchronous machine rotor speeds diverge from the COI speed, and voltage instability occurs when bus voltage magnitudes are outside their limits.We quantify these phenomena in the LTMPC objective function by defining where ωg = ω g − ω COI , and |V i | is the bus-i voltage magnitude.
The reference output vector Y ref in (12) is such that all terms associated with ω are equal to 0 (i.e., over the prediction horizon), terms associated with ω COI are equal to 1 pu, and terms associated with |V | are set to the pre-fault voltage magnitudes, determined here by an optimal power flow (OPF).The relative importance of the three stability objectives can be adjusted by tuning the weights in Q.
The dynamics of the microgrid in Fig. 1 can be expressed in the form ( 1)-( 3), based on the definitions of y, u, and the models described above.Thus, the LTMPC formulation outlined in Section II can be applied.The dynamic states of the microgrid model in the LTMPC are x = x gen,1 . . .x gen,4 x im,1 . . .x im,8 , (45) where x gen are the states of a synchronous generator, x im are the states of an induction motor, and x PI,ω and x PI,V are the states of the PI controls for the speed governor and AVR of each generator, respectively.Analytical expressions for the Jacobian matrices in (7) are straightforward to obtain, however their derivations are omitted due to space limitations.

H. Accounting for System Constraints
Several component models used in the LTMPC contain nonlinear saturation limits.These constraints must be accounted for in the LTMPC formulation.The upper and lower bounds on the PI control input u PI can be incorporated as constraints of the form (15).The impacts of all other saturation limits, such as those on fsr and V a or the turbine slew rate limiters, can be modeled exactly during the nonlinear simulation step of the LTMPC computation.However, these saturation functions are not continuously differentiable, and therefore must be approximated during the trajectory linearization step.We approximate a saturation limit on the signal z using the function, where z is the unsaturated signal value, z max and z min are the upper and lower bounds, Fig. 7 shows a comparison of the actual saturation function and the approximation (48).When evaluating the Jacobian matrices in (7), the derivative of ( 48) is given by Thus, the impact of the nonlinear saturation limits is accounted for during the LTMPC trajectory linearization process.

I. Synchronous Machine State Observer Implementation
In practice, the LTMPC requires estimates of the dynamic states of all system components (synchronous machines, turbines, exciters, induction machines, etc.).For the sake of brevity, we illustrate the implementation of the EKF-based state observer for only the synchronous machines.We assume state observers for other system components can be readily developed using the general EKF formulation in Section III. 3or each synchronous machine, we assume sensor measurements of the mechanical rotor angle, θ rm , three-phase stator currents, I abcs , and terminal voltages, V abcs , are available.As a pre-processing step, the electrical rotor angle, θ r can be calculated by multiplying θ rm by the number of pole pairs of the machine.Additionally, using θ r and Park's transformation, the abc stator currents and voltages are transformed into the qd-axis in the rotor reference frame [26].Finally, we assume that measurements (or estimates from other state observers) of the machine inputs T m and E fd are also available.
The EKF is used to estimate the synchronous machine states , and the output measurements are ỹγ = [V q V d θ r ] .The continuous-time dynamics of the machine are modeled by (34)-( 35), (37), and The algebraic model of the measured outputs, h γ (x γ ), can be found using the VBR circuit diagram of Fig. 3. Given these definitions, steps 1-3 of the discrete-time EKF calculation described in Section III can be used to estimate the dynamic states of the synchronous machine.The state δ used in the LTMPC can be found in post-processing by [26] δ = θ r − ω s t . (52) V. HIL TESTBED SETUP The real-time operation of the proposed LTMPC is validated via an HIL testbed (Fig. 8).The testbed includes a high-fidelity simulation of the microgrid and the local state observers for each synchronous machine, which runs on a real-time OPAL-RT 5600 simulator with a Xeon E5 Intel CPU and RT-LAB v2021.1 software.We emphasize that the models used to simulate the physical microgrid in the OPAL-RT differ from those used in the LTMPC and state observers.For example, the ground-truth model employs a 2.2 synchronous machine model for each generator and includes stator and network transients; meanwhile, the LTMPC and EKF use 1.1 synchronous machine models and neglect the fast dynamics of the network, machine stators, and damper windings.
The real-time microgrid simulation is interfaced via Ethernet communication (TCP/IP protocol) with a hardware implementation of the LTMPC.The rate of data exchange over the Ethernet communication channels is the same as the LTMPC sampling time T s .The LTMPC calculations are compiled as MATLAB executable (MEX) files which are run using MATLAB R2021a on a Dell OptiPlex 5050 computer with a 4-core, 3.4-GHz Intel i7-6700 processor.

VI. NUMERICAL CASE STUDY
Numerical simulations were conducted using the HIL testbed to evaluate the real-time performance of the LTMPC following various disturbances. 4The LTMPC parameters are provided in Table I.A sampling time of T s = 0.1 s is used for the controller.The weights associated with rotor speed deviation, COI speed, and voltage magnitudes in Q are 100, 1000, and 1, respectively.The weights associated with speed and voltage setpoints in V 1 and V 2 are 0.1.
In practice, the parameters of the model used in the LTMPC and EKF are not known exactly.In order to study the impact of model parameter uncertainties on the LTMPC performance, we use slightly different model parameters in the LTMPC and EKF than those used in the ground-truth model.We assume an online parameter identification procedure is used to obtain reasonable estimates of the model parameters [31], [32].Thus, we use model parameters in the LTMPC and EKF that are uniformly randomly selected within ±1% of the values used in the ground-truth model.Additionally, it is assumed that the LTMPC obtains updated network topology information from system relays after a breaker is opened or closed, with a communication delay of 0.1 s.
A sampling time of T f = 0.1 ms is used in each EKF-based state observer.A small value of T f is necessary to ensure the accuracy of the forward-Euler method used in the EKF calculation.Since T f < T s , each observer will send the most recent estimate of the synchronous machine states to the LTMPC every T s seconds.
Zero-mean Gaussian noise is added to the input and output measurements for each state observer.A standard deviation of 0.0167 pu is used for measurements of V abcs , I abcs , T m , and E fd , and a standard deviation of 1.67 deg is used for measurements of θ rm .Note that when the simulated Gaussian noise on the abc voltage and currents is propagated through the nonlinear rotor reference frame transformation, the resulting noise on the qd variables becomes non-Gaussian; however, the EKF implementation discussed above assumes Gaussian measurement noise on the qd variables.Despite this minor inaccuracy in noise modeling, the proposed EKF still performs well in the below case study.The covariance matrices Q 1,γ,κ , Q 2,γ,κ , and R γ,κ are all assigned values of I.For each EKF, the initial state estimate is x+ γ,0 is calculated based on the steady-state operating point and model parameters of the state observer.The initial estimation-error covariance is P + γ,0 = I.In practice, if the sensor measurement noise or initial machine states can be more accurately estimated, then these values can be adjusted accordingly to improve the EKF performance.

A. System Response Following a Three-Phase Fault
In this study, we consider the case when the microgrid is operating islanded from the grid.The initial operating point of the microgrid is determined by an OPF calculation, based on the initial load values and induction machine slips provided in Table II.The microgrid is operating in steady state until t = 5 s, when a three-phase fault occurs on Line 2-3b close to Bus 3. The fault is cleared after 4 cycles by opening the line.Fig. 9 shows the bus voltage magnitude and rotor speed response of each of the four generators.The LTMPC performance is compared to a baseline case, where the reference setpoints of the PI controls are fixed.Specifically, speed governor setpoints are fixed at 1 pu, and the AVR voltage setpoints are set equal to the initial bus voltage magnitudes determined by the OPF.It can be seen that the LTMPC is able to quickly stabilize the system, whereas the PI control with fixed reference setpoints results in prolonged oscillations of both the bus voltages and rotor speeds.The LTMPC performs well during the transient despite the uncertainties in the model parameters.

B. System Response to Line-Starting of a Large Induction Motor
Next, we compare the performance of the LTMPC and PI controls following the line-starting of a large induction motor.We consider the same system as in Section VI-A; however, the large induction motor on Bus 3 is initially disconnected.At t = 0.5 s, the induction motor is line-started.Fig. 10 and 11 show the response of system voltage magnitudes and generator rotor speeds following the disturbance for baseline PI control and using the proposed LTMPC, respectively.
When the induction motor is first connected and is operating at high slip, it draws a large amount of inrush current (Figs.10(c) and 11(c)), which causes system voltage magnitudes to immediately drop (Figs.10(b) and 11(b)).In response, both the conventional PI controls and LTMPC adjust the inputs to the excitation system of each generator (Fig. 12) to increase generator field winding voltages and begin to recover the bus voltages.However, as the motor nears synchronous speed (around t = 4 s in Figs.10(d) and 11(d)), the stator current magnitudes rapidly decrease causing the bus voltages to increase.The conventional PI controls yield significant overshoot in the voltages and frequency.More specifically, bus voltages reach 28% higher than nominal and the COI frequency deviates from 60 Hz by ±2%.These voltage and frequency deviations last for around 1 s or longer, which may cause overvoltage and over/underfrequency protection settings of generator relays to trip and blackout the microgrid (these relay settings are not modeled in the simulations shown) [33].In contrast, since the LTMPC has a model of the induction motor, it can predict these startup transients and begin adjusting the control setpoints of the primary controls sooner (e.g., see Fig. 12).Thus, the LTMPC is able to avoid the large voltage and frequency deviations and improve the short-term stability of the microgrid over the conventional PI controls.

C. Timing Analysis of LTMPC and EKF
The computation time required by each step of the calculation was recorded for each LTMPC sampling interval in the previous study.A breakdown of the computation time statistics is provided in Table III.In this case, the solution of the optimization problem requires the longest time, followed by the trajectory linearization process and nonlinear simulation; however, these results may depend on system size and the LTMPC parameters.The total computation time required by the LTMPC (interval (C) in Fig. 2) is around 40 ms on average and 70 ms in the worst case.The measured one-way Ethernet communication delay (intervals (B) and (D) in Fig. 2) varied between 1-10 ms for this laboratory setup.The EKF computation time required (interval (A) in Fig. 2) was estimated at around 5-20 μs, which is well within the observer sampling time of 100 μs. 6Thus, even the worst-case sum of computation times and communication delays is less than the allotted delay interval of 100 ms (N d = 1 and T s = 0.1 s).Additionally, there were zero overruns detected by the OPAL-RT, which indicates that the simulation was successfully performed in real time.These results clearly show that the proposed trajectory linearization around the nonlinear power system dynamics reduces the computational cost, and enables real-time MPC implementation.

D. Impact of Communication Quality on LTMPC Performance
In the previous section, we showed that the LTMPC is robust to communication delays.However, there may be scenarios where data sent by the state estimators to the LTMPC is lost (i.e., not received) due to communication failures.To study the LTMPC performance under such scenarios, we repeat the case study of Section VI-A, but include random communication failures.More specifically, every time the state observer for each system component (e.g., generator, induction machine, etc.) sends an updated state estimate to the LTMPC, we simulate that the data will be lost with probability p loss .We assume that when a state estimate update is lost, the LTMPC will use the most recently received data for its calculations.
To evaluate the LTMPC performance in terms of short-term stability, we introduce the metric   7We observe that for values of p loss between 5-30%, the LTMPC performs well, with no significant deterioration.For modern communication networks, data loss of 10-30% is severe (see Table IV).Thus, the LTMPC is robust to reasonable communication quality issues.
For higher levels of data loss (e.g., p loss = 90%), the LTMPC performance degrades substantially.However, a key advantage of the proposed LTMPC architecture is that it does not entirely replace the existing PI-based primary controls.Thus, it is envisioned that in situations where communication quality is poor or communication is lost completely for a long time period, the LTMPC would be disabled (using default values for PI control setpoints) until communication is re-established.

VII. CONCLUSION
The proposed LTMPC can be used to enhance the short-term rotor angle, frequency, and voltage stability of synchronous generator-based microgrids.Through the dynamic adjustment of reference setpoints, the LTMPC improves the performance of existing primary control architectures without the need to re-tune control gains.The trajectory linearization approach enables the use of high-fidelity, nonlinear component dynamic models in a real-time implementation.
The LTMPC was implemented in physical hardware, and state observers based on the EKF were introduced to estimate the dynamic states of the microgrid.The performance and real-time operation of the combined controller-observer system was validated using a HIL testbed.The robustness of the LTMPC to realistic communication delays, and uncertainties in the model parameters was also demonstrated.
As with any other MPC approach, the key ingredient for the successful implementation of the proposed LTMPC is a relatively accurate model of the power system at the timescales of interest, which depends on the availability of component parameters (e.g., reactances and time constants of generators and loads, exciter and turbine control gains, network topology and impedances, etc.).Depending on the information available in a given situation, a family of system models can be developed, ranging from low to high fidelity.Therefore, a possible line of future research is the investigation of the tradeoffs between various levels of modeling complexity and control performance in terms of stabilizing the system, as well as tradeoffs with other important implementation aspects such as computing hardware requirements.In other words, it is important to establish a systematic procedure for identifying the lowest level of modeling complexity so that the LTMPC performs reasonably well in a

Manuscript received 1
March 2023; revised 2 June 2023; accepted 10 July 2023.Date of publication 14 August 2023; date of current version 22 November 2023.Paper 2023-PSEC-0154.R1, presented at the 2022 IEEE Power and Energy Conference at Illinois, Champaign, IL, USA, Mar.10-11, and approved for publication in the IEEE TRANSACTIONS ON INDUSTRY APPLICATIONS by the Power Systems Engineering Committee of the IEEE Industry Applications Society [DOI: 10.1109/PECI54197.2022.9744014].This work was supported by Schweitzer Engineering Laboratories.(Corresponding author: Dakota Hamilton.)

Fig. 9 .
Fig. 9. Comparison of system response following a fault on Bus 3 at end of Line 2-3b for baseline PI control (gray) and LTMPC (black).(a) Rotor speeds of all generators.(b) Voltage magnitudes of all buses.

Fig. 10 .
Fig. 10.System response following the line-starting of a large induction machine on Bus 3 under baseline PI control.(a) Rotor speeds of all generators.(b) Voltage magnitudes of all buses.(c) Stator current magnitudes of all large induction machines.(d) Rotor speeds of all large induction machines.

Fig. 11 .
Fig. 11.System response following the line-starting of a large induction machine on Bus 3 under LTMPC.(a) Rotor speeds of all generators.(b) Voltage magnitudes of all buses.(c) Stator current magnitudes of all large induction machines.(d) Rotor speeds of all large induction machines.
t) − y ref ) Q (y(t) − y ref ) dt, (53)where t sim is the length of the simulation (i.e., 20 s), and the values of reference outputs in y ref and weights in Q match

TABLE III LTMPC
COMPUTATION TIME STATISTICS

TABLE IV LTMPC
PERFORMANCE WITH COMMUNICATION FAILURES those in Y ref and Q (see Section IV-G and VI), respectively.