A Comparison of DER Voltage Regulation Technologies Using Real-Time Simulations

: Grid operators are now considering using distributed energy resources (DERs) to provide distribution voltage regulation rather than installing costly voltage regulation hardware. DER devices include multiple adjustable reactive power control functions, so grid operators have the di ﬃ cult decision of selecting the best operating mode and settings for the DER. In this work, we develop a novel state estimation-based particle swarm optimization (PSO) for distribution voltage regulation using DER-reactive power setpoints and establish a methodology to validate and compare it against alternative DER control technologies (volt–VAR (VV), extremum seeking control (ESC)) in increasingly higher ﬁdelity environments. Distribution system real-time simulations with virtualized and power hardware-in-the-loop (PHIL)-interfaced DER equipment were run to evaluate the implementations and select the best voltage regulation technique. Each method improved the distribution system voltage proﬁle; VV did not reach the global optimum but the PSO and ESC methods optimized the reactive power contributions of multiple DER devices to approach the optimal solution.


Introduction
Installed variable distributed renewable energy is continuing to grow, with much of this growth being related to distribution systems [1]. In the past, utilities could assume unidirectional power flow and regulation of grid voltage within ANSI Standard C84. 1 [2] limits. The penetration of distributed energy resources (DER) presents challenges for distribution circuit voltage regulation [3] because fluctuating DER current injection on the feeder from renewable energy resources leads to voltage perturbations [4]. The voltage swings are currently corrected by load tap changing (LTC) transformers, capacitor banks, and other voltage regulation equipment-often designed or placed prior to the addition of the DER systems. Many distribution system operators (DSOs) are concerned about increasing DER deployments, and often create screening criteria to limit the deployments [5][6][7], because they cannot guarantee power quality in high DER penetration scenarios [8].
New grid code requirements in Hawaii, California, and at the national level (i.e., IEEE 1547-2018 [9]) are requiring newly installed DER grid support functions, including volt-watt, fixed power factor, and volt-VAR. When these functions are configured correctly, distribution hosting capacity can be drastically increased by mitigating thermal and voltage excursions, minimizing losses, and maintaining ANSI limits [10,11].

Voltage Regulation Methods
In a distribution circuit without generation, power flows from a substation to customers and the voltage generally drops from source to load. As previously stated, the utility will regulate the voltage with LTC transformers, capacitor banks, and other voltage regulation equipment. However, a distribution circuit with DER may create bidirectional power flows and increase voltages at the point of common coupling (PCC). The increasing voltage operating envelope and variable power output of the DER create voltage perturbations and may violate ANSI C84.1 Range A limits. The voltage perturbations increase the mechanical operations of the utility-owned voltage regulation equipment [8,35]. This equipment was not designed for such rapid switching, which increase operations and maintenance (O&M) costs and decreases the DSOs ability to provide safe and reliably power [8]. An alternative solution is to use the equipment that is causing the voltage rise to mitigate the problem by using the DER grid support functions with reactive power capabilities. In this paper, we investigate volt-VAR, ESC, and ProDROMOS controls. Volt-var is a common function required by many European and North American grid codes and interconnection standards [37][38][39] included in most DER devices. ESC and ProDROMOS are more sophisticated controllers which calculate the target reactive power of DER and then set a power factor function to produce the reactive power level. It should be noted that in new DER interconnection standards, such as IEEE Std 1547-2018, there are constant reactive power modes-but these functions were not available in the equipment at the time of these experiments.

Autonomous Volt-VAR Control
The volt-VAR (VV) function adjusts the reactive power of the DER based on the PCC voltage measurements according to an adjustable curve of (voltage and reactive power) points. This distributed autonomous control function was the first reactive power factor evaluated. In this project, the VV curve for the RT experiments was defined with non-aggressive voltage points = {95%, 99%, 101%, 105%} and Var points = {25%, 0%, 0%, −25%} because there were oscillations between the simulated DER device and the power simulation.

Extremum Seeking Control
Extremum seeking control (ESC) is a real-time optimization technique for multi-agent, nonlinear, and infinite-dimensional systems [33][34][35][36]. The decentralized, model free control algorithm operates by adjusting the reactive power to optimize measured outputs, in this case, feeder voltage. The system adjusts inputs; (u); via a sinusoidal injection, demodulates system outputs; (J(u)); to extract approximate gradients, and finally performs a gradient descent. A block diagram of this approach is shown in Figure 1. The designer sets the parameters k, l, h, a and ω. Additionally, unique probing frequencies for Energies 2020, 13, 3562 4 of 26 each controllable DER are chosen, but one must ensure that l and h << ω, ensuring proper operation of the high and low pass filters in each ESC loop. The reader is referred to [36] for more details.
Energies 2020, 13, x FOR PEER REVIEW 3 of 4 ensuring proper operation of the high and low pass filters in each ESC loop. The reader is referred to [36] for more details. Refer to [36] for description.
The optimization function for the simulations is given in Equation (1): where Vi, is the voltage measurement at Bus i and Vnon is the nominal voltage. This is a proposed new grid support function that can achieve the global optimum. We assume that each bus in the circuit feeder, can provide measurements in these experiments and more details will be provided in the RT simulation platform section.

Extremum Seeking Control Parameter Selection
These parameters were chosen based on prior experience with ESC simulations and laboratory demonstrations [33][34][35][36], and a systematic debugging process that considered communication latencies, DER power factor (PF) accuracies, and prior lessons learned. The sequence to implement ESC control follows: 1. The ProDROMOS manager constructs the objective function, J. 2. The ProDROMOS manager configures the DER with unique ωs to avoid controller conflict while producing 10 or more data points per cycle. 3. The parameters l and h are set significantly less than ω, such that the perturbation is passed through the washout filter s/(s + h) but is removed by the lowpass filter l/(s + l). 4. Parameter a is selected to produce the smallest reactive power oscillation that is observable in the objective function, J. 5. Parameter k is set based on designer experience, the stability of ESC simulation, and desired time to reach local minimum of J. 6. The ProDROMOS monitors the objective function and makes PF changes to each DER to improve the performance of the system.
The optimization of ESC parameters is a future research interest due to the time requirements necessary to tune the parameters for the application. The use of the ProDROMOS platform and the RT simulator allowed the ESC designer confidence in the parameters before real world deployment.

