State Estimation-Based Distributed Energy Resource Optimization for Distribution Voltage Regulation in Telemetry-Sparse Environments Using a Real-Time Digital Twin

: Real-time state estimation using a digital twin can overcome the lack of in-ﬁeld measurements inside an electric feeder to optimize grid services provided by distributed energy resources (DERs). Optimal reactive power control of DERs can be used to mitigate distribution system voltage violations caused by increased penetrations of photovoltaic (PV) systems. In this work, a new technology called the Programmable Distribution Resource Open Management Optimization System (ProDROMOS) issued optimized DER reactive power setpoints based-on results from a particle swarm optimization (PSO) algorithm wrapped around OpenDSS time-series feeder simulations. This paper demonstrates the use of the ProDROMOS in a RT simulated environment using a power hardware-in-the-loop PV inverter and in a ﬁeld demonstration, using a 678 kW PV system in Grafton (MA, USA). The primary contribution of the work is demonstrating a RT digital twin effectively provides state estimation pseudo-measurements that can be used to optimize DER operations for distribution voltage regulation.


Introduction
With the rapid price reductions in clean energy generation over the last several years, renewable energy installations continue to increase in number [1,2]. This trend reduces the exclusive reliance on traditional fossil fuel-based generation, helping decrease carbon dioxide (CO 2 ) emissions which are harmful to the environment [3]. One of the disadvantages of integrating these renewable energy sources into the grid is that they rely on varying natural resources (sunlight, wind, waves, etc.). Variable generation can have adverse effects on the power system, causing instabilities or voltage swings, and affecting power system reliability [4][5][6]. With recent advancements in grid-support capabilities and interoperability standardization [7], devices such as photovoltaic (PV) inverters can also be used to provide voltage regulation [8,9]. In fact, worldwide, grid codes and interconnection standards are now requiring PV inverters to provide more control options [10]. Extensive research in PV inverters grid-support functions to provide ancillary services (voltage and/or frequency regulation) utilizing techniques such as fixed power factor (PF) setpoints, volt-var (VV), and frequency-watt (FW) control [11,12]. Utilities are required to maintain customer voltages within the specific tolerances. Distributed energy resource (DER) reactive power can be adjusted to compensate for voltage deviations. A common grid-support approaches is the voltage-reactive power (volt-var, or VV) control which allows a PV inverter to provide reactive power based on the measured grid voltage at the point of common coupling from stored data [52,53]. In the area of healthcare, GE has used digital twins to monitor hospital operations and it is expected that in the next 10 years each patient will have their own digital twin that monitors health and proposes treatments [54]. Siemens employed digital twins in wastewater treatment plants to monitor water pipes with the objective of achieving energy efficiency and predicting potential risks [55]. BP used digital twin models in oil and gas plants located in remote areas to increase reliability, oil exploration and production [56].
A digital twin is a reliable method for mirroring power systems behavior [57,58]. Duplicating the power system operations provides various benefits. For example, digital twins have been used to support fault detection analysis methods [59][60][61]. Some digital twins can be used to help improve control center capabilities, by analyzing system response in RT [62]. Studies have shown that at a power system level, digital twins enable intelligent communication, support power management controls and provide inputs to optimization control algorithms [63][64][65]. With recent advances in RT communication, data handling and optimization, power systems are evolving into smarter power networks [66]. Intelligent power systems include the collection and monitoring of performance data used to understand operating conditions and ultimately prevent system failures in the power grid [67]. Due to their connectivity, data management and accessibility, digital twins have also incorporated the concept of Internet of Things (IoT), enabling more advanced services [68,69]. A drawback associated with incorporating IoT capabilities into power control applications is that these systems become vulnerable to mishandling of data as well as external threats such as cybersecurity attacks [70][71][72]. Studies have shown that communication between system control and digital twins can be achieved through implementing protection architectures such as blockchains. Blockchain provides security against cyber threats when employing IoT capabilities [73,74]. The literature describes blockchain in many energy applications, ranging from electric vehicles, PV systems and energy storage [75,76]. Moreover, recent studies illustrate how the use of blockchain can aid in securing DER data in power systems [77,78]. Some researchers have proposed using blockchain in combination with IoT power systems to help secure the exchange of information, and thus adding a layer of protection against cybersecurity threats [79,80].
Co-optimizing DER operations to meet multiple objectives is novel in distribution system management. In this work, the main contribution is the field deployment of an innovative DER management system (DERMS) as well as the development of a RT digital twin model of a simulated distribution system that can operate simultaneously with a physical distribution system, influencing physical DERs on the field is presented. This paper provides a proof of concept of the digital twin state estimator. The state estimation approach was deployed within a simulation environment and in the field; in each case the PSO voltage regulation optimization used the state estimators' outputs to control the PV inverters reactive power. The remainder of the paper is structured as follows. Section 2 discusses each of the components of the DERMS, called the Programmable Distribution Resource Open Management Optimization System (ProDROMOS) used to execute the PSO. This section summarizes the optimization approach, PV forecasting, state estimation and communication components of the ProDROMOS technology. Section 3 describers the distribution system used to test ProDROMOS. Simulation results are demonstrated in Section 4. Section 5 illustrates the experimental results obtained utilizing power hardware-in-the-loop (PHIL). Section 6 illustrates field demonstration results. Finally, the project conclusion is presented in Section 7.

