Net-zero energy communities (NZECs) are critical to assuring the sustainability and resilience of modernized power systems. System model ing help s overcome technical challenges in designing and operating NZECs. In this paper , we present an open -source NZEC virtual testbed in Modelica based on a real NZEC in Florida, USA . This testbed consists of t wo sets of models : (1) higher -fidelity physics - based models that consider the interaction between subsystems of the studied NZEC and capture fast dynamics , and (2) lowe r-fidelity data -driven models that requir e fewer resource s to establish and/or run . All models are validated against measurement s from this real NZEC . In addition , this testbed includes a simulation framework that streamlines the processes for simulation and thus allows the use of developed models to form a virtual testbed . To demonstrate the usage of the virtual testbed , a case study is conducted where a building -to-grid integration control is evaluated via simulation. The evaluation results suggest that the tested control significantly smooth s the power draw of the studied community and does not sacrifice thermal comfort to a great exte nt. Keywords: Net -Zero Energy Communities , Virtual Testbed , Modelica , Building -to- Grid Integration 2 1 Introduction From an energy perspective, a c ommunity is a set of buildings that are served by the same energy systems, e.g., electricity generation or space heating/cooling systems. Achieving energy -efficient communit ies is critical in promoting environmentally sustainable and resilient cities (Sebastiano Marasco 2012, Omidreza Saadatian 2013) . One way to do so is through net -zero energy communities (Sesil Koutra 2018) . According to Peterson, Torcellini et al. (2015) , “A net-zero energy community is an energy -efficient community where , on a source energy basis , the actual annual delivered energy is less than or equal to the on -site renewable exported energy .” A net-zero energy community ( NZEC ) usually consists of energy generation fro m renewable sources (e.g. , solar , wind) ; energy distribution through electric - and thermal -distribution systems ; energy usage by space conditioning , lighting , and a wide range of miscellaneous loads ; energy storage by electric cars and other devices ; and interaction with the electric power system. Achieving net -zero energy at the end-use side is critical to realiz ing the sustainability and resilience of modernized power system s (Hachem -Vermette, Guarino et al. 2018) . Traditionally , net-zero energy is mainly studied at the b uilding level because buildings are major consumers of electricity (Lu, Wang et al. 2015) . However , NZE Cs provide significant economic and environmental advantages over net - zero energy buildings , such as greater flexibility to accommodate large energy infrastructure s and richer load diversity to reduce the design capacity of energy systems (Hadia Awad 2018) . Because of those advantages , various efforts in designing and building NZECs have been reported . One example is the West Village , California , USA. (Gaiser and Stroeve 2014) . The West Village is an all -electric campus that consists of solar photovoltaic (PV) arrays , a bio -digester generator , and heat pumps w ith backup electric heaters. The PV arrays produce the main power , the generator provides additional electricity supply , and the heat pumps generat e the domestic hot water (DHW). Although the West Village is intended to achieve the NZEC goal , according to its annual report (UC Davis West Village 2014) , its renewable energy generation only cover ed 82% of its energy demand during 2013 –2014. Another example is the Historic Green Village (HGV) on Anna Maria Island , Florida , USA . The HGV is a small NZEC consisting of five mixed -use commercial buildings (retail , including a bakery and general store ; office ; and residences ). The community achieved its net-zero energy go al in 2014. Section 2 of th is paper provides a detailed introduction to the HGV. 3 Despite these encourag ing application s, designing and operating a n NZEC is challenging , primarily because of the coupling relationships between subsystems within NZECs. An NZEC is a complex combination of multiple subsystems , usually including renewable energy generation ; heating , ventilation , and air conditioning (HVAC) ; refrigeration ; DHW ; and other miscellaneous equipment. These subsystems interact with one another , as elaborated by Hye Soo Suh (2019) . For example , to reduce the energy use d by the DHW system , the waste heat from the HVAC system is used to preheat the water supply for the DHW system (Miranda, Huang et al. 2015) . Thus , when evaluat ing the performance of the DHW system , it is necessary to simultaneously consider the HVAC operations. In addi tion, an NZEC typically need s to achieve multiple objectives (Hirushie Karunathilake 2019) , such as reducing emissions , minimizing demand charges , and accommodating grid supply interruptions. These obje ctives pose challenges to the design process because a multi -objective optimization needs to be solved to make trade -offs among these requirement s, especially when sizing each subsystem to assure economic feasibility. These objectives also complicate the operation process. For example , to minimize demand charges with energy storage capacity , research suggests that building system control parameters should be continuously updated according to the operation conditions (Muruganantham, Gnanadass et al. 2017) . System model ing is widely used in the literature to better understand the characteristics of community scale energy systems , as summarized by Andrew Lyden (2018) . Examples of applying system m odelling to support the design and operation of NZECs are listed in Table 1. However , gaps are also identified in existing system mode ling. For example , the interaction between systems from different domains tend s to be ignored . Because an NZEC includes heterogeneous , overlapping subsystems that span multiple , distinct physical domains , it is necessary to model all the domains as a whole to understand their interactions and interdependencies (Wetter 2011) . In addition , fast dynamics are sometimes excluded from the model ing. Given the importance of controls to the system operation , modeling of the transient control processes is indispensable (Wetter 2011) . Therefore , to better evaluate how the above factors influence the design and operation of NZECs , an accurate representation of the fast dynamics is necessary . Furthermore, validation results reg arding those models are seldo m reported in the literature , limiting the usability of those models . 4 Table 1. Review of current research for NZEC Studied Community Modeling Scope Considering cross -disciplinary interactions Fast Dynamic Modeling ? A dummy community , Belgium (De Coninck, Baetens et al. 2014) DHW system (Modelica (Modelica Association) ]) Only models for the DHW system Yes West Village , California , USA (Gaiser and Stroeve 2014) - Building system (eQuest) - PV system (PolySun and SAM) - DHW system (Building America Report) Separate models for building , electrical , and DHW system s No Qingshan Lake District in Hangzhou , China (Lu, Yu et al. 2014) - Combined heat and power system - Small -scale wind turbine system - Solar photovoltaic system - Solar thermal heater system (MATLAB (MathWorks) ]) Coupled models for building and electrical systems No A dense neighborhood and a low - density suburban neighborhood in Neupré , Belgium (Marique and Reiter 2014) - Space heating , cooling , and ventilation system (in -house code) - Appliances , cooking , and DHW system (estimated based on the number of inhabitants in the building) - On-site energy production system (Townscope) Separate models for building and electrical systems No A dummy community , Calgary , Alberta, Canada (Hachem -Vermette, Cubi et al. 2015) - DHW , heating/cooling , PV system (TRNSYS) - Lighting and equipment system (EnergyPlus) Coupled models for thermal and electrical systems Yes A village in Graubünden , Switzerland (Orehounig, Mavromatidis et al. 2014) Building System (EnergyPlus) Only models for building system Yes Dummy communities located in 5 climate zones (Vázquez -Canteli, Kämpf et al. 2019) - DHW, cooling dem and (EnergyPlus ) - Heat pump, energy storage, battery (reduced -order models ) Separate models for building and electrical systems Yes Two residential buildings located in downtown Austin, TX (Vázquez -Canteli, Ulyanin et al. 2019) - Space heating, cooling syste m (CitySim) - Lighting system (CitySim) Separate models for building and electrical systems Yes In this paper, we present an open -source NZEC virtual testbed in Modelica based on a real NZEC , the HGV. This testbed consists of two sets of models : (1) higher -fidelity physics -based 5 models that consider the interaction between subsystems of the studied NZEC and captur e fast dynamics , and (2) lower -fidelity data -driven models that requir e fewer resource s to establish and/or run. All models are validate d against measurements from this real NZEC. The open -source NZEC testbed has been publicly released at https://bitbucket.org/sbslab -zuo/scc -nzec. In addition, this testbed includes a simulation framework that streamlines the simulation processes and thus allows the use of developed models to form a virtual testbed. To demonstrate the usage of the virtual testbed, a case study is conducted where a building -to-grid integration control is evaluated via simulation. This paper contributes to the literature as follows . First, we develop and validate an open - source NZEC virtual testbed. This testbed can meet different requirements for designing and operating NZECs and provides great flexibility in terms of model selection and integration for simulat ion and optimization. Specifically , This NZEC virtual testbed allows consider ing three types of buildings: office, residential, and retail buildings ; and one kind of distributed energy resource : PV panels . Second , we performed a comprehensive evaluation on a demand -response control strategy with this virtual testbed . Compared to similar evaluations in the literature, this one is more comprehensive as it employs more detailed system models . This evaluation reveals new issues that have not been reported in the literature before to our best knowledge . The rest of the paper is organized as follows . First, we provide a detailed introduction to the HGV . We then discuss how we realize the virtual testbed with different models and a software framework . After that , we present the evaluation o f a building -to-grid integration control with the proposed virtual testbed . Finally , we present our conclusions and discuss future work. 2 Historic Green Village 2.1 Overview Located on Anna Maria Island , Florida , USA, the HGV is a small community (2 ,776 m2) that consists of two retail stores , one commercial office , and two residences , as shown in Figure 1. The climate of Anna Maria Island is categorized as Zone 2A , which is hot and humid (ASHRAE 2007) . The aver age solar radiation is around 19.33 W/m2 (Wilcox 2012) , providing a significant resource for renewable energy. In addition , by employing aggressive measures in energy efficiency (e.g. , 6 cool roofs, Low-E glass , extra insulation , and waste heat recovery from water -source heat pumps) , the HGV achieved its net-zero energy goal in 2014 after operat ing for two and a half year s (beginning in June of 2011) , as shown in Figure 2. Figure 1. HGV on Anna Maria Island, FL . The HGV has three major energy systems: a solar PV electric generation subsystem , a solar thermal DHW subsystem , and a ground source heat pump (GSHP) subsystem. The solar PV system generates electricity during daylight hours. Excess generation is sold to the utility provider and purchased from the utility as needed. The GSHP subsystem provides both cooling and heati ng to all buildings in the community. During cooling -dominated periods , the GSHPs also supplement the DHW requirements. The solar thermal DHW subsystem provides domestic water. Excess hot water generated during the day is stored in a tank. The building sub system represents the energy usage of the five buildings for different activities. All five subsystems are interconnected. For instance , the buildings are conditioned by heat pumps , which are powered by electricity from either the PV panels or the utility. Figure 2. The electricity generation and demand of the HGV . 37,301102,828119,448 118,096 43,141105,133150,485 116,582 2011 2012 2013 2014[kWh] Annual Electricity Generation Annual Electricity Demand 7 2.2 Energy Subsystems 2.2.1 Solar PV Generation The solar PV system generate s on-site electric power. To achieve net-zero energy , PV panels are installed on building roofs , two solar carports , and a neighboring warehouse building. The configuration of the PV system is shown in Table 2. In total , the PV system is 612.09 m2 with a total peak capacity of 60 kW. Table 2. PV system specifications . Building Name PV Group PV Panel Area [m2] Array Tilt [o] Array Azimuth [o] Rated Module Efficiency [%] Rated DC/AC Conversion Efficiency Thelma A 15.36 30.0 -45 14.9 0.8 15.36 15.0 45 14.9 0.8 Rosedale B 46.08 5 -45 14.9 0.8 38.40 5 -45 14.9 0.8 25.60 30 45 14.9 0.8 Sears C 51.20 20.0 -45 14.9 0.8 Pickle Fish D 29.34 5.0 -45 16.0 0.8 Pillsbury E 17.93 10.0 45 16.0 0.8 19.56 5.0 45 16.0 0.8 13.04 5.0 -45 16.0% 0.8 9.78 10.0 45 16.0 0.8 North Carport NC 84.48 5.0 -45 14.9 0.8 South Carport SC 96.00 5.0 -45 14.9 0.8 Warehouse W 149.96 5.0 45 16.0 0.8 AC is alternating current; DC is direct current. 8 2.2.2 Utility -provided Power and Buildings The utility subsystem is treated as a virtual battery. If the solar panels meet the electricity load of the NZEC , the excess electricity “charge s” the utility subsystem by feeding the electricity to the power grid; if on-site production is insufficient to meet the HGV load , the utility subsystem “discharge s” the electricity by drawing electricity from the power grid to the HGV , as shown in Figure 2. By having the utility as a virtual battery with inexhaustible capacity , one avoids the high cost of installing on -site energy storage systems. In addition, each building in the HGV is equipped with miscellaneous electric equipment, including lighting, refrigeration, food preparation (dishwasher, coffee maker, and dryer), hot water, and miscellaneous equipment and plug l oads. Figure 3. Schematic of the utility distribution system . 2.2.3 Ground Source Heat Pump As shown in Figure 4, the GSHP subsystem consists of a water loop that connects nine water -to- air heat pumps in the buildings , five refrigerator racks , and a heat exchanger linked to a ground source geothermal sy stem. The nine heat pumps provide cooling and heating to the buildings. The heat pumps’ cooling capacities vary from 5.4 kW to 17.3 kW. The total cooling capacity of the nine heat pumps is 87.3 kW. Depending on the building cooling/heating loads , some buil dings 9 may have two or three heat pumps. On the waterside , each heat pump has two dedicated circulating water pumps. This water loop also supplies direct cooling for the five refrigeration units. The GSHP sub system consists of one sink well and one source w ell. Those wells penetrate a layer of limestone rock at 137 m (450 ft) deep , where the underground water has a constant temperature of 26.7 °C (80 °F). A variable speed well pump controls the flow rate of the extracted groundwater through the heat exchanger to maintain the temperature of the outlet water from the heat exchanger in the water loop. In cooling mode , when the outlet water temperature exceeds 28.3 C (83 °F), the well pump speed s up from 20% to 100% of the maximum speed to reduce the water temperature. In heating mode , when the outlet water temperature is below 23.9 C (75 F), the well pump speed s up from 20 % to 100 % of the maximum speed to increase the water temperature. 2.2.4 Domestic Hot Water The DHW system serves one retail store (Pickle Fish; Table 2). In this system , the makeup water from the municipal water is preheated by a solar collector to around 37.8 °C (100 °F) and is then stored in a tank with a capacity of 0.6 m3 (120 gallons). When makeup water is needed , the water is pumped to another tank where heat is recovered from the two GSHPs (HP8 and HP9) and Figure 4. Schematic of a GSHP subsystem in the HGV . Circulating PumpsGSHP1 GSHP2 GSHP3 GSHP4 GSHP5 GSHP6 GSHP7 GSHP8 GSHP9 Refrigerator Rack Heat ExchangerSink WellSource Well 10 absorbed by the makeup water. When the recovered heat is not available , an electric backup heater is activated to provide more heat. 3 System Modeling 3.1 Model Architecture As mentioned above, t he HGV contains different subsystem s. We package t he models within each subsystem into one subsystem model , while t he energy/mass/information connections between different subsystems are then modeled with corresponding information buses . The decomposition of the HGV results in a model architecture , as shown in Figure 5. This architecture has three advantages : (1) subsystem models can be easily isolated for calibration or test ing purposes ; (2) connections bet ween subsystems can be managed in a flexible and scalable way ; and (3) each subsystem can be assigned a dedicated solver so that each subsystem can have unique simulation time steps to capture internal dynamics at lower computational costs. Figure 5. The exchange and communication of different subsystem models in the HGV . 3.2 Model Methods It is generally preferable to have high -fidelity physics -based models when detail ed information of Domestic Hot Water SubsystemGrid Subsystem Ground -coupled Heat Pump Subsystem Buildings Subsystem BuildingSolar HeaterPV Panel Bus BusBus BusBus Bus Bus Heat Exchanger BoreholeBus Bus Heat PumpBus Electricity Domestic Hot Water Heating /Cooling Air Recovered HeatHeating /Cooling Water Renewable Energy Subsystem Wind TurbineGrid 11 the studied systems is required . Having such models , however , may not always be practical. Physics -based models require a detailed understanding o f the studied system and impose a significant computational demand (Braun and Chaturvedia 2002) . In a real -world setting , the necessary conditions may be difficult to meet. For example , occupancy has a significant influence on building energy demand and thus should be considered in load modeling. However , the stochasti c nature of occupancy in time and space (Feng, Yan et al. 2015) makes it difficult to model accurately. Computational fluid dynami cs models are effective in mode ling thermal comfort level s but are very computational ly intens ive. This limits their applicability to real -world adaptive controls (Zuo and Chen 2010) . In addition , physics -based model s may require detailed operating data to calibrate t he parameters. Such data may be inaccessible due to a lack of meters or sensors. The data -driven models , on the other hand , are often more suitable when resource s are limited . They can be built without any information regarding physics. They are usually computationally efficient because they often contain only linear equations. Further , they can be developed with limited types of data , although the accuracy of the models is directly related to the number of training data points available (Braun and Chaturvedia 2002) . Given the pros and cons of physics -based models and data -driven models , we believe it is desirable to have a hybrid model that includes both physics -based models and data -driven models. The physics -based models can be used to model the most critical subsystems , or those with more known and detailed operating data , when computational resource s permit. The data -driven models can be used to model systems that lack either a good physics -based model or detailed operating data, or where computational performance is a limiting factor. In this study , we used Modelica to implement the physics -based models. Modelica is an equation -based , object-orient ed modeling language for dynamic systems. Examples of building -related modeling with Modelica include the modeling of building envelopes , a data center cooling system , a single -zone cooling system , a DHW system , and chiller plants (Nouidui, Phalak et al. 2012, Zuo, Wetter et al. 2014, Huang, Malara et al. 2016, Zuo, Wetter et al. 2016) . 3.3 System Models Applying the model structure and the model methods from the previous two sections , we created models for the entire HGV in a hierarchical fashion . Figure 6 shows the top -level model , which mainly includes the following five subsystem models: renewable energy generation , grid, 12 buildings , water -source heat pumps , and information bus. The top -level model is implemented in Modelica. The data and signals can be exchanged using the bus connectors that exist in Modelica. The following sections discuss the lower -level models in more detail . Figure 6. Diagram of top-level models for the HGV . (1) Renewable Energy Generation Figure 7 shows the model diagram for the entire PV subsystem , which is composed of models for the eight arrays shown in Table 2. For each PV array , we developed a dedicated model for each array . The direct current (DC) output of the array is converted to alternating current (AC) through an inverter. The AC output is agg regated to one output port terminal to connect with other subsystem models. In addition , the PV energy production information is transferred to other subsystem models for further control purpose s via the data bus , which is an expandable connector for data communication. Irradiation informationWeather bus Power Grid Information bus BuildingPV GSHP DHW 13 Figure 7. Diagram of Modelica models for the PV subsystem . We create both a physics -based model and a data -driven model for the PV array. For the physics - based model , we used the PVSimpleOriented model in the Modelica Buildings Library (Wett er, Zuo et al. 2014) . This model calculates the electrical power, 𝑃, generated by the PV panel using the following equation: 𝑃=𝐴𝑓𝑎𝑐𝑡𝜂𝐺𝜂𝐷𝐶𝐴𝐶 , (1) where 𝐴 is the area of each PV panel, 𝑓𝑎𝑐𝑡 is the fraction of the aperture area, 𝜂 is the PV panel efficiency, 𝐺 is the total solar irradiation, and 𝜂𝐷𝐶𝐴𝐶 is the efficiency of the conversion between DC and AC. In this model, 𝐺 is the sum of direct and diffuse irradiation and is calculated by : 𝐺=𝐺𝐷𝑖𝑓+𝐺𝐷𝑖𝑟, (2) where 𝐺𝐷𝑖𝑓 and 𝐺𝐷𝑖𝑟 are the sum s of the diffuse and direct irradiation on a tilted surface , respectively. 14 𝐺𝐷𝑖𝑟=𝑀𝑎𝑥(0,𝑐𝑜𝑠⁡(𝑖𝑛𝑐𝐴𝑛𝑔)𝐻𝐷𝑖𝑟𝑁𝑜𝑟 , (3) where 𝑖𝑛𝑐𝐴𝑛𝑔 is the solar incidence angle on the surface and 𝐻𝐷𝑖𝑟𝑁𝑜𝑟 is the direct normal radiation. 𝐺𝐷𝑖𝑓=𝐺𝑆𝑘𝑦𝐷𝑖𝑓+𝐺𝐺𝑟𝑜𝐷𝑖𝑓 , (4) where 𝐺𝑆𝑘𝑦𝐷𝑖𝑓 and 𝐺𝐺𝑟𝑜𝐷𝑖𝑓 are the hemispherical diffuse solar irradiation on a tilted surface from the sky and hemispherical diffuse solar irradiation on a tilted surface from the ground, respectively. For the data -driven model , we employed a forward backpropagation neural network model (Hecht -Nielsen 1989) . As shown in Figure 8, there are four neurons ( 𝑥1,⁡𝑥2,⁡𝑥3,⁡𝑥4) and one output layer with one output ( 𝑦1) in this neural network model. In addition , there is also one hidden layer with 12 neurons. For the PV panel modeling , we select ⁡𝑥1,⁡𝑥2,⁡𝑥3, and⁡𝑥4 as the direct normal radiation , diffuse horizontal solar radiation , global horizontal radiation , and the power generation of the PV panel a half -hour before the current time , respectively. Variable 𝑦1 is the power generation of the PV panel at the studied time. This neural network mode l is implemented with MATLAB and more detailed settings are listed in Table 3. To enable data exchange between MATLAB and the Modelica environment during runtime , we utilize the Python interface in the Modelica Buildings Library and the MATLAB Engine API for Python (Mathworks 2017) . Table 3. Detailed settings of the neural network model for the PV panel Number of input layer nodes 4 Number of output layer nodes 1 Number of Hidden layer nodes 12 Neuron transfer function logsig, purelin Learning function Learngd（Gradient descent weight and bias learning function） Performance function Msereg（weighted mean square deviation ） Network training function Trainbr (Bayesian regularization) Weight initialization method Initnw（Nguyen -Widrow method ） Maximum training time 2000 Epochs Target error 10-3 15 Figure 8. The structure of the neural network model . (2) Power Grid In this research , we used the Electrical.AC.OnePhase.Sources.Grid model in the Modelica Buildings Library to build the physics -based model for the utility -supplied power. The input for this model is a fixed voltage signal , while the output will be the power supplied by the utility at a point of common coupling (PCC) to the HGV electric distribution system. We define a convention that the power is positive if real power is consumed from the grid and negative if it is flow ing out of the PCC. We did not build a data -driven model for the utility distribution system because we do not have sufficient historic al data to characterize the distribution system behavior. (3) Buildings The equipment in the HGV can generally be categorized as having either easy load pattern recognition or difficult load pattern recognition . As an example of an easily recognizable pattern , Figure 9 shows the load profile of the lighting equipment in a retail building for 20 days. We can see that all the days share the same pattern: from 9:00 to 17:00 , the lighting load is maintained around 300 W, and the rest of time, the lighting is close to zero. From this , we can easily build simple data -driven models for such equipment. For an example of a load pattern that is difficult to understand , Figure 10 shows the load profile of fridges in a restaurant building for 20 days. We 16 can see that the fridge power load changes randomly , mainly due to the cycling control of the compressors . For simplification , we assume the loads for this equipment are always constant. Figure 9. Electricity load profiles of lighting equipment in a retail building (each line represents one day of data) . Figure 10. Electricity load profiles of fridges in a restaurant building (each line represents one day of data) . 17 (4) GSHP Subsystem The GSHP subsystem is the major energy load type in the HGV. We create both physics -based and data -driven model s for the GSHP subsystem. When creating the physics -based model , we respect the real topology of the GSHP subsystem , which results in a top -level model for the GSHP subsystem , shown in Figure 11. In this top -level model , the major components include a weather bus, internal heat gain , a condenser water loop , a heat recovery loop , zone controllers , and GSHPs. The weather bus component provides external disturbances , while the internal heat gain component estimates internal disturbances in the GSHP operation. The condenser water loop component and the heat reco very loop component are modeled as constant flow loop s with constant speed pumps. In addition , the GSHP module contains sub -models , as shown in Figure 12, including one heat pump model and the thermal zone model. Both models are from the Modelica Buildings Library ; readers can find more information in (Wetter, Zuo et al. 2 014). Specifically, the thermal zone model calculates the power of each GSHP using 𝑃ℎ𝑝=𝑃ℎ𝑝,𝑛𝑜𝑚𝑃𝐿𝑅∙𝐶𝑂𝑃𝑛𝑜𝑚/(𝑇𝑒𝑣𝑎 𝑇𝑐𝑜𝑛−𝑇𝑒𝑣𝑎𝜀𝑐𝑎𝑟𝑛𝑜𝑡𝜀𝑃𝐿𝑅(𝑃𝐿𝑅)), (5) where 𝑃𝑐ℎ,𝑛𝑜𝑚 is the nominal power of the GSHP , 𝑃𝐿𝑅 is the partial load ratio (the ratio of the cooling load handled by the GSHP to its nominal capacity), 𝐶𝑂𝑃𝑛𝑜𝑚 is the GSHP’s coefficient of performance at the nominal condition, and 𝑇𝑒𝑣𝑎 and 𝑇𝑐𝑜𝑛 are the temperatures in the evaporator and condenser sides of the GSHP , respectively. The ⁡𝜀𝑐𝑎𝑟𝑛𝑜𝑡 is the Carnot effectiveness (assumed to be constant) , and 𝜀𝑃𝐿𝑅 is the operation effectiveness at partial loads, which is a function of the PLR: 𝜀𝑃𝐿𝑅(𝑃𝐿𝑅)=𝑐1+𝑐2𝑃𝐿𝑅+𝑐3𝑃𝐿𝑅2+(1−𝑐1−𝑐2−𝑐3)𝑃𝐿𝑅3, (6) where 𝑐1,⁡𝑐2,⁡𝑐3 are constant coefficients to mimic the internal capacity control of each GSHP . 18 Figure 11. Diagram of Modelica models for the GSHP subsystem . Figure 12. Diagram of Modelica models for the GSHP module . Weather bus Zone Controller GSHP Heat recovery loopCondenser water loopPerformance Dataset Internal heat gain Thermal zone model Heat pump 19 In addition , we employ a neural network model to establish a data-driven model for the GSHP subsystem. The structure of the neural network model for the geothermal subsystem is the same as that for the PV subsystem , as shown in Figure 8. However , 𝑥1,⁡𝑥2,⁡𝑥3, and⁡𝑥4 are different because t hey represent the day type , the hour index in one day , the outdoor dry bulb temperature , and the power of the heat pump in the previous hour, respectively. Table 3 shows the definition s of each day type. For more details regarding day type , refer to (Huang, Zuo et al. 2016) . We do not have a dedicated power meter for HP 4 , HP 5 , HP 8 , and HP 9 ; instead, there is one aggregated meter for HP 4 and HP 5 and one aggregated meter for HP 8 and HP 9 . Thus, we must model HP 4 and HP 5 as one heat pump and HP 8 and HP 9 as another heat pump. Table 4. Day type categor ies in the HGV . Index Day Category Name Description 1 Weekday in Peak Season Weekday from December 15 to May 14 2 Weekend in Peak Season Weekend from December 15 to May 14 3 Weekday in Non -Peak Season Weekday from May 15 to December 14 4 Weekend in Non -Peak Season Weekend from May 15 to December 14 (5) DHW Subsystem We create only a physics -based model for the DHW subsystem due to the lack of operation data for developing a data -driven model. We established the physics -based model based on the design specification and engineering judgment of the DHW system . Figure 13 shows the top -level model for the DHW subsystem . The major components includ e the city water module , the solar collector , the water usage profile , and the backup heater. We build the models for the above modules again with the Modelica Buildings Library. For simplicity , we assume water usage , pump power , and heater efficiency are constant. 20 Figure 13. Diagram of Modelica models for the DHW subsystem . 3.4 Validation We validate the developed model with data from the HGV . For data acquisition , we installed web servers with 12 –24 current transformers in each electric service panel. This allows us to log circuit - level power demand and generation and upload data to a cloud -based application that can aggregate the data by building and load type (HVAC , lighting , refrigeration , etc.). Solar radiation data was obtained from an authoritative source , the National Solar Radiation Database (NSRDB) (The Renewable Resource Data Center) , at a resolution of 4 km by 4 km using geostationary satellites. When validating the models, we use the developed system models to predict the power for a test period of 1 week (Aug 25 –31, 2014) for each subsystem . Then, we demonstrate the power prediction results at the community level. For the DHW system , no validation is conducted due to the lack of historical data . 3.4.1 Renewable Energy Generation Figure 14 shows the aggregated power results for the PV subsystem , which demonstrate that the data-driven models are quite close to the measurement s, except occasionally near noon . The prediction s from the physics -based models are significantly lower than the measurement s, and the deviation can be up to 17%. Water usage profileCity waterSolar collector Tank 1 Tank 2 Recovery water outlet Recovery water inletBackup heater 21 Figure 14. Comparison of the predicted and measured aggregated PV power generation . 3.4.2 Building System As shown in Figure 15, although the prediction can capture the general pattern of the measurement , relatively large deviations exist . This means that the simple method we used to model the equipment in the second category may be oversimplified . Figure 15. Comparison of the predicted and measured miscellaneous electric equipment power at the aggregated level . 22 3.4.3 GSHP Figure 16 shows the aggregated power results for the GSHP subsystem. In general , the prediction s by both the p hysic s-based models and the data -driven models match the measurement pattern well. However , for some time intervals , there are relatively large deviation s in both models . These deviations are up to 18% for the data -driven model and up to 30% for the physic s-based models . Figure 16. Comparison of the predicted and measured aggregated heat pump power . 3.4.4 Community Figure 17 shows the aggregated power results for the entire community from the data -driven models , but t he pump energy consumption in the DHW subsystem is not included in these results . Negative power means that power generation exceeds power demand , while positive power means that power generation is insufficient to meet power demand. At the community level, w e can see that the prediction closely matches the measurement. It is worth notin g that the aggregated power for the entire community turns out to be a typical “ duck curve,” i.e., after times of high solar 23 generation , the grid must rapidly increase power output around sunset to compensate for the drop in solar generation. The “ duck curve” may cause load -supply balance issues (Denholm, O’Connell et al. 2015) . An energy storage system or a peak load shaving or shifting strategy can resolve this issue. Figure 17. Comparison of the predicted and measured aggregated power for the HGV . 4 Simulation Framework As shown in Figure 18, this framework is designed to facilitate the use of the models and is composed of four modules: pre-processing, modeling, optimization, and post -processing. The pre- processing module process es the input data for the virtual testbed. For instance, gaps may exist during the data collection process. In this case, the pre-processing module is used to fill the gaps through interpolation to provide continuous input for the dynamic simulation. The modeling 24 module simulates all the models developed for the studied NZEC, while the optimization module performs model -based optimizations using the input data as well as the objective functions. Various optimization engines can integrated into this virtual testbed , such as the Python -based pyOpt framework (Perez, J ansen et al. 2012) , GenOpt (Wetter 2004) , and the MATLAB Optimizat ion Toolbox. Those optimization engines can be used to implement advanced control strategies, such as model prediction control, with arbitrary control objectives. Alternatively, those optimization engines can be also replaced with rules -based controllers. The post-processing module processes the results from simulation and optimization based on predefined metrics. This framework automates the process of identifying the optimal design and the optimal operation scheme for the HGV. Specifically, in the design phase, architects and engineers can compare the simulation results of various designs to further enhance the design. In the operation phase, facility managers can improve the operation by altering the optimal control variables or by changing demand -side management strategies. It is worth mentioning that we performed numerical tests on the proposed framework on a Desktop (IntelR CoreTM i9-10980XE CPU (3. 0GHz) and 32 GB memory, Window 10 64 -bite) and found an annual simulation (control sampling inter val: 15 min) takes 1 ,340 seconds. Figure 18. Simulation framework . 25 5 Case Study In this case study , we demonstrate the use of the developed models and the simulation framework as a virtual testbed for evaluati ng a building -to-grid integration method that adjusts the thermostat settings to smooth the power draw of NZECs from the power grid . 5.1 Testing Co ntrol Method The control method smooth s the power draw of NZECs from the grid , as defined in : 𝐽=argmin(∑|𝑃𝐷,𝑖−𝑃𝐺,𝑖|𝑛 𝑖=1), (7) where 𝑃𝐷,𝑖 is the power demand of NZECs and 𝑃𝐺,𝑖 is the power generation of NZECs at time step 𝑖. A smooth power draw helps the power grid achieve a balance between power supply and demand. As illustrated in Figure 19, this control method is a distributed control strategy and contains both campus -level con trol and building -level controls. The campus -level control takes the measured power draw as an input and generates a price of electricity at time step 𝑖, by: 𝜆𝑖=𝜆0+Δ𝜆(𝑘(𝑃𝐷,𝑖−𝑃𝐺,𝑖)+𝐼(∑(𝑃𝐷,𝑗−𝑃𝐺,𝑗)𝑖 𝑗=1)), (8) where 𝜆0 is the base price , 𝛥𝜆 is the maximum changing range for the price , and 𝑘 and 𝐼 are parameters. E quation ( 8) is essentially a proportional -integral control that uses the measured power draw as feedback to adjust the price signal. 26 Figure 19. Proposed control method . The building -level control establishes a mapping between the electricity price and the thermostat setpoint , as illustrated in Figure 20. Only the cooling mode is considered for simplicity. We can see that , in general , the thermostat setpoint is proportional to the electricity price. However , the thermostat setpoint is only allowed to vary in a range of [𝑇𝑠𝑒𝑡,0−𝑘𝛥𝑇,𝑇𝑠𝑒𝑡,0+𝑘𝛥𝑇], where 𝑇𝑠𝑒𝑡,0 is the nominal set point and 𝛥𝑇 is the maximum deviation from the nominal setpoint . Additionally , for each building , the parameter 𝑘 is determined based on occupant preference. If the occupant s highly value thermal comfort , then 𝑘 approach es the minimum value of 0; otherwise , 𝑘 is closer to the maximum value of 1. Figure 20. Mapping between the electricity price and thermostat setpoint . Campus Level Control Building Level Control 1 Building Level Control n Measured Power Draw Thermostat Setpoint for Building 1 Thermostat Setpoint for Building nElectricity Price Electricity Price Electricity PriceCooling Setpoint λ0-Δλ λ0+ ΔλTset,0 + kΔT Tset,0 -kΔT 27 5.2 Deployment of Testing Control Method in the Virtual Testbed To evaluate the proposed contr ol method , we deploy it in the virtual t estbed. Specifically , the control method is implemented in Python and coupled with the virtual t estbed via the socket connection. The virtual t estbed periodically (every 5 minutes ) sends the power information to the proposed c ontrol method , which then re turns the setpoint . Additionally , as shown in Figure 17, power generation is usually much larger than power demand in the HGV , meaning the PV subsystem is oversized and making it very difficult to achieve a smooth power draw . To provide a fair evaluation , we ignore the PV group s E, NC, SC, and W to reduce the capacity of PV generation. Furthermore , we employ the physics -based models of the renewable energy generation , power grid , GSHP , and DHW subsystems and the data -driven models for other subsystems in the evaluation. The simulation period is set to be a typical summer day , July 27 , 2014 , and measured weather data is used as the major input. The 𝑘 value for each GSHP is randomly selected between 0.6 and 0.9 to reflect diversity in occupant preference s. 5.3 Results and Discussion Figure 21 illustrates the power draw of the HGV with and without the control method. We can see that the power draw varies dramatically through out the testing day in the baseline scenario. The power dr aw reaches its minimum value of −7 kW around 11:00 and reaches its maximum value of 15 kW around 19:00. Specifically , the power draw becomes negative from 8:00 to 18:00 , which means the extra power generation must be absorbed by the power grid. In the control scenario where the proposed method is implemented , the control method achieves the control objective from approximately 8:00 to 18:00 because the power draw equal s zero during this period . To do this, the control method varies the setpoints of all GSHPs from 21.5 C to 23.5 C. Because of different 𝑘 values , however , each GSHP has a dedicated setpoint trajectory. For the remaining time, because there is very little or no power generated by the PV systems , the power grid is still larger than zero even though the setpoints are already set at the maximum values. This suggests that the HGV has little power flexibility that can be utilized during those periods . Figure 22 compares the results during the control scenario and the testing scenarios. We can see that the control scenario uses more power from 8:00 to 17:00 and less power for most of the remaining hours. 28 Figure 21. Power draw (top) and cooling temperature setpoint (bottom) under the b aseline and under control . As a result , the accumulated energy consumption in the control scenario is less than that in the baseline scenario before 11:00 and larger than that in the baseline scenario after 11:00 . It is worth mentioning that the energy consumption for the testing day in the control scenario is larger than that in the baseline scenario at the end of the testing day . This suggests that providing gird service may lead to more energy consumption by buildings . Figure 22. Power and energy consumption of the HGV under the b aseline and under control . 29 Figure 23 shows the details of how one GSHP operates under the control scenario and the baseline scenario. Under the baseline scenario , the zone temperature is quite close to the setpoint for most of the time. The speed ratio of the compressor varies from 0 to 0.9 , and the peak value occurs at around 12:00. On the other hand , in the control scenario , the setpoint is set to be 23.5 °C from 0:00 to 8:00. In that case , because the temperature is lower than the setpoint , the compressor is shut off. From 9:00 to 18:00 , the setpoint decreases , leading to a peak speed ratio of almost 1. After 17:00 , the setpoint increases again , and the compressor is shut off accordingly. Figure 23. Power and energy consumption of the HGV under the b aseline and under control . 6 Conclusion and Future Work In this paper , we proposed a virtual testbed for designing and operating NZECs. The virtual testbed consists of a series of system models and a software framework . The system models are developed with physics -based and data -driven modeling methods, respectively , and t he framework is designed to streamline simulation and optimization . We compared the prediction power results with the measurement s from a real NZEC , and the outcome suggest s that (1) in general , the system models achieve good accuracy in predicting power. For the community -level model , the model can track the trajectory of the power for the testing week while the R2 value is 0.95. However, the accuracy of the subsystem models is relatively lower. For example, The deviation of the prediction by the PV subsystem model from the measurement can be up to 17%. (2) For PV system modeling , 30 the data -driven models perform better than the physics -based models. The reason for this difference is the variabilit y of the information used in the models that should be relate d to model characteristics. (3) For miscellaneous electric equipment modeling , the load patterns of some equipment are randomly distributed. (4) The aggregated load and solar generation profiles of the HGV indicate that there is an inherent imbalance between energy demand and supply. Exploring opportunities to reduce coincident peak demand of electric and thermal loads can reduce the cost of achieving energy resilienc e at the community an d grid scale s. Beyond comparing the prediction and measurement results , we also applied the virtual testbed to evaluate building -to-grid integration control. Based on the evaluation results , apply ing this control in the HGV to address the “ duck curve” issu e is a promising avenue ; further more , the studied control respects the performance of the occupants , which can simplify control deployment. Compared to similar modeling tools in the literature, this virtual testbed has three major advantages: 1) it is the only modeling tool that is based on a real community and validated against the real measurement. Thus, this virtual testbed can better support the NZEC design and operation with more reliable test results. 2) it provides great flexibility in terms of model selection to meet the requirements of specific applications. This virtual testbed has a hierarchical model structure and includes both physic s-based and data -driven models. Users can customize the testing model based on the focus to balance the accuracy and complexity. 3) it simplifies the process for deploying design methods and control strategies to the simulation environment. This virtual testbed contains software components to streamline the process for running simulations and/or optimization. Therefore, parametric studies or model -based design/optimization can be easily set up. The research described in this paper focuses on establishing a foundation for explor ing the design and operations of a n NZEC. With additional resources , a series of studies can be conducted based on the testbed developed here . For example , control strategi es with other operational constraints, such as minimal operation/utility costs, minimal emissions, can be evaluated. Smart thermo stats and occupancy sensors would allow us to develop adaptive control strategies that maximize on -site solar energy consumptio n, demonstrate the ability to reduce coincident peak electric and thermal loads , reduce demand charges , and participate in demand -response programs. More distributed energy resources and energy storage devices can be added to the virtual tes tbed. 31 Detailed electricity models will be considered to model different operation modes of the community, including “island mode”. Occupancy -related studies should be conducted to help us better understand how occupancy affects the operation of building systems. In the future , it would be helpful to distinguish between weather and occupant effects for different load types. Additional operating data and smart controls can allow system designers to correctly size mechanical systems and energy (thermal and chemical) storage systems. Facility operators can devise operational and control strategies to pre heat and precool spaces based on typical load patterns , weather and occupancy forecasts , and available distributed energy resources. To support the studies on a large scale, there are additional: 1) to generate heterogeneous building models with a limited number of real buildings, methods need to be developed to diversify the building parameters. In addition, the building systems, such as the ground sou rce condenser water loop, need to be re -sized or even re -designed to accommodate more buildings. 2) the numerical performance of the virtual testbed may need to be re -evalu ated to understand how the required computing resources vary by the scale of the stud ied community. Based on the evaluation results, advanced solvers can be selected to meet the requirement of the supporting studies. Acknowledgment This research was partially supported by the of the U.S. Department of Energy, Energy Efficiency and Renewabl e Energy, Office of Building Technologies, Emerging Technologies Program, under Contract no. DE - AC05 -76RL01830. This research was also supported by the National Science Foundation under Awards No. IIS -1802017. This work also emerged from the IBPSA Proj ect 1, an inter - nationally collaborative project conducted under the umbrella of the International Building Performance Simulation Association (IBPSA). Project 1 aims to develop and demonstrate a BIM/GIS and Modelica Framework for building and communit y energy system design and operation.