Abstract

Latent thermal energy storage (TES) devices could enable advances in many thermal management applications, including peak load shifting for reducing energy demand and cost of HVAC or providing supplemental heat rejection in transient thermal management systems. However, real-time feedback control of such devices is currently limited by the absence of suitable state of charge estimation techniques, given the nonlinearities associated with phase change dynamics. In this paper, we design and experimentally validate a state-dependent Riccati equation (SDRE) filter for state of charge estimation in a phase change material (PCM)-based TES device integrated into a single-phase thermal-fluid loop. The advantage of the SDRE filter is that it does not require linearization of the nonlinear finite volume model; instead, it uses a linear parameter-varying system model which can be quickly derived using graph-based methods. We leverage graph-based methods to prove that the system model is uniformly detectable, guaranteeing that the state estimates are bounded. Using measurements from five thermocouples embedded in the PCM of the TES and two thermocouples measuring the fluid temperature at the inlet and outlet of the device, the state estimator uses a reduced-order finite volume model to determine the temperature distribution inside the PCM and in turn, the state of charge of the device. We demonstrate the state estimator in simulation and on experimental data collected from a thermal management system testbed to show that the state estimation error converges near zero and remains bounded.

1 Introduction

There is growing interest in the integration of thermal energy storage (TES) devices with a variety of thermal-fluid systems for improving performance. However, new control strategies are needed to take full advantage of the potential benefits of TES devices in fast-timescale transient thermal management systems (TMSs) [15]. Control strategies for many forms of TMSs with integrated latent thermal storage, or hybrid thermal management systems, have been studied in both experimental and simulation environments. In Ref. [1], Shanks et al. demonstrate a logic-based controller for state of charge management in a simulated single-phase hybrid TMS for aircraft electronics cooling. In Ref. [2], Shafiei and Alleyne develop a model predictive controller in simulation for a two-phase hybrid TMS in a transport refrigeration system. Similarly, Pangborn et al. develop a model predictive controller for a vehicle TMS with distributed TES and demonstrate it in simulation [3]. In each of these examples, state estimation and measurement are neglected because the controllers are incorporated into the simulation model, and all system states are known. However, in practice, real-time controllers for TES require accurate state estimation, including real-time estimation of the internal temperatures and the state of charge (SOC) [4,5]. Experimental testing of control strategies for fast-timescale hybrid TMS, such as vehicle TMS, is limited; therefore, development and experimental validation of real-time model-based state estimation techniques for latent TES devices will contribute to filling this gap in the literature.

For latent TES containing phase-change materials (PCMs), the SOC is typically defined as a function of either (i) the fraction of PCM in either the liquid or solid phase (phase fraction) or (ii) the amount of energy stored [5]. Methods for direct measurement of the phase fraction definition have been developed for various types of TES [5,6], but direct measurement of the stored energy is difficult when the stored energy can be in the form of both latent and sensible heat. Instead, methods for determining this SOC definition must rely on indirect measurements and model-based estimators [7]. Dynamic state estimators such as the Kalman filter or Luenberger observer are commonly used to estimate unmeasured or unmeasurable states in dynamic systems given limited measurements [8,9]. In latent TES systems, temperature and SOC estimation is particularly difficult because the temperature of the PCM is nearly constant during the phase change process when the PCM exchanges most of its energy. In addition, the nonlinear or hybrid dynamics associated with the phase change increase the complexity of the dynamic model, making observability, stability, and convergence difficult to guarantee in advance.

1.1 Related Work.

In experimental environments, state of charge measurement or estimation methods vary extensively and are often application-specific. Zsembinszki et al. present an overview of various measurement methods for determining the phase fraction [5]. These include (i) displacement or level sensors for measuring the volume change as the PCM melts and solidifies [10], (ii) pressure sensors for measuring the volume change in enclosed PCM containers [6,11], (iii) digital cameras or image sensors for locating the phase change front [12], and (iv) electrical conductivity sensors for determining the phase at selected locations in electrically conductive PCMs [6,13]. Locating the phase change front using distributed temperature sensors is possible, but calculating the stored energy from temperature measurements introduces large uncertainties because of hysteresis, undercooling, and the isothermal nature of phase change [7]. Beyne et al. discuss additional methods for estimating the stored energy; one method involves measuring the heat transfer rate between the working fluid and the PCM with temperature and mass flowrate sensors in the fluid only and then integrating the heat transfer rate over time to calculate the stored energy. This method requires no temperature sensors embedded in the PCM, but heat transfer between the TES and the surroundings must be negligible for the method to be valid [7].

Research on model-based dynamic state estimation for thermal energy storage is limited. Barz et al. design and implement an extended Kalman filter (EKF) for temperature and SOC estimation in a shell-and-tube TES device using temperature measurements. The authors find that the SOC estimated by the EKF tracks that of a high-fidelity simulation more closely than the SOC calculated directly from the temperature measurements [4]. Similarly, Pernsteiner et al. implement an EKF for state and parameter estimation of a reduced-order PCM-based TES simulation model [14]. Jaccoud et al. investigate a particle filter, a type of Monte Carlo method for state estimation, for estimating the temperature and phase change front location in a one-dimensional heat transfer problem [15]. However, a limitation of these works is the absence of guarantees on the convergence of the state estimates. Although Barz et al. [4] verify local observability, none of these authors attempt to prove uniform observability or uniform detectability of the nonlinear system model or guarantee convergence of the state estimates. One exception is Morales Sandoval et al. who utilize a nonlinear Luenberger observer for temperature estimation in a sensible TES device, a hot water thermal storage tank. The authors check observability of their five-state model but find that the system is only observable when all five states are measured [16]. Hence, the thermal storage device must be designed with observability and state estimation in mind to ensure there are a sufficient number of sensors to fully observe the system. Increasing the number of states in a thermodynamic model would improve model accuracy, but the higher order model could be unobservable, leading to divergence of the state estimate [4,16]. In many latent TES architectures, including a sensor for each state in the model is not feasible because of manufacturability or cost constraints, so the estimation model must be constructed such that the system is observable, or at least detectable, given limited measurements.

1.2 Contribution.

We fill this gap in the literature by designing and experimentally validating a state-dependent Riccati equation (SDRE) filter [1719] for a PCM-based thermal energy storage device integrated with a single-phase cooling loop. In contrast to the ubiquitous nonlinear state estimators—the extended Kalman filter [20] and the unscented Kalman filter [21,22]—the SDRE filter uses a linear parameter-varying model parameterized by the state of the system with which uniform detectability and boundedness of the estimation error can be guaranteed in advance for many types of dynamic systems [23]. A key difference between the SDRE filter and the EKF is that the SDRE filter does not require the Jacobian matrix for linearization. For latent TES, the highly nonlinear nature of phase change makes using estimation methods requiring linearization especially undesirable; derivation of the Jacobian matrix can be difficult, and its calculation for higher order models can be too computationally intensive for real-time estimators. Additionally, in highly nonlinear systems, linearization can lead to poor performance, instability, and loss of observability. Fortunately, these issues can be avoided by using the SDRE filter [24].

We leverage a finite volume model of a PCM-based TES and limited temperature measurements to estimate the temperature distribution inside the TES using the continuous-discrete SDRE formulation [19]. A state of charge metric is defined that is directly calculated from the estimated temperatures. We show that the nonlinear finite volume model, when parameterized as a linear parameter-varying model for the SDRE filter, is uniformly detectable, thereby guaranteeing boundedness of the state estimation error. Furthermore, we show that the SDRE approach can be generalized to any phase-change TES architecture by constructing the finite volume heat transfer model with a strongly connected weighted graph, which will be uniformly detectable given at least one state measurement. While there are other nonlinear state estimators (particle and Kalman-type) [25] which, like the SDRE, do not require the Jacobian, this graph-structured model for which we can directly show detectability and error boundedness with limited sensing makes the SDRE a natural choice. Through a series of simulated and experimental case studies, we demonstrate the boundedness and accuracy of the SDRE filter applied to the nonlinear estimation problem of a latent TES device.

This paper is organized as follows. In Sec. 2 we present the TES architecture as well as the system model—both a higher fidelity one used for simulation and a reduced-order model used for estimation. Section 3 describes the SDRE filter formulation and proves that the estimation model is uniformly detectable. The experimental setup is described in Sec. 4. Section 5 presents a series of case studies demonstrating the estimator's accuracy and convergence in both simulation and experiments.

2 Thermal Energy Storage Thermodynamic Model

The TES device considered in this work is designed for use in a vehicle thermal management system to provide supplementary heat rejection capability during periods of large transient heat loading. The TES design consists of a paraffin phase change material, hexadecane, with a melting point of 289.5 K [26] embedded in a rectangular fin heat sink. The working fluid enters the device and flows along a flat metal separator plate below the PCM layer, as shown in Fig. 1. The fins in the PCM layer increase the thermal conductivity of the PCM/fin composite. The entire TES module is contained in an insulating enclosure to reduce heat transfer with the surroundings.

