Computers and Chemical Engineering

Latent heat thermal energy storages (LHTES) utilize a material’s phase transition to store energy at an almost constant temperature. To fully exploit their high energy density, reliable state estimation is essential, which requires a suitable model-based observer. In previous works


Motivation
More than two-thirds of the energy consumed in the industry is generated by the conversion of fossil fuels and thus contributes to the increase in carbon dioxide emissions, see Napp et al. (2014) . Thermal energy storage systems (TES) are essential to replace these fossil fuels with renewable energy sources despite their volatile nature. TES decouple energy supply and consumption, which leads to greater flexibility and thus to higher performance of the entire energy system. TES can rely on different physical and chemical principles, depending on the required operating conditions, see Gil et al. (2010) . This work focuses on latent heat thermal energy storages (LHTES), which consist of phase change material (PCM) to store energy through phase transition processes at an almost constant temperature. PCM has high energy density but low thermal conduc-tivity, see Agyenim et al. (2010) and Zalba et al. (2003) . To increase the low thermal conductivity, heat-conducting structures are often integrated, as summarized by Tao and He (2018) and Ibrahim et al. (2017) . Such a heat conduction enhancement by aluminum fins is shown in a typical LHTES configuration in Fig. 1 . In this configuration, a section bounded by two adjacent fins is referred to as PCM cell, whose state estimation is the subject of this work. The above-mentioned characteristics of complex energytemperature dependency and low conductivity pose a significant challenge to the state estimation, which in turn is crucial for an efficient implementation of LHTES in industrial applications. The task of the state observer developed within this work is to estimate the distributed system state (temperature distribution and enthalpy content) by only a few measurements on the basis of an accurate model. Moreover, the observation should show high performance despite the expected presence of model errors.

Innovation
In this work, an observer based on the extended Kalman filter (EKF) is developed, which can accurately reconstruct the thermodynamic state of an LHTES using only a few measurements of temperature and other measures such as PCM volume. In contrast to the standard EKF, two different model complexities are employed. The observer uses a nonlinear real-time-capable PCM cell model for prediction. This high-order prediction model is comprised of the Navier-Stokes and the energy equations to represent the dominant effects of natural convection and conduction in the melting/solidification problem. In order to update the state predictions, a linear low-order observer model is incorporated in the EKF. The observer model is obtained by successively linearizing the high-order prediction model and reducing it via balanced truncation to increase the computational efficiency and improve the performance of the observer. In order to account for parameter mismatches and allow a targeted intervention of the observer, the states of the observer model are extended by selected parameters to be estimated, such as the heat transfer coefficient α in between the heat source/sink and the PCM cell. An EKF-based observer uses the system matrices of the low-order observer model as the system equations' Jacobians, resulting in accurate and wellconverging state estimates. The newly proposed observer structure for the nonlinear melting/solidification problem is tested with both, an incorrect initialization and an incorrect model parameterization, under realistic load conditions. The method is able to find the correct state of charge and accurately estimates the melting front in different test scenarios, using only a few temperature sensors and a PCM volume measurement. The observer is tested with simulated measurements generated by a more detailed, but not real-time-capable, high-fidelity model, which produces a slightly different tem perature distribution. Furthermore, the sensor selection is discussed, and a mode-shape-based placement criterion for temperature sensors is provided. The novel contributions of this work are the proposed observer concept for high-order nonlinear melting/solidification problems, its insensitivity against model errors, and the considered simulation study.
In the following literature overview, relevant observer approaches for nonlinear distributed-parameter systems are presented. Then, the state of the art regarding phase change models, which serve as the basis for the observer, is highlighted. Finally, reduction techniques to efficiently tackle high-order systems in observation and control tasks are outlined.

