Discrete-Time Kalman Filter Design for Linear Infinite-Dimensional Systems

As the optimal linear filter and estimator, the Kalman filter has been extensively utilized for state estimation and prediction in the realm of lumped parameter systems. However, the dynamics of complex industrial systems often vary in both spatial and temporal domains, which take the forms of partial differential equations (PDEs) and/or delay equations. State estimation for these systems is quite challenging due to the mathematical complexity. This work addresses discrete-time Kalman filter design and realization for linear distributed parameter systems. In particular, the structural- and energy-preserving Crank–Nicolson framework is applied for model time discretization without spatial approximation or model order reduction. In order to ensure the time instance consistency in Kalman filter design, a new discrete model configuration is derived. To verify the feasibility of the proposed design, two widely-used PDEs models are considered, i.e., a pipeline hydraulic model and a 1D boundary damped wave equation.


Introduction
State-of-the-art advancements in the realm of industrial process operations require good understanding of complex phenomena and systems. As a representative method, first principle modelling has been extensively utilized for process description, control, and operation due to its capability of being mathematically provable and physically explainable compared to data-driven methods. Among this category, lumped parameter system (or finite-dimensional systems) modelling has been widely applied to temporal dynamics modelling described by ordinary differential equations (ODEs). However, many complex industrial processes often vary along both temporal and spatial domains. For instance, material and flow processing, transporting, and reacting usually take the forms of integral and partial different equations (PDEs), which have states dependent on both time and space, leading to the realm of distributed parameter systems (DPSs) or infinite-dimensional systems [1][2][3]. Given that DPS is capable of taking both temporal and spatial variables into account, it has been widely utilized in the mathematical modelling and description of various processes and applications, including chemical [2,4], mechanical [5,6], biotechnological [7][8][9], and biomedical [10][11][12][13] engineering. More applications of distributed parameter systems can be found in [14][15][16]. Compared to lumped parameter systems, dealing with distributed parameter systems is not straightforward due to the mathematical modelling complexity induced by PDEs and/or delay equations, which pose further challenges for the ensuing estimator, controller, and regulator design.
State information is usually required for a state feedback controller and/or regulator design, whilst full state information is not often available either due to physical constraints of sensor installation or the prohibitive expense of implementing spatially-distributed sensors, if possible to be built. Hence, state estimation for distributed parameter systems attracts a lot of attention from academia and industry. water hammer dynamics is associated with the recent findings that negative pressure propagation of the transients can result in contamination of portable water systems, which implies that demand for better understanding, efficient monitoring, and regulation of this phenomenon is urgent [60][61][62]. In particular, the hydraulic transients are critical design factors in a large number of fluid flow systems spanning from automotive fuel injection processes to large-scale water supply, transmission, and distribution systems [63]. Along with the pipeline developments, the mechanically-regulated fluid control devices, including many types of pumps and valves, in synergy with sophisticated electronic digital sensors and embedded digital controls, contribute to the complex process behaviour. To better monitor the potential water hammer phenomena and ensure the pipeline integrity, this manuscript addresses the state estimation problem of the stochastic water hammer equations via a discrete Kalman filter. More specifically, there are various versions of water hammer equations; for simplicity, the basic version with the form of linear bi-directional coupled hyperbolic PDEs is discussed in detail. As for the wave equation, it has been extensively exploited in mathematical modelling of mechanical and/or light waves. The nature of the wave equation is given by a second-order-in-time linear hyperbolic PDE. Among existing contributions of the wave equation, an adaptive boundary control strategy was proposed for unstable wave equation with a backstepping method in [64]. Guo and his co-workers investigated output feedback stabilization and disturbance rejection control of the 1D anti-stable wave equation with [65,66]. In [67], a derivative-free nonlinear Kalman filter was designed for state estimation and fault diagnosis in distributed parameter systems. However, an early lumping method was utilized in the study with spatial order reduction, leading to a spatial approximation design. In our work, the novelty arises in a discrete-time infinite-dimensional Kalman filter design for the wave equation without any spatial approximation or order reduction.
Motivated by the superiority of condition monitoring with digital devices to account for the discrete nature of monitoring realization, we develop discrete DPS models of first-order coupled hyperbolic PDEs and second-order hyperbolic PDEs, which account for the infinite-dimensional nature of the pipeline hydraulic process and the wave equation. In particular, we deploy the Cayley-Tustin discretization framework so as to transfer the continuous PDEs into their discrete counterparts by avoiding model reduction and/or approximation. In this way, the discretized model preserves its structural properties (such as: stability, controllability, and observability); see [68,69]. Along the line of time discretization schemes for distributed parameter systems, we derive an appropriate discrete setting that is consistent with the standard finite-dimensional state-space setting associated with the well-known Kalman filter design. Finally, the extension of finite-dimensional Kalman filter design is provided in a closed form of a discrete infinite-dimensional distributed parameter system describing the water hammer and wave dynamics and accounts for naturally present noise in both output measurements and PDE states. The setting provides an obvious advantage since all discrete operators are ultimately bounded and there are no technical difficulties associated with the state space setting and filter realizations.
The rest of the manuscript is structured as follows: In the Mathematical Preliminaries Section, a general infinite-dimensional system framework is described and discretized with the aid of the Cayley-Tustin time discretization approach. After that, the standard, discrete-time two-step Kalman filter is realized for infinite-dimensional systems in Section 3, followed by two realistic process examples in Sections 4 and 5. In Section 4, the water hammer equation is introduced, and linearization, boundary transformations, and time discretization with explicit representation of the resolvent operators of discrete water hammer equations are elaborated, followed by an associated numerical study. In Section 5, the wave system is described, discretized, and simulated with filter implementation. Finally, some conclusions are drawn in Section 6.