ProDROMOS
The ProDROMOS system completed the DER-reactive power setpoint optimization in multiple stages. Initially, intelligent electronic devices (IEDs) from the power system would collect data from the RT, PHIL, or live feeder and send those to the WinIGS state estimation system tool. For the RT The optimization function for the simulations is given in Equation (1): where Vi, is the voltage measurement at Bus i and V non is the nominal voltage. This is a proposed new grid support function that can achieve the global optimum. We assume that each bus in the circuit feeder, can provide measurements in these experiments and more details will be provided in the RT simulation platform section.

Extremum Seeking Control Parameter Selection
These parameters were chosen based on prior experience with ESC simulations and laboratory demonstrations [33][34][35][36], and a systematic debugging process that considered communication latencies, DER power factor (PF) accuracies, and prior lessons learned. The sequence to implement ESC control follows:

1.
The ProDROMOS manager constructs the objective function, J.

2.
The ProDROMOS manager configures the DER with unique ωs to avoid controller conflict while producing 10 or more data points per cycle.

3.
The parameters l and h are set significantly less than ω, such that the perturbation is passed through the washout filter s/(s + h) but is removed by the lowpass filter l/(s + l).

4.
Parameter a is selected to produce the smallest reactive power oscillation that is observable in the objective function, J.

5.
Parameter k is set based on designer experience, the stability of ESC simulation, and desired time to reach local minimum of J.

6.
The ProDROMOS monitors the objective function and makes PF changes to each DER to improve the performance of the system.
The optimization of ESC parameters is a future research interest due to the time requirements necessary to tune the parameters for the application. The use of the ProDROMOS platform and the RT simulator allowed the ESC designer confidence in the parameters before real world deployment.

ProDROMOS
The ProDROMOS system completed the DER-reactive power setpoint optimization in multiple stages. Initially, intelligent electronic devices (IEDs) from the power system would collect data from the RT, PHIL, or live feeder and send those to the WinIGS state estimation system tool. For the RT and PHIL experiments, these data were collected from virtualized micro-phasor measurement units (µPMUs) and formatted in IEEE C37.118 format using a physical SEL-3373 Phasor Data Concentrator (PDC) in the Distributed Energy Technologies Laboratory (DETL) at Sandia National Laboratories. These data were then parsed into different sections of the WinIGS feeder. By segmenting the state estimation problem into smaller feeder pieces, WinIGS was able to solve the state estimation at 30 Hz. Details of the construction of the WinIGS models for the Public Service Company of New Mexico (PNM) and National Grid (NG) feeders in NM and MA are presented in [26]. The solution to the state estimation included active and reactive power values for different loads in the model. These values were passed to the optimization routine as another IEEE C37.118 data stream.
The Particle Swarm Optimization (PSO) was used to determine the Optimal Power Factor (OPF) settings for the DER devices because of the nonconvex fitness landscape [26]. PSO OPF uses a time-series OpenDSS simulation wrapped in a PSO loop to calculate the OPF values. The load data for OpenDSS were populated with state estimation solutions before each run of the optimization routine, which occurred every 1 or 2 min. A solar production persistence forecast was created to estimate the solar energy over the next time horizon (5 min). This forecast was used as the photovoltaic (PV) power profile in the OpenDSS simulation. The WinIGS software and the python-based PSO code were implemented in a portable Connected Energy cloud-based Windows virtual machine (VM) that allowed geographically distributed developers to access the environment. Figure 2 provides a simplified approach of the PSO OPF method.
Energies 2020, 13, x FOR PEER REVIEW  3 of 4 and PHIL experiments, these data were collected from virtualized micro-phasor measurement units (μPMUs) and formatted in IEEE C37.118 format using a physical SEL-3373 Phasor Data Concentrator (PDC) in the Distributed Energy Technologies Laboratory (DETL) at Sandia National Laboratories. These data were then parsed into different sections of the WinIGS feeder. By segmenting the state estimation problem into smaller feeder pieces, WinIGS was able to solve the state estimation at 30 Hz. Details of the construction of the WinIGS models for the Public Service Company of New Mexico (PNM) and National Grid (NG) feeders in NM and MA are presented in [26]. The solution to the state estimation included active and reactive power values for different loads in the model. These values were passed to the optimization routine as another IEEE C37.118 data stream. The Particle Swarm Optimization (PSO) was used to determine the Optimal Power Factor (OPF) settings for the DER devices because of the nonconvex fitness landscape [26]. PSO OPF uses a timeseries OpenDSS simulation wrapped in a PSO loop to calculate the OPF values. The load data for OpenDSS were populated with state estimation solutions before each run of the optimization routine, which occurred every 1 or 2 minutes. A solar production persistence forecast was created to estimate the solar energy over the next time horizon (5 minutes). This forecast was used as the photovoltaic (PV) power profile in the OpenDSS simulation. The WinIGS software and the python-based PSO code were implemented in a portable Connected Energy cloud-based Windows virtual machine (VM) that allowed geographically distributed developers to access the environment. Figure 2 provides a simplified approach of the PSO OPF method. Once populated with the anticipated PV production and the current-and assumed staticactive and reactive power levels for the loads in the OpenDSS simulations, the PSO ran the OpenDSS model with different DER PF setpoints to calculate the optimal setpoints for each of the DER devices. The optimization formulation was designed to capture the voltage regulation and economic considerations of operating PV systems with a non-unity power factor (PF): where if any |V| > Vlim (4) Once populated with the anticipated PV production and the current-and assumed static-active and reactive power levels for the loads in the OpenDSS simulations, the PSO ran the OpenDSS model with different DER PF setpoints to calculate the optimal setpoints for each of the DER devices. The optimization formulation was designed to capture the voltage regulation and economic considerations of operating PV systems with a non-unity power factor (PF): Energies 2020, 13, 3562 is the standard deviation of V − V base V is a vector of bus voltages, V base is a vector of the nominal voltages for each bus, and PF is a vector of the DER PFs. The objective function is minimized when the bus voltages are at V base and PF = 1. V lim was selected to be the ANSI C84.1 Range A limits of ±0.05 per unit (pu), so any solutions outside the limits would be highly penalized. The third term was a simplified method to discourage solutions that moved away from unity power factor, because these solutions would curtail active power (and expense the PV owner through net metering, power purchase agreements, etc.) at high irradiance times. More sophisticated methods for determining the curtailment magnitude were considered, but the simple approach shown here was implemented. For the experiments conducted in this project, w 0 = 1.0, w 1 = 2.0, and w 2 = 0.05. The optimization was configured so that if all the bus voltages were within an acceptance threshold (set to 0.2% of nominal voltage) the PSO would not run. If any of the voltages were outside ANSI Range A the PSO would run. Furthermore, if the new PF values did not change the objective function by an objective threshold (set to 1 × 10 −7 ), the new PFs would not be sent to the DER devices to minimize communications and DER memory writes.