Communications and Control Architecture
The communications and control architecture known as ProDROMOS is a software tool with the ability to tackle distributions systems with 50% or more of PV integration (available in the Supplementary Material Section). ProDROMOS provides RT insight into the power system dynamics of a distribution system and its interconnected DERs. It has the capability of calculating the optimal active and reactive power set points for the DERs in order to provide adequate voltage regulation for the distribution system. The software system provides RT protection, voltage regulation, and distribution visualization by integrating the following modules: (1) The Georgia Tech Distribution System Distributed Quasi-Dynamic State Estimator (DS-DQSE) takes IEEE C37.118 feeder telemetry from PMUs and generates a power flow estimation and validates the RT model. This information is used to populate the OpenDSS quasi-state time-series (QSTS) simulations within the optimization engine. (2) The forecasting component provides short-term (e.g., 5-min) forecasts of PV power output and load using recent system states and statistical irradiance modeling in conjunction with PV performance models. (3) An optimization engine determines the necessary reactive power settings for the DER to maintain voltage and distribution protection systems for the time horizon (e.g., 15 min). The optimization evaluates circuit performance given the state estimate loads and DER power forecasts to minimize the risk of voltage or protection violations. (4) The communications system monitors and controls multiple DER devices. PF commands were issued to the DERs using SunSpec Modbus, IEEE 1815 (DNP3), and proprietary protocols.
ProDROMOS controlled the DER reactive power output to minimize the risk of exceeding the ANSI C84.1 Range A voltage limits. Since inverter-specified PF commands can reduce active power and PV owner revenue, the control system was designed to minimize lost PV production as well. Co-optimizing DER operations to meet multiple objectives is novel in distribution system management.
The software elements shown in the block diagram of Figure 1 were assembled by Connected Energy (Philadelphia, PA, USA). The Connected Energy software ran on a separate Linux virtual machine (VM), while one or more Windows VMs hosted the GT software and the EPRI PV simulator. These were all in an Amazon Web Service (AWS) cloud environment. The state estimation was performed by Advanced Power Concepts' (Atlanta, GA, USA) Windows Integrated Grounding System (WinIGS) software. The solutions to the distribution state estimation were used to fill in the active and reactive power consumption information in an OpenDSS simulation model. Forecasted DER power values were calculated using a persistence method and used to populate the OpenDSS PV model power output settings. in order to provide adequate voltage regulation for the distribution system. The software system provides RT protection, voltage regulation, and distribution visualization by integrating the following modules: 1) The Georgia Tech Distribution System Distributed Quasi-Dynamic State Estimator (DS-DQSE) takes IEEE C37.118 feeder telemetry from PMUs and generates a power flow estimation and validates the RT model. This information is used to populate the OpenDSS quasi-state time-series (QSTS) simulations within the optimization engine.
2) The forecasting component provides short-term (e.g., 5-min) forecasts of PV power output and load using recent system states and statistical irradiance modeling in conjunction with PV performance models. 3) An optimization engine determines the necessary reactive power settings for the DER to maintain voltage and distribution protection systems for the time horizon (e.g., 15 min). The optimization evaluates circuit performance given the state estimate loads and DER power forecasts to minimize the risk of voltage or protection violations. 4) The communications system monitors and controls multiple DER devices. PF commands were issued to the DERs using SunSpec Modbus, IEEE 1815 (DNP3), and proprietary protocols.
ProDROMOS controlled the DER reactive power output to minimize the risk of exceeding the ANSI C84.1 Range A voltage limits. Since inverter-specified PF commands can reduce active power and PV owner revenue, the control system was designed to minimize lost PV production as well. Co-optimizing DER operations to meet multiple objectives is novel in distribution system management.
The software elements shown in the block diagram of Figure 1 were assembled by Connected Energy (Philadelphia, PA, USA). The Connected Energy software ran on a separate Linux virtual machine (VM), while one or more Windows VMs hosted the GT software and the EPRI PV simulator. These were all in an Amazon Web Service (AWS) cloud environment. The state estimation was performed by Advanced Power Concepts' (Atlanta, GA, USA) Windows Integrated Grounding System (WinIGS) software. The solutions to the distribution state estimation were used to fill in the active and reactive power consumption information in an OpenDSS simulation model. Forecasted DER power values were calculated using a persistence method and used to populate the OpenDSS PV model power output settings.

