Discretization

Discretization is also concerned with the transformation of continuous differential equations into discrete difference equations, suitable for numerical computing.

The following continuous-time state space model

$${\displaystyle {\begin{aligned}{\dot {\mathbf {x} }}(t)&=\mathbf {Ax} (t)+\mathbf {Bu} (t)+\mathbf {w} (t)\\[2pt]\mathbf {y} (t)&=\mathbf {Cx} (t)+\mathbf {Du} (t)+\mathbf {v} (t)\end{aligned}}}$$

where v and w are continuous zero-mean white noise sources with power spectral densities

$${\displaystyle {\begin{aligned}\mathbf {w} (t)&\sim N(0,\mathbf {Q} )\\[2pt]\mathbf {v} (t)&\sim N(0,\mathbf {R} )\end{aligned}}}$$

can be discretized, assuming zero-order hold for the input u and continuous integration for the noise v, to

$${\displaystyle {\begin{aligned}\mathbf {x} [k+1]&=\mathbf {A_{d}x} [k]+\mathbf {B_{d}u} [k]+\mathbf {w} [k]\\[2pt]\mathbf {y} [k]&=\mathbf {C_{d}x} [k]+\mathbf {D_{d}u} [k]+\mathbf {v} [k]\end{aligned}}}$$

with covariances

$${\displaystyle {\begin{aligned}\mathbf {w} [k]&\sim N(0,\mathbf {Q_{d}} )\\[2pt]\mathbf {v} [k]&\sim N(0,\mathbf {R_{d}} )\end{aligned}}}$$

where

$${\displaystyle {\begin{aligned}\mathbf {A_{d}} &=e^{\mathbf {A} T}={\mathcal {L}}^{-1}{\Bigl \{}(s\mathbf {I} -\mathbf {A} )^{-1}{\Bigr \}}_{t=T}\\[4pt]\mathbf {B_{d}} &=\left(\int _{\tau =0}^{T}e^{\mathbf {A} \tau }d\tau \right)\mathbf {B} \\[4pt]\mathbf {C_{d}} &=\mathbf {C} \\[8pt]\mathbf {D_{d}} &=\mathbf {D} \\[2pt]\mathbf {Q_{d}} &=\int _{\tau =0}^{T}e^{\mathbf {A} \tau }\mathbf {Q} e^{\mathbf {A} ^{\top }\tau }d\tau \\[2pt]\mathbf {R_{d}} &=\mathbf {R} {\frac {1}{T}}\end{aligned}}}$$

and T is the sample time. If A is nonsingular, ${\displaystyle \mathbf {B_{d}} =\mathbf {A} ^{-1}(\mathbf {A_{d}} -\mathbf {I} )\mathbf {B} .}$

The equation for the discretized measurement noise is a consequence of the continuous measurement noise being defined with a power spectral density.^[1]

A clever trick to compute A_d and B_d in one step is by utilizing the following property:^{[2]: p. 215}

$${\displaystyle e^{{\begin{bmatrix}\mathbf {A} &\mathbf {B} \\\mathbf {0} &\mathbf {0} \end{bmatrix}}T}={\begin{bmatrix}\mathbf {A_{d}} &\mathbf {B_{d}} \\\mathbf {0} &\mathbf {I} \end{bmatrix}}}$$

Where A_d and B_d are the discretized state-space matrices.

Numerical evaluation of Q_d is a bit trickier due to the matrix exponential integral. It can, however, be computed by first constructing a matrix, and computing the exponential of it^[3] $${\displaystyle {\begin{aligned}\mathbf {F} &={\begin{bmatrix}-\mathbf {A} &\mathbf {Q} \\\mathbf {0} &\mathbf {A} ^{\top }\end{bmatrix}}T\\[2pt]\mathbf {G} &=e^{\mathbf {F} }={\begin{bmatrix}\dots &\mathbf {A_{d}} ^{-1}\mathbf {Q_{d}} \\\mathbf {0} &\mathbf {A_{d}} ^{\top }\end{bmatrix}}\end{aligned}}}$$ The discretized process noise is then evaluated by multiplying the transpose of the lower-right partition of G with the upper-right partition of G: $${\displaystyle \mathbf {Q_{d}} =(\mathbf {A_{d}} ^{\top })^{\top }(\mathbf {A_{d}} ^{-1}\mathbf {Q_{d}} )=\mathbf {A_{d}} (\mathbf {A_{d}} ^{-1}\mathbf {Q_{d}} ).}$$

