Detailed Thermodynamic Modeling of Multi-Zone Buildings with Resistive-Capacitive Method

Increased use of energy in buildings and HVAC systems requires advanced control schemes like model-based control to improve energy efficiency, which in turn requires accurate thermodynamic models of buildings. The Resistive-Capacitive (RC) method is a popular and versatile approach for thermal modeling of buildings. Despite this, it is not easy to find practical solutions of implementation of the RC method. It is the goal of this paper to clarify the RC method and demonstrate simple implementation of this method, especially for multi-zone buildings, which have more potential for energy savings from use of model-based control. This paper provides two contributions. First is a detailed explanation of the RC method, focusing on its use for developing a structure of a model and first-principles approach for estimation of parameters of a model. Second is a demonstration of an algorithm that enables automatic development of the structure of a model from basic information about a building (layout, construction elements) and its combination with data-based parameter estimation. Use of the algorithm is presented with a case-study on industrial multi-zone building, for which such a grey-box model is developed and analyzed. The resulting model is rapidly developed and used in a simulation with the measured data. The outputs of the model are compared with the measured temperatures and they show good fit.


Introduction
There are several approaches to how to increase the energy efficiency in buildings: By taking care of the energy efficiency during the design phase and using contemporary technologies; by using materials (better insulation) and devices (higher energy efficiency) that increase it; by using more efficient HVAC systems, properly designed and based on renewable energy sources; by changing occupants' habits to take more care concerning conserving energy; by using building management systems for tracking and optimizing energy usage; and by using advanced control methods for HVAC and connected systems (lighting, windows, blinds, etc.) [1]. The first three approaches are usually expensive (especially concerning existing buildings), while changing occupants' habits is often hard and unreliable. Transfer to the advanced control of HVAC systems can be a relatively cheap and easy way to increase the overall energy efficiency in buildings (especially considering that HVAC systems make up over 30% of the total energy use in buildings [2]). There are strategies for increasing the energy efficiency in buildings that use alternative concepts: Papers [3][4][5] relied on occupancy-based demand controls, while paper [6] used occupant feedback (voting) to find optimal setpoints for HVAC control systems.
The main goal of HVAC control systems is to provide thermal comfort while using minimal energy. Methods based on prediction, especially model prediction control, present one practical and popular solution [7]. To have good prediction, it is necessary to have a good model of the building and HVAC system. One popular thermal modeling method

The RC Method
There are two steps in using the RC method; first it is necessary to develop the structure of the model using the first-principles approach (in [25] this is called "cognitive model development"), while parameters of the model must be estimated in the second step.

Example for a Single Wall
A popular and very simple example of the RC method is the representation of a single wall with the 2R1C model (two resistances and one capacity), shown in Figure 1, where the inside space represented by the temperature TIn is a state and the outside temperature TOut is a disturbance. The wall is represented by two resistances ROut and RIn, representing thermal processes between the wall and the outside environment and between the wall and the inside space, and by capacity CWall, which represents the thermal capacity of the wall, while TWall represents the internal temperature of the wall. The internal space has the capacity of air and elements inside the room CIn. Using the state-space representation, the resulting model is as follows: The calculation of values of R and C parameters can be done using physical properties of wall and construction materials (thickness, density, etc., see [26], for example, or [27] for more detailed analysis). In such a first-principles approach, the resulting model is of white-box type.

Example for a Simple Multi-Zone Building
A more complex task is modeling a simple multi-zone building with the RC method. In this example, a simple two-room building shown in Figure 2 will be used. For simplicity, it is presumed that ceiling and roof are perfect isolators, thus keeping the model two dimensional. Disturbances are the outside temperature and the solar radiation on the walls, while a heater is considered a controllable input. Using the state-space representation, the resulting model is as follows: ·T Out , (1) [T In ] = 0 1 · T Wall T In The calculation of values of R and C parameters can be done using physical properties of wall and construction materials (thickness, density, etc., see [26], for example, or [27] for more detailed analysis). In such a first-principles approach, the resulting model is of white-box type.