Mathematical Preliminaries
In this section, a general linear infinite-dimensional system structure is considered. Based on that, the well-known Cayley-Tustin method is utilized for model time discretization without spatial approximation or spatial order reduction. For the sake of time instance consistency associated with the model description utilized in classical Kalman filter design (x k ) → (y k ), a new Cayley-Tustin time discretization framework is derived and applied to the output equation, resulting in novel expressions of C d and D d operators.

Model Description
Let us consider a linear infinite-dimensional system given by the following general form: where ζ ∈ [0, l] and t ∈ [0, ∞) represent spatial and temporal coordinates, and the state x(·, t) ∈ X , with X = L 2 ([0, l]; R) defined as a separable Hilbert space. The input is u(t) ∈ L 2 loc ([0, ∞), U) and the output y(t) ∈ L 2 loc ([0, ∞), Y), with U and Y given as real Hilbert spaces. In this paper, let us assume that A is closed, leading to the domain D(A), also being a Hilbert space equipped with the graph norm x 2 D(A) := x 2 X + Ax 2 X , and the corresponding resolvent set ρ(A) is not empty, with the norm given by x D(A) = (α − A)x X with an arbitrary constant α ∈ ρ(A). One can denote X 1 := D(A) and apply the norm (α − A)x X , which implies (α − A) −1 projecting X isometrically to X 1 and results in the space X −1 as the completion of X with respect to the norm x X −1 = (α − A) −1 x X . Denoting a block operator node S := A&B C&D as a function X × U → X × Y, we can define the operator node on (X , U, Y) as below: (1) A is a densely-defined and closed operator on X with the resolvent set ρ(A) = (4) C&D ∈ L(D(S), Y); see [57,70,71]. The mapping S : u(·) → y(·) is to project input to output, and the transfer function G(s) = C(sI − A) −1 B + D is obtained by applying the Laplace transform. In addition, operators B and C stand for input and output operators, which can be spatially distributed or applied at the boundaries. In this manuscript, we denote A * as the adjoint operator of A where Aφ, ψ = φ, A * ψ with φ ∈ D(A) and ψ ∈ D(A * ). For the considered water hammer and wave systems, there is no feedthrough operator, i.e., D = 0.