Observers for distributed-parameter systems
Observers/estimators (used synonymously in this work) can be divided into two main categories according to their underlying system dynamics: (1) lumped-parameter systems are described by ordinary differential equations, and (2) distributed-parameter systems are described by partial differential equations and are thus, here, space-and time-dependent. Both types of systems can show linear or nonlinear behavior.
The most well-known observer techniques for linear systems are the Luenberger observer ( Luenberger (1966) ) and the Kalman filter in its discrete and continuous forms ( Kalman (1960) and Kalman and Bucy (1961) ), in which white Gaussian noise for model and measurement errors is considered. Many subsequently developed observers are modifications thereof. For nonlinear system extensions of these observers, see, e.g., Jazwinski (1970) and Wishner et al. (1969) .
In general, observers have been addressed in numerous reviews, see, e.g., Mohd Ali et al. (2015) and Dochain (2003) for chemical process systems. Daum (2005) compares nonlinear filters in a tutorial-style work, while Patwardhan et al. (2012) focussed on the development of Bayesian estimators, both involving different types of the Kalman filter. Afshar et al. (2015) analyzed different observer designs for the heat equation.
There are two main approaches for distributed-parameter systems: early-and late-lumping. For the early-lumping approach, the infinite-dimensional solutions of the system are approximated by a finite-dimensional spatial representation, such as finite elements or other discretization methods, resulting in a system of ordinary differential equations. For the late-lumping approach, the infinite-dimensional character is retained to analyze the system but also to design controllers or observers. In early-lumping, the Kalman filter is a well-known method. An unscented Kalman filter is reviewed by Julier and Uhlmann (2004) , and a broad overview of different unscented Kalman filter formulations is given by Kolås et al. (2009) . Examples for an EKF with unknown inputs, where the unknown inputs are considered as part of the states, can be found in Yang et al. (2007) and Ghahremani and Kamwa (2011) . Another well-known early-lumping observer type is based on the sliding mode approach, which provides robustness to possible discrepancies between the model and the actual system, see, for example, Drakunov and Utkin (1992) . In latelumping, on the other hand, Bitzer and Zeitz (2002) based their observer design on the injection of correction functions into the model equations and boundary conditions. Other promising latelumping designs have been found in the backstepping method, see Jadachowski et al. (2014) .
A comparison of early-lumping, in the form of an unscented Kalman filter, versus late-lumping, using a distributed-parameter observer, is conducted by Kreuzinger et al. (2008) to estimate the state of a stratified storage tank. In Marko et al. (2018) , an early-lumping EKF outperformed two late-lumping observers, in the form of a Lyapunov-based and a backstepping-based design. In this work, an early-lumping approach is employed for state estimation of the phase change problem by first discretizing the characterizing equations via finite elements and then applying an EKF to the discretized system.
In previous approaches for state estimation of phase change processes, Akkari et al. (2006) successfully applied an EKF to microwave defrosting using a finite volume model. Backi et al. (2014) and Backi et al. (2016) presented a state observer for freezing processes based on the EKF and a finite difference model. Barz et al. (2018) designed a nonlinear state observer for a liquid/solid LHTES using an EKF, a collocation method to model the energy balance equations, and temperature sensors. Zsembinszki et al. (2020) summed up the common practices to evaluate the state of charge of liquid/solid LHTES. However, the state estimators developed so far are based on modeling simplifications, such as the neglect of convection. They would therefore not fully describe the dominant system behavior in the present general phase change problem. The observer approach applied in this work can handle high-order models and thus consider all relevant heat transfer effects.
The problem associated with the EKF and high-order systems is, first, computational complexity (see, e.g., Mohd Ali et al. (2015) ) and, second, that most systems are not observable when the full state dimension is considered. Therefore, reduced-order EKFs were the subject of various research projects. Park et al. (2013) developed a reduced-order EKF, in which two EKFs with reduced-order models were operated in parallel, instead of one high-order EKF. Lee et al. (2007) combined the reduced model of the EKF with a measurement noise model and data rejection to account for model errors caused by the model simplification. Khodadadi and Jazayeri-Rad (2011) applied a dual EKF design to estimate states and parameters separately, and turned off the parameter estimator after reaching the optimal values to reduce the computational load. Potocki and Tharp (1993) reduced the model order of the EKF by balanced truncation to estimate the temperature distribution in hyperthermia treatment of cancer. In order to track storms, Farrell and Ioannou (2001) applied a Kalman filter to a high-order model by computing a Kalman gain obtained from a reduced system, which was derived from balanced truncation.
The observer approach presented in this work combines already existing techniques to estimate the states in a phase change process (such as an LHTES) while taking into account all relevant heat transfer effects, such as natural convection and conduction. Therefore, the phase change problem is first described by a real-timecapable high-order nonlinear model already existing in the literature. In order to apply the EKF to the high-order system, reduction via balanced truncation is made. The main procedure is to use the fine-grained, but real-time-capable, nonlinear simulation model to predict and a successively linearized reduced-order (dominant-states only) model to update the simulated model state via an appropriate EKF-based estimator. This approach leads to an increased computational speed and renders the estimation insensitive against model errors. Different types of measurements, such as temperature and volume, are combined in the observer to estimate the state efficiently.