Example for a Simple Multi-Zone Building
A more complex task is modeling a simple multi-zone building with the RC method. In this example, a simple two-room building shown in Figure 2 will be used. For simplicity, it is presumed that ceiling and roof are perfect isolators, thus keeping the model two dimensional. Disturbances are the outside temperature and the solar radiation on the walls, while a heater is considered a controllable input. Even for this relatively simple building, there are several approaches on how to define the structure of the model. One approach is to model every element of a building as an RC structure and then combine them in a large RC network [28]. Another approach assumes that the whole building envelope is a single RC structure (usually 3R2C, but other combinations are possible), presumes some thermal mass inside (in the form of additional thermal capacity C) and finally models the thermal flow of a heater as a current source and the outside temperature as a heat-flow source, possibly separating these two spaces by additional RC elements [16,17,25]. Another approach is to assume some "standard structure" for a zone and combine as many structures as necessary [13,20,29]. Additionally, such a structure can take into consideration some properties or features that are in the focus of authors, for example "fast and slow response" [30]. All these structures have at least two purposes: To make defining the model structure easier by providing ready-to-use structures and to keep the model relatively simple. Unfortunately, most papers completely skip explaining the process of designing this structure.
The approach that is explained in this paper is based on the first approach ("full model"), where each element of a building is represented by one RC element. This will increase the complexity of the model by adding additional states (compared to the "standard structures" approach), but it will result in better accuracy, maintain the physical representation of the model and enable isolation of some zones for detailed examination. Moreover, it will keep the process of designing the model structure straightforward (one building element equals one RC element). Added complexity can be handled later either by using methods for reduction of model complexity (e.g., [28]) or using techniques for distributed modeling and control (e.g., [15,31]).
The process of designing the model structure is explained in several steps. First, it is necessary to define two zones (or spaces) that are represented by the temperatures inside rooms and their thermal capacity (this includes the indoor air, furniture and other materials inside). These, and the outside temperature, are represented by potentials TOutside, TRoom1 and TRoom2 and capacities CRoom1 and CRoom2, as in Figure 3.  Even for this relatively simple building, there are several approaches on how to define the structure of the model. One approach is to model every element of a building as an RC structure and then combine them in a large RC network [28]. Another approach assumes that the whole building envelope is a single RC structure (usually 3R2C, but other combinations are possible), presumes some thermal mass inside (in the form of additional thermal capacity C) and finally models the thermal flow of a heater as a current source and the outside temperature as a heat-flow source, possibly separating these two spaces by additional RC elements [16,17,25]. Another approach is to assume some "standard structure" for a zone and combine as many structures as necessary [13,20,29]. Additionally, such a structure can take into consideration some properties or features that are in the focus of authors, for example "fast and slow response" [30]. All these structures have at least two purposes: To make defining the model structure easier by providing ready-to-use structures and to keep the model relatively simple. Unfortunately, most papers completely skip explaining the process of designing this structure.
The approach that is explained in this paper is based on the first approach ("full model"), where each element of a building is represented by one RC element. This will increase the complexity of the model by adding additional states (compared to the "standard structures" approach), but it will result in better accuracy, maintain the physical representation of the model and enable isolation of some zones for detailed examination. Moreover, it will keep the process of designing the model structure straightforward (one building element equals one RC element). Added complexity can be handled later either by using methods for reduction of model complexity (e.g., [28]) or using techniques for distributed modeling and control (e.g., [15,31]).
The process of designing the model structure is explained in several steps. First, it is necessary to define two zones (or spaces) that are represented by the temperatures inside rooms and their thermal capacity (this includes the indoor air, furniture and other materials inside). These, and the outside temperature, are represented by potentials T Outside , T Room1 and T Room2 and capacities C Room1 and C Room2 , as in Figure 3.
The next step is to replace each wall with an RC element. In Section 2.1, a 2R1C structure for representing walls was explained. Usually, such a structure is considered too simple, so a 3R2C element is used to define walls. This means that each wall will be represented by two temperatures, T WallIn and T WallOut , three resistances, R WallIn , R Wall and R WallOut , and two capacities, C WallIn and C WallOut , as shown in Figure 4. Designations "In" and "Out" are arbitrary, but it is easier to follow the same rule throughout the whole structure. The process of designing the model structure is explained in several steps. First, it is necessary to define two zones (or spaces) that are represented temperatures inside rooms and their thermal capacity (this includes the indoor a niture and other materials inside). These, and the outside temperature, are repre by potentials TOutside, TRoom1 and TRoom2 and capacities CRoom1 and CRoom2, as in Figure 3.   The next step is to replace each wall with an RC element. In Section 2.1., a 2R1C structure for representing walls was explained. Usually, such a structure is considered too simple, so a 3R2C element is used to define walls. This means that each wall will be represented by two temperatures, TWallIn and TWallOut, three resistances, RWallIn, RWall and RWallOut, and two capacities, CWallIn and CWallOut, as shown in Figure 4. Designations "In" and "Out" are arbitrary, but it is easier to follow the same rule throughout the whole structure. In the following steps it is necessary to add RC elements that represent doors and windows, and then to add heat flows from the sun and the heater. Doors and windows (or commonly openings) are usually represented only with resistances connected in series (due to their thinness, it is assumed that they do not have thermal capacity). So in this example, they are represented by 3R structure or three resistances (for example, RWindowIn, RWindow, RWindowOut). Heat is represented by the current input QHeater which directly affects Room 2 (TRoom2). The solar radiation is represented by current inputs QSolarNorth, QSolarEast, QSolarSouth and QSolarWest, which directly affect walls. To simplify this example, it is assumed that there is no solar radiation through the windows and that there is no internal radiation.
The final structure is represented in Figure 5. With this, the model is defined as an electrical equivalent network that represents the structure of the model of the simple building from Figure 2. There are in total 16 capacities, 33 resistances, 17 potentials and five current inputs. In the following steps it is necessary to add RC elements that represent doors and windows, and then to add heat flows from the sun and the heater. Doors and windows (or commonly openings) are usually represented only with resistances connected in series (due to their thinness, it is assumed that they do not have thermal capacity). So in this example, they are represented by 3R structure or three resistances (for example, R WindowIn , R Window , R WindowOut ). Heat is represented by the current input Q Heater which directly affects Room 2 (T Room2 ). The solar radiation is represented by current inputs Q SolarNorth , Q SolarEast , Q SolarSouth and Q SolarWest , which directly affect walls. To simplify this example, it is assumed that there is no solar radiation through the windows and that there is no internal radiation.
The final structure is represented in Figure 5. With this, the model is defined as an electrical equivalent network that represents the structure of the model of the simple building from Figure 2. There are in total 16 capacities, 33 resistances, 17 potentials and five current inputs. In the next step, the model is to be represented in a mathematical form, in the state-space representation. This is a linear model, with six inputs in total (the temperature of the outside environment, the solar radiation divided by cardinal directions and the heat flow from the heater, of which the last is a controllable input) and 16 states (all temperatures except the outside; this number is equal to the number of capacities). Vectors of the state-space representation are as follows: ]. ( The actual order of states in the vectors is arbitrary. Here, for the state vector T the rule was to put first the temperature on the "upper" or "left" side (depending on whether the wall is "horizontal" or "vertical" on a two-dimensional layout). With the input vector u, first is the outside temperature, followed by the components of the solar radiation (in "clockwise" order as was for the walls) and in the end the controllable input (heater). The output vector y is always composed of temperatures in zones, forwarding values from states. In the next step, the model is to be represented in a mathematical form, in the statespace representation. This is a linear model, with six inputs in total (the temperature of the outside environment, the solar radiation divided by cardinal directions and the heat flow from the heater, of which the last is a controllable input) and 16 states (all temperatures except the outside; this number is equal to the number of capacities). Vectors of the state-space representation are as follows: The actual order of states in the vectors is arbitrary. Here, for the state vector T the rule was to put first the temperature on the "upper" or "left" side (depending on whether the wall is "horizontal" or "vertical" on a two-dimensional layout). With the input vector u, first is the outside temperature, followed by the components of the solar radiation (in "clockwise" order as was for the walls) and in the end the controllable input (heater). The output vector y is always composed of temperatures in zones, forwarding values from states. For the first wall, there are two equations describing all relations (analogy with current): This can be written, following the order of the states as in the state vector, as .
Finally, using generalized coefficients that correspond to elements of the state and input matrices, it can be represented as: .
Such equations can be developed for all walls and finally rooms, resulting in the state-space representation of the model: .
If, for demonstration purposes, it is assumed that all parameter R and C are all equal to 1, parameter a 1,1 will be calculated as: So, the resulting state matrix A, input matrix B, output matrix C and feed-forward matrix D have values as following:

Estimation of Parameters
As mentioned beforehand, calculation of the parameters of the model can be approached in two ways: Using the physical approach for calculating the values of parameters from the nominal values of thermal properties of materials, elements and parts of the building (thus, resulting in a white-box model); or using the data-based approach to estimate the parameter values from the measured data (resulting in a grey-box model). Here, it is shown how to use the first approach, while the second approach is explained later.
The thermal transfer between the outside environment and the surface of the wall is mostly convective (represented in the 3R2C model of the wall with R WallOut and R WallIn ). Newton's laws of cooling define convective heat transfer as: stating that the heat flow is proportional to the heat transfer coefficient h (which depends on the type of convection, for example, natural or forced), surface of the wall A and the difference in temperature between the surface and the fluid ∆T (e.g., air of the environment). If this is written as: then the thermal resistance of convection can be defined as: The heat transfer through the wall is mostly conductive (represented in the 3R2C model of the wall with R Wall ). Conduction is explained by Fourier's law, which states that the heat flow is proportional to the coefficient of thermal conductivity k, the surface of the wall A and the negative difference of the temperature between the surface ∆T and inverse to the thickness of the wall l: Thus, the thermal resistance of conduction can be defined as: Radiation (for example, between two walls, which was ignored in the previous sample) is governed by the simplified Stefan-Boltzmann equation, which says that the heat flow between two walls is proportional to the radiation thermal conductance h, surface A and the temperature difference between two bodies ∆T raised to the fourth power: From this, the resistance of radiation in this case can be defined as: Finally, the thermal capacity C, a property of body (in this case, some element of a building) to hold heat, can be defined as the heat Q necessary to change the temperature of this body for some difference ∆T: The thermal capacity itself is calculated as: where c is the specific heat capacity, ρ is the density and V is the volume. For more information, several sources are available: Underwood and Yik [27], Bejan and Kraus [32], Davies [33] and Janna [34].

Detailed Modeling of Multi-Zone Buildings with Algorithm and Optimization of Parameters
As demonstrated in the previous chapter, manually developing the thermodynamic model of a building based on the RC method is relatively simple and straightforward, but a long and time-consuming process. Due to many elements in buildings, especially for larger buildings, the likelihood of some error is high.
The estimation of parameters for such model structures is theoretically simple, based only on several equations. In practice though, a large amount of information is unknown or unreliable: Thickness and composition of materials in walls, influence of ageing of material, etc. Sometimes even getting basic information about used materials can be impossible. This chapter will focus on two tools that improve developing the detailed thermodynamic models with the RC method: First is an algorithm that enables automatic development of the structure of the model; and second is the optimization of model parameters based on the measured data.

The Algorithm for Automatic Definition of the Model Structure Based on the RC Method
When analyzing matrices of the state-space representation of thermodynamic models developed using the RC method, some similarities and patterns can be recognized. First, it is evident that matrices are sparse, which can improve the speed of their operations. Next, it can be seen that the state matrix always follows the same symmetrical shape. For example, for the 3R2C representation of a wall, there will be four non-zero elements on the diagonal, representing thermal processes between inner and outer layers of the wall and two non-zero elements at the ends of the row and the column, representing thermal processes between the wall and surrounding zones. Even in such cases when states are ordered differently, elements are just accordingly moved in the matrix. Similar pattern can be recognized in the input matrix, where there will be non-zero elements for the outer layer of the walls, which represent thermal processes with the outside environment and the solar radiation (or other heat flows). The output matrix just forwards states of zones to the outputs, while the feed-forward matrix is always empty. With experience, it is possible to recognize physical interpretation intuitively or design the structure without actual equations.
This regularity enables automation of the process of model development with an algorithm that defines the model based on the layout of the building.
First, it is necessary to define structures that the algorithm will use: • Source I-Represents a heat flow to the system (radiator, the solar radiation, etc.) or input of the system. Its only parameter is type, which states whether this input is a disturbance (e.g., the solar radiation) or a manipulated variable (e.g., the influence of HVAC system).

•
Space P-Represents a temperature inside a zone. This is both an output and a state of the system. Its parameters are: Type (required), stating if this is a zone or an outside space (special case in which this temperature is a disturbance); thermal capacity C (required); source (not required), which defines if this space is affected by some (previously defined) source I.

•
Wall Z-Represents a wall between two spaces P (either interior or inside; which includes floors and ceilings) and as such adds states to the system (but not outputs). Its parameters are: Spaces P 1 and P 2 (required), defining which of two (previously defined) spaces is this wall; physical properties (required), or values of R and C parameters (exact number or R and C parameters depends on the desired complexity, e.g., 2R1C or 4R3C, etc.); sources (not required), defining if this wall is affected by some (previously defined) sources I 1 and I 2 and on which side.

•
Opening O-Represents a door or window between spaces (similar to walls). Its parameters are: Spaces (required), defining which of two (previously defined) spaces P 1 and P 2 is this opening; physical properties (required); or values of R. There are two assumptions: That openings do not have thermal capacity and that they cannot be influenced by heat sources.
The algorithm processes information about the building contained in these structures and creates a thermodynamic model of the building in the state-space representation: Vectors (states T, inputs u and outputs y) and matrices (system A, input B, output C and feed-forward D). It is assumed that these structures are ordered in indexed lists (sources, spaces, walls, opening) with the first element on position zero. Moreover, it is assumed that the outside space is first in the list of spaces, followed by the disturbances and ending with the controllable inputs. There are seven steps of the algorithm and they will be demonstrated for the 3R2C version of the RC method, as specified in Table 1. For details, see [35], where an implementation of this algorithm in application ModelBuilder is presented. Table 1. Algorithm for defining the structure of a thermodynamic model of a building.
Step Action

1.
Definition of the number of states, inputs and outputs and their vectors (T, u and y) and matrices (A, B, C and D)

2.
Definition of the structure and nominal values of elements of the state matrix A from relationships between walls Z and spaces P

3.
Modification of the nominal values of elements of the state matrix A from relationships between openings O and spaces P

4.
Definition of the structure and nominal values of elements of the input matrix B from relationships between walls Z, spaces P and sources I

5.
Modification of the nominal values of elements of the input matrix B from relationships between openings O and spaces P Step Action

6.
Modification of the nominal values of elements of the input matrix B from relationships between spaces P and sources I

7.
Definition of the structure and values of elements of the output matrix C (passing values of temperatures of spaces) and the feed-forward matrix D (always a null matrix) The first step, showed in Algorithm 1, is counting the number of states, inputs and outputs and the definition of vectors and matrices. Algorithm 1. The definition of the number of states, inputs and outputs and their vectors and matrices (pseudo-code representation) states_count = length of walls * 2 + length of spaces -1 inputs_count = length of inputs + 1 outputs_count = lengths of spaces − 1 T = concatenate walls and spaces without outside space u = concatenate outside space and sources y = spaces without outside space The second step, showed in Algorithm 2, is the definition of the structure of state matrix A and its initial elements' values from relationships between walls and space. The third step, showed in Algorithm 3, is the modification of values of elements of the state matrix A from relationships between openings and spaces. As already mentioned, openings do not increase the number of states, but they modify parameters of the model. With this, the state matrix A is fully defined. The fourth step, showed in Algorithm 4, is the definition of the structure of the input matrix B and its initial values from relationships between walls, space and sources. The fifth step, showed in Algorithm 5, is the modification of the input matrix B values from relationships between openings and spaces. The seventh and last step, showed in Algorithm 7, is the definition of the output matrix C. As already said, the feed-forward matrix D is empty. With this, the state-space representation of the appropriate thermodynamic model of a building is finished.
ModelBuilder is an implementation of the algorithm developed in the C# environment. It enables the definition of the structures (sources, spaces, walls and openings) as graphical objects. This way, it is possible to transform the building's layout to the structure of a thermodynamic model. ModelBuilders has several features: The 2R1C, 3R2C and 4R3C versions of the RC method and conversions between them; easy definition and change of parameters of the model (including values of physical properties for typical construction materials); integration with MATLAB (export of model); generation of the input data and simulation; etc. It has limitations that spaces can be affected by only one source, while walls can be affected by two sources, but both limitations can be avoided by aggregating multiple sources.

Optimization of the Model Parameters Based on the Measured Data
A model developed using the RC method previously described is a white-box model, where initial parameters are calculated from the nominal values of properties of elements: Physical dimensions and thermal properties of used materials. Of course, when dealing with real buildings, there are numerous reasons why such information is not accessible or is not reliable: • missing technical data about the building (lost documentation or documentation being not detailed enough), • suspicion about correctness of information (did the constructor really build according to the technical documentation), • parameter or properties that are unknown or hard to measure (for example, the heat transfer coefficient that depends on the speed of wind), • ageing of materials and structures, etc.
Even when all data is available, there is a possibility that the model is not correct or that some ignored properties affect it more than planned. In such cases, grey-box models, where parameters are estimated from measured data, will prove more accurate.
The estimation of parameters can be done without any initial values. In a case where some relatively correct initial values already exist (at least if they are in the right range; even thermodynamic models with slightly wrong parameters can provide adequate results [36]), this will enable much faster re-estimation or optimization of parameters.
The approach selected for optimization of parameters in this paper is the offline minimization of the difference (error) between the measured data and the output of the model. Let f R (A R ,B R ) be the referent model, where A R and B R are its state and input matrices, represented by the set of input-output data y R = f R (u), and let f T (A T ,B T ) be the testing model with A T , B T and y T . Then the goal of the minimization is to optimize the matrices A T and B T in such a way to reduce y R −y T : Before performing the minimization, it is necessary to agree on a criterion function for error between the measured data and the output of the model. There are several frequently used errors, such as absolute error, root-mean-square error and maximum error. Here, for the criterion of the minimization, but also as a rating of the accuracy of the model, Root-Mean-Square Error (RMSE) is selected. For this case (time series of several values), it is defined as: where Q is the number of outputs of the model and K is the number of measurements in the time series. During the research, several different approaches were tried for minimization: • Selection of the solver and the algorithm for optimization problem: fminsearch (Simplex/ Derivative-Free), fminunc (Quasi-Newton), fmincon (Interior-Point, SQP, Active-Set), lsqnonlin (Levenberg-Marquardt algorithm) [37]; For solving the minimization problem, the quasi-Newton algorithm based on the BFGS method was used in this paper, with the RMSE as the minimization criterion and 900 s as the sampling period. This combination of optimization settings resulted in fast (few minutes for the testing model) and reliable optimization results.
Another criterion that can be used for rating a model's accuracy is the maximum error (MaxErr) on whole range of time series: This type of error is not a good choice to use in minimization, because it is very dependent on outliers (for example, opening of some window during measuring can produce an outlier), so it will lead to minimization of outliers and damage the overall accuracy of the model. On the other hand, it is very indicative as a measure of accuracy of the model. RMSE is also sensitive to outliers (due to squaring of error), but much less than MaxErr.
As already mentioned, since the RC method maintains the physical representation, estimated parameters should be in an acceptable range so it is possible to validate them by comparing them to parameters calculated using the first-principles approach.

Case Study
The test of the proposed method was performed on a building in Croatia built in 2017 [38], shown in Figure 6. This is an industrial building with production space, offices, hot-water-based HVAC and sufficient control system and sensors for testing. The building was used in the usual production processes and the experiment did not affect its use.
where Q is the number of outputs of the model and K is the number of measurements in the time series.
For solving the minimization problem, the quasi-Newton algorithm based on the BFGS method was used in this paper, with the RMSE as the minimization criterion and 900 s as the sampling period. This combination of optimization settings resulted in fast (few minutes for the testing model) and reliable optimization results.
Another criterion that can be used for rating a model's accuracy is the maximum error (MaxErr) on whole range of time series: This type of error is not a good choice to use in minimization, because it is very dependent on outliers (for example, opening of some window during measuring can produce an outlier), so it will lead to minimization of outliers and damage the overall accuracy of the model. On the other hand, it is very indicative as a measure of accuracy of the model. RMSE is also sensitive to outliers (due to squaring of error), but much less than MaxErr.
As already mentioned, since the RC method maintains the physical representation, estimated parameters should be in an acceptable range so it is possible to validate them by comparing them to parameters calculated using the first-principles approach.

Case Study
The test of the proposed method was performed on a building in Croatia built in 2017 [38], shown in Figure 6. This is an industrial building with production space, offices, hot-water-based HVAC and sufficient control system and sensors for testing. The building was used in the usual production processes and the experiment did not affect its use.

Building, HVAC System, Measurement
The building has a concrete construction, while the outside walls are from thermal panels and triple-glazed windows. The production area of the building (614 m 2 area, with height of 6 m) was used in the experiment, consisting of four zones: Assembly, entry, mechanical and warehouse, shown in Figure 7.
The HVAC system for the production area consists of a 35 kW gas-boiler (the medium is hot water with temperature around 70 • C) and ceiling radiators (187 W/m) for distribution of heat in zones. The control system based on PLC regulates the temperature of the medium (by a three-way valve) and the flow of the medium to radiators in each zone (by an open-close valve for each zone).
The building has a concrete construction, while the outside walls are from thermal panels and triple-glazed windows. The production area of the building (614 m 2 area, with height of 6 m) was used in the experiment, consisting of four zones: Assembly, entry, mechanical and warehouse, shown in Figure 7. The HVAC system for the production area consists of a 35 kW gas-boiler (the medium is hot water with temperature around 70 °C) and ceiling radiators (187 W/m) for distribution of heat in zones. The control system based on PLC regulates the temperature of the medium (by a three-way valve) and the flow of the medium to radiators in each zone (by an open-close valve for each zone).
The measurements consist of following:  a temperature sensor in each zone (the assembly has two sensors, but an average was used),  a calorimeter for measuring the heat flow from the boiler (two temperature sensors on pipes and flow sensor),  a weather station on the outside, capable of measuring the outside temperature and the solar radiation (among other measurements).
The measurements on the site were performed in winter conditions (January and February of 2019), at a 15 min interval.
The measured data was separated into several scenarios: Three training scenarios with a duration of six days (Winter-1 with outside temperatures between −10 °C and 5 °C, Winter-2 with 0 °C to 15 °C range, and Winter-3 with −5 °C to 10 °C range) and 15 testing scenarios with a duration of two days (Winter-1-1 to Winter-1-5, Winter-2-1 to Winter-2-5 and Winter-3-1 to Winter-3-5), where some testing scenarios consisted of the same data as training scenarios (only shorter). Data used in the experiment was raw, without any preprocessing (for example, removal of outliers).
There were occupants in the building during the first shift (from 7:00 to 15:00). Temperatures in zones were set to 20 °C in assembly and 18 °C in other zones during the first shift and ignored at other times.

Initial Model of Building End Estimation of Parameters
The model structure for this building is created in the ModelBuilder application and presented in Figure 8. The 3R2C version of the RC method was selected, for giving good enough complexity to encompass thermal relationships in the building, but not adding too many states to the model (easy conversion to the 2R1C or 4R3C versions is possible). The resulting model has 38 states (13 walls, four segments of roof and four zones), six inputs (the outside temperature, the solar radiation and four controllable heat flows from ceiling radiators) and four outputs (temperatures in zones). Its state matrix A has dimension 38 × 38, input matrix B 38 × 6, output matrix C 4 × 38 and feed-forward matrix D 4 × 10 (null-matrix). The measurements consist of following: • a temperature sensor in each zone (the assembly has two sensors, but an average was used), • a calorimeter for measuring the heat flow from the boiler (two temperature sensors on pipes and flow sensor), • a weather station on the outside, capable of measuring the outside temperature and the solar radiation (among other measurements).
The measurements on the site were performed in winter conditions (January and February of 2019), at a 15 min interval.
The measured data was separated into several scenarios: Three training scenarios with a duration of six days (Winter-1 with outside temperatures between −10 • C and 5 • C, Winter-2 with 0 • C to 15 • C range, and Winter-3 with −5 • C to 10 • C range) and 15 testing scenarios with a duration of two days (Winter-1-1 to Winter-1-5, Winter-2-1 to Winter-2-5 and Winter-3-1 to Winter-3-5), where some testing scenarios consisted of the same data as training scenarios (only shorter). Data used in the experiment was raw, without any preprocessing (for example, removal of outliers).
There were occupants in the building during the first shift (from 7:00 to 15:00). Temperatures in zones were set to 20 • C in assembly and 18 • C in other zones during the first shift and ignored at other times.

Initial Model of Building End Estimation of Parameters
The model structure for this building is created in the ModelBuilder application and presented in Figure 8. The 3R2C version of the RC method was selected, for giving good enough complexity to encompass thermal relationships in the building, but not adding too many states to the model (easy conversion to the 2R1C or 4R3C versions is possible). The resulting model has 38 states (13 walls, four segments of roof and four zones), six inputs (the outside temperature, the solar radiation and four controllable heat flows from ceiling radiators) and four outputs (temperatures in zones). Its state matrix A has dimension 38 × 38, input matrix B 38 × 6, output matrix C 4 × 38 and feed-forward matrix D 4 × 10 (null-matrix).
Adjacent offices have an independent heat system. Since there is a good thermal insulation of the connecting wall and concerning the fact that the temperature in these spaces was kept constantly around 21 • C, they were ignored in this model. The thermal influence of people and equipment in zones was also ignored, as was the opening of doors (which was minimal, especially in assembly). Due to the good thermal isolation of the building's foundation (8 cm XPS bricks under the concrete base), heat losses of floors towards the ground were also ignored. Adjacent offices have an independent heat system. Since there is a good thermal insulation of the connecting wall and concerning the fact that the temperature in these spaces was kept constantly around 21 °C, they were ignored in this model. The thermal influence of people and equipment in zones was also ignored, as was the opening of doors (which was minimal, especially in assembly). Due to the good thermal isolation of the building's foundation (8 cm XPS bricks under the concrete base), heat losses of floors towards the ground were also ignored.
The initial model has parameters that are calculated using the nominal properties of typical construction materials (concrete, glass, etc.), so it is a white-box type. The used values can be found in Appendix A, Table A1.
In the next step, the initial model is converted to the grey-box model using the minimization of error between the measured data and the output of the model, with root-mean-square error (RMSE) as criterion. The minimization was performed using MATLAB and nonlinear programming solver for finding the minimum of the unconstrained multivariable function (MATLAB Optimization Toolbox function fminunc) [37].
Optimization of the parameters was concentrated only on the non-zero parameters of the initial model in the state-space representation, first because only these parameters show the relations between the states and inputs; and second, to increase the speed of the optimization. Table 2 shows the values of initial and optimized parameters for one wall (wall_618 in Figure 8.). Since the optimization of the parameters was not constrained, the optimizer reduced the values. The complete set of values of the parameters can be found in [39].  The initial model has parameters that are calculated using the nominal properties of typical construction materials (concrete, glass, etc.), so it is a white-box type. The used values can be found in Appendix A, Table A1.
In the next step, the initial model is converted to the grey-box model using the minimization of error between the measured data and the output of the model, with root-mean-square error (RMSE) as criterion. The minimization was performed using MATLAB and nonlinear programming solver for finding the minimum of the unconstrained multivariable function (MATLAB Optimization Toolbox function fminunc) [37].
Optimization of the parameters was concentrated only on the non-zero parameters of the initial model in the state-space representation, first because only these parameters show the relations between the states and inputs; and second, to increase the speed of the optimization. Table 2 shows the values of initial and optimized parameters for one wall (wall_618 in Figure 8.). Since the optimization of the parameters was not constrained, the optimizer reduced the values. The complete set of values of the parameters can be found in [39]. During the optimization, it is necessary to remember that the model is intended to be used in the model predictive control scheme and for this it should have good prediction of thermodynamic properties (good generalization). This means that the model should be able to predict future outputs for previously unseen data. Buildings will behave differently in different scenarios (winter versus summer). To test the generalization, the model was optimized (trained) with a dataset from one scenario (Winter-1) and tested with data from different scenarios (Winter-2 and Winter-3). For use in different scenarios (e.g., summer), the best course would be to re-optimize the model using a dataset with similar conditions (outside temperature, solar radiation, etc.).

Results
For the training phase, the Winter-1 scenario was selected. Figure 9 shows results of the initial model during six days of the scenario. The outside temperature (green line) varies between −5 • C and 5 • C, while the heating (yellow line, measured on the right axis) turns on when the temperature in the assembly (blue line) drops below the setpoint (20 • C during the first shift). The output of the model (red line) generally follows measured values, but it is evident that differences are large (RMSE 4. During the optimization, it is necessary to remember that the model is intended to be used in the model predictive control scheme and for this it should have good prediction of thermodynamic properties (good generalization). This means that the model should be able to predict future outputs for previously unseen data. Buildings will behave differently in different scenarios (winter versus summer). To test the generalization, the model was optimized (trained) with a dataset from one scenario (Winter-1) and tested with data from different scenarios (Winter-2 and Winter-3). For use in different scenarios (e.g., summer), the best course would be to re-optimize the model using a dataset with similar conditions (outside temperature, solar radiation, etc.).

Results
For the training phase, the Winter-1 scenario was selected. Figure 9 shows results of the initial model during six days of the scenario. The outside temperature (green line) varies between −5 °C and 5 °C, while the heating (yellow line, measured on the right axis) turns on when the temperature in the assembly (blue line) drops below the setpoint (20 °C during the first shift). The output of the model (red line) generally follows measured values, but it is evident that differences are large (   As mentioned before, to test the generalization of a model, it is necessary to test it with a different dataset than it was optimized with, but with similar working conditions (mainly, the outside temperature). In this case, the initial model is tested with scenario Winter 2-3. Results, shown in Figure 11, are RMSE of 2.74 °C and MaxErr of 5.1 °C for the temperature in the assembly. From the figure, it is evident that the outputs of the model generally follow measurements, but such large errors are not usable, especially in model predictive control. For other zones, results are similar: Mechanical (RMSE 4.12 °C and As mentioned before, to test the generalization of a model, it is necessary to test it with a different dataset than it was optimized with, but with similar working conditions  As mentioned before, to test the generalization of a model, it is necessary to test it with a different dataset than it was optimized with, but with similar working conditions (mainly, the outside temperature). In this case, the initial model is tested with scenario Winter 2-3. Results, shown in Figure    As the testing scenario moves away from the training scenario, as in Figure 13   As the testing scenario moves away from the training scenario, as in Figure 13  As the testing scenario moves away from the training scenario, as in Figure 13

Conclusions
This paper presented a grey-box approach for detailed modeling of multi-zone building thermodynamics. The approach is based on the Resistive-Capacitive equivalent method (the RC method) for structural modeling and the optimization of the initial parameters by the minimization of the difference between the measurements and the model outputs. The RC method is flexible enough to model buildings of different and complex structures and to provide models of desired complexity. This approach is implemented in the algorithm that enables simple and efficient development of models with sufficient accuracy to use in control schemes as model predictive control and is explained step by step. Using the data-based estimation of parameters ensures that parameters are close to the real values and give accurate outputs of the model, while combination with the RC method maintains physical representation and enables that parameters can be quickly compared with the expected values.
This approach was tested on a real multi-zone building with measurements acquired in January 2019. The algorithm described in the paper was used to develop a model of multi-zone building in Valpovo, Croatia, with 38 states and six inputs. In the

Conclusions
This paper presented a grey-box approach for detailed modeling of multi-zone building thermodynamics. The approach is based on the Resistive-Capacitive equivalent method (the RC method) for structural modeling and the optimization of the initial parameters by the minimization of the difference between the measurements and the model outputs. The RC method is flexible enough to model buildings of different and complex structures and to provide models of desired complexity. This approach is implemented in the algorithm that enables simple and efficient development of models with sufficient accuracy to use in control schemes as model predictive control and is explained step by step. Using the data-based estimation of parameters ensures that parameters are close to the real values and give accurate outputs of the model, while combination with the RC method maintains physical representation and enables that parameters can be quickly compared with the expected values.
This approach was tested on a real multi-zone building with measurements acquired in January 2019. The algorithm described in the paper was used to develop a model of multi-zone building in Valpovo, Croatia, with 38 states and six inputs. In the next step, the parameters of the model were optimized based on the measurements and the offline model was tested in different scenarios. Results of the tests with occupants in the building show a root-mean-square error of 0.33 • C between the model and the measurements for similar scenarios and 0.64 • C for different scenarios, with good generalization of the model, while the model was developed with minimal investment of time and effort.

R Out
Thermal resistance between the wall and the outside space K/W R In Thermal resistance between the wall and the inside (room) K/W C Wall Thermal capacity of the wall J/K T Wall Internal temperature of the wall • C C In Thermal capacity of air and element of the inside space (room) J/K Example for a simple multi-zone building T Outside Temperature of the outside space • C

T Room1
Temperature inside room 1 • C

T Room2
Temperature inside room 2 • C

C Room1
Thermal capacity of air and element inside room 1 J/K C Room2 Thermal capacity of air and element inside room 2 J/K

T WallIn
Temperature on the surface of the wall on the inside • C (considering reference zone/room)

T WallOut
Temperature on the surface of the wall on the outside • C (considering reference zone/room)

R WallIn
Thermal resistance of the wall and the inside (considering K/W reference zone/room)-convection R Wall Thermal resistance inside the wall-conduction K/W R WallOut Thermal resistance of the wall and the outside (considering K/W reference zone/room)-convection

C WallIn
Thermal capacity of the inside part of the wall J/K (considering reference zone/room)

C WallOut
Thermal capacity of the outside part of the J/K wall (considering reference zone/room)

R WindowIn
Thermal resistance of the windows and the inside K/W (considering reference zone/room)-convection R Window Thermal resistance inside the windows-conduction K/W Output vector a i,j Element of system matrix A in i-th column and j-th row b i,j Element of input matrix B in i-th column and j-th row c i,j Element of output matrix C in i-th column and j-th row d i,j Element of feed-forward matrix D in i-th column and j-th row Object representing a source (heat flow) P (Space) Object representing a space (zone, room) P.C Parameter of object Space with value of its thermal capacity P.I Parameter of object Space, pointer to object Source that is affecting this space Z (Wall) Object representing a wall (vertical wall, floor, ceiling) Z.P Parameter of object Wall, pointer to object Space that is adjacent to this wall (always two such parameters)

Z.R
Parameter of object Wall with value of its thermal resistance (one or more, ordered as points to adjacent Spaces)

Z.C
Parameter of object Wall with value of its thermal capacity (one or more, ordered as points to adjacent Spaces) Z.I Parameter of object space, pointer to object Source that is affecting this space (always two, ordered as points to adjacent Spaces) O (Opening) Object representing an opening (door, window)

O.P
Parameter of object Opening, pointer to object Space that is adjacent to this wall (always two such parameters) O.R Parameter of object Opening with value of its thermal resistance (one or more, ordered as points to adjacent Spaces) sources Ordered indexed list with all sources (objects I) spaces Ordered indexed list with all spaces (objects P) walls Ordered indexed list with all walls (objects Z) openings Ordered indexed list with all sources (objects O) states Vector with states (of the system) inputs Vector with inputs (to the system) outputs Vector with outputs (from the system) i, j, k Indexes