Distribution Systems of Study
Several distribution systems models were acquired from PNM and NG for evaluation. The PNM models were in the SynerGEE software format and NG models were in the CYME software format. To run the power hardware-in-the-loop (PHIL) experiments, a series of conversions were needed to convert the models to OPAL RT-Lab-compliant formats (MATLAB/Simulink). The CYME and SynerGEE models were first converted to OpenDSS. The objective of the project was to show voltage regulation with 50% PV penetration level (defined by feeder PV nameplate/peak feeder load). For those feeders that did not meet this threshold, additional residential PV systems were randomly added to the buses in the model to meet this goal. Time-domain RT-Lab simulations are not possible with feeder models with thousands of buses, so the full OpenDSS models were reduced in complexity to be incorporated in the RT-Lab and WinIGS platform. The circuit reduction algorithms reduced the number of buses in a model while maintaining an equivalent voltage profile throughout the model at feeder locations of interest [40,41]. Ultimately, a commercial 12.47-kV PNM feeder and a suburban 13.80 kV NG feeder were used for the evaluation of the different voltage regulation algorithms. The PNM reduced order feeder that was used in these experiments is shown in Figure 3 and the NG feeder is shown in Figure 4.

Real-Time Evaluation Platform
The evaluation platform could be run in a software-in-the-loop (SIL), a power-hardware-in-theloop (PHIL), or a combination of SIL and PHIL setup. The SIL and PHIL simulations were conducted with an Opal-RT OP5600 real-time digital simulator running distribution feeder models in RT-LAB with a 100 us timestep. To complete the PHIL experiments, a physical 3.0-kW inverter was interfaced to the power simulation via a 180 kVA Ametek RS180 grid simulator. The current waveform was measured by a Person CT110 and returned to the power simulation via an analog input channel. The DC side of the inverter was connected to an Ametek TerraSAS PV simulator configured to represent a 3.0-kWp silicon PV system.
In the RT simulations, the EPRI-developed DER Simulator [42] software, emulated smart solar inverters with several grid support functionalities and communication interfaces. Any number of DER devices, with different nameplate capacities and phasing configurations, were instantiated with independent irradiance profiles. The EPRI-developed TCP/IP protocol called Data Bus (DBus) was used as the co-simulation interface to exchange data between the OPAL-RT simulator and the DER simulator. The DER simulator determined the active and reactive power of the DER devices based on irradiance profiles, voltage, frequency, and DER settings. The active and reactive power outputs are exchanged via DBus with the RT simulation. The RT simulation passed the voltage and frequency of the points of common coupling (PCCs) where the DER were interconnected to the feeder.

Real-Time Evaluation Platform
The evaluation platform could be run in a software-in-the-loop (SIL), a power-hardware-in-theloop (PHIL), or a combination of SIL and PHIL setup. The SIL and PHIL simulations were conducted with an Opal-RT OP5600 real-time digital simulator running distribution feeder models in RT-LAB with a 100 us timestep. To complete the PHIL experiments, a physical 3.0-kW inverter was interfaced to the power simulation via a 180 kVA Ametek RS180 grid simulator. The current waveform was measured by a Person CT110 and returned to the power simulation via an analog input channel. The DC side of the inverter was connected to an Ametek TerraSAS PV simulator configured to represent a 3.0-kWp silicon PV system.
In the RT simulations, the EPRI-developed DER Simulator [42] software, emulated smart solar inverters with several grid support functionalities and communication interfaces. Any number of DER devices, with different nameplate capacities and phasing configurations, were instantiated with independent irradiance profiles. The EPRI-developed TCP/IP protocol called Data Bus (DBus) was used as the co-simulation interface to exchange data between the OPAL-RT simulator and the DER simulator. The DER simulator determined the active and reactive power of the DER devices based on irradiance profiles, voltage, frequency, and DER settings. The active and reactive power outputs are exchanged via DBus with the RT simulation. The RT simulation passed the voltage and frequency of the points of common coupling (PCCs) where the DER were interconnected to the feeder.

Real-Time Evaluation Platform
The evaluation platform could be run in a software-in-the-loop (SIL), a power-hardware-in-the-loop (PHIL), or a combination of SIL and PHIL setup. The SIL and PHIL simulations were conducted with an Opal-RT OP5600 real-time digital simulator running distribution feeder models in RT-LAB with a 100 us timestep. To complete the PHIL experiments, a physical 3.0-kW inverter was interfaced to the power simulation via a 180 kVA Ametek RS180 grid simulator. The current waveform was measured by a Person CT110 and returned to the power simulation via an analog input channel. The DC side of the inverter was connected to an Ametek TerraSAS PV simulator configured to represent a 3.0-kWp silicon PV system.
In the RT simulations, the EPRI-developed DER Simulator [42] software, emulated smart solar inverters with several grid support functionalities and communication interfaces. Any number of DER devices, with different nameplate capacities and phasing configurations, were instantiated with independent irradiance profiles. The EPRI-developed TCP/IP protocol called Data Bus (DBus) was used as the co-simulation interface to exchange data between the OPAL-RT simulator and the DER simulator. The DER simulator determined the active and reactive power of the DER devices based on irradiance profiles, voltage, frequency, and DER settings. The active and reactive power outputs are exchanged via DBus with the RT simulation. The RT simulation passed the voltage and frequency of the points of common coupling (PCCs) where the DER were interconnected to the feeder. The PHIL simulations used the RT simulation environment but replaced one of the EPRI PV systems with a 3.0-kW residential inverter that was then scaled to a device in the PNM model (inverter 3, 258 kW) and NG feeder (inverter 1, 658 kW). The 258-kW and 684-kW inverters were three-phase balanced output inverters where each phase produced the same current magnitude and phase angle. The 3-kW inverter was scaled to match its per phase output to the 684-kW inverter for the NG experiments and the 258-kW inverter for the PNM experiments. The advantage of incorporating a physical inverter into the simulations was twofold: the control algorithms had to be built with communication interfaces to interact with real equipment over the public internet networks and the abnormal operations of physical equipment (e.g., startup times, ramp rates, non-optimal maximum power point tracking, converter losses, etc.) were represented in the power simulation and needed to be managed by the control algorithms. Scaling a residential PV system to represent a utility-owned site likely produced higher reactive power errors because of the lower-cost control and power electronics components in the 3-kW inverter, but this helped approximate inverter interoperability and control functionality, with the expectation that the fielded system would perform better.

