## 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) [1–5]. 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 [17–19] 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.

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\xd7ny$ 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 $ny\u22122$ 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* = *n _{x}* from the inlet to the outlet

^{1}. 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$.

*a*) 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\u02d9jadv$ term, defined in Eq. (1

*b*), represents heat transfer in fluid control volumes due to mass transfer, or advection; $m\u02d9$ 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\u02d9i\u2192j$ term, defined in Eq. (1

*c*), represents the conductive (solid-solid) or convective (fluid-solid) heat transfer rate from

*i*to

*j*, where $i\u2208A(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. (1

*d*), which is either conductive or convective. In Eq. (1

*d*), $di,j$ is the distance between the centers of control volumes

*i*and

*j*,

*κ*is the thermal conductivity of

_{j}*j*, $ai,j$ is the area of the boundary between

*i*and

*j*, and

*U*is the convective heat transfer coefficient of fluid control volume

_{j}*j*

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. (1*a*) and (1*d*), 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].

*T*, the effective specific heat of CPCM is greatly increased so that the latent heat is modeled as sensible heat over the temperature range $Tpc\xb1\Delta Tpc2$, with $\Delta Tpc=8$ K, as shown in Fig. 3. The width of this temperature range is determined by the parameter $\alpha =8\Delta 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,

_{pc}*h*

_{fus}, and the solid and liquid specific heats, $cp,sol$ and $cp,liq$, represent properties of the CPCM.

### 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 *n _{x}* and

*n*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

_{y}*n*= 3 by

_{x}*n*= 7 control volumes for a total of

_{y}*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.*

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.▪

*L*(

*x*), in addition to a diagonal capacitance matrix,

*M*(

*x*), and an input matrix,

*B*(

*x*), where the state vector $x=[T1\cdots Tn]\u22a4\u2208\mathbb{R}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

*x*refers to the value of

_{k}*x*at time

*t*, or $xk=x(tk)$. The control input $u\u2208\mathbb{R}1$ is the mass flowrate, and $yk\u2208\mathbb{R}p$ is the (output) vector of states that are measured at time-step

_{k}*k*.

*a*)

^{2}. The matrix $M(x)\u2208\mathbb{R}n\xd7n$ is defined in Eq. (5

*b*) 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)\u2208\mathbb{R}n\xd71$ is defined in Eq. (5

*c*) and contains the advection terms. Finally,

*T*is the temperature of the fluid at the inlet, which is not a state. Recall that

_{in}*T*

_{1},

*T*

_{2}, …, $Tnx$ are the fluid control volume temperatures, which are the first

*n*states in

_{x}*x*.

Suppose measurements of the temperature of certain control volumes are available at discrete time instances *t _{k}*, $k=1,2,\u2026$; then $C\u2208\mathbb{R}p\xd7n$ is a binary output matrix where $Ci,j=1$ if $yk,i=xk,j$, and $Ci,j=0$ otherwise. Thus, the output equation (Eq. (4

*b*)) is discrete-time although the state dynamics are continuous.

### 2.3 State of Charge.

*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.

*T*. The maximum enthalpy

_{pc}*H*

_{max}and minimum enthalpy

*H*

_{min}in Eq. (6) are calculated for the PCM at a uniform temperature

*T*

_{max}= 308 K or

*T*

_{min}= 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.

## 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.

*K*is given by

_{k}Using a smaller time interval, $\Delta 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 *t _{k}*, 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 $(\Phi k,C)$ is uniformly detectable [35]. The TES model given by Eq. (4) represents a continuously connected undirected graph since $Ai,j(x)>0\u2009\u2200t$ for each existing edge (*i*, *j*) where $i\u2260j$. If a parameter-varying undirected graph is continuously connected $\u2200t$, the graph implements consensus, and is, hence, uniformly detectable from any node (see the Appendix). Furthermore, the pair $(\Phi 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.

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.

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.

The thermocouples are sampled at regular intervals to produce the measurement output vector *y _{k}*. Thermocouple TC0 measures

*T*

_{in}, 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.

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.

*a*). The process noise covariance is not as easily characterized, so we assume it is negligible; Eq. (14

*b*) defines the matrix

*W*just large enough to ensure that $Pk+1|k$ remains positive definite.

_{k}## 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.

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 $\Delta 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 *n _{x}* = 21 and

*n*= 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

_{y}*y*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.

_{k}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\u0302(tk)\u2212xsim(tk)$ where *x*_{sim} is a state variable calculated by the simulation.

Even though the temperatures of CV5 and CV6 are not included in the output vector *y _{k}*, 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.

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.

*n*control volumes at time-step

*t*, 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

_{k}*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.*

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.

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\u0302(tk)\u2212yk,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.

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.

### 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.

*N*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.

_{s}Sample rate | $erms,CV1$ (K) |
---|---|

80 S/s | 0.2436 |

10 S/s | 0.2617 |

1 S/s | 0.3811 |

0.2 S/s | 0.7526 |

Sample rate | $erms,CV1$ (K) |
---|---|

80 S/s | 0.2436 |

10 S/s | 0.2617 |

1 S/s | 0.3811 |

0.2 S/s | 0.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.

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.

Sample rate | $erms,CV3$ (K) |
---|---|

80 S/s | 1.0022 |

10 S/s | 1.1046 |

1 S/s | 1.2048 |

0.2 S/s | 1.3270 |