Fig. 1
Section view of a 3D model of the TES module. The fluid outlet is not visible, but the device is symmetric and the outlet geometry is identical to the inlet.
Fig. 1
Section view of a 3D model of the TES module. The fluid outlet is not visible, but the device is symmetric and the outlet geometry is identical to the inlet.
Close modal

Experimental testing of the TES design has been used to validate a high-fidelity nonlinear simulation model, which is briefly discussed in Sec. 2.1. More detail about the simulation model can be found in Refs. [27] and [28]. The TES model used for state estimation, discussed in Sec. 2.2, is based on this validated simulation model.

2.1 Finite Volume Method.

To model heat transfer in the TES, we derive a model based on finite volume discretization as described in Ref. [27]. Here, we summarize the model for the benefit of the reader. We discretize the cross section of the device into a grid of nx×ny control volumes, as shown in Fig. 2, and derive an energy balance equation for each control volume. The fluid channel and metal plate each comprise one layer of control volumes, and the remaining ny2 layers contain the PCM and metal fins. The set of all control volumes in the fluid channel is denoted as F, and the fluid control volumes are numbered j =1 to j = nx from the inlet to the outlet1. The set of metal plate control volumes, comprising j=nx+1 to j=2nx, is denoted as P, and the set of PCM/fin control volumes (j=2nx+1 to j=nxny) is denoted as C.