Simulations with the PNM Model
The PNM model included two scenarios. In the first scenario, there were three simulated EPRI PV inverters, and in the second scenario, the 258-kW inverter at bus 15 was replaced with a physical single phase 3.0 kW that was then scaled to the 258-kW rating in the simulation. The EPRI simulator provided three correlated, highly variable irradiance profiles from a PV site on the US east coast to feed the simulated inverters. These profiles were generated from a single Global Horizontal Irradiance (GHI) measurement and then shifted and scaled to account for geographical separation of the PV sites. The PV irradiance profiles were then smoothed based on the plant sizes using the Wavelet Variability Model [43][44][45][46] to account for spatial averaging with large PV systems. To understand the impact of spatial averaging on the irradiance profiles the reader is invited to review [47,48]. The final PV profiles are shown in Figure 5. The smoothed irradiance profiles are used because clouds do not instantaneously shade larger PV systems and it takes time for them to pass over the entire array. Additionally, variable load profiles were used on loads connected to bus 6 and bus 13 to create transient voltage variability.
Energies 2020, 13, x FOR PEER REVIEW 3 of 4 The PHIL simulations used the RT simulation environment but replaced one of the EPRI PV systems with a 3.0-kW residential inverter that was then scaled to a device in the PNM model (inverter 3, 258 kW) and NG feeder (inverter 1, 658 kW). The 258-kW and 684-kW inverters were three-phase balanced output inverters where each phase produced the same current magnitude and phase angle. The 3-kW inverter was scaled to match its per phase output to the 684-kW inverter for the NG experiments and the 258-kW inverter for the PNM experiments. The advantage of incorporating a physical inverter into the simulations was twofold: the control algorithms had to be built with communication interfaces to interact with real equipment over the public internet networks and the abnormal operations of physical equipment (e.g., startup times, ramp rates, non-optimal maximum power point tracking, converter losses, etc.) were represented in the power simulation and needed to be managed by the control algorithms. Scaling a residential PV system to represent a utility-owned site likely produced higher reactive power errors because of the lower-cost control and power electronics components in the 3-kW inverter, but this helped approximate inverter interoperability and control functionality, with the expectation that the fielded system would perform better.

Simulations with the PNM Model
The PNM model included two scenarios. In the first scenario, there were three simulated EPRI PV inverters, and in the second scenario, the 258-kW inverter at bus 15 was replaced with a physical single phase 3.0 kW that was then scaled to the 258-kW rating in the simulation. The EPRI simulator provided three correlated, highly variable irradiance profiles from a PV site on the US east coast to feed the simulated inverters. These profiles were generated from a single Global Horizontal Irradiance (GHI) measurement and then shifted and scaled to account for geographical separation of the PV sites. The PV irradiance profiles were then smoothed based on the plant sizes using the Wavelet Variability Model [43][44][45][46] to account for spatial averaging with large PV systems. To understand the impact of spatial averaging on the irradiance profiles the reader is invited to review [47,48]. The final PV profiles are shown in Figure 5. The smoothed irradiance profiles are used because clouds do not instantaneously shade larger PV systems and it takes time for them to pass over the entire array. Additionally, variable load profiles were used on loads connected to bus 6 and bus 13 to create transient voltage variability.

Baseline Simulation
The PNM model was simulated using P/Q data passed from the EPRI DER simulator to the RT power simulation using the Data Bus (DBus) exchange. The hardware PV inverter was interfaced using the Ideal Transformer Method (ITM) [48]. The simulated PV inverters connected to bus 12 and bus 14 were simulated using a P/Q controlled PV inverter model [49]. In this feeder model, controllable capacitor banks were removed. At each PV site, irradiance is converted to AC power using pvlib-python [50]; an example of the power simulation is shown in Figure 6. The voltage profile for each of the 15 buses is shown in Figure 7. As can be seen from these results, there is relatively little voltage rise when the PV systems inject power at their points of common coupling (PCCs). Overall, the feeder maximum voltage is driven by the 10-MW PV system on this feeder.
Energies 2020, 13, x FOR PEER REVIEW 3 of 4 The PNM model was simulated using P/Q data passed from the EPRI DER simulator to the RT power simulation using the Data Bus (DBus) exchange. The hardware PV inverter was interfaced using the Ideal Transformer Method (ITM) [48]. The simulated PV inverters connected to bus 12 and bus 14 were simulated using a P/Q controlled PV inverter model [49]. In this feeder model, controllable capacitor banks were removed. At each PV site, irradiance is converted to AC power using pvlib-python [50]; an example of the power simulation is shown in Figure 6. The voltage profile for each of the 15 buses is shown in Figure 7. As can be seen from these results, there is relatively little voltage rise when the PV systems inject power at their points of common coupling (PCCs). Overall, the feeder maximum voltage is driven by the 10-MW PV system on this feeder.   Energies 2020, 13, x FOR PEER REVIEW 3 of 4 The PNM model was simulated using P/Q data passed from the EPRI DER simulator to the RT power simulation using the Data Bus (DBus) exchange. The hardware PV inverter was interfaced using the Ideal Transformer Method (ITM) [48]. The simulated PV inverters connected to bus 12 and bus 14 were simulated using a P/Q controlled PV inverter model [49]. In this feeder model, controllable capacitor banks were removed. At each PV site, irradiance is converted to AC power using pvlib-python [50]; an example of the power simulation is shown in Figure 6. The voltage profile for each of the 15 buses is shown in Figure 7. As can be seen from these results, there is relatively little voltage rise when the PV systems inject power at their points of common coupling (PCCs). Overall, the feeder maximum voltage is driven by the 10-MW PV system on this feeder.