PCM modeling approach
The two most relevant heat transfer mechanisms occurring in low-viscosity PCMs are conduction and natural convection. Dutil et al. (2011) and Liu et al. (2014) review options for mathematical modeling of PCM. Analytical solutions only exist for a limited number of melting/solidification problems, such as the one-dimensional Stefan-problem, see Radhakrishnan and Balakrishnan (1992) . Fortunato et al. (2012) state that the effect of natural convection is widely neglected in modeling PCM thermal storage systems due to its complexity. However, Vogel et al. (2016) and Kasper et al. (2021) found a significant heat transfer enhancement due to convection for low-viscosity PCM depending on the geometry of its enclosure.
When modeling multidimensional heat conduction in PCMs, numerical methods such as the enthalpy method or the effective heat capacity method, summarized in Liu et al. (2014) , must be applied. Nedjar (2002) states that especially finite element methods are able to handle coupled thermomechanical problems with complex boundary conditions. The effective heat capacity method and finite elements were, for example, applied by Tenchev et al. (2005) to solve a conduction and natural convection phase change problem. In a similar approach, Kasper (2020) used the effective heat capacity method and an adaptation of the finite difference code, published by Seibold (2008) , for the twodimensional Navier-Stokes equations in convection modeling. This experimentally validated coupled finite element/difference model was adapted in Kasper et al. (2021) to study the melting and solidification process of different PCM geometries, and in Pernsteiner et al. (2020) , for a co-simulation methodology of a Ruths steam storage surrounded by PCM cells.
The computation of the relevant heat transfer effects with finite element and finite difference methods, as proposed by Kasper (2020) , requires high computational effort and leads to simulation models which typically cannot be computed in or faster than real-time as needed for observation and control tasks. Therefore, Pernsteiner et al. (2021) developed a model reduction approach to short-cut the laborious solution of the Navier-Stokes equations by a data-based stream function model with high accuracy, which is adopted and built on in this work.

Reduction methods for high-order systems
To efficiently model, analyze and simulate systems with high order, model reduction approaches are used to reduce the system dynamics' complexity while maintaining their dominant behavior. Especially model reduction in linear systems is well established, and the research topic has been covered in numerous reviews, see, e.g., Antoulas (2004) and Benner et al. (2015) . The preferred method in this work is the balanced truncation, where the system is transformed to a basis in which the hard-to-reach states are simultaneously difficult to observe and are simply eliminated by truncation to obtain the reduced model. This method was first introduced by Moore (1981) , and a short overview of implementations is given by Gugercin and Antoulas (2004) . Other model reduction methods for linear systems are the Hankelnorm approximation and the singular perturbation methods, see Glover (1984) and Liu and Anderson (1989) , respectively. It is noted that in this work the balanced realization is constructed for linearized implicit system equations with sparse coefficient matrices, for which special numerical tools as in Castagnotto et al. (2017) are available.

Main contributions
To the best of the authors' knowledge, no observer approach has been presented so far to estimate the state of a complex melting/solidification phase change problem, taking into account the effects of conduction and convection. The main contributions of this paper are as follows: • An efficient observer approach for a high-order nonlinear distributed-parameter phase change problem is developed based on a combination of existing modeling, reduction, and observer design methods. • The observer shows high performance and robust behavior, as it can update the system only according to its dominant behavior and combines different types of measurements. • A mode-shape-based sensor placement criterion by analyzing sensitivities is applied and discussed. • Simulation studies are conducted and analyzed, showing the effectiveness and accuracy of the proposed method, even despite model errors. The new approach leads to excellent convergence and is able to find both, the correct state of charge and the melting front.

Paper structure
The paper is organized as follows: In Section 2 , the PCM cell model and its fundamental equations are presented. Section 3 describes the overall observer structure as well as the reduction of the linearized high-order model and the EKF. In Section 4 , a sensor placement strategy using mode shapes is devised. In Section 5 , simulation studies to demonstrate the effectiveness of the proposed observer framework are performed.

PCM cell model
The observer approach is demonstrated using a twodimensional rectangular PCM domain with aluminum fins as heat-conducting structures. For reasons of symmetry, modeling of an LHTES configuration, as illustrated in Fig. 1 , is reduced to a single so-called PCM cell, consisting only of one fin section of the storage, see Fig. 2 . The high-fidelity PCM cell model, developed by Kasper (2020) and described in Section 2.1 , is the basis for the model reduction approach from Pernsteiner et al. (2021) . This reduced nonlinear model, described in Section 2.2 , is real-timecapable with high accuracy and is therefore implemented in the observer to predict the state.

High-fidelity PCM cell model
A detailed multi-phase thermodynamic model, including melting/solidification, heat conduction, and natural convection in two dimensions, is considered for the PCM cell, developed and thoroughly documented in Kasper (2020) . It is based on the conservation laws of energy, mass, and momentum resulting in the energy Eq. (1) , continuity Eq. (2) , and Navier-Stokes Eq. (3) as follows: Therein, the temperature field T = T (x, y, t ) is treated as the dependent variable in the energy Eq. (1) , and the velocity field v = [ u, v ] T with spatial components u (x, y, t ) and v (x, y, t ) is treated as the dependent variable in the Navier-Stokes Eqs.
(2) -(3) , see Chorin (1968) . The variable p PCM represents the pressure in the PCM. The symbols ρ, c, λ, and μ denote the parameters density, apparent heat capacity thermal conductivity, and dynamic viscosity, respectively. The phase transition takes place in a small temperature range referred to as "mushy region " [ T m − ε, T m + ε ] , which is defined by the melting temperature T m and the mushy region width parameter ε. The force density f describes the buoyancy force which is calculated via the Boussinesq approximation by the volumetric thermal expansion coefficient β, the constant (reference) density ρ 0 , a reference temperature T ref , and the gravitational standard acceleration vector g.
The time-dependent energy Eq. (1) is discretized using a standard Galerkin finite element (FE) approach with four-noded bilinear rectangular elements, and the incompressible Navier-Stokes Eqs.
(2) -(3) are solved via finite differences. The obtained velocity field is applied in the next time step of the energy equation.

