Parameter-input estimation of RC thermal models of buildings using unscented Kalman filter and nonlinear least square method

Abstract

Effective building energy management (e.g. temperature control strategies) necessitates reliable and computationally efficient building thermal models. One type of them is the resistor–capacitor (RC) model. However, estimating model parameters and inputs (e.g. solar heat gain) simultaneously is challenging, especially when some of the temperature states are missing due to instrumentation limitations and/or sensor malfunctions. The present study utilizes unscented Kalman filter (UKF) and nonlinear least squares (NLSs) methods for parameters and input estimation of RC models with possible unavailable temperature states. The estimation procedure, mathematical operations and result analysis are presented in detail. To evaluate the capability of the method, two case studies were conducted. The first case study involved a simple, made-up RC model with known parameters, inputs and states, while the second case study used monitored data from a single detached house. The capability of the method was evaluated by comparing the estimated parameters, inputs and states to the corresponding true values in both study cases. The performance evaluation shows that the proposed method can effectively estimate RC model parameters and inputs, even with certain missing states. The proposed method can be employed for timely online updating of RC model parameters to improve response prediction.

Introduction

Buildings currently consume approximately 40% of global energy and contribute to 31% of world CO₂ emissions ¹ ; therefore, reducing buildings’ energy consumption is critical to meeting global sustainability targets.^2,3 Amongst the efforts to reduce buildings’ energy consumption, building energy management with control systems that can optimize building energy consumption is being developed and deployed.^4,5 These systems are proven to have the potential to save energy by up to 28%. ⁶ Control systems require reliable thermal dynamic models in order to estimate buildings’ thermal parameters (e.g. effective thermal resistance) and behaviour (e.g. temperature response). ⁴

Different modelling strategies, including ‘white-box’, ‘black-box’ and ‘grey-box’ modelling, have been developed in the literature to predict the thermal behaviour of buildings and their systems.^7–11 White-box (physics-based) modelling, such as using Energy Plus, ¹² requires detailed construction information (e.g. dimensions) and can simulate building physics in detail. ¹³ Due to the unknowns of building conditions (e.g. material degradation, users’ behaviour and operation strategies), using this modelling strategy on existing buildings can be problematic. ¹² On the other hand, black-box modelling uses pure mathematical machine-learning techniques, such as using artificial neural network,^14,15 to create non-interpretable (hidden) relationships between input and output data without necessitating physics. ¹⁶ Black-box modelling requires a large amount of high-quality data. ¹⁷ ‘Grey-box’ modelling strategy, such as thermal resistor–capacitor (RC) models, is based on a physics-base model structure and uses mathematical optimization techniques to estimate the equivalent physical parameters of the model. The strategy incorporates the benefits of both white-box to eliminate outliers and that of black-box strategies to reduce the necessity of detailed information. The grey-box modelling approach has been widely used for thermal dynamic modelling.^18–22

RC models are based on a set of equivalent model parameters, resistors (Rs) and capacitors (Cs), to relate system inputs (e.g. heating and cooling supply) and temperature states. ²³ The model parameters can be influenced by both internal and external factors. Internal factors include people, furniture and electric consumption, while external factors encompass adjacent obstructions, air infiltration and weather-related influences. Thus, RC model parameters are not meant to be directly measured or calculated simply based on the construction details of buildings. The unknown parameters and inputs can be estimated through calibration with measured data (usually temperature states).^23–25

RC model parameter estimation is typically an inverse optimization problem, which can be solved by various algorithms (e.g. genetic algorithms, ²⁶ least squares regression system identification, ²⁷ stochastic filtering, ²⁸ linear Kalman filter, ²⁹ extended Kalman filter, ²⁸ ensemble Kalman filter ³⁰ and unscented Kalman filter (UKF) ²⁵ ). Using an extended Kalman filter and UKF with low-sampling-rate historical data, Baldi et al. ³¹ estimated the unknown states and parameters of a simple RC model. Later on, Li et al. ²⁵ demonstrated the capability of UKF to jointly estimate the unknown states and model parameters in complex RC models representing real-life buildings. Additionally, Zamani et al. ²¹ introduced a novel method integrating RC models and UKF with nonlinear least squares (NLSs) for precise estimation of heating and cooling supply, essential for effective temperature control in buildings. Through case studies, its accuracy in estimating supply was demonstrated, affirming its suitability for temperature control objectives.

Parameters can be obtained via offline estimation (i.e. batch) assuming the parameters are constant within the period of a set of collected measurements and using the measurements to estimate the parameter. This approach may result in inaccurate estimations and consequently inaccurate predictions due to time-variant building physical conditions and occupants’ behaviours.^28,32,33 To address this issue, parameter estimation can be conducted in an online manner (e.g. adaptive), in which the model parameters are updated continuously with new measurements in a computationally efficient manner.^25,34–36 There are various online parameter estimation methods in the literature, including recursive least square, ³⁷ recursive prediction error minimization, ³⁸ linear Kalman filter, extended Kalman filter and UKF. Numerous studies in the literature have employed the Kalman filter methods to demonstrate how RC models and the Kalman filter can work well together for online parameter estimations.^25,39 In particular, Radecki et al. ³⁹ demonstrated that temperature predictions are accurate when UKF is used to estimate the parameters of a multi-zone thermal RC model. Later on, a framework was proposed by Maasoumy et al. ⁴⁰ for an online building model that simultaneously estimates the states and unknown parameters and continuously tunes the model parameters. In another study, extended Kalman filter was applied to an RC model of a passive house ²⁸ to create an online thermal model that can accurately predict indoor temperature states.

The abovementioned parameter estimation work focused on cases where all system inputs (e.g. heating and cooling supply) are known; however, in real practice, system inputs and even some temperature states may be unknown for various reasons. For example, the system inputs and states can be temporarily missing due to sensor malfunction or unavailable due to limitations in practices (e.g. too difficult to obtain accurate data) and resources (e.g. too expensive to install many sensors). Meanwhile, knowing the inputs is crucial to the optimization of control strategies. For example, knowing the actual capacity of space heating/cooling supply and solar heat gain is critical to developing predictive control strategies (e.g. space preheating and cooling). This study aimed to develop a method to jointly estimate model parameters and unknown inputs. Through examples, the capability of the developed method was evaluated for estimating unknown model parameters and inputs in different scenarios with different numbers of available temperature states.

Problem statement and methodology

In practice, some inputs of a building’s thermal system can be unavailable, immeasurable, or difficult to measure for various reasons, as mentioned earlier in the Introduction section. Therefore, there is a need to estimate both unknown RC model parameters and inputs. One approach to doing so is to use available system temperature state measurements. Furthermore, it is not often that all the system states are measured and/or can be easily measured, and thus some of them may also be unknown. Accordingly, this study developed a feasible method that can be used to estimate model parameters and unknown inputs with possibly missing states, as illustrated in Figure 1. The missing states can also be simultaneously estimated by the proposed method. The development of a suitable model structure (i.e. number of temperature nodes, inputs, R, C and their connections) is beyond the scope of this study and, therefore, is not discussed in this paper. After obtaining a suitable model structure, the RC model parameters, inputs and possibly unavailable states were estimated using a UKF-based parameter-input estimation method integrated with NLS.^41–43 Specifically, the UKF was mainly used for parameter and state estimation, while the NLS was used to estimate unknown inputs. This section presents the two essential components of the developed method: (1) thermal dynamic modelling using RC models and (2) a parameter-input estimation method developed based on the integration of the UKF and NLS.

Figure 1.

General procedure of the proposed parameter-input estimation method.

Thermal dynamic modelling using RC models

In using an RC model to model the thermal dynamics of a building’s thermal zone, the temperature nodes within the zone can be represented by model nodes. Then the model nodes are connected to model parameters (i.e. Rs and Cs) and inputs (e.g. solar heat gain or mechanical heating/cooling supply). The model nodes are connected to each other through the Rs. The parameters, inputs and model nodes together create a thermal RC network. Running the RC model with inputs (i.e. simulation) can provide temperatures at the model nodes (i.e. states). The thermal dynamics at each node can be presented by ordinary differential equations (ODEs). The typical equation for a node, say $j^{\text{th}}$ , is shown in equation (1):

$\frac{{dT}_j}{dt}=\frac{1}{C_j}\left({\displaystyle \sum }_h\frac{T_h-T_j}{R_{h,j}}+{\displaystyle \sum }_j{Q}_j\right)$

(1)

