- Email: [email protected]
Fast and accurate district heating and cooling energy demand and load calculations using reduced-order modelling
Applied Energy 238 (2019) 963–971
Contents lists available at ScienceDirect
Applied Energy journal homepage: www.elsevier.com/locate/apenergy
Fast and accurate district heating and cooling energy demand and load calculations using reduced-order modelling
T
⁎
Eui-Jong Kima, , Xi Heb, Jean-Jacques Rouxc, Kévyn Johannesd, Frédéric Kuznikc a
Department of Architectural Engineering, Inha University, Incheon 22212, South Korea Graduate School of Civil Engineering, Tongji University, Shanghai, China c Université de Lyon, INSA–Lyon, CETHIL UMR5008, F–69621 Villeurbanne, France d Université de Lyon, Université Lyon 1, CETHIL UMR5008, F–69621 Villeurbanne, France b
H I GH L IG H T S
G R A P H I C A L A B S T R A C T
new fast and accurate district • Aheating and cooling loads calculation method is proposed.
reduction technics are • Mathematical applied to the set of differential equations representing physics.
solver Dymola is used to calculate • The the load profiles. 2 differential equations are re• Only quired to calculate the energy loads of 10 buildings.
A R T I C LE I N FO
A B S T R A C T
Keywords: Urban building model Reduced model Heating and cooling load calculation State-space model
Recent developments in building energy models for urban energy simulation are primarily based on bottom-up modelling (N models used for N buildings). This work aims to develop a single assembled model for multiple buildings for convenient use in detailed urban analysis. The proposed model exhibits state-space model formalism, and a state-size reduction technique is applied to maintain model accuracy, even for a low-order representation. To accelerate the calculation time and ensure numerical stability, a direct solver is proposed to eliminate the iterative calculations required in Dymola for annual load calculations. The results of the proposed reduced model are in good agreement with the reference model. For a test case of ten buildings, a 2nd order reduced model (i.e., 2 differential equations) with the proposed direct solver can predict accurately the dynamic energy behaviour, resulting in an error of about 0.43% for the annual loads.
1. Introduction Is it possible to predict the energy demand and loads of a district with less than 10 equations? As the number of people residing in cities is increasing, sustainable urban development has become an important area of interest for planners and engineers. In particular, efficient energy management is required for future cities because the energy
⁎
demand is increasing [1]. Furthermore, existing studies on future energy supply systems, such as decentralized power systems or district heating and cooling systems, emphasize the importance of energy management at the city scale [2]. This trend leads to requirements for developing detailed urban models for city-scale energy system planning and management. Recent studies reviewed various modelling approaches in this nascent field [3].
Corresponding author. E-mail address: [email protected] (E.-J. Kim).
https://doi.org/10.1016/j.apenergy.2019.01.183 Received 3 September 2018; Received in revised form 30 November 2018; Accepted 19 January 2019 0306-2619/ © 2019 Elsevier Ltd. All rights reserved.
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
X Y ∼ Y
Nomenclature Ao, Bo, Co state model matrices Ac, Bc, Cc modal-basis state model matrices A, B, C, D reduced model matrices LTI linear time-invariant MIMO multiple-input and multiple-output N number of buildings O model order P input energy gain PI proportional integral Ref reference RM reduced model S surface area (m2) SISO single-input and single-output SM simplified model SS state-space model T temperature vector T temperature (K) T time (s) U U-value (W/m2 K) U input vector V volume (m3)
basis-changed state vector output vector approximated output vector
Greek symbol
¯ Φ
average solar heat flux (W/m2)
Subscripts abs dir in k out set sky t trans uk wall win
absorbed solar flux direct solver indoor known values outdoor temperature set point temperature sky temperature total transmitted solar irradiance unknown values wall window
deduced from a comprehensive and detailed physical model [5]. In this model, a small number of selected building parameters was used to determine all required model terms through regression-based assumptions. The model represented specific building typologies, and it was difficult to generate a new grey-box model for a new building typology. CITYSIM [6] used a simple resistance and capacitance model (RC model [7]). Two capacitances were used: one for indoor air and the other for the envelope node. This simplification resulted in an error of about 6% in annual energy consumption and 5 °C in hourly indoor temperatures for a few wall types [8] compared to a detailed reference model. These models may be sufficient for evaluating annual energy demand; however, they demonstrate the necessity of a detailed model for hourly or sub-hourly information of building dynamics. In addition, commercial simulation tools have been used for urban simulation. Courchesne-Tardif et al. used DOE-2 [9] for detailed load calculations for buildings of a district [10]. Then, hourly load data were used to optimize community-scale energy systems. A similar detailed simulation was performed using EnergyPlus [11] to optimize a similar district-scale energy system [12]. While this approach is useful for a small number of buildings, fast and detailed calculations are required when optimizing urban energy systems connected to a large number of buildings. For all abovementioned models, the number of models was equal to the number of simulated buildings (N models used for N buildings). This approach implies that several numerical links are required between buildings and energy systems, such as district heating or cooling systems, which interact with all connected buildings. A highly simplified building model for a group of buildings may not be adequate for detailed simulations. The issue addressed in this work is the development of a fast and accurate model of the heating and cooling energy demand and load calculations using reduced-order modelling. This model corresponds to the lower-right domain of Fig. 1, which is a less developed domain. The scope of the work is presented in Section 2 to clarify the target area of the proposed work. The methodology developed in this work is based on the model order reduction technique described in Section 3. The results and the discussion are presented in Section 4, while the conclusions are given in Section 5.
Fig. 1 provides a schematic overview of the existing approaches and tools. Unlike previously developed urban building models, almost all newly developed models are based on bottom-up modelling. This approach implies that a thermal model represents a building and that a district- or city-level model can be constructed by combining the thermal models. As dozens to thousands of buildings are the targets of urban analysis, there can be various levels of detail in building models, and they are generally determined by the purpose and size of the target project. For instance, urban building energy simulation (UBES) shown in Fig. 1 belongs to detailed building simulation where one building is described with a large number of calculation nodes to predict hourly or sub-hourly loads. Often, these building models are combined with a detailed micro-climate model to evaluate effects of such interactions on resulting building annual energy loads. On the other hands, city energy simulation (CES) is run for a large number of buildings in a district or city. Therefore, various simplifications are required to calculate mainly the city energy demand expressed in kWh. The hourly city energy loads expressed in kW cannot be obtained with such a simplified approach. Several CES tools have been developed. One of the first CES tool was SUNtool [4], which was developed as part of a European research project. The building model was defined as a grey-box model that was
Fig. 1. Comparison of the domains of availability of different categories of energy simulation of urban buildings – extracted from [3]. 964
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
2. Scope of the work
3. Methods
Why is it required to predict heating and cooling demand and loads for a district? An answer is to provide data for decision-makers to decide various scenarios of engineering and energy solutions. For instance, new sizing of district heating or cooling networks is based on energy demand and loads. This work only focuses on district demand and loads calculations, and predicting urban energy or electricity consumption is not included in the scope of the work. In consequence, other energy uses pertaining to lighting, appliances, and HVAC operation are not discussed, neither energy systems efficiency. Recently similar works, which focus on the demand side, have been carried out. Kohler et al. developed the district energy model using a degree-day method [13], and [14] also used a similar method to evaluate the relationship between the city energy demand and building key coefficients. However, these models only account for annual heating and cooling demand (kWh). Hourly loads (kW) are useful for energy management and to test various energy conservation measures (ECM). For instance, effects of high insulated walls and windows on shaving peak loads can be evaluated with a load calculation tool. The tool is typically based on detailed modelling Such models are time-consuming as a number of buildings in the district and even all the districts in a city are simulated as shown in Fig. 2. Therefore, reduced-order modelling is proposed in this work to tackle the computational time issue. Fig. 2 shows an example of the proposed load calculation model and how it can be used. From GIS data, building models for a target district, ith district, can be defined (see also the following Section 3.1), and typical weather data or set-point temperatures are given to the models to evaluate the demand and loads. When a tested scenario is applied, e.g. harsh or future weather case, the resulting demand and loads can be calculated. Once a reduced-order model is developed, scenarios can be tested in a fast and accurate way. This calculation can be extended to test ECM scenarios. Since various optimal methods were proposed to accelerate computation time required for constructing reduced-order models [15], the ECM scenarios can also be rapidly tested. Leaving this application for future work, the current study is focusing on the modelling methodology.
3.1. Urban energy simulation framework In a previous work [16], a single-zone building model was developed for an individual building. A simplified model (SM), which we considered as a reference model (Ref) for testing purposes, was developed. The SM was suitable for model order reduction while maintaining very good accuracy. In the previous work [16], the SM model was tested with a detailed model that was defined using the EDF’s library BuildSysPro [17]. The SM gave similar results as the BuidSysPro models, showing errors of less than 1% in annual energy demand and less than 3% in hourly loads. The detailed BuildSysPro model was validated experimentally and by inter-comparison with other well-known simulation models [18]. Therefore, validation tests are reported using the SM. As the state model size reduction technique used in the previous work is applicable for only linear and time-invariant (LTI) systems, a few assumptions and modifications were made to define the SM as an LTI system. These assumptions are summarized as follows:
• Use of an equivalent wall for all walls • Use of an equivalent window for all windows • Use of constant longwave heat exchange coefficients • Externalization of the calculation of variable solar transmittance rates through windows
These assumptions are shown in Fig. 3. Based on the assumptions, we can use fewer wall component models, i.e., an equivalent wall, a window, and a floor. Orientation-dependent characteristics, such as absorbed or transmitted solar irradiance, are calculated using an external calculator. Then, mean variables are considered as inputs. There are six inputs in the SM, i.e., outdoor temperature (Tout), sky tem¯ trans ), mean absorbed perature (Tsky), mean transmitted irradiance (Φ ¯ abs_win ), mean absorbed irradiance for walls irradiance for windows (Φ ¯ abs_wall ), and heating or cooling load + internal gain (P). The resultant (Φ indoor air temperature under the given inputs is set as the output variable. The six inputs are used as boundaries; for instance, Tout is an outer boundary of walls and windows to calculate unsteady-state heat transfer rates through walls and windows. In addition, Tout is used to calculate ventilation and infiltration losses by an air exchange module. Air change rates need to be given as a model parameter and calculated by a subprogram. Then, knowing these values, the described technique can be used either for commercial or residential buildings.
Fig. 2. District loads and demand calculation for various test scenarios. 965
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
Fig. 3. Reference model (simplified linear and time-invariant model) for a single building.
as an LTI system to apply the state model size reduction technique. This technique is not well known among building simulation researchers; however, it has already been applied to common building simulation programs such as Comfie [20]. Briefly, the technique can be explained using the principle of the Fourier series, where the first dominant terms can sufficiently predict full-order dynamics. Among the reduction techniques, the state model size reduction technique is considered a mature and stable method in the field of modelling. However, its applications are limited to LTI systems, as mentioned above. The LTI system in the SS model formalism can be used to apply the technique without elaborative processes owing to existing well-developed functions. An advantage of Modelica is that the SS model can be exported easily from the user-friendly modelling interface (i.e., Dymola). In this work, the building model is defined as an LTI system, and model elements are defined using the open-source Buildings library [21] in Modelica. The SM developed in this work is similar to the SM [16] that was defined using BuildSysPro [17]. Fig. 5a shows an example of an assembly of ten building models in Dymola, and Fig. 5b shows an SM for an individual building. Global inputs and outputs are connected to individual building models through Modelica connectors. Local inputs and outputs are the same as those for the SM (see Fig. 3). The power inputs (P) that represent the energy gain
While this is a simple single-zone case, similar SMs can be defined for multi-zone cases, such as a block building or a roof including a building. These various SMs, referred to as building typology models, can be prepared using a simulation platform such as Dymola [19], which is a Modelica interface. As shown in Fig. 4, N buildings in a district are allocated to predefined typologies. Using a GIS database, this typology allocation and the instantiation of building SMs by setting parameters can be automated. The final N SMs can be used to deduce N reduced models (RMs). This work investigates the possibility of assembling these N SMs to develop a single state-space (SS) model and then obtain a single RM model for N buildings (see the bottom right part of Fig. 4). This approach may have advantages over an N RM methodology. A single SS model for N buildings with a considerably low-order representation can drastically decrease the computational resources and cost required for calculation. In addition, it is easier to export a single model for simulations using other computational tools. Finally, this low-order single model representation for a group of buildings is useful for model-based urban system management. 3.2. Linear time-invariant systems to apply a single state-space model As explained in the previous section, each building model is defined
Fig. 4. Urban energy simulation framework. 966
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
Fig. 5. (a) Assembly of ten building models in Dymola to generate a single SS model and (b) the inner structure of an SM.
are provided to all SMs. Thus, indoor temperatures (output) can be controlled by adjusting this value, which is referred to as heating or cooling load if the indoor temperature is equal to the set point temperature. Based on this assembly, a SS model for ten building models can be obtained using a Modelica linearization function.
⎧ Ẋ (t ) = AX (t ) + BU (t ) ⎨∼ ⎩ Y (t ) = CX (t ) + DU (t )
All state-space matrices (Ao, Bo, and Co) in Eq. (1) for buildings are easily obtained using a Modelica linearization function. As the outputs are indoor temperatures, which are part of state vector T, input vector U does not contribute to output values Y. Thus, only Co is used in Eq. (1). In the balanced truncation method, reduction can be achieved by eliminating the state variables that are less affected by changes in input (controllability) or have less influence on the variation in output (observability). As each state variable does not necessarily have the same degree of controllability and observability, the initial system (Ao, Bo, and Co) is transformed through basis change (T.X) so that controllability and observability are balanced. After the state coordinate transformation is complete, the system consists of Ac, Bc, and Cc (Eq. (2)). Through the basis change, the initial sparse matrices transform into full-term matrices because a full-term state-transforming matrix is applied between the initial and new state vectors [22]. As shown in Fig. 6a, the initial matrix Ao is sparse. Each subgroup relates to a building model because the buildings are not physically connected. The basis change is carried out by a reverse multiplication of the full-term state-transforming matrix, and thus, the new matrix Ac becomes fullterm, as shown in Fig. 6b. Consequently, subgroup independence is lost for all matrices (Ac, Bc, and Cc), and the energy gain input (P) for one building can affect the thermal dynamics of other buildings. As the terms of the state vector (X) in the basis-transformed state-
3.3. Balanced truncation method for reduction Among LTI system-based reduction techniques, Moore’s method [22], which is referred to as the balanced truncation method, is regarded as one of the most useful methods [23]. In this work, we use this method, and its primary steps are given in Eqs. (1)–(3). For buildings, the system of differential equations governing the calculation of heating and cooling loads (see Eq. (6)) can be deduced from Eq. (3). The reduction process starts with Eq. (1), which describes the typical temporal evolution of indoor air temperatures as shown in Fig. 3. Y(t) includes several indoor temperatures of buildings. A detailed method to replace this Y(t) with a new output of heating and cooling loads is given in Section 3.4.
⎧ Ṫ (t ) = A OT (t ) + BO U (t ) ⎨ ⎩ Y (t ) = CO T (t )
(1)
Considering the basis change:
⎧ Ẋ (t ) = A C X (t ) + BC U (t ) ⎨∼ ⎩ Y (t ) = CC X (t )
(3)
(2)
After the reduction method:
Fig. 6. Matrix transformation (Ao → Ac) through a basis change. 967
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
requires a greater number of iterations to find a set of P values that cause all the Tin values to reach the desired set Tset. Thus, the process is time consuming and may cause numerical instabilities such as convergence errors. To prevent these problems, we propose a new solver that can calculate the loads directly from the reduced model without the PI controller and consequently without iterations by substituting P with Tin in Eq. (3). The following equations (Eqs. (4)–(6)) explain these modifications. First, the input variables are gathered into known variables Uk and unknown variables Uuk, which represents heating or cooling loads.
Table 1 Parameters used for the case of 10 buildings. Parameter
Building Number
3
Vin (m ) S (m2) Windows (m2) Swin/St Uwall (W/m2 K)
Vin, S, Windows, Swin/St Uwall (W/m2 K)
N° 1
N° 2
N° 3
N° 4
N° 5
180 72 45 0.38 0.27
720 288 80 0.33 0.27
1125 450 150 0.50 0.27
1800 720 180 0.30 0.27
5625 2250 420 0.28 0.27
N° 6
N° 7
N° 8
N° 9
N° 10
″ 3.05
″ 3.05
″ 3.05
″ 3.05
″ 3.05
⎧ Ẋ (t ) = AX (t ) + B k Uk (t ) + B u Uuk (t ) ⎨∼ ⎩ Y (t ) = CX (t ) + D k Uk (t ) + Du Uuk (t )
space model have the same controllability and observability, the dominant terms whose controllability and observability are higher can be selected to establish a final reduced model. The non-selected states, referred to as truncated states, can be added to the output model to ensure static gain conservation of the system. Therefore, the final reduced model, expressed as Eq. (3), includes matrix D. The number of dominant states is referred to as model order in the reduction technique. The balred MATLAB function is used in this study for the reduction process. Previous works [16] proposed a 6-order (O (6)) RM for hourly building simulations.
(4)
The heating or cooling load variables can be rewritten using the second part of Eq. (4).
∼ Uuk (t ) = D−u1 (Y (t ) − CX (t ) − D k Uk (t ))
(5)
This equation can substitute Uuk in the first part of Eq. (4). After rearranging the equation to maintain the SS model formalism, the final direct SS model is obtained as follows:
⎧ Ẋ (t ) = (A − B u D−u1 C)X (t ) + (B k − B u D−u1 D k |B u D−u1 ) ⎪ B⏟ A⏟ dir dir ⎪ ⎨ −1 −1 −1 Uk (t ) ⎪ Uuk (t ) = −Du CX (t ) + (−Du D k |Du ) ∼Y (t ) C⏟ D⏟ ⎪ Y⏟ dir dir dir U⏟ dir ⎩
3.4. Direct load calculation solver
{ } Uk (t ) ∼ Y (t ) U⏟ dir
{ }
In Modelica, a non-causal tool, the load calculation is performed iteratively using a proportional integral (PI) controller [24]. In our initial model, before the basis change, the indoor air temperature Tin of one building is taken as the observed variable to control its input P. A carefully controlled value of P that matches Tin with Tset represents the heating or cooling load. Here, none of the values of P set for other building models are utilized in this iterative process; therefore, the initial system (Eq. (1)) is a single-input and single-output (SISO) problem. Although ten PI controllers are required in the case of ten buildings, each controller operates independently. However, the reduced model is covariant to any of the inputs since all the states are connected with each other after introducing a statetransforming matrix for the basis change. This system is called a multiple-input and multiple-output (MIMO) system. Typically, this system
(6)
From this equation, the heating or cooling loads Uuk can be directly ∼ calculated by introducing Tset as the new input term Y(t ) of Eq. (6). Using this formula, hourly heating or cooling loads are obtained under the condition that Tin is equal to Tset. To take into account scheduling such as building occupancy hours, particular attention needs to be paid. A solution is that Eq. (4) is used for non-occupancy periods by giving Uuk = 0. When a building starts to be operated, Eq. (6) is newly run. In this case, X(t) obtained using Eq. (4) must be used as the state vector in the precedent time-step for calculating Ẋ (t ) .
Fig. 7. Annual hourly load comparison of ten Ref models with (a) ten RMs (one by one) and (b) a single RM (assembled). 968
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
4. Results and discussion
model apparently affects the model accuracy. As fewer model states are used, the peak and mean errors increase. The 2-state model for ten buildings shows the highest error among the models; however, this high error is temporary and is only 0.43% for the annual load. This level of error in annual energy demand is quite accurate since the CITYSIM model [8] showed an error of 6% compared to ESP-r reference models [10]. Furthermore, the results of the 2-state RM for the peak heating and cooling periods show that the dynamics and magnitude are similar to those of the Ref model, as shown in Fig. 10. An option is to use the 10state RM. A reasonable magnitude of error is maintained with the 10state RM even for cases with a large number of buildings, e.g., up to a hundred thousand buildings [26].
4.1. Simulation conditions He [25] carried out a detailed comparison test between Buildings and BuildSysPro SM models and found that they were in very good agreement. As mentioned earlier, this outcome was expected because both models use the same physical assumptions. Therefore, the SMs developed in this work using Buildings are used as the Ref models for the following simulation results. Cases of multiple buildings were tested using the parameters presented in Table 1. Different parameters were used for comparison. The internal gains were not considered in this study; however, they can be easily taken into account by including them with the P values in the case of a PI controller or by adding it to Uk (Eq. (6)) in the case of a direct solver. To test our final single RM (assembled) using the proposed direct solver, a separate simulation with an RM or Ref model for each building was performed multiple times using a PI controller (one by one), and the results were summed. In addition, the performance of the state-size reduction technique was tested using the Ref models.
4.4. CPU time comparison The CPU time required for the simulation of different RMs was measured. A computer equipped with an Intel Core i7-4770 3.4 GHz CPU and 16 GB RAM was used as the benchmark for simulating the models in Dymola 2016 with the DASSL solver. Table 2 shows the results of the CPU time for hourly annual simulations of the ten buildings. The assembled RM model without the proposed direct solver is slower than the one-by-one RMs even if they are of similar order. As the single RM is a MIMO system, it requires a greater number of iterations for convergence. For the case of the direct solver, the assembled RM with 50 states requires approximately 4 s. To further investigate the CPU time required by the model with the proposed solver, the assembled RMs with different numbers of states were tested. Fig. 11 shows the CPU time results of the annual hourly simulations. The CPU time of the assembled RM model and its order are linearly proportional. A 4-order model requires less than 1 s for an annual hourly simulation of 10 buildings. This finding implies that a large number of buildings in urban areas can be simulated rapidly when a very-low-order RM with the proposed direct solver is used.
4.2. Validation of the reduction technique The model reduction performance was measured using the Ref model (SM model). For ten buildings, ten Ref models were prepared. Each Ref model was run separately to calculate the annual hourly load. Similarly, the ten RMs were simulated separately. Fig. 7a shows the differences in annual loads between the models (one-by-one Ref vs. one-by-one RM). For the case of ten RMs, each RM model was on the order of 6 states; thus, 60 states were used for 10 buildings. In this figure, the load ranges between +300 kW (for heating) and −300 kW (for cooling). The differences are within 2 kW, which shows good agreement between the Ref models and RMs. Then, the assembled RM model that was obtained by assembling and reducing the ten building models was compared with the Ref model. Fig. 7b shows a comparison between the annual loads of the assembled RM model and the Ref model (one by one). This single representation model (RM assembled) was our target test model. As shown in the figure, this assembled RM model shows similar accuracy to the previous one-by-one RMs. This finding implies the possibility of using a single lower-order SS model to accurately model a number of buildings. However, the model could be tested to only 50 orders. Further investigations to decrease the order of the model could not be performed because of the occurrence of a convergence problem in the MIMO system, which was mentioned previously. Thus, the current iterative method and the corresponding solver are not the best techniques for this assembled RM model.
5. Conclusions In this work, we proposed a single reduced model for hourly load calculation for multiple buildings. According to an urban energy analysis framework, the model can be deduced from the detailed parameter data of buildings. Unlike conventional approaches that require N models for N buildings, the proposed method requires a single
4.3. Single reduced model with a direct solver The proposed direct solver for the single RM was tested. The new direct representation model (Eq. (6)) was also run in Dymola to compare the CPU times. Fig. 8 shows the results of the comparison between the Ref model and the single RM with the proposed solver (RM direct (assembled)). The RM model in this case is on the order of 50 states. This new solver shows good agreement with the Ref model; however, infrequent peak errors are observed. The peak error is approximately 20 kW; however, the other periods exhibit dynamics that are very similar to those of the Ref model. When this solver is used, the convergence problems that were caused by iterative methods do not occur. Therefore, models of lower orders can be tested. Fig. 9 shows the different error levels according to the RM order. The tests were carried out down to the 2-state model. As expected, the proposed solver could run without any numerical problem even with a model of very low order. However, the order of the
Fig. 8. Performance of the direct solver for the RM assembled case. 969
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
Fig. 11. CPU time of RM direct-assembled models with different numbers of states. Fig. 9. Errors of different model orders in the RM direct (assembled).
For the test case, which included ten buildings, this assembled RM model gave similar results with the one-by-one RMs as well as with the SMs; however, the total model order could not be reduced because of convergence issues. In addition, the CPU time was increased, as the assembled RM requires more iterations in Dymola. Therefore, a direct solver was developed to calculate the loads without iteration, which enables reducing the model order further. The final assembled RM with the proposed direct solver can reduce the order down to only 2 states, with an error of 0.43% in the annual load. This means that two equations can predict the thermal dynamics of ten buildings or more. This low-order single representation of buildings could be useful for rapid computation and simulation interoperability for further detailed urban analysis. Various applications that utilize the proposed model are expected in the near future. Acknowledgment This work was supported under the framework of international cooperation program managed by the College of Engineering, Inha University. In addition, this work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIP) (No. 2016R1C1B2011097).
Fig. 10. Two-state RM direct for peak heating and cooling periods.
References Table 2 CPU time for different RMs in Dymola.
RM one by one (60 states) RM assembled (50 states) RM direct assembled (50 states)
CPU time (s)
Sys.
81 499 4
SISO MIMO Direct
[1] An J, Yan D, Hong T, Sun K. A novel stochastic modeling method to simulate cooling loads in residential districts. Appl Energy 2017;206:134–49. [2] Ma W, Fang S, Liu G, Zhou R. Modeling of district load forecasting for distributed energy system. Appl Energy 2017;204:181–205. [3] Frayssinet L, Merlier L, Kuznik F, Hubert JL, Milliez M, Roux JJ. Modeling the heating and cooling energy demand of urban buildings at city level. Renew Sustain Energy Rev 2018;81:2318–27. [4] Robinson D, Stankovic S, Morel N, Dequé F, Rylatt M, Kabele K, et al. Integrated resource flow modelling of urban neighbourhoods: project suntool. Proceedings of building simulation 2003, Einthoven, The Netherlands. 2003. [5] Déqué F, Ollivier F, Poblador A. Grey boxes used to represent buildings with a minimum number of geometric and thermal parameters. Energy Build 2000;31(1):29–35. [6] Robinson D, Haldi F, Kämpf J, Leroux P, Perez D, Rasheed A, Wilke U. CITYSIM: comprehensive micro-simulation of resource flows for sustainable urban planning. Proceedings of Building Simulation 2009, Glasgow, Scotland. 2009. [7] Roux JJ. Proposition de modèles simplifiés pour l'étude du comportement thermique des batiments Phd thesis Lyon (France): Université de Lyon; 1984 [8] Kämpf J, Robinson D. A simplified thermal model to support analysis of urban resource flows. Energy Build 2007;39:445–53. [9] DOE-2 engineers manual, version 2.1A. Technical Information Center, United States Department of Energy; 1982. [10] Courchesne-Tardif A, Kummert M, Demark S, Butler T, Pearl D, Jones S, et al.
assembled model for multiple buildings. Such method is mandatory to predict the heating and cooling loads, in kW, and then for city energy management. This is clearly a progress in comparison with common CES focusing on annual energy demand, in kWh. Similar to previous studies, the state model size reduction technique was applied to this model. The balanced truncation method is effective in reducing the computational time and ensuring the accuracy of the model. In a previous work, this reduction technique was applied to each SM model that represents a single building. In this work, the reduction technique is applied to a single state-space model for the buildings, constructed by assembling N SMs in Dymola. 970
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
[11] [12]
[13]
[14]
[15]
[16] [17]
[18] Bontemps S, Kaemmerlen A, Blatman G, Mora L. Reliability of dynamic simulation models for building energy in the context of low-energy buildings. Proceedings of building simulation 2013, Chambéry, France. 2013. [19] Dymola version 2013, (10.01.2013); 2013. Homepage: http://www.3ds.com/ products/catia/portfolio/dymola. [20] Peuportier B, Blanc Sommereux I. Simulation tool with its expert interface for the thermal design of multizone buildings. Int J Sol Energy 1990;8:109–20. [21] Wetter M, Zuo W, Nouidui TS. Modeling of heat transfer in rooms in the Modelica buildings library. 12th Conference of international building performance simulation association, Sydney, Australia. 2011. [22] Moore BC. Principal component analysis in linear system: controllability, observability, and model reduction. IEEE Trans Autom Control 1981;AC-26(1):17–32. [23] Palomo E, Bonnefous Y, Déqué F. Guidance for the selection of a reduction technique for thermal model. Proceedings of BS, Prague, Czech Republic. 1997. [24] Ghiaus C. Causality issue in the heat balance method for calculating the design heating and cooling load. Energy 2013;50:292–301. [25] He X. Simulation énergétique à l’échelle quartier: réduction des modèles de bâtiment sous Dymola Master thesis Villeurbanne (France): ENTPE; 2014 [26] Garcia Perez M, Kim E-J, Roux J-J. Reduction techniques for urban energy simulation. 2017 SAREK summer annual conference, Seoul, Republic of Korea. 2017. p. 587–90.
Assessing community-scale energy supply scenarios using TRNSYS simulations. Proceedings of building simulation 2011, Australia, Sydney. 2011. EnergyPlus Engineering Reference. The reference to EnergyPlus calculations. Ernest Orlando Lawrence Berkeley National Laboratory, US Department of Energy; 2010. Huber J, Nytsch-Geusen C. Development of modeling and simulation strategies for large scale urban districts. Proceedings of building simulation 2011, Sydney, Australia. 2011. Kohler M, Blond N, Clappier A. A city scale degree-day method to assess building space heating energy demands in Strasbourg Eurometropolis (France). Appl Energy 2016;184:40–54. Gianniou P, Reinhart C, Hsu D, Heller A, Rode C. Estimation of temperature setpoints and heat transfer coefficients among residential buildings in Denmark based on smart meter data. Build Environ 2018;139:125–33. Martin GP, Kim E-J, Roux J-J. Simulation des besoins énergétiques à l’échelle de la ville: Une méthodologie de modélisation adaptée. Proceedings of IBPSA France Conference, Bordeaux, France. 2018. Kim EJ, Plessis G, Huber JL, Roux JJ. Urban energy simulation: simplification and reduction of building envelope models. Energy Build 2014;84:193–202. Plessis G, Kaemmerlen A, Lindsay A. BuildSysPro: a Modelica library for modelling buildings and energy systems. Proceedings of the Modelica conference 2014, Lund, Sweden. 2014.
971
Contents lists available at ScienceDirect
Applied Energy journal homepage: www.elsevier.com/locate/apenergy
Fast and accurate district heating and cooling energy demand and load calculations using reduced-order modelling
T
⁎
Eui-Jong Kima, , Xi Heb, Jean-Jacques Rouxc, Kévyn Johannesd, Frédéric Kuznikc a
Department of Architectural Engineering, Inha University, Incheon 22212, South Korea Graduate School of Civil Engineering, Tongji University, Shanghai, China c Université de Lyon, INSA–Lyon, CETHIL UMR5008, F–69621 Villeurbanne, France d Université de Lyon, Université Lyon 1, CETHIL UMR5008, F–69621 Villeurbanne, France b
H I GH L IG H T S
G R A P H I C A L A B S T R A C T
new fast and accurate district • Aheating and cooling loads calculation method is proposed.
reduction technics are • Mathematical applied to the set of differential equations representing physics.
solver Dymola is used to calculate • The the load profiles. 2 differential equations are re• Only quired to calculate the energy loads of 10 buildings.
A R T I C LE I N FO
A B S T R A C T
Keywords: Urban building model Reduced model Heating and cooling load calculation State-space model
Recent developments in building energy models for urban energy simulation are primarily based on bottom-up modelling (N models used for N buildings). This work aims to develop a single assembled model for multiple buildings for convenient use in detailed urban analysis. The proposed model exhibits state-space model formalism, and a state-size reduction technique is applied to maintain model accuracy, even for a low-order representation. To accelerate the calculation time and ensure numerical stability, a direct solver is proposed to eliminate the iterative calculations required in Dymola for annual load calculations. The results of the proposed reduced model are in good agreement with the reference model. For a test case of ten buildings, a 2nd order reduced model (i.e., 2 differential equations) with the proposed direct solver can predict accurately the dynamic energy behaviour, resulting in an error of about 0.43% for the annual loads.
1. Introduction Is it possible to predict the energy demand and loads of a district with less than 10 equations? As the number of people residing in cities is increasing, sustainable urban development has become an important area of interest for planners and engineers. In particular, efficient energy management is required for future cities because the energy
⁎
demand is increasing [1]. Furthermore, existing studies on future energy supply systems, such as decentralized power systems or district heating and cooling systems, emphasize the importance of energy management at the city scale [2]. This trend leads to requirements for developing detailed urban models for city-scale energy system planning and management. Recent studies reviewed various modelling approaches in this nascent field [3].
Corresponding author. E-mail address: [email protected] (E.-J. Kim).
https://doi.org/10.1016/j.apenergy.2019.01.183 Received 3 September 2018; Received in revised form 30 November 2018; Accepted 19 January 2019 0306-2619/ © 2019 Elsevier Ltd. All rights reserved.
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
X Y ∼ Y
Nomenclature Ao, Bo, Co state model matrices Ac, Bc, Cc modal-basis state model matrices A, B, C, D reduced model matrices LTI linear time-invariant MIMO multiple-input and multiple-output N number of buildings O model order P input energy gain PI proportional integral Ref reference RM reduced model S surface area (m2) SISO single-input and single-output SM simplified model SS state-space model T temperature vector T temperature (K) T time (s) U U-value (W/m2 K) U input vector V volume (m3)
basis-changed state vector output vector approximated output vector
Greek symbol
¯ Φ
average solar heat flux (W/m2)
Subscripts abs dir in k out set sky t trans uk wall win
absorbed solar flux direct solver indoor known values outdoor temperature set point temperature sky temperature total transmitted solar irradiance unknown values wall window
deduced from a comprehensive and detailed physical model [5]. In this model, a small number of selected building parameters was used to determine all required model terms through regression-based assumptions. The model represented specific building typologies, and it was difficult to generate a new grey-box model for a new building typology. CITYSIM [6] used a simple resistance and capacitance model (RC model [7]). Two capacitances were used: one for indoor air and the other for the envelope node. This simplification resulted in an error of about 6% in annual energy consumption and 5 °C in hourly indoor temperatures for a few wall types [8] compared to a detailed reference model. These models may be sufficient for evaluating annual energy demand; however, they demonstrate the necessity of a detailed model for hourly or sub-hourly information of building dynamics. In addition, commercial simulation tools have been used for urban simulation. Courchesne-Tardif et al. used DOE-2 [9] for detailed load calculations for buildings of a district [10]. Then, hourly load data were used to optimize community-scale energy systems. A similar detailed simulation was performed using EnergyPlus [11] to optimize a similar district-scale energy system [12]. While this approach is useful for a small number of buildings, fast and detailed calculations are required when optimizing urban energy systems connected to a large number of buildings. For all abovementioned models, the number of models was equal to the number of simulated buildings (N models used for N buildings). This approach implies that several numerical links are required between buildings and energy systems, such as district heating or cooling systems, which interact with all connected buildings. A highly simplified building model for a group of buildings may not be adequate for detailed simulations. The issue addressed in this work is the development of a fast and accurate model of the heating and cooling energy demand and load calculations using reduced-order modelling. This model corresponds to the lower-right domain of Fig. 1, which is a less developed domain. The scope of the work is presented in Section 2 to clarify the target area of the proposed work. The methodology developed in this work is based on the model order reduction technique described in Section 3. The results and the discussion are presented in Section 4, while the conclusions are given in Section 5.
Fig. 1 provides a schematic overview of the existing approaches and tools. Unlike previously developed urban building models, almost all newly developed models are based on bottom-up modelling. This approach implies that a thermal model represents a building and that a district- or city-level model can be constructed by combining the thermal models. As dozens to thousands of buildings are the targets of urban analysis, there can be various levels of detail in building models, and they are generally determined by the purpose and size of the target project. For instance, urban building energy simulation (UBES) shown in Fig. 1 belongs to detailed building simulation where one building is described with a large number of calculation nodes to predict hourly or sub-hourly loads. Often, these building models are combined with a detailed micro-climate model to evaluate effects of such interactions on resulting building annual energy loads. On the other hands, city energy simulation (CES) is run for a large number of buildings in a district or city. Therefore, various simplifications are required to calculate mainly the city energy demand expressed in kWh. The hourly city energy loads expressed in kW cannot be obtained with such a simplified approach. Several CES tools have been developed. One of the first CES tool was SUNtool [4], which was developed as part of a European research project. The building model was defined as a grey-box model that was
Fig. 1. Comparison of the domains of availability of different categories of energy simulation of urban buildings – extracted from [3]. 964
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
2. Scope of the work
3. Methods
Why is it required to predict heating and cooling demand and loads for a district? An answer is to provide data for decision-makers to decide various scenarios of engineering and energy solutions. For instance, new sizing of district heating or cooling networks is based on energy demand and loads. This work only focuses on district demand and loads calculations, and predicting urban energy or electricity consumption is not included in the scope of the work. In consequence, other energy uses pertaining to lighting, appliances, and HVAC operation are not discussed, neither energy systems efficiency. Recently similar works, which focus on the demand side, have been carried out. Kohler et al. developed the district energy model using a degree-day method [13], and [14] also used a similar method to evaluate the relationship between the city energy demand and building key coefficients. However, these models only account for annual heating and cooling demand (kWh). Hourly loads (kW) are useful for energy management and to test various energy conservation measures (ECM). For instance, effects of high insulated walls and windows on shaving peak loads can be evaluated with a load calculation tool. The tool is typically based on detailed modelling Such models are time-consuming as a number of buildings in the district and even all the districts in a city are simulated as shown in Fig. 2. Therefore, reduced-order modelling is proposed in this work to tackle the computational time issue. Fig. 2 shows an example of the proposed load calculation model and how it can be used. From GIS data, building models for a target district, ith district, can be defined (see also the following Section 3.1), and typical weather data or set-point temperatures are given to the models to evaluate the demand and loads. When a tested scenario is applied, e.g. harsh or future weather case, the resulting demand and loads can be calculated. Once a reduced-order model is developed, scenarios can be tested in a fast and accurate way. This calculation can be extended to test ECM scenarios. Since various optimal methods were proposed to accelerate computation time required for constructing reduced-order models [15], the ECM scenarios can also be rapidly tested. Leaving this application for future work, the current study is focusing on the modelling methodology.
3.1. Urban energy simulation framework In a previous work [16], a single-zone building model was developed for an individual building. A simplified model (SM), which we considered as a reference model (Ref) for testing purposes, was developed. The SM was suitable for model order reduction while maintaining very good accuracy. In the previous work [16], the SM model was tested with a detailed model that was defined using the EDF’s library BuildSysPro [17]. The SM gave similar results as the BuidSysPro models, showing errors of less than 1% in annual energy demand and less than 3% in hourly loads. The detailed BuildSysPro model was validated experimentally and by inter-comparison with other well-known simulation models [18]. Therefore, validation tests are reported using the SM. As the state model size reduction technique used in the previous work is applicable for only linear and time-invariant (LTI) systems, a few assumptions and modifications were made to define the SM as an LTI system. These assumptions are summarized as follows:
• Use of an equivalent wall for all walls • Use of an equivalent window for all windows • Use of constant longwave heat exchange coefficients • Externalization of the calculation of variable solar transmittance rates through windows
These assumptions are shown in Fig. 3. Based on the assumptions, we can use fewer wall component models, i.e., an equivalent wall, a window, and a floor. Orientation-dependent characteristics, such as absorbed or transmitted solar irradiance, are calculated using an external calculator. Then, mean variables are considered as inputs. There are six inputs in the SM, i.e., outdoor temperature (Tout), sky tem¯ trans ), mean absorbed perature (Tsky), mean transmitted irradiance (Φ ¯ abs_win ), mean absorbed irradiance for walls irradiance for windows (Φ ¯ abs_wall ), and heating or cooling load + internal gain (P). The resultant (Φ indoor air temperature under the given inputs is set as the output variable. The six inputs are used as boundaries; for instance, Tout is an outer boundary of walls and windows to calculate unsteady-state heat transfer rates through walls and windows. In addition, Tout is used to calculate ventilation and infiltration losses by an air exchange module. Air change rates need to be given as a model parameter and calculated by a subprogram. Then, knowing these values, the described technique can be used either for commercial or residential buildings.
Fig. 2. District loads and demand calculation for various test scenarios. 965
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
Fig. 3. Reference model (simplified linear and time-invariant model) for a single building.
as an LTI system to apply the state model size reduction technique. This technique is not well known among building simulation researchers; however, it has already been applied to common building simulation programs such as Comfie [20]. Briefly, the technique can be explained using the principle of the Fourier series, where the first dominant terms can sufficiently predict full-order dynamics. Among the reduction techniques, the state model size reduction technique is considered a mature and stable method in the field of modelling. However, its applications are limited to LTI systems, as mentioned above. The LTI system in the SS model formalism can be used to apply the technique without elaborative processes owing to existing well-developed functions. An advantage of Modelica is that the SS model can be exported easily from the user-friendly modelling interface (i.e., Dymola). In this work, the building model is defined as an LTI system, and model elements are defined using the open-source Buildings library [21] in Modelica. The SM developed in this work is similar to the SM [16] that was defined using BuildSysPro [17]. Fig. 5a shows an example of an assembly of ten building models in Dymola, and Fig. 5b shows an SM for an individual building. Global inputs and outputs are connected to individual building models through Modelica connectors. Local inputs and outputs are the same as those for the SM (see Fig. 3). The power inputs (P) that represent the energy gain
While this is a simple single-zone case, similar SMs can be defined for multi-zone cases, such as a block building or a roof including a building. These various SMs, referred to as building typology models, can be prepared using a simulation platform such as Dymola [19], which is a Modelica interface. As shown in Fig. 4, N buildings in a district are allocated to predefined typologies. Using a GIS database, this typology allocation and the instantiation of building SMs by setting parameters can be automated. The final N SMs can be used to deduce N reduced models (RMs). This work investigates the possibility of assembling these N SMs to develop a single state-space (SS) model and then obtain a single RM model for N buildings (see the bottom right part of Fig. 4). This approach may have advantages over an N RM methodology. A single SS model for N buildings with a considerably low-order representation can drastically decrease the computational resources and cost required for calculation. In addition, it is easier to export a single model for simulations using other computational tools. Finally, this low-order single model representation for a group of buildings is useful for model-based urban system management. 3.2. Linear time-invariant systems to apply a single state-space model As explained in the previous section, each building model is defined
Fig. 4. Urban energy simulation framework. 966
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
Fig. 5. (a) Assembly of ten building models in Dymola to generate a single SS model and (b) the inner structure of an SM.
are provided to all SMs. Thus, indoor temperatures (output) can be controlled by adjusting this value, which is referred to as heating or cooling load if the indoor temperature is equal to the set point temperature. Based on this assembly, a SS model for ten building models can be obtained using a Modelica linearization function.
⎧ Ẋ (t ) = AX (t ) + BU (t ) ⎨∼ ⎩ Y (t ) = CX (t ) + DU (t )
All state-space matrices (Ao, Bo, and Co) in Eq. (1) for buildings are easily obtained using a Modelica linearization function. As the outputs are indoor temperatures, which are part of state vector T, input vector U does not contribute to output values Y. Thus, only Co is used in Eq. (1). In the balanced truncation method, reduction can be achieved by eliminating the state variables that are less affected by changes in input (controllability) or have less influence on the variation in output (observability). As each state variable does not necessarily have the same degree of controllability and observability, the initial system (Ao, Bo, and Co) is transformed through basis change (T.X) so that controllability and observability are balanced. After the state coordinate transformation is complete, the system consists of Ac, Bc, and Cc (Eq. (2)). Through the basis change, the initial sparse matrices transform into full-term matrices because a full-term state-transforming matrix is applied between the initial and new state vectors [22]. As shown in Fig. 6a, the initial matrix Ao is sparse. Each subgroup relates to a building model because the buildings are not physically connected. The basis change is carried out by a reverse multiplication of the full-term state-transforming matrix, and thus, the new matrix Ac becomes fullterm, as shown in Fig. 6b. Consequently, subgroup independence is lost for all matrices (Ac, Bc, and Cc), and the energy gain input (P) for one building can affect the thermal dynamics of other buildings. As the terms of the state vector (X) in the basis-transformed state-
3.3. Balanced truncation method for reduction Among LTI system-based reduction techniques, Moore’s method [22], which is referred to as the balanced truncation method, is regarded as one of the most useful methods [23]. In this work, we use this method, and its primary steps are given in Eqs. (1)–(3). For buildings, the system of differential equations governing the calculation of heating and cooling loads (see Eq. (6)) can be deduced from Eq. (3). The reduction process starts with Eq. (1), which describes the typical temporal evolution of indoor air temperatures as shown in Fig. 3. Y(t) includes several indoor temperatures of buildings. A detailed method to replace this Y(t) with a new output of heating and cooling loads is given in Section 3.4.
⎧ Ṫ (t ) = A OT (t ) + BO U (t ) ⎨ ⎩ Y (t ) = CO T (t )
(1)
Considering the basis change:
⎧ Ẋ (t ) = A C X (t ) + BC U (t ) ⎨∼ ⎩ Y (t ) = CC X (t )
(3)
(2)
After the reduction method:
Fig. 6. Matrix transformation (Ao → Ac) through a basis change. 967
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
requires a greater number of iterations to find a set of P values that cause all the Tin values to reach the desired set Tset. Thus, the process is time consuming and may cause numerical instabilities such as convergence errors. To prevent these problems, we propose a new solver that can calculate the loads directly from the reduced model without the PI controller and consequently without iterations by substituting P with Tin in Eq. (3). The following equations (Eqs. (4)–(6)) explain these modifications. First, the input variables are gathered into known variables Uk and unknown variables Uuk, which represents heating or cooling loads.
Table 1 Parameters used for the case of 10 buildings. Parameter
Building Number
3
Vin (m ) S (m2) Windows (m2) Swin/St Uwall (W/m2 K)
Vin, S, Windows, Swin/St Uwall (W/m2 K)
N° 1
N° 2
N° 3
N° 4
N° 5
180 72 45 0.38 0.27
720 288 80 0.33 0.27
1125 450 150 0.50 0.27
1800 720 180 0.30 0.27
5625 2250 420 0.28 0.27
N° 6
N° 7
N° 8
N° 9
N° 10
″ 3.05
″ 3.05
″ 3.05
″ 3.05
″ 3.05
⎧ Ẋ (t ) = AX (t ) + B k Uk (t ) + B u Uuk (t ) ⎨∼ ⎩ Y (t ) = CX (t ) + D k Uk (t ) + Du Uuk (t )
space model have the same controllability and observability, the dominant terms whose controllability and observability are higher can be selected to establish a final reduced model. The non-selected states, referred to as truncated states, can be added to the output model to ensure static gain conservation of the system. Therefore, the final reduced model, expressed as Eq. (3), includes matrix D. The number of dominant states is referred to as model order in the reduction technique. The balred MATLAB function is used in this study for the reduction process. Previous works [16] proposed a 6-order (O (6)) RM for hourly building simulations.
(4)
The heating or cooling load variables can be rewritten using the second part of Eq. (4).
∼ Uuk (t ) = D−u1 (Y (t ) − CX (t ) − D k Uk (t ))
(5)
This equation can substitute Uuk in the first part of Eq. (4). After rearranging the equation to maintain the SS model formalism, the final direct SS model is obtained as follows:
⎧ Ẋ (t ) = (A − B u D−u1 C)X (t ) + (B k − B u D−u1 D k |B u D−u1 ) ⎪ B⏟ A⏟ dir dir ⎪ ⎨ −1 −1 −1 Uk (t ) ⎪ Uuk (t ) = −Du CX (t ) + (−Du D k |Du ) ∼Y (t ) C⏟ D⏟ ⎪ Y⏟ dir dir dir U⏟ dir ⎩
3.4. Direct load calculation solver
{ } Uk (t ) ∼ Y (t ) U⏟ dir
{ }
In Modelica, a non-causal tool, the load calculation is performed iteratively using a proportional integral (PI) controller [24]. In our initial model, before the basis change, the indoor air temperature Tin of one building is taken as the observed variable to control its input P. A carefully controlled value of P that matches Tin with Tset represents the heating or cooling load. Here, none of the values of P set for other building models are utilized in this iterative process; therefore, the initial system (Eq. (1)) is a single-input and single-output (SISO) problem. Although ten PI controllers are required in the case of ten buildings, each controller operates independently. However, the reduced model is covariant to any of the inputs since all the states are connected with each other after introducing a statetransforming matrix for the basis change. This system is called a multiple-input and multiple-output (MIMO) system. Typically, this system
(6)
From this equation, the heating or cooling loads Uuk can be directly ∼ calculated by introducing Tset as the new input term Y(t ) of Eq. (6). Using this formula, hourly heating or cooling loads are obtained under the condition that Tin is equal to Tset. To take into account scheduling such as building occupancy hours, particular attention needs to be paid. A solution is that Eq. (4) is used for non-occupancy periods by giving Uuk = 0. When a building starts to be operated, Eq. (6) is newly run. In this case, X(t) obtained using Eq. (4) must be used as the state vector in the precedent time-step for calculating Ẋ (t ) .
Fig. 7. Annual hourly load comparison of ten Ref models with (a) ten RMs (one by one) and (b) a single RM (assembled). 968
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
4. Results and discussion
model apparently affects the model accuracy. As fewer model states are used, the peak and mean errors increase. The 2-state model for ten buildings shows the highest error among the models; however, this high error is temporary and is only 0.43% for the annual load. This level of error in annual energy demand is quite accurate since the CITYSIM model [8] showed an error of 6% compared to ESP-r reference models [10]. Furthermore, the results of the 2-state RM for the peak heating and cooling periods show that the dynamics and magnitude are similar to those of the Ref model, as shown in Fig. 10. An option is to use the 10state RM. A reasonable magnitude of error is maintained with the 10state RM even for cases with a large number of buildings, e.g., up to a hundred thousand buildings [26].
4.1. Simulation conditions He [25] carried out a detailed comparison test between Buildings and BuildSysPro SM models and found that they were in very good agreement. As mentioned earlier, this outcome was expected because both models use the same physical assumptions. Therefore, the SMs developed in this work using Buildings are used as the Ref models for the following simulation results. Cases of multiple buildings were tested using the parameters presented in Table 1. Different parameters were used for comparison. The internal gains were not considered in this study; however, they can be easily taken into account by including them with the P values in the case of a PI controller or by adding it to Uk (Eq. (6)) in the case of a direct solver. To test our final single RM (assembled) using the proposed direct solver, a separate simulation with an RM or Ref model for each building was performed multiple times using a PI controller (one by one), and the results were summed. In addition, the performance of the state-size reduction technique was tested using the Ref models.
4.4. CPU time comparison The CPU time required for the simulation of different RMs was measured. A computer equipped with an Intel Core i7-4770 3.4 GHz CPU and 16 GB RAM was used as the benchmark for simulating the models in Dymola 2016 with the DASSL solver. Table 2 shows the results of the CPU time for hourly annual simulations of the ten buildings. The assembled RM model without the proposed direct solver is slower than the one-by-one RMs even if they are of similar order. As the single RM is a MIMO system, it requires a greater number of iterations for convergence. For the case of the direct solver, the assembled RM with 50 states requires approximately 4 s. To further investigate the CPU time required by the model with the proposed solver, the assembled RMs with different numbers of states were tested. Fig. 11 shows the CPU time results of the annual hourly simulations. The CPU time of the assembled RM model and its order are linearly proportional. A 4-order model requires less than 1 s for an annual hourly simulation of 10 buildings. This finding implies that a large number of buildings in urban areas can be simulated rapidly when a very-low-order RM with the proposed direct solver is used.
4.2. Validation of the reduction technique The model reduction performance was measured using the Ref model (SM model). For ten buildings, ten Ref models were prepared. Each Ref model was run separately to calculate the annual hourly load. Similarly, the ten RMs were simulated separately. Fig. 7a shows the differences in annual loads between the models (one-by-one Ref vs. one-by-one RM). For the case of ten RMs, each RM model was on the order of 6 states; thus, 60 states were used for 10 buildings. In this figure, the load ranges between +300 kW (for heating) and −300 kW (for cooling). The differences are within 2 kW, which shows good agreement between the Ref models and RMs. Then, the assembled RM model that was obtained by assembling and reducing the ten building models was compared with the Ref model. Fig. 7b shows a comparison between the annual loads of the assembled RM model and the Ref model (one by one). This single representation model (RM assembled) was our target test model. As shown in the figure, this assembled RM model shows similar accuracy to the previous one-by-one RMs. This finding implies the possibility of using a single lower-order SS model to accurately model a number of buildings. However, the model could be tested to only 50 orders. Further investigations to decrease the order of the model could not be performed because of the occurrence of a convergence problem in the MIMO system, which was mentioned previously. Thus, the current iterative method and the corresponding solver are not the best techniques for this assembled RM model.
5. Conclusions In this work, we proposed a single reduced model for hourly load calculation for multiple buildings. According to an urban energy analysis framework, the model can be deduced from the detailed parameter data of buildings. Unlike conventional approaches that require N models for N buildings, the proposed method requires a single
4.3. Single reduced model with a direct solver The proposed direct solver for the single RM was tested. The new direct representation model (Eq. (6)) was also run in Dymola to compare the CPU times. Fig. 8 shows the results of the comparison between the Ref model and the single RM with the proposed solver (RM direct (assembled)). The RM model in this case is on the order of 50 states. This new solver shows good agreement with the Ref model; however, infrequent peak errors are observed. The peak error is approximately 20 kW; however, the other periods exhibit dynamics that are very similar to those of the Ref model. When this solver is used, the convergence problems that were caused by iterative methods do not occur. Therefore, models of lower orders can be tested. Fig. 9 shows the different error levels according to the RM order. The tests were carried out down to the 2-state model. As expected, the proposed solver could run without any numerical problem even with a model of very low order. However, the order of the
Fig. 8. Performance of the direct solver for the RM assembled case. 969
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
Fig. 11. CPU time of RM direct-assembled models with different numbers of states. Fig. 9. Errors of different model orders in the RM direct (assembled).
For the test case, which included ten buildings, this assembled RM model gave similar results with the one-by-one RMs as well as with the SMs; however, the total model order could not be reduced because of convergence issues. In addition, the CPU time was increased, as the assembled RM requires more iterations in Dymola. Therefore, a direct solver was developed to calculate the loads without iteration, which enables reducing the model order further. The final assembled RM with the proposed direct solver can reduce the order down to only 2 states, with an error of 0.43% in the annual load. This means that two equations can predict the thermal dynamics of ten buildings or more. This low-order single representation of buildings could be useful for rapid computation and simulation interoperability for further detailed urban analysis. Various applications that utilize the proposed model are expected in the near future. Acknowledgment This work was supported under the framework of international cooperation program managed by the College of Engineering, Inha University. In addition, this work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIP) (No. 2016R1C1B2011097).
Fig. 10. Two-state RM direct for peak heating and cooling periods.
References Table 2 CPU time for different RMs in Dymola.
RM one by one (60 states) RM assembled (50 states) RM direct assembled (50 states)
CPU time (s)
Sys.
81 499 4
SISO MIMO Direct
[1] An J, Yan D, Hong T, Sun K. A novel stochastic modeling method to simulate cooling loads in residential districts. Appl Energy 2017;206:134–49. [2] Ma W, Fang S, Liu G, Zhou R. Modeling of district load forecasting for distributed energy system. Appl Energy 2017;204:181–205. [3] Frayssinet L, Merlier L, Kuznik F, Hubert JL, Milliez M, Roux JJ. Modeling the heating and cooling energy demand of urban buildings at city level. Renew Sustain Energy Rev 2018;81:2318–27. [4] Robinson D, Stankovic S, Morel N, Dequé F, Rylatt M, Kabele K, et al. Integrated resource flow modelling of urban neighbourhoods: project suntool. Proceedings of building simulation 2003, Einthoven, The Netherlands. 2003. [5] Déqué F, Ollivier F, Poblador A. Grey boxes used to represent buildings with a minimum number of geometric and thermal parameters. Energy Build 2000;31(1):29–35. [6] Robinson D, Haldi F, Kämpf J, Leroux P, Perez D, Rasheed A, Wilke U. CITYSIM: comprehensive micro-simulation of resource flows for sustainable urban planning. Proceedings of Building Simulation 2009, Glasgow, Scotland. 2009. [7] Roux JJ. Proposition de modèles simplifiés pour l'étude du comportement thermique des batiments Phd thesis Lyon (France): Université de Lyon; 1984 [8] Kämpf J, Robinson D. A simplified thermal model to support analysis of urban resource flows. Energy Build 2007;39:445–53. [9] DOE-2 engineers manual, version 2.1A. Technical Information Center, United States Department of Energy; 1982. [10] Courchesne-Tardif A, Kummert M, Demark S, Butler T, Pearl D, Jones S, et al.
assembled model for multiple buildings. Such method is mandatory to predict the heating and cooling loads, in kW, and then for city energy management. This is clearly a progress in comparison with common CES focusing on annual energy demand, in kWh. Similar to previous studies, the state model size reduction technique was applied to this model. The balanced truncation method is effective in reducing the computational time and ensuring the accuracy of the model. In a previous work, this reduction technique was applied to each SM model that represents a single building. In this work, the reduction technique is applied to a single state-space model for the buildings, constructed by assembling N SMs in Dymola. 970
Applied Energy 238 (2019) 963–971
E.-J. Kim et al.
[11] [12]
[13]
[14]
[15]
[16] [17]
[18] Bontemps S, Kaemmerlen A, Blatman G, Mora L. Reliability of dynamic simulation models for building energy in the context of low-energy buildings. Proceedings of building simulation 2013, Chambéry, France. 2013. [19] Dymola version 2013, (10.01.2013); 2013. Homepage: http://www.3ds.com/ products/catia/portfolio/dymola. [20] Peuportier B, Blanc Sommereux I. Simulation tool with its expert interface for the thermal design of multizone buildings. Int J Sol Energy 1990;8:109–20. [21] Wetter M, Zuo W, Nouidui TS. Modeling of heat transfer in rooms in the Modelica buildings library. 12th Conference of international building performance simulation association, Sydney, Australia. 2011. [22] Moore BC. Principal component analysis in linear system: controllability, observability, and model reduction. IEEE Trans Autom Control 1981;AC-26(1):17–32. [23] Palomo E, Bonnefous Y, Déqué F. Guidance for the selection of a reduction technique for thermal model. Proceedings of BS, Prague, Czech Republic. 1997. [24] Ghiaus C. Causality issue in the heat balance method for calculating the design heating and cooling load. Energy 2013;50:292–301. [25] He X. Simulation énergétique à l’échelle quartier: réduction des modèles de bâtiment sous Dymola Master thesis Villeurbanne (France): ENTPE; 2014 [26] Garcia Perez M, Kim E-J, Roux J-J. Reduction techniques for urban energy simulation. 2017 SAREK summer annual conference, Seoul, Republic of Korea. 2017. p. 587–90.
Assessing community-scale energy supply scenarios using TRNSYS simulations. Proceedings of building simulation 2011, Australia, Sydney. 2011. EnergyPlus Engineering Reference. The reference to EnergyPlus calculations. Ernest Orlando Lawrence Berkeley National Laboratory, US Department of Energy; 2010. Huber J, Nytsch-Geusen C. Development of modeling and simulation strategies for large scale urban districts. Proceedings of building simulation 2011, Sydney, Australia. 2011. Kohler M, Blond N, Clappier A. A city scale degree-day method to assess building space heating energy demands in Strasbourg Eurometropolis (France). Appl Energy 2016;184:40–54. Gianniou P, Reinhart C, Hsu D, Heller A, Rode C. Estimation of temperature setpoints and heat transfer coefficients among residential buildings in Denmark based on smart meter data. Build Environ 2018;139:125–33. Martin GP, Kim E-J, Roux J-J. Simulation des besoins énergétiques à l’échelle de la ville: Une méthodologie de modélisation adaptée. Proceedings of IBPSA France Conference, Bordeaux, France. 2018. Kim EJ, Plessis G, Huber JL, Roux JJ. Urban energy simulation: simplification and reduction of building envelope models. Energy Build 2014;84:193–202. Plessis G, Kaemmerlen A, Lindsay A. BuildSysPro: a Modelica library for modelling buildings and energy systems. Proceedings of the Modelica conference 2014, Lund, Sweden. 2014.
971