Boundary conditions
On the upper and lower sides of the PCM cell, the domain boundaries are modeled as adiabatic, For the left and right boundaries, type-3 boundary conditions, also known as Robin or heat-transfer boundary conditions, are prescribed following Newton's law of cooling, where α in , α out are heat transfer coefficients, q q q is the specific heat flux across the boundary, and T in , T out present boundary temperatures at the left and right wall surfaces, respectively.

Regarding the velocity field
Furthermore, the velocity is set to zero wherever the PCM temperature lies below its liquidus temperature, v | T <T m + ε = 0 0 0 . (10)

Reduced nonlinear and real-time-capable model
The model presented above is able to simulate heat transfer by convection and conduction in PCM with high accuracy. Still, it cannot be computed in real-time for adequate problem sizes and on standard computing units. To achieve real-time capability, the model reduction method developed by Pernsteiner et al. (2021) is applied. This reduction method is applicable when the system exhibits simple convective dominant flow patterns. Then, the Navier-Stokes and continuity equations are replaced by a simplified databased stream function model, but the energy equation is fully solved via the FE model, see Fig. 3 . Therein, the state vector x x x k contains the temperatures at the nodes of the FE model for the time step k . The input vector u u u is composed of the temperatures at the boundary surfaces T in ,k and T out ,k , and the approximated velocity field ˆ v v v k −1 is derived from the stream function model.
In order to develop the data-based stream function model, are derived from the velocity field obtained in reference simulation studies using the high-fidelity model from Section 2.1 . In order to consider the solution-dependent flow region due to melting and solidification processes, the stream function is mapped from the original flow domain to a solution-independent unit domain. Then, the stream function snapshots are decomposed into modes of space and time using singular value decomposition (SVD). The temporal behavior V V V red of the dominant spatial modes U U U red can The stream function model can replace the Navier-Stokes and continuity equations. The reduced model is no longer limited by the solver requirements of the Navier-Stokes equations and is therefore robust with respect to mesh and time step enlargement. Thus, a real-time capable model is available that serves as the prediction model in the observer framework and is suitable for reliable state estimation.

Observer framework
In contrast to the standard EKF, two different model complexities are employed. Therefore, the presented observer approach uses a high-order nonlinear, but real-time-capable, prediction model to predict (simulate) the states of the PCM cell one time step ahead. Then, a low-order observer model is derived from the prediction model by linearizing it around the current trajectory (state estimate), subsequently reducing it to its dominant behavior via balanced truncation, and finally augmenting its states by selected parameters to be additionally estimated. The observer model system matrices serve as Jacobians in an EKF. The EKF computes a state update of the prediction model according to the difference between the measurements taken from reality and the predicted model outputs, see Fig. 4 . This structure leads to a reliable state estimation with increased computational speed and improved performance because the observer intervention is restricted to the dominant behavior of the system.
In the following Section, predicted quantities are denoted by the superscript " − " and estimated quantities by "ˆ ", so that the predicted state estimates are denoted by ˆ x x x − , for example. To estimate the state of a PCM cell as described above, the following steps are carried out in the observer structure as illustrated in Fig. 4 .

Measurements
Measurements are taken from an LHTES system (reality) and provided to the observer along with the predicted (simulated) model outputs to compute updated state estimates.

Prediction model
The discrete-time prediction model is simulated forward in time to evolve the model states x x x and outputs y y y , where k is the time step index with a chosen time step length t , so t k = k t , and signals with subscript k denote their values updates at the next time step based on the first-order differential equations of the model, and h ( x x x k ) is the nonlinear function mapping the states to the model outputs y y y k , which correspond to the available measurements. The underlying equations of (13) are based on the reduced nonlinear model described in Section 2.2 , which are solved with a backward Euler scheme. Using the finite-element representation and applying Galerkin's method of weighted residuals (see Pepper and Heinrich (2005) ), the phase change problem (1) -(13) results in a continous-time implicit system of high-order ordinary differential equations of the form with the large but sparse system matrices J J J , K K K and the right-handside vector R R R .