where $T_j$ represents the temperature of the $j^{\text{th}}$ node, $C_j$ is the thermal capacity associated with the $j^{\text{th}}$ node, $T_h$ shows the temperature of the $h^{\text{th}}$ node adjacent to the $j^{\text{th}}$ node, $R_{h,j}$ denotes the thermal resistance between the $j^{\text{th}}$ and $h^{\text{th}}$ nodes, and $Q_j$ is the heat input to the $j^{\text{th}}$ node, if any. An RC mode can be described by a system of continuous-time linear state ODEs, which can be represented in a matrix form as a shown in equation (2):

$\dot{\mathit{s}}\left(t\right)={\mathit{A}}^c\mathit{s}\left(t\right)+{\mathit{B}}^c\mathit{u}\left(t\right)$

(2)

where $\mathit{s}$ represents a vector containing the nodal temperatures (i.e. temperature state $T_j$ ) of all nodes considered, $\mathit{u}$ denotes for system inputs vector, $\dot{\mathit{s}}$ is the first-order time derivative of the state vector $\mathit{s}$ , and ${\mathit{A}}^{\mathit{c}}$ and ${\mathit{B}}^{\mathit{c}}$ are the state and input matrices, respectively, defined by the model parameters. Note that the superscript c denotes a continuous-time model.

Parameter-input estimation method

This section describes the method for estimating unknown model parameters, inputs and unmeasured states, if any, using the UKF integrated with NLS. The model parameters were grouped to form a vector of parameters, denoted by θ . The model parameters were assumed to be time-invariant in this study, but for estimation purposes, these were modelled as a random process.^23,25 In order to adjust the method for dealing with cases in which information regarding one or more states is also unavailable/unmeasured in addition to the model parameters, the state-parameter vector $\mathit{x}$ at time t is defined by augmenting state vector s (t) (which includes all states of the thermal system) with the model parameter vector $\mathit{\theta }\left(t\right)$ , at time t, as shown in equation (3). ⁴⁴

$\mathit{x}\left(t\right)=\left[\mathit{s}\left(t\right)\hspace{0.17em}\mathit{\theta }\left(t\right)\right]$

(3)

After considering additive process noise, which represents the estimation error, equation (3) can be transformed to a state-transition equation in the state-space format as shown in equation (4):

$\begin{array}{l}\dot{\mathit{x}}\left(t\right)=\left[\dot{\mathit{s}}\left(t\right)\hspace{0.17em}\dot{\mathit{\theta }}\left(t\right)\right]\\ \dot{\mathit{x}}\left(t\right)=f^c\left(\mathit{x}\left(t\right),{\mathit{u}}^{\mathit{k}\mathit{n}}\left(t\right),{\mathit{u}}^{\mathit{u}\mathit{n}}\left(t\right)\right)+{\mathit{w}}^c\left(t\right)\end{array}$

(4)

where $f^c\left(\mathit{x}\left(t\right),{\mathit{u}}^{\mathit{k}\mathit{n}}\left(t\right),{\mathit{u}}^{\mathit{u}\mathit{n}}\left(t\right)\right)={\mathit{A}}^{\mathit{c}\mathit{x}}\mathit{x}\left(t\right)+{\mathit{B}}^{\mathit{c}\mathit{x}}\left[\begin{array}{l}{\mathit{u}}^{\mathit{k}\mathit{n}}\left(t\right)\\ {\mathit{u}}^{\mathit{u}\mathit{n}}\left(t\right)\end{array}\right]$ which $f^c$ represents the continuous-time nonlinear state function, and ${\mathit{u}}^{\mathit{k}\mathit{n}}$ and ${\mathit{u}}^{\mathit{u}\mathit{n}}$ are the known and unknown input vectors, respectively. In addition, ${\mathit{w}}^{\mathit{c}}$ represents the process noise, which was assumed to be Gaussian white noise with zero mean and covariance matrix G and ${\mathit{A}}^{\mathit{c}\mathit{x}}$ and ${\mathit{B}}^{\mathit{c}\mathit{x}}$ are the state and input matrices adjusted to the state-parameter vector as shown in equation (5):

${\mathit{A}}^{cx}={\left[\begin{array}{ll}{{\mathit{A}}^c}_{n_s\times n_s}& 0_{n_s\times n_{\theta }}\\ 0_{n_{\theta }\times n_s}& 0_{n_{\theta }\times n_{\theta }}\end{array}\right]}_{n_x\times n_x}{\mathit{B}}^{cx}={\left[\begin{array}{l}{{\mathit{B}}^c}_{n_s\times n_u}\\ 0_{n_{\theta }\times n_u}\end{array}\right]}_{n_x\times n_u}$

(5)

where $n_s$ , $n_{\theta }$ , $n_u$ and $n_x$ represent, respectively, the number of the states, model parameters, model inputs (including known and unknown) and state-parameter vector components. Equation (4) can be discretized by using the matrix exponential method, ⁴⁴ and become equation (6):

${\mathit{x}}_{k+1}=f\left({\mathit{x}}_k,{\mathit{u}}_k^{kn},{\mathit{u}}_k^{un}\right)+{\mathit{w}}_k$

(6)

where $f\left({\mathit{x}}_k,{\mathit{u}}_k^{kn},{\mathit{u}}_k^{un}\right)={\mathit{A}}^{\mathit{d}\mathit{x}}{\mathit{x}}_k+{\mathit{B}}^{\mathit{d}\mathit{x}}\left[\begin{array}{l}{\mathit{u}}_{\mathit{k}}^{\mathit{k}\mathit{n}}\\ {\mathit{u}}_{\mathit{k}}^{\mathit{u}\mathit{n}}\end{array}\right]$ , ${\mathit{x}}_{k+1}$ is the state-parameter parameter vector at time step ( $k+1$ ), and ${\mathit{u}}_{\mathit{k}}^{\mathit{k}\mathit{n}}$ and ${\mathit{u}}_{\mathit{k}}^{\mathit{u}\mathit{n}}$ are the known and unknown input vectors at time step k. ${\mathit{A}}^{dx}$ and ${\mathit{B}}^{dx}$ in this equation are, respectively, derived from ${\mathit{A}}^{cx}$ and ${\mathit{B}}^{cx}$ to produce the discrete-time state and input matrices. ${\mathit{A}}^{dx}=e^{\mathit{A}cx·∆t}$ and ${\mathit{B}}^{\mathit{d}\mathit{x}}=∆t{\mathit{A}}^{dx}{\mathit{B}}^{cx}$ , with $∆t$ being the time duration of a time step. The discretization process involved nonlinear operations amongst state, parameter and input matrices. Therefore, while the RC model retained its inherent linearity, its representation in state-space format incorporated essential nonlinear elements to effectively capture the dynamics of real-world systems.

Corresponding measurements of the system states can be linked to predictions from the thermal dynamic model through the measurement equation as shown in equation (7):

${\mathit{y}}_{k+1}=h\left({\mathit{x}}_{k+1},{\mathit{u}}_{k+1}^{kn},{\mathit{u}}_{k+1}^{un}\right)+{\mathit{v}}_{k+1}$

(7)

where $h\left({\mathit{x}}_{k+1},{\mathit{u}}_{k+1}^{kn},{\mathit{u}}_{k+1}^{un}\right)={\mathit{x}}_{m,k+1}$ , in which $\mathit{y}$ is the measurement vector, $h$ is the nonlinear measurement function whose output is the predicted state vector corresponding to available measurements, $v$ is the measurement noise modelled by Gaussian white noise with zero mean and covariance matrix Z , ${\mathit{x}}_{m,k+1}$ corresponds to the available measurements in the state-parameter vector with m representing the index of the measured states and subscript $k$ denotes the $k^{\text{th}}$ step in the discrete-time ( $k=0,1,\dots ,N-1,$ where $N$ is the total number of data samples). ⁴²

Using discretely monitored system states (i.e. node temperatures) to estimate unknown thermal dynamic model parameters and inputs necessitates an estimation method dealing with the nonlinear discrete-time state-space model. In this regard, the UKF method was selected for parameter and state estimation because it has been proven to be robust and effective for highly nonlinear problems.^43,45 Furthermore, to efficiently estimate the unknown model inputs, the NLS algorithm was integrated into UKF. The combination of UKF and NLS allowed the estimation of unknown model inputs and parameters, together with unmeasured states.^41–43

UKF is a Kalman-based technique, employing a prediction–correction two-step strategy based on unscented transformation (UT). ⁴⁶ In the prediction step, UT estimates the mean and covariance matrix of the nonlinear state function by evaluating the function using a set of deterministically selected points, known as sigma points (SPs). The unscented transformation requires a selection of $2n_x+1$ SPs, where ${\hspace{0.17em}n}_x$ is the dimension of the state-parameter vector. ⁴⁷ Accordingly, SPs were calculated using equation (8):

