Discretization






A solution to a discretized partial differential equation, obtained with the finite element method.


In applied mathematics, discretization is the process of transferring continuous functions, models, variables, and equations into discrete counterparts. This process is usually carried out as a first step toward making them suitable for numerical evaluation and implementation on digital computers. Dichotomization is the special case of discretization in which the number of discrete classes is 2, which can approximate a continuous variable as a binary variable (creating a dichotomy for modeling purposes, as in binary classification).


Discretization is also related to discrete mathematics, and is an important component of granular computing. In this context, discretization may also refer to modification of variable or category granularity, as when multiple discrete variables are aggregated or multiple discrete categories fused.


Whenever continuous data is discretized, there is always some amount of discretization error. The goal is to reduce the amount to a level considered negligible for the modeling purposes at hand.


The terms discretization and quantization often have the same denotation but not always identical connotations. (Specifically, the two terms share a semantic field.) The same is true of discretization error and quantization error.


Mathematical methods relating to discretization include the Euler–Maruyama method and the zero-order hold.




Contents






  • 1 Discretization of linear state space models


    • 1.1 Discretization of process noise


    • 1.2 Derivation


    • 1.3 Approximations




  • 2 Discretization of continuous features


  • 3 See also


  • 4 References


  • 5 Further reading


  • 6 External links





Discretization of linear state space models


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



x˙(t)=Ax(t)+Bu(t)+w(t){displaystyle {dot {mathbf {x} }}(t)=mathbf {A} mathbf {x} (t)+mathbf {B} mathbf {u} (t)+mathbf {w} (t)}{dot {{mathbf {x}}}}(t)={mathbf A}{mathbf {x}}(t)+{mathbf B}{mathbf {u}}(t)+{mathbf {w}}(t)

y(t)=Cx(t)+Du(t)+v(t){displaystyle mathbf {y} (t)=mathbf {C} mathbf {x} (t)+mathbf {D} mathbf {u} (t)+mathbf {v} (t)}{mathbf {y}}(t)={mathbf C}{mathbf {x}}(t)+{mathbf D}{mathbf {u}}(t)+{mathbf {v}}(t)


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



w(t)∼N(0,Q){displaystyle mathbf {w} (t)sim N(0,mathbf {Q} )}{mathbf {w}}(t)sim N(0,{mathbf Q})

v(t)∼N(0,R){displaystyle mathbf {v} (t)sim N(0,mathbf {R} )}{mathbf {v}}(t)sim N(0,{mathbf R})


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



x[k+1]=Adx[k]+Bdu[k]+w[k]{displaystyle mathbf {x} [k+1]=mathbf {A} _{d}mathbf {x} [k]+mathbf {B} _{d}mathbf {u} [k]+mathbf {w} [k]}{mathbf {x}}[k+1]={mathbf A}_{d}{mathbf {x}}[k]+{mathbf B}_{d}{mathbf {u}}[k]+{mathbf {w}}[k]

y[k]=Cdx[k]+Ddu[k]+v[k]{displaystyle mathbf {y} [k]=mathbf {C} _{d}mathbf {x} [k]+mathbf {D} _{d}mathbf {u} [k]+mathbf {v} [k]}{mathbf {y}}[k]={mathbf C}_{d}{mathbf {x}}[k]+{mathbf D}_{d}{mathbf {u}}[k]+{mathbf {v}}[k]


with covariances



w[k]∼N(0,Qd){displaystyle mathbf {w} [k]sim N(0,mathbf {Q} _{d})}{mathbf {w}}[k]sim N(0,{mathbf Q}_{d})

v[k]∼N(0,Rd){displaystyle mathbf {v} [k]sim N(0,mathbf {R} _{d})}{mathbf {v}}[k]sim N(0,{mathbf R}_{d})


where



Ad=eAT=L−1{(sI−A)−1}t=T{displaystyle mathbf {A} _{d}=e^{mathbf {A} T}={mathcal {L}}^{-1}{(smathbf {I} -mathbf {A} )^{-1}}_{t=T}}{mathbf A}_{d}=e^{{{mathbf A}T}}={mathcal {L}}^{{-1}}{(s{mathbf I}-{mathbf A})^{{-1}}}_{{t=T}}