Starting with the continuous model $${\displaystyle \mathbf {\dot {x}} (t)=\mathbf {Ax} (t)+\mathbf {Bu} (t)}$$ we know that the matrix exponential is $${\displaystyle {\frac {d}{dt}}e^{\mathbf {A} t}=\mathbf {A} e^{\mathbf {A} t}=e^{\mathbf {A} t}\mathbf {A} }$$ and by premultiplying the model we get $${\displaystyle e^{-\mathbf {A} t}\mathbf {\dot {x}} (t)=e^{-\mathbf {A} t}\mathbf {Ax} (t)+e^{-\mathbf {A} t}\mathbf {Bu} (t)}$$ which we recognize as $${\displaystyle {\frac {d}{dt}}{\Bigl [}e^{-\mathbf {A} t}\mathbf {x} (t){\Bigr ]}=e^{-\mathbf {A} t}\mathbf {Bu} (t)}$$ and by integrating, $${\displaystyle {\begin{aligned}e^{-\mathbf {A} t}\mathbf {x} (t)-e^{0}\mathbf {x} (0)&=\int _{0}^{t}e^{-\mathbf {A} \tau }\mathbf {Bu} (\tau )d\tau \\[2pt]\mathbf {x} (t)&=e^{\mathbf {A} t}\mathbf {x} (0)+\int _{0}^{t}e^{\mathbf {A} (t-\tau )}\mathbf {Bu} (\tau )d\tau \end{aligned}}}$$ which is an analytical solution to the continuous model.

Now we want to discretise the above expression. We assume that u is constant during each timestep. $${\displaystyle {\begin{aligned}\mathbf {x} [k]&\,{\stackrel {\mathrm {def} }{=}}\ \mathbf {x} (kT)\\[6pt]\mathbf {x} [k]&=e^{\mathbf {A} kT}\mathbf {x} (0)+\int _{0}^{kT}e^{\mathbf {A} (kT-\tau )}\mathbf {Bu} (\tau )d\tau \\[4pt]\mathbf {x} [k+1]&=e^{\mathbf {A} (k+1)T}\mathbf {x} (0)+\int _{0}^{(k+1)T}e^{\mathbf {A} [(k+1)T-\tau ]}\mathbf {Bu} (\tau )d\tau \\[2pt]\mathbf {x} [k+1]&=e^{\mathbf {A} T}\left[e^{\mathbf {A} kT}\mathbf {x} (0)+\int _{0}^{kT}e^{\mathbf {A} (kT-\tau )}\mathbf {Bu} (\tau )d\tau \right]+\int _{kT}^{(k+1)T}e^{\mathbf {A} (kT+T-\tau )}\mathbf {B} \mathbf {u} (\tau )d\tau \end{aligned}}}$$ We recognize the bracketed expression as ${\displaystyle \mathbf {x} [k]}$, and the second term can be simplified by substituting with the function ${\displaystyle v(\tau )=kT+T-\tau }$. Note that ${\displaystyle d\tau =-dv}$. We also assume that u is constant during the integral, which in turn yields

$${\displaystyle {\begin{aligned}\mathbf {x} [k+1]&=e^{\mathbf {A} T}\mathbf {x} [k]-\left(\int _{v(kT)}^{v((k+1)T)}e^{\mathbf {A} v}dv\right)\mathbf {Bu} [k]\\[2pt]&=e^{\mathbf {A} T}\mathbf {x} [k]-\left(\int _{T}^{0}e^{\mathbf {A} v}dv\right)\mathbf {Bu} [k]\\[2pt]&=e^{\mathbf {A} T}\mathbf {x} [k]+\left(\int _{0}^{T}e^{\mathbf {A} v}dv\right)\mathbf {Bu} [k]\\[4pt]&=e^{\mathbf {A} T}\mathbf {x} [k]+\mathbf {A} ^{-1}\left(e^{\mathbf {A} T}-\mathbf {I} \right)\mathbf {Bu} [k]\end{aligned}}}$$