Model Discretization
Continuous-time distributed parameter systems can be taken as good representations of complex industrial processes, while it is more practical and preferable for designers to come up with discrete-time models for observer, controller, and regulator design when it comes to actual implementation in digital devices.
Along this line of late lumping, this work utilizes the Crank-Nicolson time discretization framework for late lumping of DPSs, which can ensure a symmetric and symplectic numerical integration of linear systems without spatial approximation or model reduction [69,[72][73][74].
Let us revisit the considered linear continuous-time infinite-dimensional system (1) and (2). By using the Crank-Nicolson discretization framework for a time discretization interval ∆t, we discretize the process (1) and (2) as follows: where the discrete-and continuous-time input signals are linked by: ∆t k∆t (k−1)∆t u(t)dt, similar to output approximation. This discretization framework preserves input-to-output mapping in the continuous time setting [69,73]. After simple manipulation, the associated discrete time system is described as below: where the corresponding discrete-time spatial operators are denoted as: and its transfer function is given as In addition, we denote the resolvent operator of A as R(s, A)(·) := (sI − A) −1 (·) with s ∈ ρ(A) and the identity operator I, and R is a mapping from the initial state to the state solution in the Laplace domain, i.e., x(ζ, s) = R(ζ, s)x 0 . Based on that, the discrete operators can be further expressed as: where G(δ) = C(δI − A) −1 B = CR(ζ, δ)B is a transfer function of the associated continuous system (1)-(2) with s evaluated at δ, where δ = 2/∆t, and ∆t is the time discretization interval. Moreover, one can evaluate s = δ in the resolvent operator to get the specific expressions for A d , B d , C d , and D d in Equation (8). During the Cayley-Tustin discretization, it is apparent that the resulting discrete operators preserve the spatial characteristics, and the discretized state x k and x k−1 are spatial functions, which implies that the integration setting is applied in the time domain without spatial approximation and/or model reduction.

Remark 1.
For the sake of consistency in time instances, a new discretization structure is derived by manipulating the original state and output Equations (3) and (4) as below: with the associated discrete operators given as: whereC d andD d can be determined by evaluating s = −δ in the resolvent and transfer function. With Remark 1, one can map the current-step state x(k∆t) into the current-step output y(k∆t), which is suitable for the standard Kalman filter model structure. For the sake of simplicity, we denote x(k∆t), x((k − 1)∆t) as x k and x k−1 , similar for u(k∆t) and y(k∆t) in the ensuing sections.

Stochastic Discrete Model
In this section, the discrete-time infinite-dimensional system is extended to a stochastic framework by adding the measurement and state noises in an affine manner. To ensure the time step consistency, we introduce process noise w k and output measurement noise v k to the novel discrete model (9) and (10) as follows: with the same notations given in Equation (11). In particular, the process noise and output noise are given as independent Gaussian noises, i.e., where δ k,j is the Dirac delta function, i.e., δ k,j = 1, if k = j and δ k,j = 0 otherwise. The process and output noises are added in the affine manner to the discrete-time model since it is more convenient to form a stochastic discrete-time infinite-dimensional model setting compared to the case of adding continuous-time noise signal and then applying discretization. In addition, we define a spatial function Φ(ζ) for the description of the state noise distribution in the spatial domain. For simplicity, we take Φ(ζ) = I(ζ), where I denotes an identity function.

Discrete-Time Kalman Filter Design
In this section, a discrete-time two-step infinite-dimensional Kalman filter is constructed for state estimation of the stochastic discrete model (12) and (13), including a prediction step and a updating step, also referred as a priori estimation and a posteriori estimation [75]. In order to proceed to the design of the two-step Kalman filter, we denote that for a priori estimation, one has measurement up to k − 1 steps, i.e.,x − k = E(x k |y 1 , y 2 , ..., y k−1 ),P − k = E(P k |y 1 , y 2 , ..., y k−1 ); in a posteriori estimation, one has measurement up to k steps, i.e.,x + k = E(x k |y 1 , y 2 , ..., y k ),P + k = E(P k |y 1 , y 2 , ..., y k ). In addition, we assume that all input information {u k } is available at any given time instant. As the first step, we estimate the initial condition of statex + 0 and covariance P + 0 as follows: where the estimated initial statex + 0 is independent of process and output noises. Based on the output measurement up to time step k − 1, the prediction step is given as below: With the output measurement y k accumulated to time step k, applying the recursive least squares estimation leads to the updating step as follows [75]: This basic framework (14)- (20) extends a standard discrete-time finite-dimensional Kalman filter design algorithm [75]. Notably, some novelties arise in the following aspects. Compared to the scalar matrix (A, B, C, D) for finite-dimensional systems, this design introduces spatial operators (A d , In addition, the covariance P − k or P + k is two-dimensional and self-adjoint with spatial characteristics, and the evolution given by Equations (16)-(18) preserves its self-adjoint and semi-positive definiteness because of the positive definiteness of Q k and R k . Similarly, Kalman filter gain K k is endowed with spatial characteristics compared to its counterpart in lumped parameter system cases. Moreover, the treatment in the input time instance is also different. The current input u k (instead of u k−1 in the standard Kalman filter design) is mapped to the formulation of the prior estimated statex − k induced by state evolution dynamics given by Equation (17).

Single-Phase Pipeline Hydraulic Model
In this section, the so-called water hammer PDEs are used for pipeline hydraulic modelling in the continuous-time setting. In order to discretize the water hammer model, the Crank-Nicolson framework is utilized for time discretization without spatial approximation or spatial model order reduction. To determine the associated discrete operators, the resolvent operator is solved in a closed analytical form. Based on that, a discrete water hammer system is obtained in the discrete-time setting, and it is augmented with the measurement and state noise signals, which is amenable for discrete Kalman filter design [76].

Model Description
In this section, a one-dimensional liquid flow is considered, and the flow is assumed to be viscous, locally compressible, and isothermal, which implies that the energy loss due to the friction effects and pressure changes can be neglected. In addition, the pipe is assumed to be a buried straight pipe without any elbow. By applying continuity balance law, momentum balance law, and the speed of sound of liquid flow, the following water hammer system is formulated [77]: with associated boundary conditions given as follows: where v(ζ, t), p(ζ, t), and ρ(ζ, t) represent mass flow velocity, pressure, and liquid density, respectively. t and ζ are time and the spatial coordinate along the pipeline. u(t) stands for the boundary control action on mass flow velocity at the upstream end. The multiplication term gsinα represents the ζ-component of the gravity acceleration g; a stands for the velocity of sound in the liquid flow; λ is a dimensionless friction coefficient given by the Colebrook equation (25) for turbulent flow (with a Reynolds number larger than 2320) and Equation (26) for laminar flow (with a Reynolds number smaller than 2320) [59,78].
where D represents the pipe inner diameter, k R D is defined as the relative pipe roughness, and Re stands for the Reynolds number. The system of Equations (21)-(23) is often referred as the water hammer equations [77]. In the following sections, the numerical value of λ is adopted from [77], and it is illustrated in Table 1 with other pipeline parameters, as well.  (21)-(23), one can easily obtain the following nonlinear infinite-dimensional pipeline system:

Model Linearization
In this section, a positive mass flow velocity v is firstly considered. In order to linearize the given nonlinear system of Equations (27) By using the pipeline parameters in Table 1 and assuming the incoming liquid flow with boundary conditions of v 0 = 2.05 m/s, p 0 = 1.68 × 10 6 Pa, and ρ 0 = 680 kg/m 3 , one can determine the steady states profiles as illustrated in Figure 2. Then, a linearized water hammer system is obtained by applying the change of variables as Θ(ζ, t) = Θ(ζ, t) − Θ ss (ζ, t) (where Θ := v, p, or ρ) and the proceeding Taylor series expansion as follows: with associated boundary conditions described as: In the above derivation, a two-state system (pressure and velocity) is obtained by consideration of p =ρa 2 . For simplicity, a matrix form is further obtained as follows: where all notations are defined as follows: As shown in Figure 2, the spatial derivative of steady velocity ∂v ss ∂ζ is close to zero (spatial variation across the space is small enough), so T can be treated as zero. Compared to the large numerical value of N = a 2 ρ ss , the other terms L and R can be approximately neglected as zero. Along this line, a simplified matrix form is yielded and an output operator C is formulated to complete the infinite-dimensional water hammer system representation: , and φ 2 (ζ) being absolutely continuous (abs. con.)}, and continuous-time spatial operators are defined as: (40) where δ(ζ) and δ(ζ − L d ) are Dirac delta functions, while the corresponding notations are simplified as below: (27)- (29), the same linearization steps (30)-(36) can be replicated, and the linear water hammer model structure is identical to Equations (37) and (39) with the same notations of N and G given in Equation (41), but different expressions of Q and F, described as follows:

Remark 2. As for negative velocity in Equations
Up to now, a linearized infinite-dimensional continuous-time water hammer system with boundary control has been derived. To ensure the modelling consistency with the stochastic discrete model (12) and (13), the boundary control water hammer system is converted into an abstract in-domain control model in the next section.

In-Domain Control Problem
In order to realize Kalman filter design, an in-domain control model is obtained from the boundary control model (37)- (40). In order to construct the spatial distribution function B that achieves the boundary to in-domain model transformation, the following theorem is utilized [3].
Then, to convert the boundary control model into the associated in-domain control counterpart with respect to new states [x(ζ, t),ȳ(ζ, t)] , we consider the following system: It is apparent thatĀ has the same expression asL, but with a different domain defined as , and φ 2 (ζ) being abs. con.}.
By solving the inner product equation Aφ 1 , φ 2 = φ 1 , A * φ 2 , one can readily get: with the associated domain defined as D( and φ 2 (ζ) are abs. con.}. Along this line, one can utilize Theorem 1 to yield the following expression: where the boundary control operator is defined as: Up to now, a standard linearized continuous-time infinite-dimensional water hammer system (45)- (48) with in-domain control has been obtained.