Bd=(∫τ=0TeAτdτ)B=A−1(Ad−I)B{displaystyle mathbf {B} _{d}=left(int _{tau =0}^{T}e^{mathbf {A} tau }dtau right)mathbf {B} =mathbf {A} ^{-1}(mathbf {A} _{d}-I)mathbf {B} }{mathbf B}_{d}=left(int _{{tau =0}}^{{T}}e^{{{mathbf A}tau }}dtau right){mathbf B}={mathbf A}^{{-1}}({mathbf A}_{d}-I){mathbf B}, if A{displaystyle mathbf {A} }mathbf {A} is nonsingular

Cd=C{displaystyle mathbf {C} _{d}=mathbf {C} }{mathbf C}_{d}={mathbf C}

Dd=D{displaystyle mathbf {D} _{d}=mathbf {D} }{mathbf D}_{d}={mathbf D}

Qd=∫τ=0TeAτQeA⊤τdτ{displaystyle mathbf {Q} _{d}=int _{tau =0}^{T}e^{mathbf {A} tau }mathbf {Q} e^{mathbf {A} ^{top }tau }dtau }{displaystyle mathbf {Q} _{d}=int _{tau =0}^{T}e^{mathbf {A} tau }mathbf {Q} e^{mathbf {A} ^{top }tau }dtau }

Rd=RT{displaystyle mathbf {R} _{d}=mathbf {R} T}{displaystyle mathbf {R} _{d}=mathbf {R} T}


and T{displaystyle T}T is the sample time, although A⊤{displaystyle mathbf {A} ^{top }}{displaystyle mathbf {A} ^{top }} is the transposed matrix of A{displaystyle mathbf {A} }mathbf {A} .


A clever trick to compute Ad and Bd in one step is by utilizing the following property:[1]:p. 215


e[AB00]T=[M11M120I]{displaystyle e^{{begin{bmatrix}mathbf {A} &mathbf {B} \mathbf {0} &mathbf {0} end{bmatrix}}T}={begin{bmatrix}mathbf {M_{11}} &mathbf {M_{12}} \mathbf {0} &mathbf {I} end{bmatrix}}}e^{{{begin{bmatrix}{mathbf {A}}&{mathbf {B}}\{mathbf {0}}&{mathbf {0}}end{bmatrix}}T}}={begin{bmatrix}{mathbf {M_{{11}}}}&{mathbf {M_{{12}}}}\{mathbf {0}}&{mathbf {I}}end{bmatrix}}

and then having



Ad=M11{displaystyle mathbf {A} _{d}=mathbf {M} _{11}}{mathbf A}_{d}={mathbf M}_{{11}}

Bd=M12{displaystyle mathbf {B} _{d}=mathbf {M} _{12}}{mathbf B}_{d}={mathbf M}_{{12}}



Discretization of process noise


Numerical evaluation of Qd{displaystyle mathbf {Q} _{d}}{mathbf {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 [2]



F=[−AQ0A⊤]T{displaystyle mathbf {F} ={begin{bmatrix}-mathbf {A} &mathbf {Q} \mathbf {0} &mathbf {A} ^{top }end{bmatrix}}T}{displaystyle mathbf {F} ={begin{bmatrix}-mathbf {A} &mathbf {Q} \mathbf {0} &mathbf {A} ^{top }end{bmatrix}}T}

G=eF=[…Ad−1Qd0Ad⊤].{displaystyle mathbf {G} =e^{mathbf {F} }={begin{bmatrix}dots &mathbf {A} _{d}^{-1}mathbf {Q} _{d}\mathbf {0} &mathbf {A} _{d}^{top }end{bmatrix}}.}{displaystyle mathbf {G} =e^{mathbf {F} }={begin{bmatrix}dots &mathbf {A} _{d}^{-1}mathbf {Q} _{d}\mathbf {0} &mathbf {A} _{d}^{top }end{bmatrix}}.}


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:


Qd=(Ad⊤)⊤(Ad−1Qd)=Ad(Ad−1Qd).{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}).}{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}).}