Sample rate | $erms,CV3$ (K) |
---|---|

80 S/s | 1.0022 |

10 S/s | 1.1046 |

1 S/s | 1.2048 |

0.2 S/s | 1.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.

## 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 (m

^{2})*A*=state-space coefficient matrix

- $A$ =
set of adjacent control volumes

*B*=state-space input matrix

*c*=_{p}specific heat ($Jkg\u2212K$)

*C*=state-space output matrix

- $C$ =
set of PCM control volumes

*d*=distance (m)

*e*=estimation error

*e*=_{rms}root-mean-square error

*f*=_{m}melt fraction of a PCM control volume

- $F$ =
set of fluid control volumes

*h*=specific enthalpy $(Jkg)$

*h*=_{fus}specific enthalpy of fusion $(Jkg)$

*H*=enthalpy (J)

*i*,*j*,*l*=control volume indices

*I*=_{n}*n*×*n*identity matrix*k*=discrete time step index

*K*=Kalman gain

*L*=graph Laplacian matrix

*m*=mass (kg)

- $m\u02d9$ =
mass flow rate $(kgs)$

*M*=capacitance matrix

*n*=number of control volumes

*N*=_{s}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\u02d9$ =
heat transfer rate (w)

*R*=thermal resistance $(KW)$

*t*=time (s)

*T*=temperature (K)

*T*=_{pc}phase-change temperature (K)

*u*=input vector

*U*=convective heat transfer coefficient $(Wm2\u2212K)$

*v*=measurement noise vector

*V*=measurement noise covariance

*w*=process noise vector

*W*=process noise covariance

*x*=state vector

- $x\u0302$ =
estimated state vector

*x*=_{soc}state of charge

*y*=output or measurement vector

*α*=specific heat function width parameter

- Γ =
discrete-time input matrix

*κ*=thermal conductivity $(Wm\u2212K)$

- $\Phi $ =
discrete-time state transition matrix

### Appendix

Definition 1 (Connected Undirected Graph). *A graph is a pair*$G=(V,E)$*where*$V$*is is a finite set of n nodes or vertices*${v1,\u2026,vn}$*, and*$E\u2282{1,\u2026,n}\xd7{1,\u2026,n}$*is a set of edges. The pair (i, j) denotes the edge that links vertex v _{i} to v_{j}. If*$(i,j)\u2208E\u21d4(j,i)\u2208E$

*, 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 with*$Ec$.

*A (weighted) graph is fully described by the adjacency matrix,*$\Lambda \u2208\mathbb{R}n\xd7n$

*where*

*where*$\lambda i,j$*is the weight assigned to any of its edges (i, j). Consider the diagonal matrix D whose diagonal entry*$Di,i=\u2211j=1n\lambda 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 as*$L=D\u2212\Lambda $.

*Given a parameter-varying graph*$G(\rho )=(V,E(\rho ))$

*for some time-dependent parameter ρ, the integral graph of*$G(\rho )$

*on*$[0,\u221e)$

*is a constant graph*$G\xaf[0,\u221e):=(V,E\xaf)$

*where*$V$

*is the same vertex set of*$G(\rho )$

*, and the adjacency matrix is defined by*

*x*is given by

_{i}where *L* is the Laplacian. Lemma 1 holds for the system [38].

Lemma 1 ([38]). *The dynamics of (A3) implements consensus, that is,*$limt\u2192+\u221ex(t)\u2208span{1n}$*if and only if*$G\xaf[0,\u221e)$*is connected*.

Lemma 2 follows directly from Definition 3.

Lemma 2. *Equation (A3) achieves consensus if its corresponding parameter-varying graph*$G(\rho )$*is connected*$\u2200t$*, that is,*$\lambda i,j(\rho )\u2265\lambda l>0$*for*$(i,j)\u2208Ec(\rho )\u2009\u2200t$.

where $xk\u2208\mathbb{R}n$ is the state vector $\Phi k\u2208\mathbb{R}n\xd7n$ and $C\u2208\mathbb{R}m\xd7n$ are the system state and output matrices, respectively. Further, for two positive successive integers $k,l$ ($l\u2265k$), let the *state transition matrix* be defined as $\Phi l|k=\Phi l|l\u22121\Phi l\u22121|k$ where $\Phi k|k=In$ and $\Phi k+1|k=\Phi k$. Then the system is said to be uniformly detectable if Definition 4 holds.

*The pair*$(Ak,C)$

*is uniformly detectable if there exist integers q,*$r\u22650$

*, and some*$0\u2264a<1$

*, b > 0 such that whenever*

*for some*$\zeta \u2208\mathbb{R}n$

*and k, then*

*where*$Wk+r|k$

*is the discrete-time observability gramian given by*

Lemma 3. $(\Phi 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*

*v*,

*r*$\u22650$. Since $G(\rho )$ implements consensus, that is, $limt\u2192+\u221ex(t)=x*\u2208span{1n}$; for any $tk\u22650$ ($tk=k\Delta t$) where $x(tk)\u2209span{1n},\u2009\u2203a\u2208[0,1)$ and $tk+q$ with

*q*>

*0 such that*

*C*. Since

*C*has at least one nonzero row sum, $\u2203$ some $b>0$ such that

▪

## Footnotes

In Eq. (1*b*), when j = 1, the term $Tj\u22121=T0$ is defined as the temperature of the fluid entering the TES, or *T*_{in}.

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