Observer model
The prediction model is nonlinear and of high order m , m ∈ O(10 3 ... 4 ) . To obtain Jacobians for an efficient EKF implementation, a linear low-order observer model is derived from the prediction model by first linearizing it around the current trajectory. Then, balanced truncation reduces the observer model to r remaining dominant modes, r ∈ O(10 1 ) . Finally, the reduced model states are extended by a selected parameter to be estimated. The system matrices of this observer model serve as the Jacobians of the EKF. Since they evolve slowly in the present phase change problem, they are only recomputed every n lin th time steps, n lin ∈ N + . For the observer model system matrices, the time index is omitted for reasons of brevity.

Linearization
The successive linearization (see Zhakatayev et al. (2017) ) of (15) at the current states and inputs x x x (t k −1 ) , u u u (t k −1 ) leads to the continuous-time system in descriptor form, where indicates the deviations from the current trajectory, x x x is the state vector consisting of the temperature field represented by the FE nodal temperatures, and y y y the output vector modeling the measurements. The input vector u u u consists of the temperature at the left wall T in , the temperature at the right wall T out , and a selected parameter to be estimated Θ, such as the potentially unknown heat transfer coefficient α in . In (16) -(17) , the matrices are and the observation matrix C C C c , which maps the states to the outputs. The input matrix B B B c is comprised of The obtained successively linearized model is of high order, x x x ∈ R m ×1 , m ∈ O(10 3 ... 4 ) , with the sparse but large system matrices A A A c ,

Reduction using balanced truncation
The matrices of the linearized system A A A c , B B B c , C C C c , and E E E c from Section 3.3.1 are of high order and thus, as also pointed out by Mohd Ali et al. (2015) , cannot be implemented in an EKF for realtime applications due to its computational complexity. Therefore, the system is reduced via balanced truncation (see Moore (1981) ) to r dominant modes.
In balanced truncation, the system states are first transformed to an energy-balanced basis, in which the states are equally controllable and observable as evident in the corresponding Gramians respectively. In such a transformed representation, the controllability and observability Gramians are equal and correspond to a diagonal matrix, see Brunton and Kutz (2019) . The balanced Gramians' eigenvalues are the system's Hankel singular values and provide an energy measure for each state, so that the high-energy states can be preserved in the reduced system while the low-energy states are eliminated by truncating them entirely. In order to perform the balanced truncation of the descriptorform high-order system (16) -(17) efficiently, it is necessary to exploit the sparsity of the matrices. Therefore, Castagnotto et al. (2017) developed a MATLAB® toolbox to analyze and reduce large-scale linear dynamic systems in sparse state-space representations, see also the sparse algorithms developed by Davis (2006) , Duff et al. (1986) , Saad (2003) , and Gilbert et al. (1992) . In this work, the toolbox is used to preserve the sparsity of the high-order system and exploit it for the computations of the balanced truncation. This is achieved by applying the Petrov-Galerkin projections of the form with the projection matrices ˜ V V V and ˜ W W W T . These projection matrices result from the balanced realization, see Castagnotto et al. (2017) . Subsequently, in the vastly reduced form, the descriptor form is resolved to an explicit state space representation and discretized in time, resulting in the state-space representation can be obtained through the (r × m ) projection matrix ˜ W W W T , which retains the r most controllable and observable modes, ˜ x x x r ∈ R r×1 with r ∈ O(10 1 ) .

Augmentation of the system matrix with unknown inputs
The reduced system states are augmented by a selected parameter Θ that should be estimated along with the states in order to deal with parameter mismatch or uncertainty as well as to introduce a targeted correction option for the observer. Previously, in Section 3.3.1 , the selected parameter to be estimated, such as a potentially unknown heat transfer coefficient α in , was treated as an input to the system. Now the reduced system matrices are reassembled so that the selected parameter to be estimated is subsequently part of the system vector Then, the augmented system matrices have the following form and serve as the Jacobians of the EKF to determine an update for the augmented state vector ˜ x x x re ,k .

Observer
First, in Section 3.4.1 , the formulation of the EKF based on Simon (2006) is recapitulated. Then, aspects of the EKF's temporal evolution across different state-space representations are highlighted in Section 3.4.2 .