Derivation


Starting with the continuous model


x˙(t)=Ax(t)+Bu(t){displaystyle mathbf {dot {x}} (t)=mathbf {A} mathbf {x} (t)+mathbf {B} mathbf {u} (t)}{mathbf {{dot {x}}}}(t)={mathbf A}{mathbf x}(t)+{mathbf B}{mathbf u}(t)

we know that the matrix exponential is


ddteAt=AeAt=eAtA{displaystyle {frac {d}{dt}}e^{mathbf {A} t}=mathbf {A} e^{mathbf {A} t}=e^{mathbf {A} t}mathbf {A} }{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


e−Atx˙(t)=e−AtAx(t)+e−AtBu(t){displaystyle e^{-mathbf {A} t}mathbf {dot {x}} (t)=e^{-mathbf {A} t}mathbf {A} mathbf {x} (t)+e^{-mathbf {A} t}mathbf {B} mathbf {u} (t)}e^{{-{mathbf A}t}}{mathbf {{dot {x}}}}(t)=e^{{-{mathbf A}t}}{mathbf A}{mathbf x}(t)+e^{{-{mathbf A}t}}{mathbf B}{mathbf u}(t)

which we recognize as


ddt(e−Atx(t))=e−AtBu(t){displaystyle {frac {d}{dt}}(e^{-mathbf {A} t}mathbf {x} (t))=e^{-mathbf {A} t}mathbf {B} mathbf {u} (t)}{frac {d}{dt}}(e^{{-{mathbf A}t}}{mathbf x}(t))=e^{{-{mathbf A}t}}{mathbf B}{mathbf u}(t)

and by integrating..



e−Atx(t)−e0x(0)=∫0te−AτBu(τ)dτ{displaystyle e^{-mathbf {A} t}mathbf {x} (t)-e^{0}mathbf {x} (0)=int _{0}^{t}e^{-mathbf {A} tau }mathbf {B} mathbf {u} (tau )dtau }e^{{-{mathbf A}t}}{mathbf x}(t)-e^{0}{mathbf x}(0)=int _{0}^{t}e^{{-{mathbf A}tau }}{mathbf B}{mathbf u}(tau )dtau

x(t)=eAtx(0)+∫0teA(t−τ)Bu(τ)dτ{displaystyle mathbf {x} (t)=e^{mathbf {A} t}mathbf {x} (0)+int _{0}^{t}e^{mathbf {A} (t-tau )}mathbf {B} mathbf {u} (tau )dtau }{mathbf x}(t)=e^{{{mathbf A}t}}{mathbf x}(0)+int _{0}^{t}e^{{{mathbf A}(t-tau )}}{mathbf B}{mathbf u}(tau )dtau


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.



x[k] =def x(kT){displaystyle mathbf {x} [k] {stackrel {mathrm {def} }{=}} mathbf {x} (kT)}{mathbf x}[k] {stackrel {{mathrm {def}}}{=}} {mathbf x}(kT)

x[k]=eAkTx(0)+∫0kTeA(kT−τ)Bu(τ)dτ{displaystyle mathbf {x} [k]=e^{mathbf {A} kT}mathbf {x} (0)+int _{0}^{kT}e^{mathbf {A} (kT-tau )}mathbf {B} mathbf {u} (tau )dtau }{mathbf x}[k]=e^{{{mathbf A}kT}}{mathbf x}(0)+int _{0}^{{kT}}e^{{{mathbf A}(kT-tau )}}{mathbf B}{mathbf u}(tau )dtau

x[k+1]=eA(k+1)Tx(0)+∫0(k+1)TeA((k+1)T−τ)Bu(τ)dτ{displaystyle 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 {B} mathbf {u} (tau )dtau }{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 B}{mathbf u}(tau )dtau

x[k+1]=eAT[eAkTx(0)+∫0kTeA(kT−τ)Bu(τ)dτ]+∫kT(k+1)TeA(kT+T−τ)Bu(τ)dτ{displaystyle 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 {B} mathbf {u} (tau )dtau right]+int _{kT}^{(k+1)T}e^{mathbf {A} (kT+T-tau )}mathbf {B} mathbf {u} (tau )dtau }{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 B}{mathbf u}(tau )dtau right]+int _{{kT}}^{{(k+1)T}}e^{{{mathbf A}(kT+T-tau )}}{mathbf B}{mathbf u}(tau )dtau