Model Discretization
In order to realize a corresponding discrete-time pipeline system in the form of Equations (9) and (10) with the notations given in Equation (11), the resolvent operator needs to be determined beforehand for the considered infinite-dimensional pipeline system representation, as an important link to connect continuous-and discrete-time models. Since the resolvent operator is only related to the A 1 operator in Equation (45), one can directly apply Laplace transformation to Equation (45) without consideration of the input injection operator B as follows: where: Furthermore, one can define e Mζ in the following form: It is simple to solve for e Mζ in the following analytic expression: (59) where ϑ and are denoted as follows: Finally, the distributed parameter pipeline system can be manipulated in the following compact form: By substituting boundary conditionsp(L d , s) = 0 andv(0, s) = 0, one can solve for the expression of p(0, s) such that the resolvent operators are generated as follows: where: (63) with all notations given by Equations (59) and (60). Up to now, a linear infinite-dimensional discrete-time water hammer system has been derived. In the same manner as Equations (12) and (13), the corresponding stochastic linear infinite-dimensional discrete-time water hammer system can be yielded.

Simulation Results of the Water Hammer Equation
In this section, the developed discrete-time Kalman filter design procedures (14)- (20) are applied to the resulting stochastic linear infinite-dimensional discrete-time water hammer system. In addition, the temporal and spatial intervals are 8 s and 500 m. As for the plant and measurement noises, we take with R k = 0.01 and Q k = 0.001. In particular, the initial conditions and the control action are given as follows: Based on that, the stochastic model dynamic and the associated filtering performance are simulated, as shown in Figures 3-7. The transient velocity and pressure are considered here, and it is apparent that the profiles roughly follow the periodic wave trend, induced by the input action, but they are quite noisy in Figure 4. After applying the developed discrete Kalman filter, one can easily see that the noise has been filtered out, and the original dynamic evolution is revealed, which verifies the effectiveness of the proposed discrete Kalman filter, as shown in Figures 3 and 5. As shown in Figure 6, the estimation errors are relatively small compared to the real state magnitudes. From the output filtering perspective in Figure 7, one can readily observe that the noisy pressure and velocity at boundary ends are quite noisy, as shown in the blue dashed lines, and the proposed Kalman filter can smooth out the noisy output profiles and make the filtered output converge to real output profiles, as shown in red dashed lines and green solid lines. In addition, the overall simulation time is 160 s. In this manuscript, the simulations are carried out based on MATLAB R2017a on a MacBook Pro 2017 (processor: 3.1 GHz, IntelCore i5' memory: 8 GB 2133 MHz LPDDR3; operation system: Mac OS Sierra 10.12.06).
To account for the potential model mismatch, we consider the effect of process noises on output filtering performance. More specifically, we simulate Q k = [0.0005 : 0.0005 : 0.005] under the same measurement noise with a fixed covariance R k = 0.01. To quantity the estimation accuracy, we calculate the root-mean-squared-error (RMSE) based on the real (noise-free) outputs. To account for the randomness of the generated noises, we simulate ten times and calculate the average value of RMSE. As shown in Figure 8, we can observe that the estimation RMSE increases with the process noise Q k increase. A proper selection of Q k is in the range of [0.0005,0.0015] where the RMSE of pressure and velocity is relatively small.

Model Description
In this section, we consider a damped wave equation varying within a 1D spatial domain z with a unit length, boundary control, and observation as follows: where κ, ρ, and T represent the model parameter, mass density, and Young's modulus. For the sake of simplicity, we assume these parameters are constants.
Let us revisit Theorem 1 and generate the following abstract in-domain state space model: where: B(·) := −δ(z) 0 (·) (77) C(·) := 1 In addition, the feedthrough operator from the input to the output D is zero.
Proof. It is apparent to show that A has the same expression as L, but a different domain as: Based on Aφ, ψ = φ, A * ψ , one can readily obtain:

Model Discretization
As a key of the realization of discrete operators by using the Cayley-Tustin discretization method, we apply the Laplace transform to the plant (74) as follows: After simple manipulation, this leads to: where: Then, we have the following solution by simple integration: where: Using the boundary conditions of the model (74) and (75) given by D(A), we can solve for x 2 (0, s) = 0 and: Plugging the above Equation (84) into Equation (83) leads to: Then, the resolvent operator is given as below: Finally, one can substitute Equation (87) into Equations (9)-(11) to generate the discrete wave equation system with associated discrete-time operators.

Simulation Results for the Wave Equation
In this section, we proceed with the discrete Kalman filter design for the resulting discretized wave equation system. More specifically, we take the numerical parameters as: ρ = 0.01, T = 255, κ = 0.85. In addition, the temporal and spatial intervals are 1 s and 0.02 m. As for the plant and measurement noises, we adopt w k ∼ N (0, j ] = 0, with R k = 100 and Q k = 0.1. Here, a relatively large value of R is set because of the amplified output induced by the small numerical value of ρ. The initial conditions and input signal are described as below: x 0 (z) = 0.1cos(zπ) 0.1sin(0.5zπ) (88) To elaborate the performance of the developed Kalman filter of the waver system, we present three simulation scenarios here, including: noise-free system, stochastic system, and filtered case. After 40 s of simulation, we emulate the noise-free state evolution, as shown in Figure 9. It is apparent that the two states x 1 (z, t) = ρ ∂w(z,t) ∂t and x 2 (z, t) = ∂w(z,t) ∂z follow sinusoidal wave propagation naturally. In the stochastic state profiles, the actual states are masked by noise, as in Figure 10. By using the discrete-time Kalman filter design, the noise has been filtered out from the stochastic system, leading to original state periodic propagation profiles, which can be seen in Figure 11. The effectiveness of the proposed Kalman filter can also be verified through the small state estimation errors in Figure 12 and output filtering performance in Figure 13 as well. As shown in Figure 13, we can readily see that the output is perfectly smoothed and reconstructed with the existence of process and measurement noises.      To consider the potential model mismatch, we simulate Q k = [0.05 : 0.05 : 0.5] under the same measurement noise with a fixed covariance R k = 100 and utilize RMSE to evaluate the estimation accuracy. As shown in Figure 14, it is apparent that with the process noise Q k increase, the RMSE of estimated output goes up rapidly, where each RMSE value is given based on the averaging of ten simulations. In this case, the estimation RMSE is small enough in the range of Q k = [0.05 : 0.05 : 0.5], which can be applied in real-world applications.  Figure 14. RMSE of the estimated output of the wave system.

Conclusions
In this work, we developed the discrete-time infinite-dimensional Kalman filter design for state estimation of linear distributed parameter systems with boundary controls and point observations. Using the Crank-Nicolson framework, a novel discrete-time distributed parameter system structure was derived and realized without spatial approximation or model order reduction. Based on that, a two-step infinite-dimensional discrete-time Kalman filter was developed with consideration of the spatial distribution of the process noise. To elaborate the filtering performance, two practical models with boundary control action applied were considered in this work, including a linear coupled hyperbolic PDE pipeline system (as a spectral system) and a 1D boundary damped wave model (as a non-spectral system). To account for the effect of process noise on estimation accuracy, we calculated the RMSE of estimated outputs under various simulation scenarios with different process noises and found that under the same measurement noise, with the process noise increase, the estimation RMSE went up. In the future, we will consider the constraint and stability issue and further extend the proposed design into nonlinear DPS models and applications.