Editing Discretization (section)

== Discretization of linear state space models {{anchor|discrete function}} ==
Discretization is also concerned with the transformation of continuous [[differential equation]]s into discrete [[difference equations]], suitable for [[Numerical analysis|numerical computing]].

The following continuous-time [[State space (controls)|state space model]]

<math display=block>\begin{align}
  \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{align}</math>

where {{mvar|v}} and {{mvar|w}} are continuous zero-mean [[white noise]] sources with [[power spectral density|power spectral densities]]

<math display=block>\begin{align}
  \mathbf{w}(t) &\sim N(0,\mathbf Q) \\[2pt]
  \mathbf{v}(t) &\sim N(0,\mathbf R)
\end{align}</math>

can be discretized, assuming [[zero-order hold]] for the input {{mvar|u}} and continuous integration for the noise {{mvar|v}}, to

<math display=block>\begin{align}
  \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{align}</math>

with covariances

<math display=block>\begin{align}
  \mathbf{w}[k] &\sim N(0,\mathbf{Q_d}) \\[2pt]
  \mathbf{v}[k] &\sim N(0,\mathbf{R_d})
\end{align}</math>

where

<math display=block>\begin{align}
\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{align}</math>

and {{mvar|T}} is the [[sample time]]. If {{math|'''A'''}} is [[Invertible matrix|nonsingular]], <math>\mathbf{B_d} = \mathbf A^{-1}(\mathbf{A_d} - \mathbf{I})\mathbf B.</math>

The equation for the discretized measurement noise is a consequence of the continuous measurement noise being defined with a power spectral density.<ref>{{Cite book|last=Analytic Sciences Corporation. Technical Staff.|title=Applied optimal estimation.|url=https://archive.org/details/appliedoptimales00agel|url-access=limited|date=1974|publisher=M.I.T. Press|others=Gelb, Arthur, 1937-|isbn=0-262-20027-9|location=Cambridge, Mass.|pages=[https://archive.org/details/appliedoptimales00agel/page/n128 121]|oclc=960061}}</ref>

A clever trick to compute {{math|'''A<sub>d</sub>'''}} and {{math|'''B<sub>d</sub>'''}} in one step is by utilizing the following property:<ref>Raymond DeCarlo: ''Linear Systems: A State Variable Approach with Numerical Implementation'', Prentice Hall, NJ, 1989</ref>{{rp|p. 215}}

<math display=block>
  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}
</math>

Where {{math|'''A<sub>d</sub>'''}} and {{math|'''B<sub>d</sub>'''}} are the discretized state-space matrices.

=== Discretization of process noise ===
Numerical evaluation of {{math|'''Q{{sub|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<ref>Charles Van Loan: ''Computing integrals involving the matrix exponential'', IEEE Transactions on Automatic Control. 23 (3): 395–404, 1978</ref>
<math display=block>\begin{align}
  \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{align}</math>
The discretized process noise is then evaluated by multiplying the transpose of the lower-right partition of {{math|'''G'''}} with the upper-right partition of {{math|'''G'''}}:
<math display=block>\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}). </math>

=== Derivation ===
Starting with the continuous model
<math display=block>\mathbf{\dot{x}}(t) = \mathbf{Ax}(t) + \mathbf{Bu}(t)</math>
we know that the [[matrix exponential]] is
<math display=block>\frac{d}{dt}e^{\mathbf{A}t} = \mathbf{A}e^{\mathbf{A}t} = e^{\mathbf{A}t} \mathbf A</math>
and by premultiplying the model we get
<math display=block>e^{-\mathbf{A}t} \mathbf{\dot{x}}(t) = e^{-\mathbf{A}t} \mathbf{Ax}(t) + e^{-\mathbf{A}t} \mathbf{Bu}(t)</math>
which we recognize as
<math display=block>\frac{d}{dt}\Bigl[e^{-\mathbf{A}t}\mathbf x(t) \Bigr] = e^{-\mathbf{A}t} \mathbf{Bu}(t)</math>
and by integrating,
<math display=block>\begin{align}
  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{align}</math>
which is an analytical solution to the continuous model.

Now we want to discretise the above expression. We assume that {{mvar|u}} is [[mathematical constant|constant]] during each timestep.
<math display=block>\begin{align}
  \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{align}</math>
We recognize the bracketed expression as <math>\mathbf x[k]</math>, and the second term can be simplified by substituting with the function <math>v(\tau) = kT + T - \tau</math>. Note that <math>d\tau=-dv</math>. We also assume that {{math|'''u'''}} is constant during the [[integral]], which in turn yields

<math display=block>\begin{align}
  \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{align}</math>

which is an exact solution to the discretization problem.

When {{math|'''A'''}} is singular, the latter expression can still be used by replacing <math> e^{\mathbf{A}T} </math> by its [[Taylor series|Taylor expansion]],
<math display=block> e^{\mathbf{A}T} = \sum_{k=0}^{\infty} \frac{1}{k!} (\mathbf{A}T)^k  .</math>
This yields
<math display=block>\begin{align} 
  \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{align}</math>
which is the form used in practice.

=== Approximations ===
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 <math>e^{\mathbf{A}T} \approx \mathbf I + \mathbf A T</math>. The approximate solution then becomes:
<math display=block>\mathbf x[k+1] \approx (\mathbf I + \mathbf{A}T) \mathbf x[k] + T \mathbf{Bu}[k] </math>

This is also known as the [[Euler method]], which is also known as the forward Euler method. Other possible approximations are <math>e^{\mathbf{A}T} \approx (\mathbf I - \mathbf{A}T)^{-1}</math>, otherwise known as the backward Euler method and <math>e^{\mathbf{A}T} \approx (\mathbf I +\tfrac{1}{2} \mathbf{A}T) (\mathbf I - \tfrac{1}{2} \mathbf{A}T)^{-1}</math>, 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.