We recognize the bracketed expression as x[k]{displaystyle mathbf {x} [k]}{mathbf x}[k], and the second term can be simplified by substituting with the function v(τ)=kT+T−τ{displaystyle v(tau )=kT+T-tau }{displaystyle v(tau )=kT+T-tau }. Note that dτ=−dv{displaystyle dtau =-dv}{displaystyle dtau =-dv}. We also assume that u{displaystyle mathbf {u} }mathbf u is constant during the integral, which in turn yields


x[k+1]=eATx[k]−(∫v(kT)v((k+1)T)eAvdv)Bu[k]=eATx[k]−(∫T0eAvdv)Bu[k]=eATx[k]+(∫0TeAvdv)Bu[k]=eATx[k]+A−1(eAT−I)Bu[k]{displaystyle {begin{matrix}mathbf {x} [k+1]&=&e^{mathbf {A} T}mathbf {x} [k]-left(int _{v(kT)}^{v((k+1)T)}e^{mathbf {A} v}dvright)mathbf {B} mathbf {u} [k]\&=&e^{mathbf {A} T}mathbf {x} [k]-left(int _{T}^{0}e^{mathbf {A} v}dvright)mathbf {B} mathbf {u} [k]\&=&e^{mathbf {A} T}mathbf {x} [k]+left(int _{0}^{T}e^{mathbf {A} v}dvright)mathbf {B} mathbf {u} [k]\&=&e^{mathbf {A} T}mathbf {x} [k]+mathbf {A} ^{-1}left(e^{mathbf {A} T}-mathbf {I} right)mathbf {B} mathbf {u} [k]end{matrix}}}{displaystyle {begin{matrix}mathbf {x} [k+1]&=&e^{mathbf {A} T}mathbf {x} [k]-left(int _{v(kT)}^{v((k+1)T)}e^{mathbf {A} v}dvright)mathbf {B} mathbf {u} [k]\&=&e^{mathbf {A} T}mathbf {x} [k]-left(int _{T}^{0}e^{mathbf {A} v}dvright)mathbf {B} mathbf {u} [k]\&=&e^{mathbf {A} T}mathbf {x} [k]+left(int _{0}^{T}e^{mathbf {A} v}dvright)mathbf {B} mathbf {u} [k]\&=&e^{mathbf {A} T}mathbf {x} [k]+mathbf {A} ^{-1}left(e^{mathbf {A} T}-mathbf {I} right)mathbf {B} mathbf {u} [k]end{matrix}}}

which is an exact solution to the discretization problem.



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 eAT≈I+AT{displaystyle e^{mathbf {A} T}approx mathbf {I} +mathbf {A} T}e^{{{mathbf A}T}}approx {mathbf I}+{mathbf A}T. The approximate solution then becomes:


x[k+1]≈(I+AT)x[k]+TBu[k]{displaystyle mathbf {x} [k+1]approx (mathbf {I} +mathbf {A} T)mathbf {x} [k]+Tmathbf {B} mathbf {u} [k]}{mathbf x}[k+1]approx ({mathbf I}+{mathbf A}T){mathbf x}[k]+T{mathbf B}{mathbf u}[k]