State Estimation
Distribution system state estimation has recently gained substantial for its ability to enable a range of distribution services [81,82]. The WinIGS software tool was used to complete the state estimation [83]. This software generated the best estimate of the distribution power flows based on a set of field measurements. However, before running the state estimation, the power system topology, locations of DER and other end-devices, and the models of distribution circuits were required. This information was provided by the utility and reconstructed into state and control algebraic quadratic companion form (SCAQCF). Given the measurements and the device SCAQCF models in a feeder section, the domain-specific design space (DS-DSE) created the measurement mathematical model at the device-level. Then, with the help of network formulation techniques, the measurement mathematical models from device-level were converted to network level measurement models. The state estimation algorithm worked directly with the measurement mathematical model at the network level. The results of a DS-DSE are the best estimate of the states and the validated model of that feeder section. Finally, the output of each DS-DSE for each section was sent to the distribution management system where the states and RT model of the distribution feeder. The distribution feeder model was updated as WinIGS models with the DER for the PHIL simulations and the field demonstration. To generate the state estimation, measured values from the PV devices and other field data must be passed to the state estimator via C37.118 streams or COMTRADE files [84,85]. For this project, only virtualized PMU devices where used; for each of the buses of the feeder model and passed voltage and/or current phasor data as C37.118 streams to the WinIGS state estimation tool. The state estimator passed the system state to the optimization tools, which determined the best set points for the DER devices.
The network model was created using the same object-oriented device modeling approach to create the example feeder in [86]. For the NG feeder, a compact component model consisting of set of algebraic and differential nonlinear equations were converted to a second-order quadratized model and then into SCAQCF, represented in Equations (1) and (2): . .
x T F i f eq,xx x . . .
where, I(t) was the through variables of the device model, x was the external and internal state variables of the device model, u were the control variables of the device model, Y eqx was the matrix defining the linear portion for the state variables, Y equ was the matrix defining the linear part for control variables; nonlinear matrices F i eq,xx , F i eq,uu , F i eq,ux defined the quadratic part for the state variables, control variables, and state and control variables for the ith device, respectively; B eq was a history-dependent vector of the device model, and C f eqc was the constraint history dependent vector of the device model. The system was constrained by operating limits h min ≤ h(x, u) ≤ h max and control variable bounds u min ≤ u ≤ u max . Measurements at the device-level were expressed as a vector function in Equation (3): where z were the measurement variables at time t, Y devm,x and Y devm,u were matrices defining the linear part for the state and control variables of the device-level measurement model; F i devm,xx , F i devm,uu , and F i devm,ux were matrices defining the quadratic part for the state variables, control variables, and both for the ith measurement; C devm was the history dependent vector of the device-level measurement model, and η was a vector of measurement noise errors. An Unconstrained Weighted Least Square (UWLS) Dynamic Distribution System State Estimation (DSSE) [87] was used to minimize the sum of weighted squares of a residual vector: where W was the weight matrix defined as the inverse of the squared standard deviations, , and σ i was the standard deviation corresponding to each measurement z i . The unknown state vector x was obtained by the optimal condition dJ/dx = 0, as indicated in [88]. The results of the DSSE were the best estimate of the states and the validated model of that feeder section, that were merged into the full feeder RT model and states [88]. The DSSE implemented in WinIGS performed an observability test to determine that there are enough measurements to observe/compute the state. If PMU measurements were not present, it would generate an error and not solve for system states. If all the PMU streams were present, WinIGS performed the dynamic state estimation and completed the chi-squared (χ 2 ) test to validate consistency between the estimated state and the network model. In the event of low χ 2 results-indicating bad measurements-WinIGS could initiate bad data identification and removed those bad data. For the NG experiments, once the PMU stream from the Opal-RT simulation were correctly aligned with the SCAQCF models, no data were identified as bad or removed. If one of the PMUs failed, the DSSE would not have enough data to make the model observable, but further revisions to WinIGS to persist the measurements or fallback to pseudo-measurements could be used to solve this issue and make the model more resilient to IED failures.

Forecasting
Short-term forecasts of PV output power can be made using several different models and techniques- [89,90]. Short-term forecasts were made using a persistence method that required the DER location, PV system AC and DC capacity, and historical power [91]. The ability to map forecasts to other DER devices was also established so that if power data was not collected by some of the DER equipment, forecasts could still be created by scaling the production forecasts based on the capacities of each system.