Extended Kalman filter formulation
The discrete-time nonlinear system model (13) -(14) with process noise w w w k and measurement noise v v v k , both assumed to be Gaussian uncorrelated and white, is given as It is further assumed that the zero-mean Gaussian noises, have the known covariance matrices Q Q Q k and R R R k . In general, the EKF consists of two steps, the prediction and update step. In the prediction step of the EKF, the predicted state vector and error covariance matrix are determined a priori. Thus, the predicted state vector is computed using the high-order nonlinear prediction model. The predicted error covariance matrix is obtained in the low-order observer model state space using the observer model's system matrix ˜ A A A re as Jacobian of the EKF with the corresponding matrix ˜ Q Q Q . In the update step of the EKF, an estimated state vector and error covariance matrix in the low-order observer model state space are computed a posteriori based on the innovation ˆ y y y k = y y y real , which is defined by the difference between the measurements y y y real taken from reality and the nonlinear prediction model output. The innovation covariance is computed using the matrix ˜ C C C re of the observer model as Jacobian of the EKF and the corresponding matrix ˜ R R R . Together with the Kalman gain the updated state vector estimation of the observer model is obtained. The state estimates of the observer model are used to update the prediction model states after back-transformation into the high-order state representation of the prediction model, see Section 3.5 . The estimation error covariance matrix update of the observer model can be formulated as where I I I is the identity matrix of appropriate size. The updated estimation error covariance evolves over time, which requires special considerations of the different model state-space representations, see Section 3.4.2 .

Consistent covariance mapping over different observer model representations
Since the Jacobians of the EKF change only slowly in the present phase change problem, the successive linearization and thus the observer model matrices are not recomputed in each time step but are considered constant for n lin time steps. With each new linearization step, the observer model and its balanced realization are newly built. Thus, not only the choice of states changes, but also the representation of the corresponding error covariance matrix of the EKF. Therefore, in order to correctly evolve the error covariance matrix across different state representations, it must be transformed between two observer model representations (reduced state choice "1" → "2"). To do so, the covariance matrix ˜ P P P 1 of the representation "1" is first transformed into the common high-order prediction model state space via the projection matrix ˜ V V V 1 , and then brought into the low-order state space of the observer model representation "2" using the projection matrix ˜ W W W 2 . The resulting covariance matrix reads This mapping avoids the actual computation of a large, dense ( m × m ) matrix.

State update of the prediction model
The observer in Section 3.4 computes updated state estimates of the observer model. This updated augmented state vector estimate ˜ x x x re ,k is first decomposed into the reduced states ˜ x x x r ,k and the parameter ˆ Θ k (30) . Then, the reduced state vector estimate ˜ x x x r ,k is transformed from the balanced realization of the observer model back into the original state representation of the prediction model using the (m × r) projection matrix ˜ In this original state representation, both, the states as well as the selected parameter, are updated accordingly, The performance of the developed observer structure is shown in Section 5 . It is seen that for the present system, an observer of this form is robust and converges well. The reduction of states restricts the intervention options of the observer. The observer computes the update of the nonlinear prediction model based on the reduced states and thus can only correct the system according to its dominant behavior.

Sensor selection
In many applications, only a few sensors can be realized for space and cost reasons, so an efficient choice of few sensor locations to maximize observability measures is required.
In the following, a data-based approach is proposed to select s efficient temperature sensor locations based on simulated temperature field time-series data for representative operation scenarios. The sensors and their corresponding measurements, are defined by the selection matrix C, where the rows of C ∈ R s ×m are unit vectors picking out rows of the state vector x x x ∈ R m ×1 . The basic idea of the presented model-free, and thus conservative, approach is to reconstruct the dominant modal coordinates by only measuring few selected states. The sensor positions that explain the dominant modes best over time are obtained by analyzing the corresponding sensitivity matrix.
First, a data matrix containing the temperature field time-series data from representative reference simulations is constructed and decomposed via singular value decomposition (SVD), whereby the state vector x x x k contains the temperatures at the m nodes of the FE model at timestep k , = diag ( σ i ) is the diagonal matrix of singular values σ 1 ≥ σ 2 ≥ . . . ≥ 0 , and U and V contain the left and right singular vectors, respectively. The i th mode shape of the temperature field is found in the i th column of U , starting with the most dominant mode (for i = 1 ). Its normalized modal coordinate time-series is contained in the i th row of V T . Truncating the SVD representation to the first q dominant modes is indicated by U q , q , and V q , respectively.
The basic idea is to find sensors whose information allows a direct reconstruction of the q most dominant modes' modal coordinate values. For a given sensor selection defined by the index vector l = { l 1 , ., l s } , the output selection matrix C reads The resulting sensitivity matrix expresses the gains from the considered modal coordinates to the measured outputs as Finally, a suitable matrix measure on the sensitivity matrix F l is formulated to evaluate the studied sensor combination. To discriminate well between the different modes, the choice that is, the minimal eigenvalue of the sensitivity matrix F l , is an appropriate utility function to be maximized over a considered (large) set of candidate sensor combinations L . The optimal sensor choice of s sensors reads arg max l∈L η(l) .