${\mathit{\chi }}_{k|k}^{\left(i\right)}=\{\begin{array}{ll}{\hspace{0.17em}\widehat{\mathit{x}}}_{k|k}\hspace{0.17em}& if\hspace{0.17em}i=0\\ {\widehat{\mathit{x}}}_{k|k}+{\left[{\left(\gamma \sqrt{{\widehat{\mathit{P}}}_{k|k}^{xx}}\right)}_i\right]}^T& \hspace{0.17em}where\hspace{0.17em}i=1,\dots \dots ,n_x\hspace{0.17em}\\ {\widehat{\mathit{x}}}_{k|k}-{\left[{\left(\gamma \sqrt{{\widehat{\mathit{P}}}_{k|k}^{xx}}\right)}_i\right]}^T& where\hspace{0.17em}i=n_x+1,\dots \dots ,2n_x\end{array}$

(8)

where $\widehat{\mathit{x}}$ and ${\widehat{\mathit{P}}}^{xx}$ are the mean and covariance of the state-parameter vector, respectively, $i$ denotes the ith row of the corresponding matrix and k is the kth time step. The value in front of the vertical bar $‘|$ ’ indicates the time step of the prior estimate (prediction step) while the value after the ‘ $|$ ’ indicates the time step of the posterior estimate (correction step). Furthermore, $\gamma =\sqrt{n_x+\lambda }$ and $\lambda ={\alpha }^2\left(n_x+\kappa \right)-n_x$ , in which $\alpha$ is a constant parameter related to the distribution of the SPs around the mean and $\kappa$ is a secondary scaling parameter. In this paper, α and κ were assigned with 0.01 and 0, respectively. ⁴⁶

By using the initially determined SPs (i.e. ${\mathit{\chi }}_{k|k}$ ) in, a new set of SPs (i.e. ${\mathit{\chi }}_{k+1|k}$ ) can be obtained by applying the ${\mathit{\chi }}_{k|k}$ to the nonlinear state function as given by equation (9):

${\mathit{\chi }}_{k+1|k}^{\left(i\right)}=f\left({\mathit{\chi }}_{k|k}^{\left(i\right)},{\mathit{u}}_k^{kn},{\mathit{u}}_k^{un}\right)$

(9)

With these propagated sigma points, prior estimates of the mean vector and the covariance matrix of the nonlinear state function at time step (k + 1) can be obtained. However, in cases with unknown inputs in the nonlinear state function, model input estimation becomes a prerequisite for computing a new set of SPs (i.e. ${\mathit{\chi }}_{k+1|k}$ ) based on equation (9). Therefore, it is essential to first estimate the unknown input. Using UKF alone to simultaneously estimate states, parameters and inputs can be computationally demanding. ⁴⁵ In this study, NLS was integrated into UKF to estimate the unknown model inputs. This integration reduced the number of SPs needed if UKF is used alone, thereby resulting in a more efficient estimation process. This integration not only adds efficiency to the estimation process but also improves the performance of the model.^41,43 NLS algorithm numerically solves the resulting least squares problem using gradient-based methods, such as the Levenberg–Marquardt method, also known as the damped least squares method. ⁴⁸ The Levenberg–Marquardt method interpolates between the Gauss–Newton algorithm and the gradient descent method, which results in increasing the robustness of the approach in locating the global minimum in optimization problems. ⁴⁹

In order to estimate the unknown model inputs, the developed method works as an optimization tool to minimize the estimation error (i.e. ${∆}_{k+1}$ ) which is defined as the difference between the predicted mean of the measurement vector (i.e. ${\widehat{\mathit{y}}}_{k+1|k}$ ), or in other words, the output of a measurement function $h$ (i.e. SPs ${\mathit{\chi }}_{m,k+1}^{\left(i\right)}$ ) and measurements ( ${\mathit{y}}_{k+1})$ , as shown in equation (10):

${∆}_{k+1}={\mathit{y}}_{k+1}-\hspace{0.17em}{\mathit{\chi }}_{m,k+1}^{\left(i\right)}$

(10)

Therefore, by adjusting the unknown inputs for minimizing estimation error (i.e. ${∆}_{k+1}$ ), which is the difference between the output of a measurement function and measurements, the estimated model input (i.e. ${\mathit{u}}_{k+1}^{un\left(i\right)}$ ) was calculated as given by equation (11). Accordingly, the final values of the estimated inputs equal

${\mathit{u}}_{k+1}^{un}=\frac{{\displaystyle \sum }_{i=0}^{2n_x}{\mathit{u}}_{k+1}^{un\left(i\right)}}{2n_x+1}$

(11)

SPs are defined to represent a nonlinear function. These weighting coefficients can be classified as $W_m^{\left(i\right)}$ and $W_c^{\left(i\right)}$ , where $m$ and $c$ correspond to the mean and covariance, respectively. These coefficients were determined using equation (12). Weighting coefficients were applied so they can be used to calculate the mean and covariance of the nonlinear function as shown in equations (13) and (14).

$\begin{array}{l}{W}_m^{\left(0\right)}=\frac{\lambda }{n_x+\lambda }\\ W_c^{\left(0\right)}=\frac{\lambda }{n_x+\lambda }+\left(1-{\alpha }^2+\beta \right)\\ W_m^{\left(i\right)}=W_c^{\left(i\right)}=\frac{1}{2\left(n_x+\lambda \right)},if\hspace{0.17em}i=1,\dots \dots ,2n_x\end{array}$

(12)

where $\beta$ is a factor used to emphasize the relative weight of each SPs. In this paper, β equals to 2. ⁴⁶

Once the unknown model input was estimated and SPs weighting coefficients were calculated, the prior estimates for the mean vector $\widehat{\mathit{x}}$ and the covariance matrix ${\widehat{\mathit{P}}}^{\mathit{x}\mathit{x}}$ at time step $\left(k+1\right)$ can be determined using equations (13) and (14):

${\widehat{\mathit{x}}}_{k+1|k}={\displaystyle \sum }_{i=0}^{2n_x}{W}_m^{\left(i\right)}{\mathit{\chi }}_{k+1|k}^{\left(i\right)}$

(13)

${\widehat{\mathit{P}}}_{k+1|k}^{xx}={\displaystyle \sum }_{i=0}^{2n_x}{W}_c^{\left(i\right)}\left[{\mathit{\chi }}_{k+1|k}^{\left(i\right)}\hspace{0.17em}-{\widehat{\mathit{x}}}_{k+1|k}\right]{\left[{\mathit{\chi }}_{k+1|k}^{\left(i\right)} – {\widehat{\mathit{x}}}_{k+1|k}\right]}^T+{\mathit{G}}_k$

(14)

Then, the predicted mean and covariance matrices of the measurement vector at time step k + 1 can be determined by equations (15) and (16):

${\widehat{\mathit{y}}}_{k+1|k}={\displaystyle \sum }_{i=0}^{2n_x}{W}_m^{\left(i\right)}{\widehat{\mathit{y}}}_{k+1|k}^{\left(i\right)}$

(15)

${\widehat{\mathit{P}}}_{k+1|k}^{yy}={\displaystyle \sum }_{i=0}^{2n_x}{W}_c^{\left(i\right)}\left[{\widehat{\mathit{y}}}_{k+1|k}^{\left(i\right)}-{\widehat{\mathit{y}}}_{k+1|k}\right]{\left[{\widehat{\mathit{y}}}_{k+1|k}^{\left(i\right)}\mathsf{‐}{\widehat{\mathit{y}}}_{k+1|k}\right]}^T+{\mathit{Z}}_{k+1}$

(16)

At time step $\left(k+1\right)$ , the correction step was performed by integrating the measured response $y_{k+1}$ with the posterior mean vector and covariance matrix of vector $x_{k+1}$ , as illustrated in equations (17) to (20):

${\widehat{\mathit{x}}}_{k+1|k+1}={\widehat{\mathit{x}}}_{k+1|k}+{\mathit{K}}_{k+1}\left({{\mathit{y}}_{k+1}-\widehat{\mathit{y}}}_{k+1|k}\right)$

(17)

${\widehat{\mathit{P}}}_{k+1|k+1}^{xx}={\widehat{\mathit{P}}}_{k+1|k}^{xx}-{\mathit{K}}_{k+1}\hspace{0.17em}{\widehat{\mathit{P}}}_{k+1|k}^{yy}{\mathit{K}}_{k+1}^T$

(18)

where

${\mathit{K}}_{k+1}={\widehat{\mathit{P}}}_{k+1|k}^{xy}{\left({\widehat{\mathit{P}}}_{k+1|k}^{yy}\right)}^{-1}$

(19)

${\widehat{\mathit{P}}}_{k+1|k}^{xy}={\displaystyle \sum }_{i=0}^{2n_x}{W}_c^{\left(i\right)}\left[{\mathit{\chi }}_{k+1}^{\left(i\right)}-\hspace{0.17em}{\widehat{\mathit{x}}}_{k+1|k}\right]{\left[{\widehat{\mathit{y}}}_{k+1|k}^{\left(i\right)}‐{\widehat{\mathit{y}}}_{k+1|k}\right]}^T$

(20)

in which ${\widehat{\mathit{P}}}_{k+1|k}^{xy}$ is the cross-covariance matrix between the measurement vector and the state-parameter vector, and ${\mathit{K}}_{k+1}$ is the Kalman gain.

To summarize, Figure 2 presents a flowchart for the developed method that can be used for parameter and input estimations in cases with unknown model inputs and unmeasured states.

Figure 2.

Parameter-input estimation method based on UKF integrated with NLS.

Application examples

Two case studies are presented to illustrate and evaluate the developed method: (1) one simple made-up RC model consisting of two thermal resistors and two thermal capacitances (i.e. 2R2C), for which simulated responses were considered as measurements, and (2) one real-world single-family house with noticeably more complexity (10R6C) and real-world recorded data. The capability of the method was evaluated by comparing the estimated parameters, inputs and states to the corresponding values in both study cases.

A case study with a made-up RC model

The made-up RC model, labelled as 2R2C, is shown in Figure 3. This model consists of four model parameters, two thermal resistors (i.e. R₂ and R₃ in $\raisebox{1ex}{$°\mathrm{C}$}\!\left/ \!\raisebox{-1ex}{$\text{kW}$}\right.$ ) and two thermal capacitors (i.e. C₂ and C₃ in $\raisebox{1ex}{$\text{kWh}$}\!\left/ \!\raisebox{-1ex}{$°\mathrm{C}$}\right.$ ), as well as three model inputs which is shown in Figure 4 including outdoor temperature T₁ in $℃$ , thermal load Q₁ in kW, and space heating and cooling supply Q₂ in kW.

Figure 3.

2R2C model.

Figure 4.

2R2C made-up model inputs.

The ordinary differential equations and state-space form of this RC model are shown in the equations (21) to (23), and the made-up (true) model parameters are shown in equation (25). To generate responses (i.e. states), the discrete-time linear state equations of the model were derived based on the true model parameters and then run with made-up inputs. The simulated responses (i.e. temperature states) were considered the true states in this study, and then they were contaminated by adding white Gaussian noise to generate measured states. The noise has a mean of zero and a standard deviation of 0.16 $°\mathrm{C}$ , which was based on the expected accuracy of ±0.5°C of the thermocouples used in real practice. ⁵⁰

$\begin{array}{l}{\hspace{0.17em}C}_2\frac{d{T}_2}{dt}=\frac{T_1-T_2}{R_2}+\frac{T_3-T_2}{R_3}\\ C_3\frac{d{T}_3}{dt}=\frac{T_2-T_3}{R_3}+Q_1+Q_2\end{array}$

(21)

Accordingly, the state-space equation of the 2R2C model, when Q₂ was assumed to be unknown, can be written as equation (22) and the discretized form is described by equations (23) and (24):

$\dot{\mathit{s}}\left(t\right)={\mathit{A}}^c\mathit{s}\left(t\right)+B^c\mathit{u}\left(t\right)$

$\frac{d}{dt}\left[\begin{array}{l}{T}_2\\ T_3\end{array}\right]=\left[\begin{array}{ll}\frac{-1}{R_2{C}_2}-\frac{1}{R_3{C}_2}& \frac{1}{R_3{C}_2}\\ \frac{1}{R_3{C}_3}& -\frac{1}{R_3{C}_3}\end{array}\right]\left[\begin{array}{l}{T}_2\\ T_3\end{array}\right]+\left[\begin{array}{lll}\frac{1}{R_2{C}_2}& 0& 0\\ 0& \frac{1}{C_3}& \frac{1}{C_3}\end{array}\right]\left[\begin{array}{l}{T}_1\\ Q_1\\ Q_2\end{array}\right]$

(22)

where $A^c=\left[\begin{array}{ll}\frac{-1}{R_2{C}_2}-\frac{1}{R_3{C}_2}& \frac{1}{R_3{C}_2}\\ \frac{1}{R_3{C}_3}& -\frac{1}{R_3{C}_3}\end{array}\right]$ and ${\mathit{B}}^c=\left[\begin{array}{lll}\frac{1}{R_2{C}_2}& 0& 0\\ 0& \frac{1}{C_3}& \frac{1}{C_3}\end{array}\right]$ .

$\begin{array}{l}\dot{\mathit{x}}\left(t\right)=\left[\dot{\mathit{s}}\left(t\right)\hspace{0.17em}\dot{\theta }\left(t\right)\right]\\ {\mathit{x}}_{k+1}=f\left({\mathit{x}}_k,{\mathit{u}}_k^{kn},{\mathit{u}}_k^{un}\right)+{\mathit{w}}_k\end{array}$

(23)

$f\left({\mathit{x}}_k,{\mathit{u}}_k^{kn},{\mathit{u}}_k^{un}\right)={\mathit{A}}^{\mathit{d}\mathit{x}}{\mathit{x}}_k+{\mathit{B}}^{\mathit{d}\mathit{x}}\left[\begin{array}{l}{\mathit{u}}_{\mathit{k}}^{\mathit{k}\mathit{n}}\\ {\mathit{u}}_{\mathit{k}}^{\mathit{u}\mathit{n}}\end{array}\right]$

(24)

where ${\mathit{A}}^{dx}={\left[\begin{array}{ll}{A^d}_{2\times 2}& 0_{2\times 4}\\ 0_{4\times 2}& {diag\left(1\right)}_{4\times 4}\end{array}\right]}_{6\times 6}$ and ${\mathit{B}}^{dx}={\left[\begin{array}{l}{{\mathit{B}}^d}_{2\times 3}\\ 0_{4\times 3}\end{array}\right]}_{6\times 3}$ .

Hourly data ⁵¹ (true inputs and states) were generated for 1 month (i.e. 720 h). The first 540 h of data were used for estimation and evaluation purposes, while the final 180 h of data were employed to evaluate the prediction accuracy of the calibrated model, which used the model parameters estimated at the end of the 540 h estimation period. The allocation of hours in each dataset was substantiated by prior studies^23,25 indicating the sufficiency of 540 h for capturing crucial thermal dynamics, ensuring the methodology is aligned with established practices and empirical findings. Furthermore, dedicating 180 h to prediction evaluation is in accordance with the requirements to demonstrate the model’s predictive capability effectively.

$\begin{array}{l}{\mathit{\theta }}_{\text{true}}={\left[\begin{array}{ll}\begin{array}{ll}{R}_{2\hspace{0.17em}\text{true}}& R_{3\hspace{0.17em}\text{true}}\end{array}& \begin{array}{ll}{C}_{2\hspace{0.17em}\text{true}}\hspace{0.17em}& C_{3\hspace{0.17em}\text{true}}\end{array}\end{array}\right]}^T\\ {\mathit{\theta }}_{\text{true}}={\left[\begin{array}{ll}\begin{array}{ll}3.1& 28.5\end{array}& \begin{array}{ll}2.06\hspace{0.17em}& 1.04\end{array}\end{array}\right]}^T\in {\mathbb{R}}^{4\times 1}\end{array}$

(25)

Considering the contaminated responses as measurements, one parameter-input estimation scenario (2R2C-1) and two state-parameter-input estimation scenarios (2R2C-2 and 2R2C-3) were investigated, as summarized in Table 1. All model parameters and space heating and cooling supply (Q₂) are unknown in all scenarios. In Scenarios 2 and 3, one of the states (T₂ or T₃) is unknown in addition to the unknown model parameters and input which would need to be estimated.

Table 1.

Different scenarios defined for the 2R2C model.

Scenario(C)	Available states	Unavailable states	Unknown input
2R2-1	$T_2,T_3$	None	$Q_2$
2R2-2	$T_2$	$T_3$	$Q_2$
2R2-3	$T_3$	$T_2$	$Q_2$

The state vector of this RC model is defined by $\mathit{s}={\left[\begin{array}{ll}{T}_2& T_3\end{array}\right]}^T\in {\mathbb{R}}^{2\times 1}$ , model parameters vector is $\mathit{\theta }={\left[\begin{array}{ll}\begin{array}{ll}{R}_2& R_3\end{array}& \begin{array}{ll}{C}_2& C_3\end{array}\end{array}\right]}^T\in {\mathbb{R}}^{4\times 1}$ and the model inputs vector is $\mathit{u}={\left[\begin{array}{ll}{\mathit{u}}^{\mathit{k}\mathit{n}}& {\mathit{u}}^{\mathit{u}\mathit{n}}\end{array}\right]}^T\in {\mathbb{R}}^{3\times 1}$ , in which ${\mathit{u}}^{\mathit{k}\mathit{n}}={\left[\begin{array}{ll}{T}_1& Q_1\end{array}\right]}^T$ and ${\mathit{u}}^{\mathit{u}\mathit{n}}=\left[Q_2\right]$ .The initial mean estimation of the state-parameter vector, $\widehat{\mathit{x}}$ , is presented in equation (26):

$\begin{array}{l}{\widehat{\mathit{x}}\hspace{0.17em}}_{0|0}={\left[\begin{array}{ll}{\widehat{\mathit{s}}}_{0|0}& {\widehat{\mathit{\theta }}}_{0|0}\end{array}\right]}^T\in {\mathbb{R}}^{6\times 1}\\ {\widehat{\mathit{s}}}_{0|0}={\left[\begin{array}{ll}21& 30\end{array}\right]}^T\\ {\widehat{\mathit{\theta }}}_{0|0}={\left[\begin{array}{ll}\begin{array}{ll}4.45& 26.35\end{array}& \begin{array}{ll}2.64& 1.2\end{array}\end{array}\right]}^T\end{array}$

(26)

Since the states and parameters of the RC model were assumed to be uncorrelated, the initial covariance matrix ${\widehat{\mathit{P}}}_{\mathbf{0}|\mathbf{0}}^{\mathit{x}\mathit{x}}$ is diagonal. The diagonal terms of the initial covariance matrix were determined by ${{(\widehat{\mathit{P}}}_{0|0}^{xx})}_{\text{states}}={\left[1.0\times \left({\widehat{\mathit{s}}}_{0|0}\right)\right]}^2$ for the states $\left(\mathit{s}\right)$ . In addition, the diagonal terms of the initial covariance matrix for the model parameters $\left(\mathit{\theta }\right)$ were defined using ${{(\widehat{\mathit{P}}}_{0|0}^{xx})}_{\text{parameters}}={\left[0.05\times \left({\widehat{\mathit{\theta }}}_{0|0}\right)\right]}^2$ . The measurement noise covariance matrix is diagonal, time-invariant and equal to $\mathit{Z}={\left(0.3°\mathrm{C}\right)}^2$ . Similarly, the process noise covariance matrix is diagonal, time-invariant and equal to $\mathit{G}={\left[1\times {10}^{-5}\times {\widehat{\mathit{x}}}_{0|0}\right]}^2.$ The values for initial covariance matrices and process noise covariance matrix were calculated based on Refs. 44 and 45.

Figure 5 presents the time histories of model parameter estimation for Scenario 2R2C-1, where all the states were measured. In these figures, the dashed black line shows true model parameter values, and the red line shows the mean estimate. The shaded grey area shows the mean +/− two standard deviations obtained from the covariance matrix ${\widehat{\mathit{P}}}_{k|k}^{xx}$ . As shown in Figure 5(a), the thermal resistance parameters R₂ and R₃ were estimated accurately and converged near their true values. Similar estimation results for the thermal capacitances C₂ and C₃ are shown in Figure 5(b). The estimated model parameters’ standard deviation is higher at the beginning of the estimation history, indicating a higher level of uncertainty in these values initially. While the estimated model parameters may be sufficiently accurate, they may not reach their true values. This difference is due to the presence of nonlinearity in the system and measurement noise. This nonlinearity can introduce complexities in the input–output relationship, leading to potential numerical errors in the estimation results.^52,53 It should be emphasized that the model parameter estimation results for Scenarios 2R2C-2 and 2R2C-3 are also converged but are not shown here for brevity. The final estimated model parameters for these scenarios can be found in Table 2.

Figure 5.

Time histories of model parameter estimations for scenario 2R2C-1: (a) thermal resistance and (b) thermal capacitances.

Table 2.

Last estimated model parameters for scenarios 2R2C-2 and 2R2C-3.

Model parameters	True values	2R2C-2	2R2C-3
R2	$4.45$	2.85	3.45
R3	$26.35$	25.1	25.23
C2	$2.64$	2.3	2.63
C3	1.2	1.2	1.14

Figure 6 presents the estimation history of the unknown model input $Q_2$ compared to its true value in Scenario 2R2C-1. To better illustrate the estimation results, the 540 h of time histories were divided into three 180-h segments, which is shown in Figure 6. These results reaffirm that the developed method can accurately estimate the unknown model input. Furthermore, to provide a quantitative comparison, mean absolute percentage error (MAPE), a measure of precision for constructing fitted time series values when trend estimation is involved, ⁵⁴ was used in this study, as defined by equation (27):

$\text{MAPE}=\left(\frac{{\displaystyle \sum }_{i=1}^n\left|\frac{{True}_i\hspace{0.17em}-{Calculated}_i\hspace{0.17em}}{{True}_i}\right|}{n}\right)\times 100(\%)$

(27)

where True represents the true value, Calculated is the predicted or estimated value and n shows the number of data points. MAPE of model input estimation in Scenario 2R2C-1 was 1.1%, indicating that the developed method can accurately estimate the model input.

Figure 6.

Comparison between the true and estimated model input for scenario 2R2C-1.

Using the developed method, parameters, input and states were estimated at each time step of the estimation period. Figure 7 compares the estimated and the true profiles of the two states of the 2R2C model (i.e. T₂ and T₃) in Scenario 2R2C-1. The MAPE of the estimated state and the true one for both temperatures were 0.4% and 0.82% for T₂ and T₃, respectively, which indicates excellent agreement between the estimated results and the corresponding true states. The accuracy of the estimated T₃ is not as high as that for T₂. This is due to the unknown input, Q₂, in the estimation period being linked to T₃, as shown in Figure 3. Figure 7 also shows that the developed method is able to estimate the parameters, inputs and states relatively accurately within the first few time steps. This observation indicates the suitability of the method for online estimations.

Figure 7.

Comparison between true and estimated temperature responses for Scenario 2R2C-1.

The MAPE value for the estimated model input and states for all scenarios is shown in Table 3. The MAPE for the estimated model input in Scenarios 2R2C-2 and 2R2C-3 was slightly higher than that of Scenario 2R2C-1. Additionally, the MAPE of the estimated states for Scenarios 2R2C-2 and 2R2C-3 was also higher but less than 2.5%. These observations reaffirm that more known measurements can result in a more accurate estimation.

Table 3.

MAPE of estimated input and states for all scenarios of 2R2C.

MAPE	2R2C-1	2R2C-2 ( $T_2$ available)	2R2C-3 ( $T_3$ available)
Q 2	1.1%	2.12%	1.36%
T 2	0.4%	1.83%	2.16%
T 3	0.82%	2.32%	1.26%

Using the parameters obtained at the end of the 540-h estimation period, the RC models were used to simulate the states in the prediction period (with known inputs). The predicted temperature responses were generated at each time step, hourly from hours 541 to 720. Figure 8 compares the predicted temperature responses for Scenario 2R2C-1. Both temperature responses T₂ and T₃ were predicted with good accuracy, with the MAPE values being 0.92% and 2.55%, respectively. Again, the accuracy of the predicted response for T₃ was not as high as that for T₂, due to the unknown input, Q₂, in the estimation period linked to T₃, as shown in Figure 3. This comparison indicates that the inaccuracy in the parameter estimation would not result in a significant error in state prediction; therefore, the developed method is sufficient for the analysis of building thermal dynamics when all states are measured.

Figure 8.

Comparison between true and predicted temperature responses for Scenario 2R2C-1.

Furthermore, the MAPE values calculated for both states in Scenarios 2R2C-2 and 2R2C-3 were found to be less than 5%, as shown in Table 4. These values indicate that the prediction results are acceptable, even when certain state measurements are unavailable. Figures 9–12 show the estimated profile of Scenarios 2R2C-2 and 2R2C-3, scenarios with a limited number of temperature responses.

Table 4.

MAPE of predicted states for all scenarios of the 2R2C model.

MAPE	2R2C-1 (%)	2R2C-2 (%)	2R2C-3 (%)
T 2	0.92	3.91	4.41
T 3	2.55	4.92	3.05

Figure 9.

Comparison between true and estimated model input for scenario 2R2C-2.

Figure 10.

Comparison between true and estimated temperature responses for scenario 2R2C-2.

Figure 11.

Comparison between true and estimated model input for scenario 2R2C-3.

Figure 12.

Comparison between true and estimated temperature responses for scenario 2R2C-3.

A case study with a real-world house

The developed method was tested in a real-world low-energy wooden-frame house in Eastman, Quebec, Canada, as shown in Figure 13. ²⁶ The study used 720 h of measured data. The first 540 h of data, organized on an hourly basis, were designated as the estimation dataset, and the remainder as the prediction dataset. The allocation of hours in each dataset was supported by previous studies,^23,25 which demonstrated the adequacy of 540 h for capturing essential thermal dynamics. This ensures that our methodology is consistent with established practices and empirical findings. Various scenarios were defined based on the availability of different state measurements. Accordingly, the effectiveness of the developed method was evaluated by (1) comparing the estimated values of states and input with their measured (true) values in the estimation dataset and (2) comparing the predicted responses obtained using the final estimated model parameters with the measured responses in the prediction dataset.

Figure 13.

A single detached house as a case study. ²⁶

The RC model structure of the detached single-family house shown in Figure 14 was developed by Wang et al. ²⁶ In this RC model, the basement, second floor and main floor were labelled as b, u and m, respectively. This RC model was called 10R6C, as it consists of 10 thermal resistances and 6 thermal capacitances, T _m , and T _b in °C were the zone temperatures of the second floor, main floor and basement, respectively. T _um , T _mm and T _bm were the internal thermal mass temperature responses, which were not measured. The ordinary differential equations and state-space form of this RC model are shown in equations (28) and (29).

$\begin{array}{l}{C}_{ui}\frac{d{T}_u}{dt}=\frac{T_0-T_u}{R_3}+\frac{T_{um}-T_u}{R_1}+\frac{T_{um}-T_u}{R_2}+{\alpha }_{hpu}{Q}_{hp}+{\alpha }_{su}{Q}_s\\ C_{um}\frac{d{T}_{um}}{dt}=\frac{T_u-T_{um}}{R_1}\\ C_{mi}\frac{d{T}_m}{dt}=\frac{T_0-T_m}{R_6}+\frac{T_{mm}-T_m}{R_4}+\frac{T_b-T_m}{R_5}+\frac{T_u-T_m}{R_2}+{\alpha }_{hpm}{Q}_{hp}+{\alpha }_{sm}{Q}_s+{\alpha }_{em}{Q}_{elec}\\ C_{mm}\frac{d{T}_{mm}}{dt}=\frac{T_m-T_{mm}}{R_4}+\frac{T_b-T_{mm}}{R_7}\\ C_{bi}\frac{d{T}_b}{dt}=\frac{T_0-T_b}{R_9}+\frac{T_{bm}-T_b}{R_8}+\frac{T_{mm}-T_b}{R_7}+\frac{T_m-T_b}{R_5}+{\alpha }_{hpb}{Q}_{hp}+{\alpha }_{sb}{Q}_s+{\alpha }_{eb}{Q}_{elec}\\ C_{bm}\frac{d{T}_{bm}}{dt}=\frac{T_b-T_{bm}}{R_8}+\frac{T_g-T_{bm}}{R_{10}}+Q_{vcs}\end{array}$

(28)

Figure 14.

RC model structure for the single detached house.

Therefore, the state-space form of the 10R6C model can be written as follows:

$\frac{d}{dt}\left[\begin{array}{l}{T}_u\\ T_{um}\\ T_m\\ T_{mm}\\ T_b\\ T_{bm}\end{array}\right]={\mathit{A}}^{\mathit{c}}\left[\begin{array}{l}{T}_u\\ T_{um}\\ T_m\\ T_{mm}\\ T_b\\ T_{bm}\end{array}\right]+{\mathit{B}}^{\mathit{c}}\left[\begin{array}{l}{T}_0\\ T_{\mathrm{g}}\\ \begin{array}{l}{Q}_{\mathrm{s}}\\ Q_{\text{hp}}\end{array}\\ Q_{\text{elec}}\\ Q_{\text{vcs}}\end{array}\right]$

(29)

where

${\mathit{A}}^{\mathit{c}}=\left[\begin{array}{cccccc}\frac{-1}{C_{ui}}\left(\frac{1}{R_1}+\frac{1}{R_2}+\frac{1}{R_3}\right)& \frac{1}{R_1*C_{ui}}& \frac{1}{R_2*C_{ui}}& 0& 0& 0\\ \frac{1}{R_1*C_{um}}& \frac{-1}{R_1*C_{um}}& 0& 0& 0& 0\\ \frac{1}{R_2*C_{mi}}& 0& \frac{-1}{C_{mi}}\left(\frac{1}{R_2}+\frac{1}{R_4}+\frac{1}{R_5}+\frac{1}{R_6}\right)& \frac{1}{R_4*C_{mi}}& \frac{1}{R_5*C_{mi}}& 0\\ 0& 0& \frac{1}{R_4*C_{mm}}& \frac{-1}{C_{mm}}\left(\frac{1}{R_4}+\frac{1}{R_7}\right)& \frac{1}{R_7*C_{mm}}& 0\\ 0& 0& \frac{1}{R_5*C_{bi}}& \frac{1}{R_7*C_{bi}}& \frac{-1}{C_{bi}}\left(\frac{1}{R_5}+\frac{1}{R_7}+\frac{1}{R_8}+\frac{1}{R_9}\right)& \frac{1}{R_8*C_{bi}}\\ 0& 0& 0& 0& \frac{1}{R_8*C_{bm}}& \frac{-1}{C_{bm}}\left(\frac{1}{R_8}+\frac{1}{R_{10}}\right)\end{array}\right]$

${\mathit{B}}^{\mathit{c}}=\left[\begin{array}{llllll}\frac{1}{R_3*C_{\text{ui}}}& 0& \frac{{\alpha }_{\text{su}}}{C_{\text{ui}}}& \frac{{\alpha }_{\text{hpu}}}{C_{\text{ui}}}& 0& 0\\ 0& 0& 0& 0& 0& 0\\ \frac{1}{R_6*C_{\text{mi}}}& 0& \frac{{\alpha }_{\text{sm}}}{C_{\text{mi}}}& \frac{{\alpha }_{\text{hpm}}}{C_{\text{mi}}}& \frac{{\alpha }_{\text{em}}}{C_{\text{mi}}}& 0\\ 0& 0& 0& 0& 0& 0\\ \frac{1}{R_9*C_{\text{bi}}}& 0& \frac{{\alpha }_{\text{sb}}}{C_{\text{bi}}}& \frac{{\alpha }_{\text{hpb}}}{C_{\text{bi}}}& \frac{{\alpha }_{\text{eb}}}{C_{\text{bi}}}& 0\\ 0& \frac{1}{R_{10}*C_{\text{bm}}}& 0& 0& 0& \frac{1}{C_{bm}}\end{array}\right]$

${\mathit{A}}^{dx}={\left[\begin{array}{ll}{{\mathit{A}}^d}_{6\times 6}& 0_{6\times 24}\\ 0_{24\times 6}& {diag\left(1\right)}_{24\times 24}\end{array}\right]}_{30\times 30}{\mathit{B}}^{dx}={\left[\begin{array}{l}{{\mathit{B}}^d}_{6\times 6}\\ 0_{24\times 6}\end{array}\right]}_{30\times 6}$

The hourly model inputs for this RC thermal model are shown in Figure 15. The model inputs include (1) global irradiation Q _s in kW per unit area on the south facade, which was used to obtain approximate effective solar heat gains weighted by solar gain factors ( $\alpha$ _s ), and (2) heating supply from the geothermal heat pump Q _hp in kW, which was, considered as the unknown model input. This heating supply was distributed across three levels using the distribution factor ${\alpha }_{hp}$ . (3) T _g was the ground temperature, which was approximately 13°C constant. (4) T₀ was the outdoor temperature in °C. (5) Q _elec was the gross electricity demand in kW, which was used to calculate internal heat gains weighted according to internal gain factors ( $\alpha$ _e ). (6) Q _vcs in kW was the thermal energy charged to the slab located in the basement.

Figure 15.

10R6C model inputs.

The state vector of this RC model can be defined by $\mathit{s}={\left[T_u,T_{um},T_m,T_{mm},T_b,T_{bm}\right]}^T\in {\mathbb{R}}^{6\times 1}$ , model inputs vector $\mathit{u}={\left[\begin{array}{ll}{\mathit{u}}^{kn}& {\mathit{u}}^{un}\end{array}\right]}^T\in {\mathbb{R}}^{6\times 1}$ , known input vector ${\mathit{u}}^{kn}={\left[T_0,T_{\mathrm{g}},Q_{\mathrm{s}},Q_{\text{elec}},Q_{\text{vcs}}\right]}^T$ and unknown input vector ${\mathit{u}}^{un}$ = $\left[Q_{hp}\right]$ . The state-parameter vector can also be expressed as equation (30), where s is the previously described state vector and θ is the unknown model parameters defined by equation (31):

$\mathit{x}={\left[\begin{array}{ll}\mathit{s}& \mathit{\theta }\end{array}\right]}^T\in {\mathbb{R}}^{30\times 1},$

(30)

$\mathit{\theta }={\left[\begin{array}{lllll}{R}^T& C^T& {\alpha }_s^T& {\alpha }_e^T& {\alpha }_{hp}^T\end{array}\right]}^T\in {\mathbb{R}}^{24\times 1},$

(31)

where $\mathit{R}={\left[\begin{array}{llllllllll}{R}_1& R_2& R_3& R_4& R_5& R_6& R_7& R_8& R_9& R_{10}\end{array}\right]}^T\in {\mathbb{R}}^{10\times 1}$ , in $\raisebox{1ex}{$°\mathrm{C}$}\!\left/ \!\raisebox{-1ex}{$\text{kW}$}\right.$ , and $\mathit{C}={\left[\begin{array}{llllll}{C}_{ui}& C_{um}& C_{mi}& C_{mm}& C_{bi}& C_{bm}\end{array}\right]}^T\in {\mathbb{R}}^{6\times 1}$ , in $\raisebox{1ex}{$\text{kWh}$}\!\left/ \!\raisebox{-1ex}{$°\mathrm{C}$}\right.$ . Additionally, the solar gain factors are ${\mathit{\alpha }}_{\mathit{s}}={\left[\begin{array}{lll}{\alpha }_{su}& {\alpha }_{sm}& {\alpha }_{sb}\end{array}\right]}^T\in {\mathbb{R}}^{3\times 1}$ , in ${\mathrm{m}}^2$ , and the internal gain factors are ${\mathit{\alpha }}_{\mathit{e}}={\left[\begin{array}{ll}{\alpha }_{em}& {\alpha }_{eb}\end{array}\right]}^T\in {\mathbb{R}}^{2\times 1}$ . The heating power distribution factor for the second floor, main floor and basement is ${\mathit{\alpha }}_{\mathit{h}\mathit{p}}={\left[\begin{array}{lll}{\alpha }_{hpu}& {\alpha }_{hpm}& {\alpha }_{hpb}\end{array}\right]}^T\in {\mathbb{R}}^{3\times 1}$ , respectively, which should satisfy equation (32):

${\displaystyle \sum }_{i=u,m,b}{\alpha }_{hpi}=1$

(32)

To evaluate the performance of the developed method, this study considered various scenarios with different available measured temperature responses, as shown in Table 5. Scenario 10R6C-1 produced all three room temperature state measurements ( $T_u,T_m,T_b$ ). Interior thermal mass temperature states ( $T_{um},T_{mm},T_{bm}$ ) were not available in all scenarios.

Table 5.

Different scenarios for 10R6C model.

Scenarios (C)	Available room temperature states	Unavailable room temperature states	Unknown input
10R6-1	$T_u,T_m,T_b$	None	$Q_{hp}$
10R6-2	$T_u,T_m$	$T_b$	$Q_{hp}$
10R6-3	$T_u,T_b$	$T_m$	$Q_{hp}$
10R6-4	$T_m,T_b$	$T_u$	$Q_{hp}$
10R6-5	$T_u$	$T_m,T_b$	$Q_{hp}$
10R6-6	$T_m$	$T_u,T_b$	$Q_{hp}$
10R6-7	$T_b$	$T_u,T_m$	$Q_{hp}$

In the estimation phase, the initial mean estimates for the state-parameter vector can be expressed as equation (33):

${\widehat{\mathit{x}}}_{0|0}={\left[\begin{array}{ll}{\widehat{\mathit{s}}}_{0|0}& {\widehat{\mathit{\theta }}}_{0|0}\end{array}\right]}^T$

(33)

where

${\widehat{\mathit{s}}}_{0|0}=\left[\begin{array}{llllll}22& 21.3& 21.1& 22.5& 21.2& 20.5\end{array}\right]\left[\text{unit}:℃\right]$

(34)

and equation(34) defined ${\widehat{\mathit{\theta }}}_{0|0}$ :

${\widehat{\mathit{\theta }}}_{0|0}={\left[\begin{array}{lllll}{\widehat{\mathit{R}}}_{0|0}^T& {\widehat{\mathit{C}}}_{0|0}^T& {\widehat{\alpha }}_{s\hspace{0.17em}0|0}^T& {\widehat{\alpha }}_{e\hspace{0.17em}0|0}^T& {\widehat{\mathit{\alpha }}}_{hp\hspace{0.17em}0|0}^T\end{array}\right]}^T$

${\widehat{\mathit{R}}}_{0|0}=\left[\begin{array}{llllllllll}1.36& 1.69& 20.7& 2.4& 4.98& 35.25& 12.75& 2.49& 29.59& 16.13\end{array}\right]\left[\text{unit}:\raisebox{1ex}{$°\mathrm{C}$}\!\left/ \!\raisebox{-1ex}{$\text{kW}$}\right.\right]$

${\widehat{\mathit{C}}}_{0|0}=\left[\begin{array}{llllll}1.42& 2.67& 2.60& 25.67& 3.74& 9.06\end{array}\right]\hspace{0.17em}\left[\text{unit}:\raisebox{1ex}{$\text{kWh}$}\!\left/ \!\raisebox{-1ex}{$°\mathrm{C}$}\right.\right]$

${\widehat{\mathit{\alpha }}}_{s\hspace{0.17em}0|0}={\left[\begin{array}{lll}1.83& 2.93& 1.51\end{array}\right]}^T\hspace{0.17em}\left[\text{unit}:{\mathrm{m}}^2\right]$

${\widehat{\mathit{\alpha }}}_{e\hspace{0.17em}0|0}={\left[\begin{array}{ll}0.50& 0.14\end{array}\right]}^T\hspace{0.17em}\left[\text{unit}:--\right]$

${\widehat{\mathit{a}}}_{hp\hspace{0.17em}0|0}={\left[\begin{array}{lll}0.25& 0.38& 0.35\end{array}\right]}^T\hspace{0.17em}\left[\text{unit}:--\right]$

As explained in the case study of a simple RC model, the diagonal terms for the states $\left(\mathit{s}\right)$ of the initial covariance matrix were selected as ${{(\widehat{\mathit{P}}}_{0|0}^{\text{xx}})}_{states}={\left[1.0\times \left({\widehat{\mathit{s}}}_{0|0}\right)\right]}^2$ , ${{(\widehat{\mathit{P}}}_{0|0}^{\text{xx}})}_{\text{parameters}}={\left[0.05\times \left({\widehat{\mathit{\theta }}}_{0|0}\right)\right]}^2$ , $\mathit{R}={\left(0.3°\mathrm{C}\right)}^2$ and $\mathit{G}={\left[1\times {10}^{-5}\times {\widehat{\mathit{x}}}_{0|0}\right]}^2$ , respectively.

The results of the estimation of model parameters, including thermal resistances ( R ), thermal capacitances ( C ), solar gain factors ( ${\mathit{\alpha }}_s$ ), internal gain factors ( ${\mathit{\alpha }}_e$ ) and heating power distribution factors ( ${\mathit{\alpha }}_{hp}$ ), for Scenario 10R6C-1, are presented in Figures 16–18. These results show convergence in some of the thermal resistances and capacitances. No clear convergence trend was observed for the solar gain factors, internal gain factors and heating power distribution factors. In general, the variations in the values were small. The variations over the estimated time were due to the constant changes in building conditions, such as closed-open windows or blinds, on-off ventilation systems and moving furniture.

Figure 16.

10R6C model parameters estimation: thermal resistances for scenario 10R6C-1.

Figure 17.

10R6C model parameters estimation: thermal capacitances for scenario 10R6C-1.

Figure 18.

10R6C model parameters estimation: solar and internal gain factors for scenario 10R6C-1.

The performance of the developed method in estimating the model input is demonstrated in Figure 19, which compares the estimated model input Q _hp to its true value in Scenario 10R6C-1. To illustrate the estimation results more clearly, the 540 h of time histories were divided into three 180-h segments, which is shown in Figure 19. The MAPE of estimated input was about 10%. Notably, there were specific time periods where Q_hp experienced both overestimations and underestimations. A significant contributor to these fluctuations in input estimation can be the global irradiation on the south facade of the house. Elevated levels of irradiation, particularly during the early morning and late afternoon when the sun was not positioned directly in front of the house, tend to result in overestimation. Consequently, these overestimations may introduce errors in the model’s parameter estimation during these specific time periods.

Figure 19.

Comparison between true and estimated input for Scenario 10R6C-1.

Figure 20 compares the estimated and true temperature response histories for Scenario 10R6C-1. The MAPE of the estimated states (T _u , T _m and T _b ) and the corresponding measured (true) values was less than 2% for the three available measurements, indicating a strong agreement between the estimated and the true states. The maximum difference was approximately 1°C. These comparisons demonstrate that the developed method can be effectively applied to real-time temperature prediction for real-world practices.

Figure 20.

Comparison between true and estimated temperature responses for Scenario 10R6C-1.

Table 6 displays estimation results for the unknown input for all scenarios. Since the maximum MAPE difference between Scenario 10R6C-1 and the other scenarios was only about 4%, it can be concluded that the developed method can estimate the unknown model input with acceptable accuracy even when there are limited available measurements. Additionally, Table 6 displays the MAPEs of state estimation of temperature responses and their true values for the seven scenarios. The difference between Scenario 10R6C-1 with three available measurements and scenarios with only one measurement was approximately 2%, which proves that when sparse measurements (one or two measurements) are available and included in the estimation process, the state estimation can still be accurate. This observation illustrated the robustness and excellent agreement of both measured and unmeasured responses. Additionally, this verified that the developed method is suitable for application in buildings where only a limited number of temperature responses can be measured.

Table 6.

MAPE of the estimated model input and states for all scenarios of 10R6C.

MAPE of estimation	Scenarios
	10R6C-1	10R6C-2	10R6C-3	10R6C-4	10R6C-5	10R6C-6	10R6C-7
Qhp	10.38%	12.13%	11.74%	12.36%	14.32%	13.42%	14.03%
T u	0.82%	1.54%	1.62%	2.12%	3.26%	2.75%	2.95%
T m	1.13%	1.32%	1.73%	1.45%	3.54%	2.68%	2.74%
T b	1.58%	1.91%	1.23%	2.01%	3.38%	2.94%	2.65%

To show the performance of the calibrated RC model for temperature state prediction, Figure 21 compares the measured temperature responses from the prediction dataset to those predicted responses using the parameters obtained at the end of the 540-h estimation period. The RC models were used to simulate the states in the prediction period. Predicted temperature responses were computed at each hourly time step within the hour 541 to 720. The MAPE in the predicted state and measured temperature responses were around 2.5% for each available measurement, indicating that the measured and predicted temperatures are in excellent agreement. This also indicates that the error in estimating the input is acceptable because it does not cause a significant error in the parameter estimation and state prediction. Table 7 presents the MAPEs of the measured temperatures and predicted states based on the final estimated model parameters derived from the estimation dataset for all scenarios. All MAPE values were less than 5.2%. The accurate predictions of future states reflected in Figure 21 and Table 7 indicate that the model parameters were effectively estimated even with a limited number of missing state measurements.

Figure 21.

Comparison between true and predicted temperature responses for Scenario 10R6C-1.

Table 7.

MAPE of the predicted responses for the 10R6C model in different scenarios.

MAPE of prediction	Scenarios
	10R6C-1	10R6C-2	10R6C-3	10R6C-4	10R6C-5	10R6C-6	10R6C-7
T u	2.74%	4.14%	4.11%	4.46%	4.69%	4.87%	4.79%
T m	2.18%	3.92%	4.21%	4.12%	4.76%	4.36%	5.2%
T b	2.53%	4.18%	3.82%	4.25%	5.12%	5.03%	4.55%

RC model parameters play crucial roles, as they can be used to interpret a building’s thermal properties ²⁶ and be used in RC thermal models to predict future states. Table 8 summarizes the estimated model parameters at the final time step for all seven scenarios. It can be seen that even with certain missing state measurements, the final estimated parameter values for each scenario are relatively consistent. The estimation profiles of Scenario 10R6C-3 selected as an example of a scenario with a limited number of temperature responses are shown in Figures 22–27, other scenarios are not shown for brevity.

Table 8.

Final estimated parameters in different scenarios for the 10R6C model.

	Scenarios							Max error (%) compared to 10R6C-1
Parameters	10R6C-1	10R6C-2	10R6C-3	10R6C-4	10R6C-5	10R6C-6	10R6C-7
$R_1$	1.32	1.27	1.30	1.38	1.39	1.38	1.36	5
$R_2$	1.59	1.57	1.76	1.63	1.72	1.65	1.68	11
$R_3$	19.06	18.52	17.72	14.48	14.74	12.87	18.13	32
$R_4$	2.78	2.71	2.67	2.79	2.53	2.66	2.42	13
$R_5$	5.26	5.06	5.44	5.51	4.98	5.02	5.06	5
$R_6$	26.35	18.96	30.41	29.07	30.31	27.51	31.51	28
$R_7$	12.49	12.69	12.34	12.49	12.76	12.77	12.9	3
$R_8$	2.42	2.48	2.39	2.42	2.48	2.49	2.47	3
$R_9$	20.62	23.18	20.55	20.73	26.92	25.62	26.31	31
$R_{10}$	13.72	14.37	13.44	13.70	15.20	14.78	14.81	11
$C_{ui}$	1.39	1.40	1.39	1.35	1.35	1.36	1.42	3
$C_{um}$	2.37	2.46	2.21	2.47	2.44	2.57	2.67	13
$C_{mi}$	2.16	1.95	2.52	2.35	2.46	2.21	2.60	20
$C_{mm}$	25.96	25.32	25.99	24.87	26.56	25.86	26.47	4
$C_{bi}$	3.74	3.66	3.67	3.83	3.71	3.66	3.52	6
$C_{bm}$	9.62	9.36	9.92	9.68	9.11	9.03	9.08	6
${\alpha }_{su}$	1.97	1.78	2.06	2.01	1.80	1.80	1.75	11
${\alpha }_{sm}$	3.05	2.83	3.05	3.26	2.81	2.83	2.69	12
${\alpha }_{sb}$	1.22	1.43	1.13	1.14	1.46	1.46	1.38	20
${\alpha }_{em}$	0.37	0.39	0.43	0.4	0.45	0.42	0.45	22
${\alpha }_{eb}$	0.13	0.13	0.14	0.14	0.14	0.14	0.13	8
${\alpha }_{hpb}$	0.35	0.32	0.34	0.35	0.34	0.33	0.36	9
${\alpha }_{hpm}$	0.40	0.40	0.37	0.39	0.38	0.40	0.36	10
${\alpha }_{hpu}$	0.23	0.25	0.26	0.24	0.26	0.25	0.25	13

Figure 22.

10R6C model parameters estimation: thermal resistances for scenario 10R6C-3.

Figure 23.

10R6C model parameters estimation: thermal capacitances for scenario 10R6C-3.

Figure 24.

10R6C model parameters estimation: solar and internal gain factors for scenario 10R6C-3.

Figure 25.

Comparison between true and estimated model input for scenario 10R6C-3.

Figure 26.

Comparison between true and estimated temperature responses for scenario 10R6C-3.

Figure 27.

Comparison between true and predicted temperature responses for scenario 10R6C-3.

Conclusion

Previous research focused on developing methods for model parameter estimation with all model inputs available. This study developed an effective method for estimating RC model parameters and inputs with possible unmeasured temperature states. In this regard, the thermal resistor–capacitor network (RC) was used for thermal dynamic modelling and the unscented Kalman filter (UKF) was integrated with the nonlinear least squares (NLS) method to estimate unknown model parameters and inputs. In order to show the applicability of the developed method, two case studies were conducted: one with made-up data and the other with real-world data. Scenarios with different numbers of available temperature measurements were created to evaluate the capability and performance of the developed method in different circumstances.

The estimated model input and states were compared with their true values. The estimated model inputs in all scenarios of the case study with made-up data have a MAPE (mean absolute percentage error) less than 2.5%, while less than 14.5% in the scenarios of the case study with real-world data. Furthermore, comparing model predictions with the true measurements in all scenarios of both case studies shows a maximum MAPE of approximately 5%. These low MAPE values prove that the developed method can effectively estimate unknown model parameters, input and possible missing states.