Optimization Engine
The optimization engine used the PSO method to determine the optimal PFs of controllable DER systems. This method determined the optimal PV PF setpoints and associated optimal power flow (OPF) by wrapping an OpenDSS time series simulation of the reduced-order feeder model inside a PSO. The active and reactive components of the loads in the OpenDSS model were populated using live state estimation results. PV forecasts for each of the PV systems were populated in the OpenDSS model to optimize operations over a future time horizon. The OpenDSS load data was changed by the WinIGS state estimation solution and the PV production was populated by the PV forecasts. A simplified representation of the PSO approach is shown in Figure 2. A Python interface was created to capture the state estimation IEEE C37.118 data streams from WinIGS. These phasor data for each of the buses and PV systems were used to calculate the active and reactive power levels for the dynamic loads in the OpenDSS model. Then using the communication interface to OpenDSS, the active and reactive multipliers, p_mult and q_mult, were updated in the OpenDSS environment. The optimization was completed every minute over a 15-min horizon using 3 periods with a 5-min step size. The forecast PV production for each of the epochs was calculated using the forecasting code. In cases where there was no PV production data available, scaled surrogate PV system forecasts were used. The p_mult and q_mult values persisted for the entire time-domain simulation. PSO was selected to locate the optimal DER PF settings because the fitness landscape was nonconvex due to the voltage regulators and other binary components. The PSO engine determined DER PF settings by evaluating circuit performance to minimize the risk of voltage or protection violations while also maximizing economic value. The formulation optimization problem was designed to acquire the voltage regulation values of operating DERs off the unity PF value, as shown in Equation (5). where: C(PF) = ∑ 1 − |PF| (8) where, the variable V represents the vector of bus voltages of the distribution system, while the variable V base is a vector the nominal voltages at each one of the buses. The variable PF is a vector representing the PFs of each one of the DERs in the system. For the experiment, the objective function shown in Equation (5) was minimized when the bus voltages of the distribution system were at the desired V base and PF values. The variable V lim represents the maximum and minimum voltage limit. These values were selected to comply with voltage range A of ANSI C84.1, which specifies voltage limits of ±0.05 pu. During the execution of the optimization algorithm, calculated solutions that deviate from the specified V lim would be discouraged. Equation (8)  lated values that deviate from unity PF, due to active power curtailment during high irradiance condition. For the experiments performed in this paper, the weights used for variables 0 , 1 and 2 were 1.0, 2.0 and 0.05, respectively. The PSO algorithm was configured so that it would not execute if all the bus voltages of the distribution feeder were within the selected nominal voltage threshold of 0.20%. Alternatively, bus voltage values outside the ANSI C84.

Communications System
Communications to/from the Connected Energy system used the DNP3 Application Note AN2013-001 information model to change the grid-support functions [92,93]. PF functions were used to change the active and reactive power behaviors of the PV systems. Data Bus (DBus) is a TCP/IP protocol developed by EPRI used to enable communication between the EPRI PV simulator and the RT MATLAB/Simulink power system simulation [94]. DBus is a co-simulation interface that facilitates communication between co-simulation components running in different threads enabling them to exchange data in an asynchronous or synchronous interaction depending on the use case of the simulation [95]. This interface allows the components of the co-simulation to request information and publish data. The Opal-RT system exchanged voltage, and frequency data points from various nodes through DBus interface with the EPRI simulator. PV Simulator determines the output powers of the devices based on the irradiance, voltage, frequency, and DER settings and then exchanges the data with the Opal-RT system every simulation cycle.

Distribution Systems Under Study
A National Grid (NG) feeder with one large utility-owned PV system was selected for the PHIL simulations and field demonstration (more information on the NG feeder model is available on the Supplementary Section). To generate the state estimation, measured field values from the RT feeder simulation on an Opal-RT OP5600 at Sandia National Laboratories' Distributed Energy Technologies Laboratory (DETL) were sent to the state estimator on the AWS Windows machine using a phasor data concentrator (PDC) IEEE C37.118 stream. Virtualized PMUs were placed on buses in the Opal-RT NG distribution feeder model and sent to a physical SEL-3373 PDC which sent the voltage and current phasor data to the state estimator. The distribution state estimator passed the system state to the optimization tools, which determined the best set points for the DER devices. Figure 3 illustrates the one-line diagram and voltage profile generated by the NG distribution feeder OpenDSS model. Laboratories' Distributed Energy Technologies Laboratory (DETL) were sent to the state estimator on the AWS Windows machine using a phasor data concentrator (PDC) IEEE C37.118 stream. Virtualized PMUs were placed on buses in the Opal-RT NG distribution feeder model and sent to a physical SEL-3373 PDC which sent the voltage and current phasor data to the state estimator. The distribution state estimator passed the system state to the optimization tools, which determined the best set points for the DER devices. Figure  3 illustrates the one-line diagram and voltage profile generated by the NG distribution feeder OpenDSS model.  Circuit reduction was performed on the full NG distribution feeder model to reduce the number of buses in model [96] while maintaining an equivalent model at locations of interest on the feeder [97]. This approach was necessary due to bus limitations when implementing the RT simulations. The NG distribution feeder reduced-order OpenDSS models were built as a MATLAB/Simulink/RT-Lab model and simulated in RT with the Opal-RT 5600 [98,99]. This reduced-order model was also converted to WinIGS format to complete the state estimation and incorporated in the PSO optimization. Figure 4 illustrates the PV inverter distribution along the feeder for the reduced-order model-a 15-bus feeder circuit consisting of a 13.8 kV source and a 13.8/13.8 kV transformer rated at 9990 kVA. The baseline reactive power and active power load demand was 5354.31 kVar and 9494.76 kW. The system consisted of thirty single-phase and 1 three-phase PV inverter for a total PV penetration of 5495.36 kVA. PV inverter 1 was a 684 kVA controllable PV inverter located in Old Upton Rd. All other PV inverters in Figure 4 were uncontrolled. There was also a load tapping changer transformer at the secondary of the substation.  Circuit reduction was performed on the full NG distribution feeder model to reduce the number of buses in model [96] while maintaining an equivalent model at locations of interest on the feeder [97]. This approach was necessary due to bus limitations when implementing the RT simulations. The NG distribution feeder reduced-order OpenDSS models were built as a MATLAB/Simulink/RT-Lab model and simulated in RT with the Opal-RT 5600 [98,99]. This reduced-order model was also converted to WinIGS format to complete the state estimation and incorporated in the PSO optimization. Figure 4 illustrates the PV inverter distribution along the feeder for the reduced-order model-a 15bus feeder circuit consisting of a 13.8 kV source and a 13.8/13.8 kV transformer rated at 9990 kVA. The baseline reactive power and active power load demand was 5,354.31 kVar and 9,494.76 kW. The system consisted of thirty single-phase and 1 three-phase PV inverter for a total PV penetration of 5,495.36 kVA. PV inverter 1 was a 684 kVA controllable PV inverter located in Old Upton Rd. All other PV inverters in Figure 4 were uncontrolled. There was also a load tapping changer transformer at the secondary of the substation.