Volt-VAR Simulation
The volt-VAR function was implemented by programming the VV function in the EPRI DER devices using a DNP3 command. The deadband was set to ±0.01 pu so the DER would support the feeder voltage. Aggressive VV parameters produced reactive power oscillations because the EPRI simulator-to-Opal-RT communication rates were limited to once/second-e.g., DBus read the DER bus voltage once a second and then DBus updated the Opal-RT DER-reactive power set point on the next DBus write. As a result, a stable VV curve defined as V = {92, 99, 101, 108} pu, Q = {25, 0, 0, −25}% nameplate VA capacity were used for the simulations. The reactive power contribution from the 10-MW plant was modest, but it reduced the maximum and average feeder voltages, as shown in Figure 8-where the minimum, maximum, and average voltages were plotted using for all buses and all phases. The line represents the average voltage for all the buses and the colored patch represents the range of voltages over time for the feeder. Notice from the simulation results that the VV function can bring the average bus voltage closer to the desired nominal value. These effects can also be observed from the maximum voltage, reducing the voltage band created from the maximum and minimum values.

Volt-VAR Simulation
The volt-VAR function was implemented by programming the VV function in the EPRI DER devices using a DNP3 command. The deadband was set to ±0.01 pu so the DER would support the feeder voltage. Aggressive VV parameters produced reactive power oscillations because the EPRI simulator-to-Opal-RT communication rates were limited to once/second-e.g., DBus read the DER bus voltage once a second and then DBus updated the Opal-RT DER-reactive power set point on the next DBus write. As a result, a stable VV curve defined as V = {92, 99, 101, 108} pu, Q = {25, 0, 0, −25}% nameplate VA capacity were used for the simulations. The reactive power contribution from the 10-MW plant was modest, but it reduced the maximum and average feeder voltages, as shown in Figure  8-where the minimum, maximum, and average voltages were plotted using for all buses and all phases. The line represents the average voltage for all the buses and the colored patch represents the range of voltages over time for the feeder. Notice from the simulation results that the VV function can bring the average bus voltage closer to the desired nominal value. These effects can also be observed from the maximum voltage, reducing the voltage band created from the maximum and minimum values.

Extremum Seeking Control Simulation
ESC was implemented with a RT model of the PNM feeder using multiple parameter sets. When selecting ESC parameters, there are multiple tradeoffs that make it challenging. For instance, large probing magnitudes help the objective function signal stand out from feeder noise, but it causes larger voltage deviations once a solution is reached. The normalized voltage deviation objective function was used with probing signal frequencies that could be independently demodulated. Ultimately, the parameters shown in Table 1 were used. As shown in Figure 9, the reactive power from the DER devices to move the average bus voltage toward the nominal set point, i.e., 1 pu as compared to the baseline simulation. The maximum and minimum bus voltage values were adjusted closer to nominal as well. The reactive power probing signal's impact on the feeder voltage is clearly seen in the results.

Extremum Seeking Control Simulation
ESC was implemented with a RT model of the PNM feeder using multiple parameter sets. When selecting ESC parameters, there are multiple tradeoffs that make it challenging. For instance, large probing magnitudes help the objective function signal stand out from feeder noise, but it causes larger voltage deviations once a solution is reached. The normalized voltage deviation objective function was used with probing signal frequencies that could be independently demodulated. Ultimately, the parameters shown in Table 1 were used. As shown in Figure 9, the reactive power from the DER devices to move the average bus voltage toward the nominal set point, i.e., 1 pu as compared to the baseline simulation. The maximum and minimum bus voltage values were adjusted closer to nominal as well. The reactive power probing signal's impact on the feeder voltage is clearly seen in the results. Energies 2020, 13, x FOR PEER REVIEW 3 of 4

State-Estimation-Based Particle Swarm Optimization
The setup for the PSO simulations is shown in Figure 2. In this case, the state estimation was used to update the load data in the OpenDSS time series simulation and forecasts of three-minute average DER power were used to update the OpenDSS power levels. Figure 10 illustrates the average power forecasts. The forecasts track the average irradiance reasonably well with some lag as is common for persistence forecasts. Note that it takes 1 hour for the forecast to begin because it uses prior production data and the clear sky index to generate the power prediction.
The OpenDSS time series simulation was run multiple times using PSO to determine the optimal PF set points for each of the DER devices. It was found that absolute care must be taken to use the same component and settings for the Simulink/RT-Lab and OpenDSS models; otherwise erroneous solutions were found. Binary variables (e.g., capacitor banks) and discrete variables (e.g., tap changes) must also be re-initialized (i.e. reset) with each new run of the OpenDSS time series simulation to prevent initial conditions from being set by the prior simulation-an unfortunate byproduct of using the OpenDSS communication interface implementation. As shown in Figure 11, the reactive power from the DER devices to move the average bus voltage toward the nominal set point, i.e., 1 pu as compared to the baseline simulation. Interestingly, with the constraints placed on updating a new PF level (e.g., the solution must be a certain amount better than the previous one), there are only 3 times the DER devices had their PF setting updated as shown in Figure 12. This threshold can be tuned, but it is desirable because it reduces the number of PF write commands and shows the solution is robust to changing PV irradiance. As seen in Figure 12, the initial PF setpoints for DER 1 and 2 absorbed reactive power while DER 3 injected reactive power. For this feeder topology and DER locations, the global optimum required reactive power absorption or injection to minimize the voltage deviation from nominal at all buses in the feeder.