(53)
This combinatorial optimization problem is solved with appropriate search algorithms, usually in an approximated or heuristic way.
This data-based sensor selection approach does not utilize modeled system dynamics, so the number of sensors s must be at least as large as the number of reconstructed modes q , i.e., q ≤ s . In this work, the proposed method will be applied for the most important temperature sensor only ( s = 1 ) to ensure high observability of the most relevant mode ( q = 1 ). Moreover, besides temperature measurements, any sensor type with a linear output Eq. (48) can be considered directly in the sensor selection process. In the present application, a PCM volume measurement is found to improve estimation robustness in the presence of model errors. Its output function is, however, nonlinear in x x x and it is not considered in the selection process of temperature sensors.

Demonstration of the observer & discussion
The performance of the developed observer is demonstrated in this section via simulation studies. First, the simulation setup is defined. Then, the observer is tested under two scenarios: an incorrect initialization and an incorrect parameterization of the model, respectively.

Simulation setup
The setup for the simulation studies and the parameters of the observer are specified below. The geometry of the PCM cell, the material parameters, and the discretization correspond to the values in Pernsteiner et al. (2021) .

Geometry
The geometry of the PCM cell is defined in Fig. 2 and Table 1 . The heated wall (left) side of the PCM cell has good heat transfer properties for charging and discharging via heat flows. On the outer wall (right) side, the PCM cell is isolated to reduce heat loss to the environment, see Fig. 2 . The assumed heat transfer coefficients are listed in Table 2 .

Material parameters
As in Vogel et al. (2016) , the PCM is chosen as a eutectic mixture of potassium nitrate and sodium nitrate KNO 3 -NaNO 3 , enclosed in an aluminum encapsulation. The PCM is either liquid (L), solid (S), or mushy (M). The material properties of the aluminum encapsulation and the PCM are listed in Table 3 . The density of solid PCM is used for both phases. The thermal expansion coefficient is only used to compute the buoyancy force (Boussinesq approximation) (5) and, together with the assumed expansion due to phase change, to determine the PCM volume.

Mesh and time step size
The geometry of the high-fidelity model is discretized by a mesh of square elements with a side length of x = y = 0 . 5 mm , and the time step is set to t = 0 . 1 s . These mesh resolutions resulted as optimal values from a mesh convergence study, see Kasper et al. (2021) . The reduction method described in Section 2.2 enables a coarser mesh and larger time step. Therefore, the time step and mesh sizes are doubled to x = y = 1 mm and t = 0 . 2 s for the prediction model that is implemented in the observer, see Table 4 .

Observer parameters
The observer is realized as described in Section 3 . The linearization and thus the creation of the observer model is executed every t lin = 60 s . The balanced truncation leads to a model with r = 5 reduced states. This model is then extended by the selected parameter to be estimated, the heat transfer coefficient α in . The volume measurement is scaled by 10 7 . The matrices Q Q Q and R R R in the observer are selected to be diagonal and constant. All diagonal components of the matrix Q Q Q are set to 10 −1 . The diagonal components of the matrix R R R are 10 3 for the temperature measurements and 10 −1 for the volume measurement. The values of Q Q Q and R R R have been adjusted to balance convergence speed and performance in the presence of model errors.

Load profile
The performance of the observer is demonstrated using a 5hour load profile with varying amplitude. Consequently, the PCM cell is charged/discharged by a heat flow resulting from the temperature difference between its heated (left) wall and the domain interior (7) . The load profile's input temperature on the left side T in of (7) is set as shown in Fig. 5 , while the temperature responsible for the heat loss on the outer (right) wall T out of (8) is fixed at 20 • C. The initial temperature in the PCM cell is T 0 = 218 • C , just below the solidus temperature of the PCM.

Optimized sensor selection
The optimized sensor selection is derived from the method presented in Section 4 . The sensor candidate positions are located in the symmetry plane of the PCM cell. One ( s = 1 ) temperature sensor is selected optimally for q = 1 dominant mode. Fig. 6 a shows the temperature sensitivity for the considered load case, and Fig. 6 b depicts the optimized sensor placement (left-most sensor is optimized, three additional temperature sensors are distributed equidistantly due to application constraints).
Simulation studies show that the optimized sensor placement is preferable to the conventional one. During the melting process, a sensor located close to the left heated wall quickly detects the beginning of the melting process and thus provides information to the observer. If the left-most sensor is too far away from this heated wall, however, as in the case of the conventional sensor placement, it remains much longer in the solid region, where it only measures an almost constant temperature (the melting temperature). Therefore the sensor does not provide the observer with significant information until the melting front reaches the first sensor and the temperature changes. These relations are also depicted in the data matrix outlined in Section 4 . The following additional observations were made: • The total enthalpy of the PCM cell and the melting front characteristics are estimated more accurately and with faster convergence with the optimized sensor placement compared to the conventional one.
• The parameter α in can be determined correctly with the optimized sensor placement but not with the conventional one when no model error is assumed.
Therefore, in the following simulation studies, the optimized temperature sensor placement is applied. Additionally, a measurement of the total PCM volume, with total PCM mass M, is assumed and considered in the observer. The additional measurement of PCM volume helps to make the observer insensitive to model errors and extends the field of appli-  Table 3 Material properties of the LHTES, phases in solid (S) and liquid (L) state.   cation to pure substances as well as eutectics, where the melting process takes place without temperature change.