Power Hardware-in-the-Loop Experimental Results
To further validate the development of the digital twin, as well as the operational effectiveness of the PSO voltage regulation methods, a DER PHIL simulation was conducted to provide better fidelity and to demonstrate that the ProDROMOS can operate with real PV communications systems, ramp rates, and other transients not represented in the EPRI PV simulator. PHIL allowed testing the behavior a physical hardware device in a RT simulated environment as well as identifying any communications challenges before field deployment [100]. PHIL in combination with digital twins has been proposed in the past, but not with the intent of co-simulating a RT distribution model to help manage DERs located in a physical live feeder. [101,102]. A physical residential hardware 3-kW PV inverter located in DETL was scaled using a proportional gain in the simulation to represent the per phase output of the controllable 684 kVA inverter Old Upton Rd PV system for the NG distribution feeder experiments. The PSO was implemented on the PHIL with only the Old Upton Rd PV inverter site being controlled, as shown in Figure 5.

Power Hardware-in-the-Loop Experimental Results
To further validate the development of the digital twin, as well as the operational effectiveness of the PSO voltage regulation methods, a DER PHIL simulation was conducted to provide better fidelity and to demonstrate that the ProDROMOS can operate with real PV communications systems, ramp rates, and other transients not represented in the EPRI PV simulator. PHIL allowed testing the behavior a physical hardware device in a RT simulated environment as well as identifying any communications challenges before field deployment [100]. PHIL in combination with digital twins has been proposed in the past, but not with the intent of co-simulating a RT distribution model to help manage DERs located in a physical live feeder. [101,102]. A physical residential hardware 3-kW PV inverter located in DETL was scaled using a proportional gain in the simulation to represent the per phase output of the controllable 684 kVA inverter Old Upton Rd PV system for the NG distribution feeder experiments. The PSO was implemented on the PHIL with only the Old Upton Rd PV inverter site being controlled, as shown in Figure 5.
The PSO OPF optimization technique was deployed using the NG distribution feeder WinIGS state estimation code and the NG distribution feeder OpenDSS time-series simulations as shown in Figure 5. While minor, some differences between these two simulation models included voltage regulation equipment at the substation and differences in noncontrollable DER devices in the PSO OPF code. For the experimental PHIL results, a Baseline case scenario for the DERs operating without control was obtained and compared with the PSO OPF results. The baseline represented a scenario where no control was applied to PV inverter 1. Four-hour irradiance profiles were created for the 31 DER devices in this model (1 three-phase and 30 single phase profiles). Figure 6 illustrates the baseline and PSO active power profile used for PV inverter 1. These experimental results illustrate that there was close agreement between the two results when the DER devices were operating at unity PF.