which is an exact solution to the discretization problem.

When A is singular, the latter expression can still be used by replacing ${\displaystyle e^{\mathbf {A} T}}$ by its Taylor expansion, $${\displaystyle e^{\mathbf {A} T}=\sum _{k=0}^{\infty }{\frac {1}{k!}}(\mathbf {A} T)^{k}.}$$ This yields $${\displaystyle {\begin{aligned}\mathbf {x} [k+1]&=e^{\mathbf {A} T}\mathbf {x} [k]+\left(\int _{0}^{T}e^{\mathbf {A} v}dv\right)\mathbf {Bu} [k]\\[2pt]&=\left(\sum _{k=0}^{\infty }{\frac {1}{k!}}(\mathbf {A} T)^{k}\right)\mathbf {x} [k]+\left(\sum _{k=1}^{\infty }{\frac {1}{k!}}\mathbf {A} ^{k-1}T^{k}\right)\mathbf {Bu} [k],\end{aligned}}}$$ which is the form used in practice.

Exact discretization may sometimes be intractable due to the heavy matrix exponential and integral operations involved. It is much easier to calculate an approximate discrete model, based on that for small timesteps ${\displaystyle e^{\mathbf {A} T}\approx \mathbf {I} +\mathbf {A} T}$. The approximate solution then becomes: $${\displaystyle \mathbf {x} [k+1]\approx (\mathbf {I} +\mathbf {A} T)\mathbf {x} [k]+T\mathbf {Bu} [k]}$$

This is also known as the Euler method, which is also known as the forward Euler method. Other possible approximations are ${\displaystyle e^{\mathbf {A} T}\approx (\mathbf {I} -\mathbf {A} T)^{-1}}$, otherwise known as the backward Euler method and ${\displaystyle e^{\mathbf {A} T}\approx (\mathbf {I} +{\tfrac {1}{2}}\mathbf {A} T)(\mathbf {I} -{\tfrac {1}{2}}\mathbf {A} T)^{-1}}$, which is known as the bilinear transform, or Tustin transform. Each of these approximations has different stability properties. The bilinear transform preserves the instability of the continuous-time system.

In statistics and machine learning, discretization refers to the process of converting continuous features or variables to discretized or nominal features. This can be useful when creating probability mass functions.

In generalized functions theory, discretization arises as a particular case of the Convolution Theorem on tempered distributions

${\displaystyle {\mathcal {F}}\{f*\operatorname {III} \}={\mathcal {F}}\{f\}\cdot \operatorname {III} }$

${\displaystyle {\mathcal {F}}\{\alpha \cdot \operatorname {III} \}={\mathcal {F}}\{\alpha \}*\operatorname {III} }$

where ${\displaystyle \operatorname {III} }$ is the Dirac comb, ${\displaystyle \cdot \operatorname {III} }$ is discretization, ${\displaystyle *\operatorname {III} }$ is periodization, ${\displaystyle f}$ is a rapidly decreasing tempered distribution (e.g. a Dirac delta function ${\displaystyle \delta }$ or any other compactly supported function), ${\displaystyle \alpha }$ is a smooth, slowly growing ordinary function (e.g. the function that is constantly ${\displaystyle 1}$ or any other band-limited function) and ${\displaystyle {\mathcal {F}}}$ is the (unitary, ordinary frequency) Fourier transform. Functions ${\displaystyle \alpha }$ which are not smooth can be made smooth using a mollifier prior to discretization.

As an example, discretization of the function that is constantly ${\displaystyle 1}$ yields the sequence ${\displaystyle [..,1,1,1,..]}$ which, interpreted as the coefficients of a linear combination of Dirac delta functions, forms a Dirac comb. If additionally truncation is applied, one obtains finite sequences, e.g. ${\displaystyle [1,1,1,1]}$. They are discrete in both, time and frequency.

• Discrete event simulation
• Discrete space
• Discrete time and continuous time
• Finite difference method
• Finite volume method for unsteady flow
• Interpolation
• Smoothing
• Stochastic simulation
• Time-scale calculus