This is also known as Euler's method. Other possible approximations are eAT≈(I−AT)−1{displaystyle e^{mathbf {A} T}approx left(mathbf {I} -mathbf {A} Tright)^{-1}}e^{{{mathbf A}T}}approx left({mathbf I}-{mathbf A}Tright)^{{-1}} and eAT≈(I+12AT)(I−12AT)−1{displaystyle e^{mathbf {A} T}approx left(mathbf {I} +{frac {1}{2}}mathbf {A} Tright)left(mathbf {I} -{frac {1}{2}}mathbf {A} Tright)^{-1}}e^{{{mathbf A}T}}approx left({mathbf I}+{frac {1}{2}}{mathbf A}Tright)left({mathbf I}-{frac {1}{2}}{mathbf A}Tright)^{{-1}}. Each of them have different stability properties. The last one is known as the bilinear transform, or Tustin transform, and preserves the (in)stability of the continuous-time system.



Discretization of continuous features


Main article: Discretization of continuous features

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.



See also



  • Discrete event simulation

  • Discrete space

  • Discrete time and continuous time

  • Finite difference method

  • Finite volume method for unsteady flow

  • Stochastic simulation

  • Time-scale calculus



References




  1. ^ Raymond DeCarlo: Linear Systems: A State Variable Approach with Numerical Implementation, Prentice Hall, NJ, 1989


  2. ^ Charles Van Loan: Computing integrals involving the matrix exponential, IEEE Transactions on Automatic Control. 23 (3): 395–404, 1978



Further reading




  • Robert Grover Brown & Patrick Y. C. Hwang. Introduction to random signals and applied Kalman filtering (3rd ed.). ISBN 978-0471128397..mw-parser-output cite.citation{font-style:inherit}.mw-parser-output .citation q{quotes:"""""""'""'"}.mw-parser-output .citation .cs1-lock-free a{background:url("//upload.wikimedia.org/wikipedia/commons/thumb/6/65/Lock-green.svg/9px-Lock-green.svg.png")no-repeat;background-position:right .1em center}.mw-parser-output .citation .cs1-lock-limited a,.mw-parser-output .citation .cs1-lock-registration a{background:url("//upload.wikimedia.org/wikipedia/commons/thumb/d/d6/Lock-gray-alt-2.svg/9px-Lock-gray-alt-2.svg.png")no-repeat;background-position:right .1em center}.mw-parser-output .citation .cs1-lock-subscription a{background:url("//upload.wikimedia.org/wikipedia/commons/thumb/a/aa/Lock-red-alt-2.svg/9px-Lock-red-alt-2.svg.png")no-repeat;background-position:right .1em center}.mw-parser-output .cs1-subscription,.mw-parser-output .cs1-registration{color:#555}.mw-parser-output .cs1-subscription span,.mw-parser-output .cs1-registration span{border-bottom:1px dotted;cursor:help}.mw-parser-output .cs1-ws-icon a{background:url("//upload.wikimedia.org/wikipedia/commons/thumb/4/4c/Wikisource-logo.svg/12px-Wikisource-logo.svg.png")no-repeat;background-position:right .1em center}.mw-parser-output code.cs1-code{color:inherit;background:inherit;border:inherit;padding:inherit}.mw-parser-output .cs1-hidden-error{display:none;font-size:100%}.mw-parser-output .cs1-visible-error{font-size:100%}.mw-parser-output .cs1-maint{display:none;color:#33aa33;margin-left:0.3em}.mw-parser-output .cs1-subscription,.mw-parser-output .cs1-registration,.mw-parser-output .cs1-format{font-size:95%}.mw-parser-output .cs1-kern-left,.mw-parser-output .cs1-kern-wl-left{padding-left:0.2em}.mw-parser-output .cs1-kern-right,.mw-parser-output .cs1-kern-wl-right{padding-right:0.2em}


  • Chi-Tsong Chen (1984). Linear System Theory and Design. Philadelphia, PA, USA: Saunders College Publishing. ISBN 978-0030716911.


  • C. Van Loan (Jun 1978). "Computing integrals involving the matrix exponential". IEEE Transactions on Automatic Control. 23 (3): 395–404. doi:10.1109/TAC.1978.1101743.


  • R.H. Middleton & G.C. Goodwin (1990). Digital control and estimation: a unified approach. p. 33f. ISBN 978-0132116657.



External links









-