AC Current Measurement
Ametek RS180 Grid Simulator The PSO OPF optimization technique was deployed using the NG distribution feeder WinIGS state estimation code and the NG distribution feeder OpenDSS time-series simulations as shown in Figure 5. While minor, some differences between these two simulation models included voltage regulation equipment at the substation and differences in noncontrollable DER devices in the PSO OPF code. For the experimental PHIL results, a Baseline case scenario for the DERs operating without control was obtained and compared with the PSO OPF results. The baseline represented a scenario where no control was applied to PV inverter 1. Four-hour irradiance profiles were created for the 31 DER devices in this model (1 three-phase and 30 single phase profiles). Figure 6 illustrates the baseline and PSO active power profile used for PV inverter 1. These experimental results illustrate that there was close agreement between the two results when the DER devices were operating at unity PF.
The PSO OPF algorithm was programmed to execute every 60 s in order to calculate the optimum active and reactive settings of the Old Upton Rd PV site in the PHIL simulation environment. Figure 7 illustrates the reactive power profile used for the baseline and PSO OPF grid-support functions for PV inverter 1. The NG distribution feeder was a highly unbalanced feeder with the average of the phase voltages close to nominal, so there was little that the PSO OPF control algorithm could do to substantially improve the voltage profile. As a result, for ~16 min of the PHIL simulation, the PSO-derived reactive power settings configured the PV inverter to absorb reactive power in order to pull down the local feeder voltage before injecting smaller reactive power levels for the remainder of the simulation.     The PSO OPF algorithm was programmed to execute every 60 s in order to calculate the optimum active and reactive settings of the Old Upton Rd PV site in the PHIL simulation environment. Figure 7 illustrates the reactive power profile used for the baseline and PSO OPF grid-support functions for PV inverter 1. The NG distribution feeder was a highly unbalanced feeder with the average of the phase voltages close to nominal, so there was little that the PSO OPF control algorithm could do to substantially improve the voltage profile. As a result, for~16 min of the PHIL simulation, the PSO-derived reactive power settings configured the PV inverter to absorb reactive power in order to pull down the local feeder voltage before injecting smaller reactive power levels for the remainder of the simulation.   Figure 8 illustrates a comparison of the experimental PHIL results from PV inverter 1 phase voltages measured at the PCC for the baseline and the PSO OPF regulation method. These experimental results illustrate that for the baseline, the bus voltage of phase A was oscillating at its nominal voltage value. For the case of phase B, experimental results illustrate that the bus voltage was above nominal voltage value. Alternatively, the experimental results for phase C, the bus voltage was below nominal voltage value. Due to the unbalance nature of the NG distribution feeder, the PSO OPF was not able to significantly improve the overall bus voltage, because three-phase PV inverters are not currently capable of independently controlling reactive power levels on each phase independently. Regardless, the PSO OPF did have an impact on the NG distribution feeder voltage. For phase B there was slight improvement, as the PSO OPF attempts to decrease the voltage, moving it closer to the desired target voltage value. For the beginning of the simulation, as the voltage value is below 0.96 pu, the PSO OPF attempts to increase the voltage value of phase C, for the remainder of the experiment there is no significant voltage improvement. In addition, phase B also demonstrated voltage improvement for the first half of the PHIL experimental results.
To draw a conclusion on the performance of the PSO OPF approach, the minimum, maximum, and average bus voltages of the distribution feeder were plotted in Figure 9. The PHIL experimental results show there was slight improvement in the system voltage range for the first 6000 s-constricting the maximum and minimum voltage band. For the remainder of the experimental PHIL results, the PSO OPF only managed to slightly reduce the maximum voltage. However, in terms of the average system voltage, after 6000 s, the average voltage is regulated closer to the desired target voltage of 1.00 pu. A detailed comparison using a formulated metric that evaluates the performance of different voltage regulation methods has been performed in previous research [103,104]. It was determined the PSO OPF method was operating effectively, and, more importantly, these results indicated the PSO OPF would provide safe results when controlling the fielded PV system.  These experimental results illustrate that for the baseline, the bus voltage of phase A was oscillating at its nominal voltage value. For the case of phase B, experimental results illustrate that the bus voltage was above nominal voltage value. Alternatively, the experimental results for phase C, the bus voltage was below nominal voltage value. Due to the unbalance nature of the NG distribution feeder, the PSO OPF was not able to significantly improve the overall bus voltage, because three-phase PV inverters are not currently capable of independently controlling reactive power levels on each phase independently. Regardless, the PSO OPF did have an impact on the NG distribution feeder voltage. For phase B there was slight improvement, as the PSO OPF attempts to decrease the voltage, moving it closer to the desired target voltage value. For the beginning of the simulation, as the voltage value is below 0.96 pu, the PSO OPF attempts to increase the voltage value of phase C, for the remainder of the experiment there is no significant voltage improvement. In addition, phase B also demonstrated voltage improvement for the first half of the PHIL experimental results.   To draw a conclusion on the performance of the PSO OPF approach, the minimum, maximum, and average bus voltages of the distribution feeder were plotted in Figure 9. The PHIL experimental results show there was slight improvement in the system voltage range for the first 6000 s-constricting the maximum and minimum voltage band. For the remainder of the experimental PHIL results, the PSO OPF only managed to slightly reduce the maximum voltage. However, in terms of the average system voltage, after 6000 s, the average voltage is regulated closer to the desired target voltage of 1.00 pu. A detailed comparison using a formulated metric that evaluates the performance of different voltage regulation methods has been performed in previous research [103,104]. It was determined Energies 2021, 14, 774 13 of 21 the PSO OPF method was operating effectively, and, more importantly, these results indicated the PSO OPF would provide safe results when controlling the fielded PV system.