Simulation scenarios
Two different scenarios are evaluated in the simulation study: • wrong initialization -The initial temperature field of the simulated reality is set to a complex temperature distribution with parts in liquid and solid state, see Fig. 7 . • wrong parameterization -The heat transfer coefficient of the simulated reality is set to 900 Wm −2 K −1 .
The measurements used by the observer come from the perturbed high-fidelity model, which hence serves as simulated reality. This model is different from the prediction model implemented in the observer and exhibits a moderately different temperature distribution, especially in the liquid phase of the PCM. This model error occurs in the prediction model due to the substitution of the Navier-Stokes Eqs.
(2) -(3) by a stream function model, which is optimized to accurately reconstruct the correct total enthalpy content as well as position of the melting front but not the exact temperature distribution.
The measurements from the simulated reality are subject to Gaussian noise. For this purpose, the standard deviation of the temperature sensors is assumed to be 1 K and of the volume measurement 1% of the PCM domain size.
The simulated reality allows an analysis of the entire distributed state of the PCM cell to evaluate the performance of the observer.

Results
The results are evaluated under two aspects: the state of charge as well as the location and shape of the melting front. These criteria are compared between the simulated reality, the observed model, and an uncorrected model (i.e., the prediction model simulation only).

State of charge
Knowing the state of charge and its distribution in the PCM cell is essential to control the LHTES. One representative formulation of the state of charge is the stored normalized enthalpy, In (55) , H(t ) is the total enthalpy at time t, H 0 the initial total enthalpy, and H latent the latent heat of the PCM cell. Figure 8 compares the stored normalized enthalpy of the different models for the wrong-initialization and wrong-parameterization scenarios, respectively. The total enthalpy can be estimated with high accuracy for both simulation cases, and only small offsets are seen. Therefore the developed observer is considered suitable to estimate the state of charge.

Melting front
As a second criterion, the shape and progression of the melting fronts in the PCM are evaluated qualitatively, see Fig. 9 for an illustration in the wrong-initialization scenario. The location of the melting front is accurately estimated and its shape is being retained correctly. Such a good estimate of the melting front is essential as it determines, for example, how quickly the stored enthalpy can be retrieved from the PCM cell.
These results are highly promising in terms of robustness and performance, especially considering the large initialization error of the first test scenario and the model error in the observer leading to temperature deviations in the liquid phase between the prediction model and the simulated reality. With such a state observer, it is also possible to dynamically operate the LHTES and integrate it into complex storage systems using modern control methods.

Joint state and parameter estimation
The augmented parameter Θ, in this case the heat transfer coefficient α in , was correctly estimated when using measurements from the real-time-capable nonlinear prediction model. However, when using measurements from the high-fidelity model (the simulated reality), that is, in the presence of a model error, the parameter to be estimated showed oscillating behavior rather than smooth convergence. Still, the augmented parameter can be used by the observer as an additional correction variable and is therefore considered a valuable design tool.

Conclusions
The developed observer framework for melting/solidification phase change problems with conduction/convection is able to handle nonlinear distributed-parameter systems. The method is based on a high-order but real-time-capable finite element model, which is used to predict future states. This prediction model is linearized around the current trajectory and then efficiently reduced via balanced truncation to obtain an observer model with dominant modes only. The low-order observer model's states are augmented by selected parameters that should be estimated along with the states in order to deal with parameter mismatch or uncertainty. The observer model system matrices serve as Jacobians in an extended Kalman filter, and a state update of the prediction model is computed according to the difference between the measurements and the predicted model outputs. Thereby the observer combines temperature and volume measurements.
It is seen that an observer of this form converges well and is insensitive against model errors. The reduction of states restricts the intervention options of the observer. The observer computes the update of the nonlinear prediction model based on the reduced states and thus can only correct the system according to its dominant behavior. The temperature sensor placement is optimized using a proposed mode-shape-based selection criterion.
The observer can be computed in real-time and its performance is demonstrated in simulation studies of two scenarios: a wrong-initialization and a wrong-parameterization case. It is tested with simulated measurements generated by a detailed (but not real-time-capable) high-fidelity reference model, which exhibits a slightly different temperature distribution. The results are satisfactory, especially considering the large initialization error of the first test scenario and the modeling error of the temperature distribution between the prediction model and the simulated reality. The total enthalpy and thus the state of charge can be estimated with high accuracy. Also, the location and shape of the melting front are estimated accurately.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.