State-Estimation-Based Particle Swarm Optimization
The setup for the PSO simulations is shown in Figure 2. In this case, the state estimation was used to update the load data in the OpenDSS time series simulation and forecasts of three-minute average DER power were used to update the OpenDSS power levels. Figure 10 illustrates the average power forecasts. The forecasts track the average irradiance reasonably well with some lag as is common for persistence forecasts. Note that it takes 1 h for the forecast to begin because it uses prior production data and the clear sky index to generate the power prediction.
The OpenDSS time series simulation was run multiple times using PSO to determine the optimal PF set points for each of the DER devices. It was found that absolute care must be taken to use the same component and settings for the Simulink/RT-Lab and OpenDSS models; otherwise erroneous solutions were found. Binary variables (e.g., capacitor banks) and discrete variables (e.g., tap changes) must also be re-initialized (i.e., reset) with each new run of the OpenDSS time series simulation to prevent initial conditions from being set by the prior simulation-an unfortunate byproduct of using the OpenDSS communication interface implementation. As shown in Figure 11, the reactive power from the DER devices to move the average bus voltage toward the nominal set point, i.e., 1 pu as compared to the baseline simulation. Interestingly, with the constraints placed on updating a new PF level (e.g., the solution must be a certain amount better than the previous one), there are only 3 times the DER devices had their PF setting updated as shown in Figure 12. This threshold can be tuned, but it is desirable because it reduces the number of PF write commands and shows the solution is robust to changing PV irradiance. As seen in Figure 12, the initial PF setpoints for DER 1 and 2 absorbed reactive power while DER 3 injected reactive power. For this feeder topology and DER locations, the global optimum required reactive power absorption or injection to minimize the voltage deviation from nominal at all buses in the feeder.

Simulations with the National Grid Model
Like the PNM model, VV, ESC, and PSO voltage regulation approaches were compared on the NG feeder model. In this model, there is a 684-kW utility-owned three-phase PV site and 30 single phase PV systems. The NG system was unbalanced, with phase B above nominal voltage, phase C below nominal voltage, and phase A operating approximately at nominal voltage. This prevented many of the voltage regulation methods from making significant improvements to the feeder voltage profile when only the three-phase utility-owned inverter was controlled. Therefore, we proposed controlling all 31 DER on the NG model for further testing.

Baseline Simulation with National Grid Model
Four-hour irradiance profiles, at one second resolution were created for the 31 DER devices in this model. Figure 13 shows the voltage profiles for each of the phases.

Simulations with the National Grid Model
Like the PNM model, VV, ESC, and PSO voltage regulation approaches were compared on the NG feeder model. In this model, there is a 684-kW utility-owned three-phase PV site and 30 single phase PV systems. The NG system was unbalanced, with phase B above nominal voltage, phase C below nominal voltage, and phase A operating approximately at nominal voltage. This prevented many of the voltage regulation methods from making significant improvements to the feeder voltage profile when only the three-phase utility-owned inverter was controlled. Therefore, we proposed controlling all 31 DER on the NG model for further testing.

Baseline Simulation with National Grid Model
Four-hour irradiance profiles, at one second resolution were created for the 31 DER devices in this model. Figure 13 shows the voltage profiles for each of the phases.

Simulations with the National Grid Model
Like the PNM model, VV, ESC, and PSO voltage regulation approaches were compared on the NG feeder model. In this model, there is a 684-kW utility-owned three-phase PV site and 30 single phase PV systems. The NG system was unbalanced, with phase B above nominal voltage, phase C below nominal voltage, and phase A operating approximately at nominal voltage. This prevented many of the voltage regulation methods from making significant improvements to the feeder voltage profile when only the three-phase utility-owned inverter was controlled. Therefore, we proposed controlling all 31 DER on the NG model for further testing.

Baseline Simulation with National Grid Model
Four-hour irradiance profiles, at one second resolution were created for the 31 DER devices in this model. Figure 13 shows the voltage profiles for each of the phases.

Volt-VAR
There were negligible reactive power contributions from the utility-owned 684-kW PV site because the average voltage at the PCC was within the VV dead band. However, the single-phase inverters were able to shrink the feeder voltage envelope, as shown in Figure 14.

Volt-VAR
There were negligible reactive power contributions from the utility-owned 684-kW PV site because the average voltage at the PCC was within the VV dead band. However, the single-phase inverters were able to shrink the feeder voltage envelope, as shown in Figure 14.

Extremum Seeking Control
ESC was conducted on the NG feeder model. In this feeder, there was significant phase imbalance and at the PV system: phase A voltages are close to nominal (~0.99-1.01 pu), phase B voltages are significantly above nominal (at times >1.05 pu), and phase C voltages are <0.95 pu-as shown in Figure 15. This presented a challenge for the ESC approach (like all the other voltage regulation methods) because the DER can only inject symmetric (positive sequence) reactive power. This means all the phase voltages are shifted in the same direction with a change in the DER power factor. Grouping the DER based on what phase they were interconnected was a good method of limiting the number of probing frequencies but allowed each of the phases to adjust their reactive power contributions independently. Unfortunately, this also produced a sizable voltage ripple on the power system because all the inverters on each phase had the same probing frequency. Ultimately, the parameters shown in Table 2 were used for ESC testing on the NG model, and the feeder voltage profile is shown in Figure 16.

Extremum Seeking Control
ESC was conducted on the NG feeder model. In this feeder, there was significant phase imbalance and at the PV system: phase A voltages are close to nominal (~0.99-1.01 pu), phase B voltages are significantly above nominal (at times >1.05 pu), and phase C voltages are <0.95 pu-as shown in Figure 15. This presented a challenge for the ESC approach (like all the other voltage regulation methods) because the DER can only inject symmetric (positive sequence) reactive power. This means all the phase voltages are shifted in the same direction with a change in the DER power factor. Grouping the DER based on what phase they were interconnected was a good method of limiting the number of probing frequencies but allowed each of the phases to adjust their reactive power contributions independently. Unfortunately, this also produced a sizable voltage ripple on the power system because all the inverters on each phase had the same probing frequency. Ultimately, the parameters shown in Table 2 were used for ESC testing on the NG model, and the feeder voltage profile is shown in Figure 16.       Figure 16. Minimum, maximum, and average bus voltages vs time for the ESC test, compared to baseline results controlling all PV inverters.