Field Demonstration Results
The ProDROMOS optimization was implemented on a live distribution power system with one utility-scale PV site. Inverters at the NG-owned Old Upton Rd site were controlled for the field demonstrations. Trimark Associates provided a DNP3 API to Connected Energy to read average bus voltage, PF, and active and reactive power. The API also included the ability to enable the PF command and set the PF to three decimal places. Writing the DNP3 point to the Trimark system, issued commands to the SMA Cluster Controller which wrote Modbus registers in the 28 24-kW SMA Tripower PV inverters at the site. Figure 10a illustrates an aerial view of the Old Upton Rd PV Site, while Figure  10b illustrates five of the 28 24-kW SMA Tripower PV inverters located at the Old Upton Rd site.

Field Demonstration Results
The ProDROMOS optimization was implemented on a live distribution power system with one utility-scale PV site. Inverters at the NG-owned Old Upton Rd site were controlled for the field demonstrations. Trimark Associates provided a DNP3 API to Connected Energy to read average bus voltage, PF, and active and reactive power. The API also included the ability to enable the PF command and set the PF to three decimal places. Writing the DNP3 point to the Trimark system, issued commands to the SMA Cluster Controller which wrote Modbus registers in the 28 24-kW SMA Tripower PV inverters at the site. Figure 10a illustrates an aerial view of the Old Upton Rd PV Site, while Figure 10b illustrates five of the 28 24-kW SMA Tripower PV inverters located at the Old Upton Rd site.

Field Demonstration Results
The ProDROMOS optimization was implemented on a live distribution power system with one utility-scale PV site. Inverters at the NG-owned Old Upton Rd site were controlled for the field demonstrations. Trimark Associates provided a DNP3 API to Connected Energy to read average bus voltage, PF, and active and reactive power. The API also included the ability to enable the PF command and set the PF to three decimal places. Writing the DNP3 point to the Trimark system, issued commands to the SMA Cluster Controller which wrote Modbus registers in the 28 24-kW SMA Tripower PV inverters at the site. Figure 10a illustrates an aerial view of the Old Upton Rd PV Site, while Figure  10b illustrates five of the 28 24-kW SMA Tripower PV inverters located at the Old Upton Rd site.  In a similar manner to the PHIL simulations, in order to issue optimal PF setpoints to the Old Upton Rd PV system, a PSO approach wrapped multiple OpenDSS time series simulations to find the optimal setpoints over the specified time horizon. Unfortunately, there were not enough field measurements to make the state estimation observable, so a digital twin model was created that represented the NG distribution feeder power system as well as any associated telemetry. The NG distribution feeder digital twin consisted of a RT-Lab simulation model coupled with the EPRI PV simulator. Figure 11 illustrates the implementation of the digital twin concept for the NG distribution feeder used for the PSO OPF field demonstration.
to the Old Upton Rd PV system, a PSO approach wrapped multiple OpenDSS time series simulations to find the optimal setpoints over the specified time horizon. Unfortunately, there were not enough field measurements to make the state estimation observable, so a digital twin model was created that represented the NG distribution feeder power system as well as any associated telemetry. The NG distribution feeder digital twin consisted of a RT-Lab simulation model coupled with the EPRI PV simulator. Figure 11 illustrates the implementation of the digital twin concept for the NG distribution feeder used for the PSO OPF field demonstration.  Figure 11. Digital Twin Concept for the PSO OPF Field Demonstration Experimental Setup Employed on the National Grid Distribution Feeder. Simulated PMU data from RT-Labs was sent to a SEL-3373 PDC and transmitted to the WinIGS state estimation platform. The results of the state estimation, e.g., voltage and current phasor data at each bus were sent as a C37.118 stream to the PSO OPF algorithm. OpenDSS active and reactive power load values were updated based on the phasor data to represent current conditions. The Old Upton Rd PV forecast was used to update the expected power production levels for all 31 PV devices in the OpenDSS time-series simulation to match the local irradiance. The simulated PV systems were configured to provide full power output (the irradiance was set to 1000 W/m 2 ) and the curtailment function of the PV inverter was used to adjust the output of all the devices to match the production at the Old Upton Rd PV site. Spatial variability was not included in the study because there was no information about cloud fields or speed. When the PSO OPF was solved for the optimal PF for the Old Upton Rd PV site, this PF setting was issued to the physical site as well as the simulated digital twin PV system connected to Opal-RT.