Fig. 2
Cross section of the TES module considered in this work (not to scale). A finite volume heat transfer model is constructed by discretizing the cross section into a grid of rectangular control volumes. In this figure, the dashed lines demarcate the grid of control volumes for the model used for state estimation.
Fig. 2
Cross section of the TES module considered in this work (not to scale). A finite volume heat transfer model is constructed by discretizing the cross section into a grid of rectangular control volumes. In this figure, the dashed lines demarcate the grid of control volumes for the model used for state estimation.
Close modal
Equation (1a) defines an energy balance for control volume j. Heat transfer along the width of the TES (z direction) is assumed to be negligible, and temperatures across the width are assumed to be uniform. The Q˙jadv term, defined in Eq. (1b), represents heat transfer in fluid control volumes due to mass transfer, or advection; m˙ is the mass flowrate, and cp,f is the specific heat of the working fluid. Conductive heat transfer in the fluid channel is neglected because it is negligible compared to the advective and convective heat transfer rates. The Q˙ij term, defined in Eq. (1c), represents the conductive (solid-solid) or convective (fluid-solid) heat transfer rate from i to j, where iA(j), and A(j) is the set of up to four control volumes adjacent to j. Heat transfer between adjacent control volumes is modeled using the thermal resistance, given in Eq. (1d), which is either conductive or convective. In Eq. (1d), di,j is the distance between the centers of control volumes i and j, κj is the thermal conductivity of j, ai,j is the area of the boundary between i and j, and Uj is the convective heat transfer coefficient of fluid control volume j
mjcp,jdTjdt=Q˙jadv+iA(j)Q˙ij
(1a)
Q˙jadv={m˙cp,f(Tj1Tj)jF0jPC
(1b)
Q˙ij=TiTjRj,i+Ri,j
(1c)
Rj,i={1Ujai,jiP;jFdi,j2κjai,jjPC
(1d)

The PCM/fin layer is modeled as a composite material composed of metal fins and PCM, which we call a composite PCM or CPCM. In Eqs. (1a) and (1d), the thermal conductivity κj and specific heat cp,j are temperature-dependent composite properties for control volumes in the CPCM layer. Thermal properties of the metal plate and fluid layers are assumed constant. Equations for the composite properties can be found in Refs. [29] and [30]. Validation of the composite assumption for closely-spaced parallel rectangular fins is given in Ref. [30].

The phase change dynamics are modeled using the effective specific heat function given in Eq. (2), which is validated in Refs. [31] and [32]. Near the phase change temperature Tpc, the effective specific heat of CPCM is greatly increased so that the latent heat is modeled as sensible heat over the temperature range Tpc±ΔTpc2, with ΔTpc=8 K, as shown in Fig. 3. The width of this temperature range is determined by the parameter α=8ΔTpc. Differential scanning calorimetry (DSC) of hexadecane has shown that the phase change occurs over this range of temperatures, although hexadecane also exhibits undercooling during the solidification process and thermal hysteresis between the melting and solidifying temperatures [33]. These additional phase-change phenomena are neglected in the simulation model. Outside the latent temperature range, the specific heat approaches the liquid specific heat cp,liq for Tj>Tpc or the solid specific heat cp,sol for Tj<Tpc. The specific enthalpy of fusion, hfus, and the solid and liquid specific heats, cp,sol and cp,liq, represent properties of the CPCM.
cp,j(Tj)=cp,sol+(cp,liqcp,sol)1+eα(TjTpc)+hfusα2+eα(TjTpc)+eα(TjTpc)jC
(2)
Fig. 3
Effective specific heat function for simulating phase change in the TES device. Shading represents the approximate range of temperatures over which the phase change occurs. The central vertical line is the true melting point, 289.5 K.
Fig. 3
Effective specific heat function for simulating phase change in the TES device. Shading represents the approximate range of temperatures over which the phase change occurs. The central vertical line is the true melting point, 289.5 K.
Close modal
During phase change, a control volume will contain both solid and liquid CPCM. The partially melted CPCM is then treated as an isotropic composite material [32]; the volume fraction of liquid CPCM is given in Eq. (3) [31].
fm,j=11+eα(TjTpc)jC
(3)

2.2 Model Used for State Estimation.

The number of states in the finite volume model is a function of the level of discretization chosen by the user. For a higher fidelity model, the user may choose nx and ny to be of the order of tens of control volumes. However, this would result in a model that is not suitable for state estimation. Instead, for the purpose of state estimation, we choose a grid of nx = 3 by ny = 7 control volumes for a total of n =21 states. Defining an energy balance for each control volume yields a system of differential equations that can be represented using a linear parameter-varying state-space model. This model can be quickly derived using graph-based methods by representing the finite volume model as the thermal resistance network in Fig. 4.

Fig. 4
The finite volume model forms a resistance network and can be analyzed as a connected graph. Each edge consists of two series resistances, one for each control volume connected by the edge. The fluid control volumes are not connected by edges because conduction in the fluid is neglected.
Fig. 4
The finite volume model forms a resistance network and can be analyzed as a connected graph. Each edge consists of two series resistances, one for each control volume connected by the edge. The fluid control volumes are not connected by edges because conduction in the fluid is neglected.
Close modal

Proposition 1. The thermal resistance network in Fig. 4 forms a connected undirected weighted graph. The control volumes constitute the vertices of the graph, and the thermal resistances between control volumes are the graph's edges.

Proof. If a path, or a set of contiguous edges, exists between any two vertices, an undirected graph is said to be connected. Proposition 1 follows from this definition.▪

The state-space model in Eq. (4) is derived by calculating the graph's Laplacian matrix, L(x), in addition to a diagonal capacitance matrix, M(x), and an input matrix, B(x), where the state vector x=[T1Tn]n contains the temperatures of all n control volumes; this includes the fluid control volumes, metal plate control volumes, and CPCM control volumes. The notation xk refers to the value of x at time tk, or xk=x(tk). The control input u1 is the mass flowrate, and ykp is the (output) vector of states that are measured at time-step k.
20x˙=M1(x)L(x)x+B(x)u=A(x)x+B(x)u
(4a)
yk=Cxk
(4b)
The Laplacian matrix L(x)n×n for the weighted graph is defined in Eq. (5a)2. The matrix M(x)n×n is defined in Eq. (5b) and contains the thermal capacitance mjcp,j of each control volume; for CPCM control volumes, the specific heat is defined using Eq. (2). The input matrix B(x)n×1 is defined in Eq. (5c) and contains the advection terms. Finally, Tin is the temperature of the fluid at the inlet, which is not a state. Recall that T1, T2, …, Tnx are the fluid control volume temperatures, which are the first nx states in x.
Li,j={1Rj,i+Ri,jijandedge(i,j)0ijandedge(i,j)lA(j)(Ll,j)i=j
(5a)
M(x)=[m1cp,1000m2cp,2000mncp,n]
(5b)
B(x)=cp,fM1[TinT1T1T2Tnx1Tnx00]n×1
(5c)

Suppose measurements of the temperature of certain control volumes are available at discrete time instances tk, k=1,2,; then Cp×n is a binary output matrix where Ci,j=1 if yk,i=xk,j, and Ci,j=0 otherwise. Thus, the output equation (Eq. (4b)) is discrete-time although the state dynamics are continuous.

2.3 State of Charge.

We define the state of charge (SOC) of the TES in Eq. (6) as a piecewise linear function of the energy stored in the device, or the enthalpy of the CPCM, H [1]. Following the convention described in Ref. [5], minimum stored energy corresponds to a maximum SOC of 1, and maximum stored energy corresponds to a minimum SOC of 0. Note that this definition of SOC considers both latent and sensible energy storage within the TES.
xsoc={1H<HminHmaxHHmaxHminHminHHmax0H>Hmax
(6)
With the finite volume model, the total stored energy, given in Eq. (7), is the sum of the enthalpies of each CPCM control volume. The specific enthalpy of the PCM/fin composite in Eq. (8) is found by integrating the effective specific heat function, Eq. (2), with respect to temperature. Figure 5 shows the shape of the specific enthalpy curve; the specific enthalpy is defined such that it is zero at the phase change temperature, Tpc. The maximum enthalpy Hmax and minimum enthalpy Hmin in Eq. (6) are calculated for the PCM at a uniform temperature Tmax = 308 K or Tmin = 278 K, which are user-defined parameters corresponding to the maximum and minimum expected temperatures, respectively. Thus, the SOC of the TES can be calculated given the temperature of each CPCM control volume.
H=jCmjh(Tj)
(7)
h(TjC)=TpcTjcp,j(x)dx=hfus2tanh[α2(TjTpc)]+(TjTpc)cp,sol+cp,liqcp,solαln[1+eα(TjTpc)2]
(8)
Fig. 5
Specific enthalpy function for the CPCM derived by integrating the effective specific heat function. Shading represents the approximate range of temperatures over which the phase change occurs. The central vertical line is the phase change temperature, 289.5 K.
Fig. 5
Specific enthalpy function for the CPCM derived by integrating the effective specific heat function. Shading represents the approximate range of temperatures over which the phase change occurs. The central vertical line is the phase change temperature, 289.5 K.
Close modal

3 State-Dependent Riccati Equation Filter

In Ref. [17], Mracek et al. discuss the continuous-time formulation of the SDRE filter in which the observer gain is obtained by solving the continuous-time state-dependent Riccati equation [17]. In Ref. [18], Jaganath et al. present the discrete-time SDRE filter in which the discrete-time Riccati equation is solved recursively with prediction and update steps resembling the EKF algorithm. For systems with continuous dynamics and discrete sampling, Berman et al. [19] describe the continuous-discrete SDRE filter in which the state estimates and covariance matrix are propagated in discrete time steps using a state-dependent state transition matrix. In this section we outline the continuous-discrete SDRE filter algorithm and describe its application to the latent TES state estimation problem. As we will show in Sec. 5.4, the continuous-discrete formulation can be particularly advantageous when the sample rate is limited.

Let the input u=m˙(t) be constant over each interval Δtk=tk+1tk. Further assume some process noise wkN(0,Wk) and measurement noise vkN(0,Vk) are present, where Wkn×n,Vkp×p are covariance matrices. If the system dynamics change little during the interval Δtk, the linear parameter-varying state-space model can be “frozen-in-time” (see Ref. [34]) and discretized using the exact solution for the state transition matrix of a linear dynamic system. The discrete-time system is given by
xk+1=Φkxk+Γkuk+wk
(9a)
yk=Cxk+vk
(9b)
where
Φk=eA(xk)Δtk
(10a)
Γk=0ΔtkeA(xk)τB(xk)dτ
(10b)
To estimate the unmeasured temperature states of the TES, we use the two-step recursive formulation of the SDRE filter given in Refs. [18] and [19] where the predicted state x̂k+1|k and estimation error covariance Pk+1|k are given by
x̂k+1|k=Φkx̂k+Γkuk
(11a)
Pk+1|k=ΦkPk|kΦk+Wk
(11b)
the Kalman gain Kk is given by
Kk=Pk+1|kC(CPk+1|kC+Vk)1
(12)
and the updated state, output and covariance estimates are
x̂k+1|k+1=x̂k+1|k+Kk(ykŷk)
(13a)
ŷk=Cx̂k
(13b)
Pk+1|k+1=(IKkC)Pk+1|k
(13c)

Using a smaller time interval, Δtk in the discrete-time model results in a better approximation of the continuous nonlinear dynamics. If measurements are not available at every time-step tk, the prediction step in Eq. (11) can be performed multiple times between each update step; we call this method the continuous-discrete SDRE filter [19] although it still uses a discrete approximation of the continuous dynamics.

Now we establish the boundedness of the SDRE filter when applied to the latent TES under consideration.

proposition 2. The estimation error covariance of the SDRE filter for the TES model given by Eq. (4) will remain bounded.

Proof. For the discrete-time system of Eq. (9), the SDRE error covariance is bounded if the pair (Φk,C) is uniformly detectable [35]. The TES model given by Eq. (4) represents a continuously connected undirected graph since Ai,j(x)>0t for each existing edge (i, j) where ij. If a parameter-varying undirected graph is continuously connected t, the graph implements consensus, and is, hence, uniformly detectable from any node (see the Appendix). Furthermore, the pair (Φk,C) is uniformly detectable for any C with at least one nonzero row sum (see Lemma 3 in the Appendix). Therefore, it follows that the error covariance of the state estimate will remain bounded.▪

Remark. Proposition 2 can be generalized to show that any thermal energy storage device that can be spatially discretized into a contiguous lattice of control volumes (that is, the thermal resistance between any two control volumes is finite) is uniformly detectable given that at least one state measurement exists. In other words, when the SDRE filter is used for state estimation, boundedness of the state estimate error covariance is guaranteed.

4 Experimental Testbed

The TES is integrated into the experimental single-phase thermal-fluid loop testbed shown in Fig. 6. The testbed has four TES modules arranged in series, each with a capacity of approximately 25.6 kJ. Each module uses a separate state estimator with a model tuned to the specific module's parameters, but in this work, we will show validation of the state estimator on only the first module in the series.

Fig. 6
Images of the experimental single-phase thermal-fluid testbed with thermal energy storage
Fig. 6
Images of the experimental single-phase thermal-fluid testbed with thermal energy storage
Close modal

Figure 7 shows a simplified diagram of the thermal management system; a pump circulates the working fluid (water) though the primary fluid loop with a variable flowrate up to 0.183 kg/s. Heat is added through a 6 kW 600 VDC electrical resistive heater mounted to an Advanced Thermal Solutions ATS-CP-1002 cold plate. The primary mode of heat rejection is a shell-and-tube heat exchanger cooled by a secondary chilled water loop from a Neslab HX300-DD 10 kW chiller. A pair of Burkert Type 2875 variable position solenoid valves control the mass flowrate through the TES. The four modules are custom designed and fabricated for the testbed, and together they provide up to 2 kW of heat rejection.

Fig. 7
Simplified schematic of the experimental thermal management system
Fig. 7
Simplified schematic of the experimental thermal management system
Close modal

Temperature measurements of the PCM layer in the TES are provided by five type T thermocouples bonded to the PCM side of the plate, labeled TC2a, TC2b, TC3, TC4a, and TC4b, in the pattern shown in Fig. 8. TC2a and TC2b are placed at the same distance along the length of the module, and likewise for TC4a and TC4b. Type T thermocouples labeled TC0 and TC1 in the fluid channel measure the inlet and outlet fluid temperature, respectively. Not shown in Fig. 8 is an Omega FTB-1313 turbine flowmeter upstream of the TES that measures the mass flowrate of the working fluid.

Fig. 8
Section view of a 3D model of the TES module. Thermocouples TC0 and TC1 measure the fluid inlet and outlet temperatures, respectively, and TC2a, TC2b, TC3, TC4a, and TC4b are embedded in the PCM.
Fig. 8
Section view of a 3D model of the TES module. Thermocouples TC0 and TC1 measure the fluid inlet and outlet temperatures, respectively, and TC2a, TC2b, TC3, TC4a, and TC4b are embedded in the PCM.
Close modal

The thermocouples are sampled at regular intervals to produce the measurement output vector yk. Thermocouple TC0 measures Tin, which is not a state but a parameter in the input matrix B(x). Thermocouple TC1 measures the output from control volume CV1; the fluid outlet temperature is assumed to be equal to the temperature of the fluid in CV1. The five thermocouples in the PCM layer lie within three control volumes (see Fig. 9). Thermocouples TC2a and TC2b lie within the same control volume, so measurements from these thermocouples are averaged to provide one output measurement, TC2, for the control volume labeled CV2. Thermocouple TC3 measures the output for control volume CV3. Measurements from thermocouples TC4a and TC4b are also averaged to provide a single output measurement, TC4, for control volume CV4.

Fig. 9
Locations of thermocouples TC0-4 relative to the control volume grid, denoted by × marks. Selected control volumes used for validation are also labeled CV1-6.
Fig. 9
Locations of thermocouples TC0-4 relative to the control volume grid, denoted by × marks. Selected control volumes used for validation are also labeled CV1-6.
Close modal

During experimental tests, the SDRE filter is deployed on a National Instruments PXIe-8820 embedded controller with an Intel Celeron 1020E 2.20 GHz dual-core processor. Data aquisition peripherals include a NI PXI-6225 16-bit analog input module for reading the flowmeter and a NI PXIe-4353 24-bit thermocouple module. Additionally, experimental data is recorded so the state estimator can be simulated offline for validation. A custom LabVIEW VI manages the data acquisition, executes the state estimator algorithm, and translates user input into voltage signals to operate the valves, heaters, and pump.

The thermocouple sensor noise variance is determined by repeatedly sampling a nearly constant temperature, fitting a linear regression model to the sampled data, and calculating the residual variance. Each thermocouple is assumed to have independent Gaussian white noise. The sample variance of each thermocouple is calculated from a dataset containing 4000 samples; the median value of the variances is σ2=0.0068K2. Given the similarity in noise variance across the thermocouples (which is expected given that each is identical and purchased from the same manufacturer), a single value of σ2=0.007K2 (to the nearest thousandth) is assumed for each thermocouple. Averaging the measurements from laterally adjacent thermocouples reduces the variance of these outputs to Var[ya+yb2]=σ22=0.0035K2. Assuming the output is yk=[yk,TC1yk,TC2yk,TC3yk,TC4], the measurement noise covariance matrix is given in Eq. (14a). The process noise covariance is not as easily characterized, so we assume it is negligible; Eq. (14b) defines the matrix Wk just large enough to ensure that Pk+1|k remains positive definite.
Vk=[0.00700000.003500000.00700000.0035]
(14a)
Wk=107In
(14b)

5 Testing and Validation

The objective of the following case studies is to validate the SDRE filter by demonstrating that the errors in the estimated temperatures and state of charge converge near zero and remain bounded. As it is not possible to validate the SOC estimate experimentally with only five thermocouples embedded in PCM layer of the TES module, we first test the estimator using the high-fidelity simulation model validated in Refs. [27] and [28]. However, in order to enable a comparison with the experimental validation and draw conclusions about those estimation errors that are due to modeling error, we use experimental data as inputs to the simulation model, as discussed in Sec. 5.1. The simulation-based case study in Sec. 5.2 is followed by validation of the estimator's convergence and boundedness on the experimental testbed as discussed in Sec. 5.3. In Sec. 5.4, we investigate the effect of increased and reduced sample rates by comparing 80 samples per second (S/s), 10 S/s, 1 S/s, and 0.2 S/s. Although the SDRE filter is capable of real-time estimation on the testbed, we perform experimental validation offline so that multiple filter configurations can be tested using the same dataset.

5.1 Experimental Data Collection.

We first collect a dataset containing temperature measurements from the seven thermocouples and a measurement of the mass flowrate on the experimental testbed at the maximum sample rate of 80 S/s. The experimentally-measured mass flowrate and fluid inlet temperature are used as inputs for the simulation-based testing described in Sec. 5.2, whereas the entire experimental dataset is used for experimental validation of the estimator in Sec. 5.3. The dataset is downsampled when testing the SDRE filter with slower sample rates. The mass flowrate through the TES is controlled by manually opening and closing the flow control valves through the LabVIEW VI. Figure 10 shows the measured mass flowrate through the TES. Figure 11 shows the temperatures of the fluid entering and leaving the TES measured by TC0 and TC1, respectively. Note that periods of zero flowrate were intentionally introduced as it would be reasonable for the TES to only be charged or discharged intermittently. Conductive heat transfer still occurs inside the TES even when the mass flowrate is zero, so the state estimator must still work during these periods to track the changing temperature distribution. In the subsequent figures, blue shading denotes periods of zero mass flowrate when this information is relevant to the result presented in the figure.

Fig. 10
Measured mass flowrate through the TES. Blue shading in this figure and subsequent figures highlights periods of zero mass flowrate.
Fig. 10
Measured mass flowrate through the TES. Blue shading in this figure and subsequent figures highlights periods of zero mass flowrate.
Close modal
Fig. 11
Inlet and outlet fluid temperature measurements
Fig. 11
Inlet and outlet fluid temperature measurements
Close modal

In Sec. 5.4, we investigate the effect of sample rate on the accuracy of the estimator. To make an unbiased comparison between different sample rates, the prediction equations given in Eq. (11) are executed with the time-step Δtk=0.0125 s so that the state estimator model is identical for all sample rates. The update step given in Eqs. (12) and (13) is performed once every 1, 8, 80, or 400 time steps depending on whether the sample rate is 80 S/s, 10 S/s, 1 S/s or 0.2 S/s, respectively. Note that the filters receiving measurements at 10 S/s, 1 S/s, and 0.2 S/s are continuous-discrete filters, but the filter receiving measurements at 80 S/s is a discrete filter since the prediction and update steps are executed at the same rate.

5.2 Simulation-Based Testing.

Using the experimentally-measured mass flowrate and fluid inlet temperature (TC0) as inputs, we simulate the high-fidelity simulation model [27,28] to generate a simulated dataset containing the SOC and temperatures of all control volumes. Note that this model has the same finite volume structure as the model used by the estimator albeit with a finer spatial discretization of nx = 21 and ny = 22; the finer discretization was chosen so that the simulation model is an accurate representation of the experimental system. In other words, differences in accuracy between the simulation and the model used for state estimation are representative of the differences between the experimental TES and the model used for state estimation. In the model used for state estimation, low computation time is prioritized over model accuracy. The simulation is executed in MATLAB using the ode23 solver [36]. The simulation calculates the temperature distribution inside the TES, converts the high-fidelity temperature distribution to a set of temperatures corresponding to the state estimator model's control volumes, and generates “thermocouple measurements” with a sample rate of 10 S/s by adding Gaussian noise to the temperatures of those control volumes that map to the measurement locations on the experimental TES module at discrete time steps. These simulated thermocouple measurements comprise the measurement vector yk for the state estimator; the thermocouple locations and the labels of the corresponding control volumes are given in Fig. 9. We then test the state estimator on the simulated dataset, using simulated values of the thermocouple measurements, to establish convergence of the SOC estimate.

First, we verify that the temperature estimates of unmeasured control volumes converge toward the temperatures calculated by the simulation, and that the estimation errors remain bounded. In this case study, the state estimator receives the full measurement output vector at a rate of 10 S/s. Figure 12(a) compares the estimated temperatures of the control volumes labeled CV5 and CV6 (see Fig. 9 for the locations of these control volumes) to the corresponding simulated temperatures. The gray shading represents the range of temperatures over which the latent heat of fusion is approximated by a large increase in specific heat. Since it is difficult to distinguish the estimates from the simulated values in Fig. 12(a), we show the estimation errors in Fig. 12(b). Note that all estimation errors are calculated as e(tk)=x̂(tk)xsim(tk) where xsim is a state variable calculated by the simulation.

Fig. 12
(a) Estimated temperatures (red) and simulated temperatures (black) at CV5 and CV6. Gray shading denotes the range of temperatures corresponding to the latent region, and the thin horizontal line is the melting point of the PCM. (b) Estimation errors for CV5 and CV6.
Fig. 12
(a) Estimated temperatures (red) and simulated temperatures (black) at CV5 and CV6. Gray shading denotes the range of temperatures corresponding to the latent region, and the thin horizontal line is the melting point of the PCM. (b) Estimation errors for CV5 and CV6.
Close modal

Even though the temperatures of CV5 and CV6 are not included in the output vector yk, the SDRE filter is able to track these temperatures with an estimation error magnitude less than 0.6 K. During periods when the temperatures are changing slowly, the estimation error magnitude is less than 0.05 K, which is less than one standard deviation of the measurement noise (0.0035K2=0.059K). In fact, the largest estimation errors only occur during periods of phase change; the larger negative transient estimation errors correspond to periods when the PCM is melting, suggesting that the SDRE filter introduces a small phase lag between the estimates and the true values. This claim is supported by the observation that the estimation error has a larger negative value when the temperatures change rapidly, such as at t =1050 s. In turn, a faster change in temperature results in more phase lag, which increases the transient estimation error. There are fewer periods of large positive error because the PCM only fully solidifies at the start (t =50 s) and end (t =1650 s) of the dataset; this means that the largest transient errors are expected when the experimental TES completes a phase change and begins to experience sensible heating while the state estimates lag behind in the latent temperature range.

The primary motivation for using the SDRE as compared to other nonlinear filters is its suitability for the application at hand, namely, SOC estimation for latent TES systems with graph-based models. Nevertheless, comparison between the SDRE filter, the extended Kalman filter (EKF), and the unscented Kalman filter (UKF) suggests that the SDRE filter achieves the lowest estimation error, as shown in Fig. 13, which compares the estimation errors at CV5 for the three state estimators. Additionally, the two Kalman filters have large estimation errors during the same periods when the SDRE filter has a large estimation error, so we can conclude that these errors result from limitations in the model used for state estimation and not from the choice of state estimator.

Fig. 13
Comparison of estimation error at CV5 for the EKF, the UKF, and the SDRE filter
Fig. 13
Comparison of estimation error at CV5 for the EKF, the UKF, and the SDRE filter
Close modal

Next, we quantify the temperature estimation error when the simulated measurement from (i) TC1 or (ii) TC3 is withheld from the state estimator. Note that when testing the state estimator on experimental data (see Sec. 5.3), these thermocouples are withheld to be used for validation. Figure 14(a) shows the temperatures at CV1, calculated by the simulation and estimated by the SDRE filter, when TC1 is withheld, and Fig. 14(b) shows the corresponding estimation error. As expected, the estimation error is more substantial when an output is removed; the error magnitude tends to be less than 0.5 K when the control volume temperature is changing slowly, but it increases to nearly 2 K when the temperature changes suddenly, as it does at t =1050 s or t =1200 s. The steady-state estimation errors from t =850 to t =1000 s and t =1350 to t =1550 s are due to the zero mass flowrate through the TES. Since advection is the dominant form of heat transfer in the fluid channel, when advection is removed, heat transfer in the fluid control volumes is greatly reduced. This causes the temperature estimate to converge more slowly.

Fig. 14
(a) Simulated temperature at CV1 and the estimated temperature for CV1. (b) Estimation error.
Fig. 14
(a) Simulated temperature at CV1 and the estimated temperature for CV1. (b) Estimation error.
Close modal
Figure 15(a) shows that the estimated temperature for control volume CV3 still tracks the simulated temperature for this control volume closely when TC3 is withheld from the estimator; Fig. 15(b) shows that the error tends to be largest near the melting point and during periods of zero mass flowrate. The estimation error magnitude stays within 1 K, but it only converges to zero outside the latent temperature range. Figure 16 shows the root-mean-square error (RMSE) averaged over all n control volumes at time-step tk, given in Eq. (15). In this figure, the RMSE is shown for the state estimator with full access to all measurements, the state estimator with access to all measurements except TC1, and the state estimator with access to all measurements except TC3. The RMS errors of the estimator with all measurements and the estimator without TC1 remain below 0.4 K (except at t =0 s, when the state estimates are initialized) and are nearly indistinguishable; this indicates that the measurement from TC1 has a negligible effect on the estimates of the other control volumes in the model. Removing TC3 slightly increases the RMSE during phase changes and periods of zero mass flowrate. Given that this sensitivity to TC3 is not a result of the system model itself, we verify that it occurs when testing the estimator on the experimental dataset in Sec. 5.3.
erms(tk)=[i=1n(x̂i(tk)xi(tk))2n]12
(15)
Fig. 15
(a) Simulated temperature at CV3 and the estimated temperature for CV3. (b) Estimation error.
Fig. 15
(a) Simulated temperature at CV3 and the estimated temperature for CV3. (b) Estimation error.
Close modal
Fig. 16
RMS errors of all control volumes for the state estimator with all available measurements, all measurements except TC1, and all measurements except TC3
Fig. 16
RMS errors of all control volumes for the state estimator with all available measurements, all measurements except TC1, and all measurements except TC3
Close modal

Although it is important to analyze the error of individual temperature estimates, the overarching goal of the proposed estimator is to estimate the SOC of the TES for online control decision-making. Figure 17(a) compares the simulated and estimated SOC values calculated from the temperature estimates of the SDRE filter with access to all outputs. The shading in Fig. 17(a) represents the range of SOC in which all or part of the PCM is undergoing a phase change. Figure 17(a) shows simulated PCM temperatures for three control volumes over the duration of the simulation; refer to Fig. 9 for the location of these control volumes. Comparing Fig. 17(a) to Fig. 18 shows that the SOC is near 1 when the PCM temperatures are low and near 0 when the PCM temperatures are high. In Fig. 18, the horizontal line is the melting point of the PCM near which small changes in temperature result in large changes in SOC.

Fig. 17
(a) Simulated and estimated SOC; gray shading represents the latent SOC region where errors are expected to be large. (b) SOC estimation error; blue shading represents periods when the mass flowrate through the TES is zero.
Fig. 17
(a) Simulated and estimated SOC; gray shading represents the latent SOC region where errors are expected to be large. (b) SOC estimation error; blue shading represents periods when the mass flowrate through the TES is zero.
Close modal
Fig. 18
Simulated temperature values from select PCM control volumes
Fig. 18
Simulated temperature values from select PCM control volumes
Close modal

This sensitivity near the melting point is also evident in Fig. 17(b), which shows the SOC estimation error. The error is maximized when the PCM temperatures are near the melting point and the SOC is within the latent heat range. Despite the nonlinear relationship between SOC and temperature, the estimated SOC tracks the simulation well, and the estimate error remains within ±0.02. Outside the latent range when the system dynamics are approximately linear (between t =50 and t =380 s, for example), the SOC error converges to zero (although there is some noise with amplitude less than 0.001). The SOC error tends to be positive when the SOC is decreasing and negative when the SOC is increasing, which indicates that there is a small phase lag between the SOC estimate and the true SOC. Additionally, the error does not converge to zero when the mass flowrate through the TES is zero and the PCM temperatures are within the latent range (for example, between t =850 and t =1000 s). This is again because heat transfer along the length of the module is dominated by advection. When the mass flowrate is zero, it takes longer for the PCM to reach an equilibrium temperature, which means it takes longer for temperature estimation errors to converge to zero. Small temperature errors near the phase change temperature result in large errors in SOC. Outside the latent temperature range, the SOC error converges to zero when the mass flowrate is zero—see the period between t =1050 and t =1200 s for an example.

5.3 Experimental Validation.

In this case study, we test the SDRE filter on experimental data. First, we investigate convergence of the SDRE filter estimates when (i) thermocouple TC1 is withheld for validation and (ii) thermocouple TC3 is withheld for validation. We also investigate how the SDRE filter performs with different sample rates. All estimation errors are calculated as e(tk)=x̂(tk)yk,j where yk,j is a measurement from an excluded thermocouple.

When TC1, the fluid outlet thermocouple measurement, is excluded from the measurement set, the state estimator tracks the fluid outlet temperature closely over long time periods, as shown in Fig. 19(a). However, the estimation error increases when the temperature changes rapidly, as shown in Fig. 19(b) at t =1000 s, t =1200 s, and t =1550 s. For most of the duration of the experiment, the magnitude of the estimation error at CV1 is less than 1 K. This error behavior is similar to that of the simulated results shown in Fig. 14(b), but the larger magnitude of the error when implementing the estimator on the experimental testbed is due to modeling errors. These include inaccuracies in the reduced-order finite volume model due to spatial discretization [27,28], simplifications of the phase-change dynamics such as neglecting undercooling and hysteresis [33,37], and exogenous disturbances in the physical system like heat transfer with the surroundings.

Fig. 19
(a) Estimated temperature of CV1 compared to the fluid outlet temperature measurement TC1. (b) Estimate error for CV1.
Fig. 19
(a) Estimated temperature of CV1 compared to the fluid outlet temperature measurement TC1. (b) Estimate error for CV1.
Close modal

Excluding TC3 from the measurement set results in degradation of the estimator's performance, as expected based on the simulated results. The estimate of CV3 still tracks TC3 closely over long time periods, as Fig. 20(a) shows. When the TES is experiencing rapid changes in temperature, Fig. 20(b) shows that the magnitude of the error tends to increase; the convergence rate of the estimate error of CV3 is slower than the transient dynamics, resulting in larger errors up to 6 K at t =1050 s and 3 K at t =1550 s. During the period between t =450 and t =900 s and the period between t =1200 and t =1400 s, the PCM undergoes partial phase change. Estimation errors during these periods of time are due to the phase change hysteresis and undercooling effects which delay solidification when the PCM is cooled rapidly [37]. The system model does not account for these effects, so the estimator incorrectly predicts that CV3 (and other unmeasured control volumes) begins to solidify at 293.5 K (the top boundary of the gray shaded region in Fig. 20(a)) whereas the PCM in the experimental TES begins to solidify at a lower temperature. Over sufficiently long time periods during which the PCM temperatures reach equilibrium, such as t =750 to t =1000 s, errors caused by hysteresis and undercooling are corrected and the estimate error returns to zero.

Fig. 20
(a) Estimated temperature of CV3 compared to the measurement of TC3. (b) Estimate error for CV3.
Fig. 20
(a) Estimated temperature of CV3 compared to the measurement of TC3. (b) Estimate error for CV3.
Close modal

5.4 Effect of Sample Rate.

In this section, we investigate how the SDRE filter performs when measurements are provided at four different rates: 80 samples per second (S/s), 10 S/s, 1 S/s, and 0.2 S/s. All previous experimental and simulation results used a sample rate of 10 S/s. The continuous-discrete SDRE filter formulation allows for more flexibility with respect to sample rates because the state estimates and the estimate covariance matrix are propagated with the state transition matrix of the continuous-time system; estimation errors introduced by discretizing the nonlinear system are mitigated by propagating the state estimates and error covariance with a time-step smaller than the interval between samples.

Figure 21 shows (a) the estimated temperatures at CV1 and (b) the estimation errors when TC1 is excluded from the measurement set during a validation test with the experimental dataset. As discussed in the previous section, the state estimator can track the temperature at CV1 well; this is also the case at the faster and slower sample rates. In Fig. 21(b), 80 S/s, 10 S/s, and 1 S/s all appear to have similar estimation errors, but 0.2 S/s produces larger error during periods of transient temperature changes. To more easily compare the performance of the four sample rates, we calculate the root-mean-square errors for the temperature at CV1 averaged over all time steps using Eq. (16) where Ns is the total number of time steps in the dataset for the given sample rate. Results for the four sample rates are included in Table 1. Note that this RMSE definition represents an average over all time steps, but the previous RMSE definition in Eq. (15) averages over all control volumes. As expected, the 80 S/s rate results in the smallest estimation error, but only by a small margin; 10 S/s results in a similar RMSE. The slowest sample rate, 0.2 S/s, results in almost twice the RMSE of the next slowest sample rate, 1 S/s. In the fluid channel, the control volume temperatures can change rapidly due to the high advective heat transfer rates. The larger error in the 0.2 S/s sample rate is due to the fast dynamics of control volume CV1; the estimation error can increase quickly when update steps are not performed fast enough.
erms,j=[k=1Ns(x̂j(tk)yk,j)2Ns]12
(16)
Fig. 21
(a) Estimated temperatures of CV1 compared to the fluid outlet temperature measurement TC1. (b) Estimation errors for CV1.
Fig. 21
(a) Estimated temperatures of CV1 compared to the fluid outlet temperature measurement TC1. (b) Estimation errors for CV1.
Close modal
Table 1

RMS errors at CV1 for four sample rates

Sample rateerms,CV1 (K)
80 S/s0.2436
10 S/s0.2617
1 S/s0.3811
0.2 S/s0.7526
Sample rateerms,CV1 (K)
80 S/s0.2436
10 S/s0.2617
1 S/s0.3811
0.2 S/s0.7526

Figure 22 shows the estimated SOC for the four sample rates. All four estimates are similar with almost no difference outside the latent range. During the periods of partial phase change, a noticeable difference in the peak of the SOC estimate can be seen at t =500 s, t =650 s, and t =1250 s. The faster sample rates reach a sharper and higher peak SOC, suggesting that the use of a faster sample rate results in a faster response time to transient changes.

Fig. 22
Estimated SOC for the four sample rates when TC1 is withheld
Fig. 22
Estimated SOC for the four sample rates when TC1 is withheld
Close modal

Finally, we compare the performance of the four sample rates when TC3 is excluded from the measurement set. Figure 23 shows (a) the estimated temperatures at CV3 and (b) the estimation errors. In Fig. 23(b), it is clear that the state estimator receiving 80 S/s has a slightly smaller error most of the time. This is noticeable at t =500 s, t =1100 s, and t =1600 s. The RMSE of the estimate at CV3 averaged over all time steps is given in Table 2. The state estimator with a sample rate of 80 S/s again performs the best, but by a small margin; this time, there is only a small difference in RMSE between the fastest and slowest sample rates. Since the thermocouple TC3 is embedded in the PCM layer, the dynamics of control volume CV3 are slower than the dynamics of CV1, which explains why sample rate has a lesser effect on RMSE at CV3 compared to the RMSE at CV1 in Table 1. Note that the RMS errors in Table 2 are larger than those in Table 1; this indicates that the measurement from TC3 is more critical to estimate the SOC accurately than that of TC1. This makes sense given that TC3 is located near the center of the TES; without the measurement from TC3, only measurements of the temperatures at the periphery of the PCM layer are available. Estimation errors in CV3 and the control volumes above CV3 can be large for long time periods (greater than 2 K for more than 20 s, according to Fig. 23(b)) before the state estimator compensates for the errors.

Fig. 23
(a) Estimated temperature of CV3 compared to the measurement of TC3. (b) Estimation errors for CV3.
Fig. 23
(a) Estimated temperature of CV3 compared to the measurement of TC3. (b) Estimation errors for CV3.
Close modal
Table 2

RMS errors at CV3 for four sample rates

Sample rateerms,CV3 (K)
80 S/s1.0022
10 S/s1.1046
1 S/s1.2048
0.2 S/s1.3270
Sample rateerms,CV3 (K)
80 S/s1.0022
10 S/s1.1046
1 S/s1.2048
0.2 S/s1.3270

Figure 24 shows the estimated SOC for the four sample rates when the measurement from TC3 is excluded. The results shown in this figure are almost identical to those in Fig. 22, except there is less difference between the 80 S/s results and the 10 S/s results. These results demonstrate an important advantage of the continuous-discrete SDRE filter over other discrete-time state estimators. Since the state estimates are propagated with a small time-step, accuracy of the prediction step is not degraded as the time interval between update steps is increased. Errors due to disturbances or model inaccuracies still affect the accuracy of the estimate, so update steps must be performed sufficiently often to avoid large errors during transient events.

Fig. 24
Estimated SOC for the three sample rates when TC3 is withheld
Fig. 24
Estimated SOC for the three sample rates when TC3 is withheld
Close modal

6 Conclusion

In this paper we designed and experimentally demonstrated a state-dependent Riccati equation (SDRE) filter for nonlinear state and state of charge estimation in a phase-change thermal energy storage device. Unlike other state estimators for nonlinear systems, the SDRE filter uses a linear parameter-varying system model which offers several advantages. First, the LPV model can be computed quickly using graph-based methods, and linearization is not required. Second, the SDRE filter is generalizable to many TES architectures because the finite volume heat transfer model can be applied to any heat conduction problem, even those with highly nonlinear dynamics and material properties. Finally, the structure of the LPV model remains the same, so boundedness and detectability of the SDRE filter can be guaranteed in advance.

The SDRE filter is first tested in simulation to verify that the SOC estimate converges to the true SOC. The SOC estimation error remains bounded within ±0.02, even when the TES experiences rapid temperature changes or the PCM undergoes phase change. In experimental tests, the SDRE filter is shown to accurately estimate temperatures at two locations in the TES. The largest sources of error are model inaccuracies such as (i) the assumption that there is no heat transfer with the surroundings and (ii) a simplified phase-change model which does not account for hysteresis or undercooling during solidification. Additionally, when the mass flowrate through the TES is zero, estimate accuracy is reduced because heat transfer along the length of the TES is typically dominated by mass transfer. Despite the limitations of the finite volume model, the accuracy of the estimator is not significantly affected by the sample rate of the measurements; sample rates of 80 samples/s and 1 sample/s yield similar root-mean-square errors.

Experimental validation of model-based state estimation techniques for fast-timescale thermal storage devices advances the continuing research on control strategies for hybrid transient thermal management systems. In future work, the boundedness and convergence of other Kalman-type and particle filters suitable for SOC estimation may be explored. Additionally, investigation of the effects of model uncertainty and the previously discussed model inaccuracies has the potential to decrease estimation error.

Funding Data

  • U.S. Office of Naval Research Thermal Science and Engineering Program under (Award No. N00014-21-1-2352; doi: 10.13039/100000006).

Nomenclature

a =

area (m2)

A =

state-space coefficient matrix

A =

set of adjacent control volumes

B =

state-space input matrix

cp =

specific heat (JkgK)

C =

state-space output matrix

C =

set of PCM control volumes

d =

distance (m)

e =

estimation error

erms =

root-mean-square error

fm =

melt fraction of a PCM control volume

F =

set of fluid control volumes

h =

specific enthalpy (Jkg)

hfus =

specific enthalpy of fusion (Jkg)

H =

enthalpy (J)

i, j, l =

control volume indices

In =

n × n identity matrix

k =

discrete time step index

K =

Kalman gain

L =

graph Laplacian matrix

m =

mass (kg)

m˙ =

mass flow rate (kgs)

M =

capacitance matrix

n =

number of control volumes

Ns =

number of time steps

N =

Gaussian random variable

p =

number of outputs

P =

state estimate error covariance

P =

set of metal plate control volumes

Q˙ =

heat transfer rate (w)

R =

thermal resistance (KW)

t =

time (s)

T =

temperature (K)

Tpc =

phase-change temperature (K)

u =

input vector

U =

convective heat transfer coefficient (Wm2K)

v =

measurement noise vector

V =

measurement noise covariance

w =

process noise vector

W =

process noise covariance

x =

state vector

x̂ =

estimated state vector

xsoc =

state of charge

y =

output or measurement vector

α =

specific heat function width parameter

Γ =

discrete-time input matrix

κ =

thermal conductivity (WmK)

Φ =

discrete-time state transition matrix

Appendix

Definition 1 (Connected Undirected Graph). A graph is a pairG=(V,E)whereVis is a finite set of n nodes or vertices{v1,,vn}, andE{1,,n}×{1,,n}is a set of edges. The pair (i, j) denotes the edge that links vertex vi to vj. If(i,j)E(j,i)E, the graph is said to be undirected. An undirected graph is connected if there exists a path (set of edges) between any two nodes; the connected edge set is denoted withEc.

Definition 2 (Adjacency and Laplacian Matrices). A (weighted) graph is fully described by the adjacency matrix,Λn×nwhere
Λi,j={λi,jif(i,j)E0otherwise
(A1)

whereλi,jis the weight assigned to any of its edges (i, j). Consider the diagonal matrix D whose diagonal entryDi,i=j=1nλi,j. D contains the total edges' weight that is incident to each vertex, and is termed the degree matrix. The graph Laplacian is defined asL=DΛ.

Definition 3 (Integral Graph). Given a parameter-varying graphG(ρ)=(V,E(ρ))for some time-dependent parameter ρ, the integral graph ofG(ρ)on[0,)is a constant graphG¯[0,):=(V,E¯)whereVis the same vertex set ofG(ρ), and the adjacency matrix is defined by
Λ¯i,j={1,if0+λi,j(ρ)dt=0,if0+λi,j(ρ)dt<
Consider the system where the dynamics of each node xi is given by
x˙i=jN(i)λi,j(ρ)(xjxi)i{1,,n}
(A2)
and λi,j(ρ) is the weight on the edge if (i,j)E (otherwise, λi,j=0). This system corresponds to a weighted graph and can be represented in state-space form as
x˙=L(ρ)x
(A3)

where L is the Laplacian. Lemma 1 holds for the system [38].

Lemma 1 ([38]). The dynamics of (A3) implements consensus, that is,limt+x(t)span{1n}if and only ifG¯[0,)is connected.

Lemma 2 follows directly from Definition 3.

Lemma 2. Equation (A3) achieves consensus if its corresponding parameter-varying graphG(ρ)is connectedt, that is,λi,j(ρ)λl>0for(i,j)Ec(ρ)t.

Suppose the system in Eq. (A3) has a piecewise constant output ykm over time interval Δt=tk+1tk, the discrete-time formulation of the system may be written as
xk+1=Φkxk
(A4a)
yk=Cxk
(A4b)

where xkn is the state vector Φkn×n and Cm×n are the system state and output matrices, respectively. Further, for two positive successive integers k,l (lk), let the state transition matrix be defined as Φl|k=Φl|l1Φl1|k where Φk|k=In and Φk+1|k=Φk. Then the system is said to be uniformly detectable if Definition 4 holds.

Definition 4 (Uniform Detectability [35]). The pair(Ak,C)is uniformly detectable if there exist integers q,r0, and some0a<1, b > 0 such that whenever
||Φk+q|kζ||a||ζ||
(A5)
for someζnand k, then
ζTWk+r|kζbζTζ
(A6)
whereWk+r|kis the discrete-time observability gramian given by
Wk+r|k:=i=kk+rΦi|kTCTCΦi|k
(A7)

Lemma 3. (Φk,C)is uniformly detectable if C has at least one nonzero row sum.

Proof. Following Definition 4, we need to evaluate for what v
||Φk+q|kv||||v||a,a[0,1)
(A8)
and show that for such v,
vTWk+r|kvbvTv(b>0)
(A9)
for some integer r0. Since G(ρ) implements consensus, that is, limt+x(t)=x*span{1n}; for any tk0 (tk=kΔt) where x(tk)span{1n},a[0,1) and tk+q with q >0 such that
||x(tk+q)x*||<a||x(tk)x*||
(A10)
Without loss of generality, let v=x(tk)x*. Equation (A10) implies Φk+q|k such that
||Φk+q|kv||||v||<a,vspan{1n}
(A11)
Thus, we only need to show that Eq. (A9) is satisfied with v=1n: In this case, the left-hand side of Eq. (A9)
vTWk+r|kv=vTi=kk+rΦi|kTCTCΦi|kv={1nTCTC1nforr=0r1nTCTC1nelse
(A12)
is, at least, equal to the 2-norm of the row sum(s) of C. Since C has at least one nonzero row sum, some b>0 such that
vTWk+r|kvbn>0
(A13)

Footnotes

1

In Eq. (1b), when j = 1, the term Tj1=T0 is defined as the temperature of the fluid entering the TES, or Tin.

2

The Laplacian of a graph is more formally defined as the difference between the degree matrix and the adjacency matrix; see Ref. [34]

References

1.
Shanks
,
M.
, and
Jain
,
N.
,
2022
, “
Control of a Hybrid Thermal Management System: A Heuristic Strategy for Charging and Discharging a Latent Thermal Energy Storage Device
,”
21st IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (iTherm)
, San Diego, CA, May 31–June 3, pp.
1
10
.10.1109/iTherm54085.2022.9899546
2.
Shafiei
,
S. E.
, and
Alleyne
,
A.
,
2015
, “
Model Predictive Control of Hybrid Thermal Energy Systems in Transport Refrigeration
,”
Appl. Therm. Eng.
,
82
, pp.
264
280
.10.1016/j.applthermaleng.2015.02.053
3.
Pangborn
,
H. C.
,
Laird
,
C. E.
, and
Alleyne
,
A. G.
,
2020
, “
Hierarchical Hybrid MPC for Management of Distributed Phase Change Thermal Energy Storage
,”.
Proceedings of the 2020 American Control Conference (ACC)
, Denver, CO, July 1–3, pp.
4147
4153
.10.23919/ACC45564.2020.9147698
4.
Barz
,
T.
,
Seliger
,
D.
,
Marx
,
K.
,
Sommer
,
A.
,
Walter
,
S. F.
,
Bock
,
H. G.
, and
Körkel
,
S.
,
2018
, “
State and State of Charge Estimation for a Latent Heat Storage
,”
Control Eng. Pract.
,
72
, pp.
151
166
.10.1016/j.conengprac.2017.11.006
5.
Zsembinszki
,
G.
,
Orozco
,
C.
,
Gasia
,
J.
,
Barz
,
T.
,
Emhofer
,
J.
, and
Cabeza
,
L. F.
,
2020
, “
Evaluation of the State of Charge of a Solid/Liquid Phase Change Material in a Thermal Energy Storage Tank
,”
Energies
,
13
(
6
), p.
1425
.10.3390/en13061425
6.
Paberit
,
R.
, and
Öjerborn
,
J.
,
2016
, “
Detecting State of Charge in PCMs - Experimental Investigation of Changes in Chemical and Physical Properties During Phase Transitions
,” Master's thesis,
Chalmers University of Technology
,
Göteborg, Sweden
.
7.
Beyne
,
W.
,
Couvreur
,
K.
,
T'Jollyn
,
I.
,
Lecompte
,
S.
, and
De Paepe
,
M.
,
2022
, “
Estimating the State of Charge in a Latent Thermal Energy Storage Heat Exchanger Based on Inlet/Outlet and Surface Measurements
,”
Appl. Therm. Eng.
,
201
, p.
117806
.10.1016/j.applthermaleng.2021.117806
8.
Kalman
,
R. E.
,
1960
, “
A New Approach to Linear Filtering and Prediction Problems
,”
ASME J. Basic Eng.
,
82
(
1
), pp.
35
45
.10.1115/1.3662552
9.
Luenberger
,
D.
,
1966
, “
Observers for Multivariable Systems
,”
IEEE Trans. Autom. Control
,
11
(
2
), pp.
190
197
.10.1109/TAC.1966.1098323
10.
Henze
,
G.
,
Kalz
,
D.
,
Liu
,
S.
, and
Felsmann
,
C.
,
2005
, “
Experimental Analysis of Model-Based Predictive Optimal Control for Active and Passive Building Thermal Storage Inventory
,”
HVACR Res.
,
11
(
2
), pp.
189
213
.10.1080/10789669.2005.10391134
11.
Steinmaurer
,
G.
,
Krupa
,
M.
, and
Kefer
,
P.
,
2014
, “
Development of Sensors for Measuring the Enthalpy of PCM Storage Systems
,”
Energy Procedia
,
48
pp.
440
446
.10.1016/j.egypro.2014.02.052
12.
Charvát
,
P.
,
Štětina
,
J.
,
Mauder
,
T.
, and
Klimeš
,
L.
,
2017
, “
Visual Monitoring of the Melting Front Propagation in a Paraffin-Based PCM
,”
EPJ Web Conf.
,
143
, p.
02042
.10.1051/epjconf/201714302042
13.
Ezan
,
M. A.
,
Çetin
,
L.
, and
Erek
,
A.
,
2011
, “
Ice Thickness Measurement Method for Thermal Energy Storage Unit
,”
J. Therm. Sci. Technol.
,
31
(
1
), pp.
1
10
.10.5072/ZENODO.31894
14.
Pernsteiner
,
D.
,
Schirrer
,
A.
,
Kasper
,
L.
,
Hofmann
,
R.
, and
Jakubek
,
S.
,
2021
, “
State Estimation Concept for a Nonlinear Melting/Solidification Problem of a Latent Heat Thermal Energy Storage
,”
Comput. Chem. Eng.
,
153
, p.
107444
.10.1016/j.compchemeng.2021.107444
15.
Jaccoud
,
B. R.
,
Orlande
,
H. R. B.
,
Colaço
,
M. J.
,
Fudym
,
O.
, and
Caldeira
,
A. B.
,
2018
, “
State Estimation for the Thermal Storage in Phase Change Materials Containing Nanoparticles
,”
High Temp.-High Pressures
,
47
(
2
), pp.
117
137
.https://www.webofscience.com/wos/woscc/fullrecord/WOS:000430789900002
16.
Morales Sandoval
,
D. A.
,
De La Cruz Loredo
,
I.
,
Bastida
,
H.
,
Badman
,
J. J. R.
, and
Ugalde-Loo
,
C. E.
,
2021
, “
Design and Verification of an Effective State-of-Charge Estimator for Thermal Energy Storage
,”
IET Smart Grid
,
4
(
2
), pp.
202
214
.10.1049/stg2.12024
17.
Mracek
,
C.
,
Clontier
,
J.
, and
D'Souza
,
C.
,
1996
, “
A New Technique for Nonlinear Estimation
,”
Proceeding of the 1996 IEEE International Conference on Control Applications
, Dearborn, MI, Sept. 15–18, pp.
338
343
.10.1109/CCA.1996.558760
18.
Jaganath
,
C.
,
Ridley
,
A.
, and
Bernstein
,
D.
,
2005
, “
A SDRE-Based Asymptotic Observer for Nonlinear Discrete-Time Systems
,”
Proceedings of the 2005 American Control Conference
, Vol.
5
, Portland, OR, June 8–10, pp.
3630
3635
.10.1109/ACC.2005.1470537
19.
Berman
,
A.
,
Zarchan
,
P.
, and
Lewis
,
B.
,
2014
, “
Comparisons Between the Extended Kalman Filter and the State-Dependent Riccati Estimator
,”
J. Guid., Control, Dyn.
,
37
(
5
), pp.
1556
1567
.10.2514/1.G000332
20.
Gelb
,
A.
,
1974
,
Applied Optimal Estimation
,
The Analytic Sciences Corporation, The MIT Press
,
London
.
21.
Julier
,
S. J.
, and
Uhlmann
,
J. K.
,
1997
, “
A New Extension of the Kalman Filter to Nonlinear Systems
,”
Signal Processing, Sensor Fusion, and Target Recognition VI, Proc. SPIE
,
3068
, pp.
182
193
.10.1117/12.280797
22.
Wan
,
E. A.
, and
van der Merwe
,
R.
,
2001
, “
The Unscented Kalman Filter
,”
Kalman Filtering and Neural Networks
,
Wiley Ltd.
,
New York
, pp.
221
280
.
23.
Beikzadeh
,
H.
, and
Taghirad
,
H. D.
,
2012
, “
Exponential Nonlinear Observer Based on the Differential State-Dependent Riccati Equation
,”
Int. J. Autom. Comput.
,
9
(
4
), pp.
358
368
.10.1007/s11633-012-0656-y
24.
Ewing
,
C.
,
2000
, “
An Analysis of the State Dependent Riccati Equation Method Nonlinear Estimation Technique
,” AIAA Paper No.
2000
4275
.10.2514/6.2000-4275
25.
Simon
,
D.
,
2006
,
Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches
,
Wiley
,
Hoboken, NJ
.
26.
Hale
,
D. V.
,
Hoover
,
M. J.
, and
O'Neill
,
M. J.
,
1971
, “
Phase Change Materials Handbook
,” Report No. NASA-CR-61363.
27.
Gohil
,
K. N.
,
Deckard
,
M.
,
Shamberger
,
P. J.
, and
Jain
,
N.
,
2020
, “
A Reduced-Order Model for Analyzing Heat Transfer in a Thermal Energy Storage Module
,”
19th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm)
, Orlando, FL, July 21–23, pp.
681
689
.10.1109/ITherm45881.2020.9190427
28.
Shanks
,
M.
,
Shoalmire
,
C. M.
,
Deckard
,
M.
,
Gohil
,
K. N.
,
Lewis
,
H.
,
Lin
,
D.
,
Shamberger
,
P. J.
, and
Jain
,
N.
,
2022
, “
Design of Spatial Variability in Thermal Energy Storage Modules for Enhanced Power Density
,”
Appl. Energy
,
314
, p.
118966
.10.1016/j.apenergy.2022.118966
29.
Shamberger
,
P. J.
, and
Fisher
,
T. S.
,
2018
, “
Cooling Power and Characteristic Times of Composite Heatsinks and Insulants
,”
Int. J. Heat Mass Transfer
,
117
, pp.
1205
1215
.10.1016/j.ijheatmasstransfer.2017.10.085
30.
Tamraparni
,
A.
,
Hoe
,
A.
,
Deckard
,
M.
,
Zhang
,
C.
,
Elwany
,
A.
,
Shamberger
,
P. J.
, and
Felts
,
J. R.
,
2021
, “
Design and Optimization of Lamellar Phase Change Composites for Thermal Energy Storage
,”
Adv. Eng. Mater.
,
23
(
1
), p.
2001052
.10.1002/adem.202001052
31.
Gillis
,
B.
, and
Jain
,
N.
,
2021
, “
Numerical Validation of Effective Specific Heat Functions for Simulating Melting Dynamics in Latent Heat Thermal Energy Storage Modules
,”
20th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm)
, San Diego, CA, June 1–4, pp.
544
550
.10.1109/ITherm51669.2021.9503136
32.
Yang
,
H.
, and
He
,
Y.
,
2010
, “
Solving Heat Transfer Problems With Phase Change Via Smoothed Effective Heat Capacity and Element-Free Galerkin Methods
,”
Int. Commun. Heat Mass Transfer
,
37
(
4
), pp.
385
392
.10.1016/j.icheatmasstransfer.2009.12.002
33.
Sgreva
,
N. R.
,
Noel
,
J.
,
Métivier
,
C.
,
Marchal
,
P.
,
Chaynes
,
H.
,
Isaiev
,
M.
, and
Jannot
,
Y.
,
2022
, “
Thermo-Physical Characterization of Hexadecane During the Solid/Liquid Phase Change
,”
Thermochim. Acta
,
710
, p.
179180
.10.1016/j.tca.2022.179180
34.
Inyang-Udoh
,
U.
,
Shanks
,
M.
, and
Jain
,
N.
,
2022
, “
A (Strongly) Connected Weighted Graph is Uniformly Detectable Based on Any Output Node
,” e-print arXiv:2209.13119.10.48550/arXiv.2209.13119 Focus to learn more
35.
Anderson
,
B. D. O.
, and
Moore
,
J. B.
,
1981
, “
Detectability and Stabilizability of Time-Varying Discrete-Time Linear Systems
,”
SIAM J. Control Optim.
,
19
(
1
), pp.
20
32
.10.1137/0319002
36.
Shampine
,
L. F.
, and
Reichelt
,
M. W.
,
1997
, “
The MATLAB ODE Suite
,”
SIAM J. Sci. Comput.
,
18
(
1
), pp.
1
22
.10.1137/S1064827594276424
37.
Barz
,
T.
, and
Emhofer
,
J.
,
2021
, “
Paraffins as Phase Change Material in a Compact Plate-Fin Heat exchanger - Part I: Experimental Analysis and Modeling of Complete Phase Transitions
,”
J. Energy Storage
,
33
, p.
102128
.10.1016/j.est.2020.102128
38.
Cao
,
L.
,
Zheng
,
Y.
, and
Zhou
,
Q.
,
2011
, “
A Necessary and Sufficient Condition for Consensus of Continuous-Time Agents Over Undirected Time-Varying Networks
,”
IEEE Trans. Autom. Control
,
56
(
8
), pp.
1915
1920
.10.1109/TAC.2011.2157393