State-Estimation-Based Particle Swarm Optimization
The DER PF PSO optimization technique wrapped NG OpenDSS time-series simulations. The code was configured to execute once a minute to determine the optimal set point for all the DER devices. Using the same swarm size of 60, the solution often took 10 iterations, where it was capped. In those cases, it took longer than 60 s, so the optimization was configured to run every two minutes. As shown in Figure 17, optimizing 31 DER devices were difficult with a swarm size of 60 because the first two solutions produced results that pushed the minimum or maximum voltages outside the baseline envelope. However, after a few more solutions, the PSO significantly improved the voltage profile of the feeder for the remainder of the simulation. It should be noted that we did not attempt to optimize the PSO code to minimize run time.

State-Estimation-Based Particle Swarm Optimization
The DER PF PSO optimization technique wrapped NG OpenDSS time-series simulations. The code was configured to execute once a minute to determine the optimal set point for all the DER devices. Using the same swarm size of 60, the solution often took 10 iterations, where it was capped. In those cases, it took longer than 60 seconds, so the optimization was configured to run every two minutes. As shown in Figure 17, optimizing 31 DER devices were difficult with a swarm size of 60 because the first two solutions produced results that pushed the minimum or maximum voltages outside the baseline envelope. However, after a few more solutions, the PSO significantly improved the voltage profile of the feeder for the remainder of the simulation. It should be noted that we did not attempt to optimize the PSO code to minimize run time.

PHIL Results
To further validate the operational effectiveness of the PSO voltage regulation method, realistic DER power hardware-in-the-loop (PHIL) simulations were conducted with both the PNM and NG feeder models. These simulations provide better fidelity because they show ProDROMOS can work with real PV communications systems and ramp rates.

PNM Baseline
To validate the PHIL setup for the PNM feeder model, a comparison was made between active power, reactive power, and bus voltage for the RT and PHIL simulations. Figure 18 shows the voltage profile for the RT and PHIL simulations. The positive and negative reactive power changes around 1000 seconds in the baseline case were due to setting the PF setpoint on the 10-MW system to +0.85 and −0.85 to verify communications and that the PHIL simulation was closed-loop. Overall, the RT and PHIL simulations are closely matched.

PHIL Results
To further validate the operational effectiveness of the PSO voltage regulation method, realistic DER power hardware-in-the-loop (PHIL) simulations were conducted with both the PNM and NG feeder models. These simulations provide better fidelity because they show ProDROMOS can work with real PV communications systems and ramp rates.

PNM Baseline
To validate the PHIL setup for the PNM feeder model, a comparison was made between active power, reactive power, and bus voltage for the RT and PHIL simulations. Figure 18 shows the voltage profile for the RT and PHIL simulations. The positive and negative reactive power changes around 1000 s in the baseline case were due to setting the PF setpoint on the 10-MW system to +0.85 and −0.85 to verify communications and that the PHIL simulation was closed-loop. Overall, the RT and PHIL simulations are closely matched.

PNM Particle Swarm Optimixation PHIL
The only change from the RT simulations was that commands were communicated with the physical device to change the power factor and the DER was connected to the power simulation using a PHIL interface, as opposed to DBUS. The 258-kW system was recreated by scaling a 3-kW DER device. It was found that the physical equipment recreated the power profile reasonably well. Slight differences in the active power are due to efficiencies of the devices and slightly oversizing the simulated PV system in the Ametek PV simulator. The PSO solutions were repeatable for multiple runs. There were some deviations in the reactive power contributions from the PSO solutions, but overall, they matched well, as shown in Figure 19. Do note the PHIL PSO simulation was not started for approximately 10 minutes at the beginning of the RT-Lab simulation, which is why it takes some time for the reactive power on the physical DER to change. Similar results for the NG RT and PHIL simulations are described in [26].

PNM Particle Swarm Optimixation PHIL
The only change from the RT simulations was that commands were communicated with the physical device to change the power factor and the DER was connected to the power simulation using a PHIL interface, as opposed to DBUS. The 258-kW system was recreated by scaling a 3-kW DER device. It was found that the physical equipment recreated the power profile reasonably well. Slight differences in the active power are due to efficiencies of the devices and slightly oversizing the simulated PV system in the Ametek PV simulator. The PSO solutions were repeatable for multiple runs. There were some deviations in the reactive power contributions from the PSO solutions, but overall, they matched well, as shown in Figure 19. Do note the PHIL PSO simulation was not started for approximately 10 min at the beginning of the RT-Lab simulation, which is why it takes some time for the reactive power on the physical DER to change. Similar results for the NG RT and PHIL simulations are described in [26].
Energies 2020, 13, x FOR PEER REVIEW 3 of 4 Figure 19. Comparison between simulated and PHIL baseline minimum, maximum, and average bus voltages for the PNM model.

Discussion
Each of the voltage regulation methods have their advantages and disadvantages in terms of communications and computational overhead, implementation complexity, and availability in DER equipment. Overall, each of them can reduce the voltage deviation from nominal and maintain the feeder well within ANSI Range A voltage limits. For the PNM model, the DER PCC voltages are barely outside the dead-band, so the DER does not absorb much reactive power. The ESC method produces a DER-reactive power probing signal that causes a voltage ripple on the feeder, but it allows the DER to track the optimal reactive power setpoint well. The PSO method is the most complicated to implement, but issues optimal set points to each DER every minute. A comparison of each of the voltage regulation methods is shown in Figure 20. The simulation results indicate ESC and PSO can regulate the voltage closer to nominal, compared to VV.
To better understand the differences in these approaches, an analytical score was developed to summarize the effectiveness of each voltage regulation method, and a best score was calculated where the voltage regulation approach drove the solution to 1 pu: where vbl is the baseline voltage, vnom is the nominal voltage (1 pu), vreg is the voltage from the voltage regulation method, T is the time period of the simulation, b is the bus, and t is the simulation time.
The scores representing the average voltage improvement for all buses averaged over a four-hour simulation period in units of pu. Table 3 summarizes the effectiveness for the PNM model of each approach per phase as well as the average of each phase, calculated with: Figure 19. Comparison between simulated and PHIL baseline minimum, maximum, and average bus voltages for the PNM model.