Opal-RT Real-Time
An example of the power production forecast and active power profile is shown in Figure 12. Notice from these results that there is a lag in the energy profile from the PV inverter. Regardless, these results demonstrate that the forecasting approach can track the DER power production reasonably well.
full power output (the irradiance was set to 1000 W/m ) and the curtailment function of the PV inverter was used to adjust the output of all the devices to match the production at the Old Upton Rd PV site. Spatial variability was not included in the study because there was no information about cloud fields or speed. When the PSO OPF was solved for the optimal PF for the Old Upton Rd PV site, this PF setting was issued to the physical site as well as the simulated digital twin PV system connected to Opal-RT.
An example of the power production forecast and active power profile is shown in Figure 12. Notice from these results that there is a lag in the energy profile from the PV inverter. Regardless, these results demonstrate that the forecasting approach can track the DER power production reasonably well.  Figure 13 illustrates the calculated PF and reactive power target solution generated by the PSO OPF algorithm for the NG distribution feeder field demonstration. Notice that the optimal calculated PF values were very close to unity, as shown in Figure 13, and only produced significant reactive power (>50 kVar) during the spike around 5500 s. Looking into the internal optimization states during that period, the PSO OPF found a PF solution which improved the objective function but on the next loop that solution was no longer optimal and returned to the lower reactive power level. It is not clear what change in the OpenDSS initial conditions or simulation environment caused this deviation.   Figure 13 illustrates the calculated PF and reactive power target solution generated by the PSO OPF algorithm for the NG distribution feeder field demonstration. Notice that the optimal calculated PF values were very close to unity, as shown in Figure 13, and only produced significant reactive power (>50 kVar) during the spike around 5500 s. Looking into the internal optimization states during that period, the PSO OPF found a PF solution which improved the objective function but on the next loop that solution was no longer optimal and returned to the lower reactive power level. It is not clear what change in the OpenDSS initial conditions or simulation environment caused this deviation.
full power output (the irradiance was set to 1000 W/m 2 ) and the curtailment function of the PV inverter was used to adjust the output of all the devices to match the production at the Old Upton Rd PV site. Spatial variability was not included in the study because there was no information about cloud fields or speed. When the PSO OPF was solved for the optimal PF for the Old Upton Rd PV site, this PF setting was issued to the physical site as well as the simulated digital twin PV system connected to Opal-RT.
An example of the power production forecast and active power profile is shown in Figure 12. Notice from these results that there is a lag in the energy profile from the PV inverter. Regardless, these results demonstrate that the forecasting approach can track the DER power production reasonably well.  Figure 13 illustrates the calculated PF and reactive power target solution generated by the PSO OPF algorithm for the NG distribution feeder field demonstration. Notice that the optimal calculated PF values were very close to unity, as shown in Figure 13, and only produced significant reactive power (>50 kVar) during the spike around 5500 s. Looking into the internal optimization states during that period, the PSO OPF found a PF solution which improved the objective function but on the next loop that solution was no longer optimal and returned to the lower reactive power level. It is not clear what change in the OpenDSS initial conditions or simulation environment caused this deviation.   Figures 14 and 15 illustrates a comparison for the measurements obtained from the digital twin and the field measurements for the active and reactive power at the PCC of the PV inverter, respectively. These results demonstrate that the developed PHIL platform for the NG distribution feeder provides an adequate representation of the distribution system in the field. The PV forecast is seen lagging the irradiance changes but reasonably approximates the energy production. As the PV power changes, the reactive power produced by the site changes significantly. This produced the swings shown in Figure 15 and caused some of the voltage variability at the PCC. Overall, the PSO operated near unity and could do little to help the voltage imbalance of the distribution feeder.
the PV inverter, respectively. These results demonstrate that the developed PHIL platform for the NG distribution feeder provides an adequate representation of the distribution system in the field. The PV forecast is seen lagging the irradiance changes but reasonably approximates the energy production. As the PV power changes, the reactive power produced by the site changes significantly. This produced the swings shown in Figure 15 and caused some of the voltage variability at the PCC. Overall, the PSO operated near unity and could do little to help the voltage imbalance of the distribution feeder.     the PV inverter, respectively. These results demonstrate that the developed PHIL platform for the NG distribution feeder provides an adequate representation of the distribution system in the field. The PV forecast is seen lagging the irradiance changes but reasonably approximates the energy production. As the PV power changes, the reactive power produced by the site changes significantly. This produced the swings shown in Figure 15 and caused some of the voltage variability at the PCC. Overall, the PSO operated near unity and could do little to help the voltage imbalance of the distribution feeder.      the PV inverter, respectively. These results demonstrate that the developed PHIL platform for the NG distribution feeder provides an adequate representation of the distribution system in the field. The PV forecast is seen lagging the irradiance changes but reasonably approximates the energy production. As the PV power changes, the reactive power produced by the site changes significantly. This produced the swings shown in Figure 15 and caused some of the voltage variability at the PCC. Overall, the PSO operated near unity and could do little to help the voltage imbalance of the distribution feeder.

Conclusions
Some distribution systems may not have the necessary measurement infrastructure required to provide centralized DER control. With more measurement devices and DERs being connected to the grid, understanding the power system dynamics could help in calculating the optimum active and reactive settings of these DERs, guarantee proper voltage stability and system reliability. Herein, this paper demonstrates these concepts utilizing PHIL simulations as well as field demonstrations. A dynamic state estimation solution was demonstrated in which state estimation data is compared to live power measurements. Field demonstrations of the PSO OPF using a 684 kVA PV site were conducted on a live NG distribution feeder located in Grafton (MA, USA). The PSO OPF