Discussion
Each of the voltage regulation methods have their advantages and disadvantages in terms of communications and computational overhead, implementation complexity, and availability in DER equipment. Overall, each of them can reduce the voltage deviation from nominal and maintain the feeder well within ANSI Range A voltage limits. For the PNM model, the DER PCC voltages are barely outside the dead-band, so the DER does not absorb much reactive power. The ESC method produces a DER-reactive power probing signal that causes a voltage ripple on the feeder, but it allows the DER to track the optimal reactive power setpoint well. The PSO method is the most complicated to implement, but issues optimal set points to each DER every minute. A comparison of each of the voltage regulation methods is shown in Figure 20. The simulation results indicate ESC and PSO can regulate the voltage closer to nominal, compared to VV.
To better understand the differences in these approaches, an analytical score was developed to summarize the effectiveness of each voltage regulation method, and a best score was calculated where the voltage regulation approach drove the solution to 1 pu: where v bl is the baseline voltage, v nom is the nominal voltage (1 pu), v reg is the voltage from the voltage regulation method, T is the time period of the simulation, b is the bus, and t is the simulation time.
The scores representing the average voltage improvement for all buses averaged over a four-hour simulation period in units of pu. Table 3 summarizes the effectiveness for the PNM model of each approach per phase as well as the average of each phase, calculated with: Impact = score best score (9) VV slightly improved the average bus voltage, with an impact percentage of 12.3%. Implementing ESC and PSO demonstrated a larger impact on the feeder, with an average impact percentage of 74.5% and 73.7% respectively. A more detailed look at the impact of this regulation approaches per phase, shows that, since the system is relatively balanced, each phase is affected equally.
Energies 2020, 13, x FOR PEER REVIEW 3 of 4 = (9) VV slightly improved the average bus voltage, with an impact percentage of 12.3%. Implementing ESC and PSO demonstrated a larger impact on the feeder, with an average impact percentage of 74.5% and 73.7% respectively. A more detailed look at the impact of this regulation approaches per phase, shows that, since the system is relatively balanced, each phase is affected equally.  A comparison of the bus voltages for the NG model is illustrated in Figure 21 for each of the control methods. VV, ESC, and PSO all collapse the voltage envelope toward nominal voltage. In this case, ESC slightly outperforms the other two methods, shown in Figure 21. Table 4 summarizes the effectiveness for the NG model of each method. Interestingly, all the methods caused phase A (that was close to nominal to start) to deviate from the nominal voltage. The phases of the distribution system are coupled; in order to reach the global minimum for the system, and substantially improve the voltages on phases B and C, there is an adverse impact on phase A. This results in phase A deviating further from the nominal voltage.  A comparison of the bus voltages for the NG model is illustrated in Figure 21 for each of the control methods. VV, ESC, and PSO all collapse the voltage envelope toward nominal voltage. In this case, ESC slightly outperforms the other two methods, shown in Figure 21. Table 4 summarizes the effectiveness for the NG model of each method. Interestingly, all the methods caused phase A (that was close to nominal to start) to deviate from the nominal voltage. The phases of the distribution system are coupled; in order to reach the global minimum for the system, and substantially improve the voltages on phases B and C, there is an adverse impact on phase A. This results in phase A deviating further from the nominal voltage.

Conclusions
The novel RT and PHIL platform and scoring criterion introduced in this paper allow for the development of a standardized methodology to tune, evaluate, and compare different inverter-based distribution voltage regulation approaches on different feeder models. Utilities and distribution system operators generally do not have the sensor infrastructure or DERMS communication network to execute centralized control of DER. However, as the number of measurement devices and interoperable DER increases, it will become possible to calculate power system states and calculate optimal DER setpoints to provide voltage regulation and provide protection assurance on the system. These capabilities were demonstrated in this project with both SIL and PHIL simulations.
This project investigated and compared three voltage regulation approaches: volt-VAR, extremum seeking control, and particle swarm optimization. Each of the approaches were shown to help provide voltage support on a feeder with symmetrically elevated phase voltages. In the case of the imbalanced NG feeder voltages, the approaches were ineffective when employed with only the three-phase inverter, but an aggregate of single-phase devices were shown to improve the NG feeder voltage profiles. This was demonstrated using RT simulations with simulated PV devices and power hardware-in-the-loop simulations with a 3-kW PV inverter. The extremum seeking control voltage regulation technique was shown to be effective at controlling groups of DER devices to improve feeder voltages, even in cases of phase imbalance. This is the first reported demonstration of this ESC

Conclusions
The novel RT and PHIL platform and scoring criterion introduced in this paper allow for the development of a standardized methodology to tune, evaluate, and compare different inverter-based distribution voltage regulation approaches on different feeder models. Utilities and distribution system operators generally do not have the sensor infrastructure or DERMS communication network to execute centralized control of DER. However, as the number of measurement devices and interoperable DER increases, it will become possible to calculate power system states and calculate optimal DER setpoints to provide voltage regulation and provide protection assurance on the system. These capabilities were demonstrated in this project with both SIL and PHIL simulations.
This project investigated and compared three voltage regulation approaches: volt-VAR, extremum seeking control, and particle swarm optimization. Each of the approaches were shown to help provide voltage support on a feeder with symmetrically elevated phase voltages. In the case of the imbalanced NG feeder voltages, the approaches were ineffective when employed with only the three-phase inverter, but an aggregate of single-phase devices were shown to improve the NG feeder voltage profiles. This was demonstrated using RT simulations with simulated PV devices and power hardware-in-the-loop simulations with a 3-kW PV inverter. The extremum seeking control voltage regulation technique was shown to be effective at controlling groups of DER devices to improve feeder voltages, even in cases of phase imbalance. This is the first reported demonstration of this ESC application. The particle swarm optimization approach worked well when the OpenDSS feeder matched the RT/PHIL simulation environment. The technique incorporated live telemetry and PV forecast data to generate the optimal PF setpoints for a collection of PV systems over a fixed time horizon.

Patents
A U.S. Provisional Patent titled "Digital Twin Advanced Distribution Management System (ADMS) and Methods" was filed on 9 March 2020 based on this work.