# How to find appropriate discretization method for a continuous time domain state space model?

I have a dsp algorithm which is based on the below given state space model in the continuous-time domain

$$\begin{bmatrix} \frac{\mathrm{d}\hat{\psi}_{r_{\alpha}}}{\mathrm{d}t} \\ \frac{\mathrm{d}\hat{\psi}_{r_{\beta}}}{\mathrm{d}t} \end{bmatrix} = \begin{bmatrix} -\frac{R_R}{L_L + L_M} & -p_p\cdot\omega_m \\ p_p\cdot\omega_m & -\frac{R_R}{L_L + L_M} \end{bmatrix} \cdot \begin{bmatrix} \hat{\psi}_{r_{\alpha}} \\ \hat{\psi}_{r_{\beta}} \end{bmatrix} + \begin{bmatrix} \frac{L_M\cdot R_R}{L_L + L_M} & 0 \\ 0 & \frac{L_M\cdot R_R}{L_L + L_M} \end{bmatrix} \cdot \begin{bmatrix} i_{s_{\alpha}} \\ i_{s_{\beta}} \end{bmatrix}$$

It is basically a model of a dynamic system, where the variables $$i_{s_{\alpha}}$$, $$i_{s_{\beta}}$$ and $$\omega_m$$ are the inputs of the dynamic system and the $$\hat{\psi}_{r_{\alpha}}$$, $$\hat{\psi}_{r_{\beta}}$$ are its outputs (the unmeasurable state variables of the system). This model is intended to be used as a state observer. I know that there are much more robust approaches for estimation of the unmeasurable state variables of a dynamic system but I would like to do a comparison between several methods. As far as the parameters of the state space model: $$R_S = 7.400\cdot 10^{-3}\,\Omega$$, $$R_R = 7.548\cdot 10^{-3}\,\Omega$$, $$L_M = 4.265\cdot 10^{-3}\,\mathrm{H}$$, $$L_L = 0.231\cdot 10^{-3}\,\mathrm{H}$$, $$p_p = 3.0$$.

For the software implementation purposes it is necessary to discretize the continuous-time domain model (sampling period is $$T = 100\cdot 10^{-6}\,\mathrm{s}$$). I have been thinking about the simple Euler method (i.e. the state transition matrix $$\Psi = \mathbf{I} + \frac{(\mathbf{A}\cdot T)}{2!} + \frac{(\mathbf{A}\cdot T)^2}{3!} + \ldots$$ approximated only via the first two terms of the series expansion) usage. But I have noticed that the above given system is probably much more complex for computation than it seems at first glance. I have calculated the eigenvalues of the system matrix (i.e. poles of the system) in respect to the input $$\omega_m$$ and I have found that the poles are pretty near to the stability boundary in the s-plane.

If I transformed the system into the discrete form via the ZOH method I found that the system poles (in respect to the input $$\omega_m$$) are even nearly on the stability boundary.

My question is how I can implement this system in a software (in the single precision floating point) to be sure that the solution (the system state variables) will be stable?

The idea behind my question is that I have noticed that the poles of the discretized system (zero-order hold discretization) are very close to the unit circle boundary (the "reserve" is only a couple of ten thousandths). Source of my doubt is that in case I implement the discrete system in a control software it is from my point of view very good chance that due to the limited precision the poles can move outside the unit circle and the system can become very easily unstable.

• Why not diagonalize the state transition matrix? That should make the pole positions and eigenvectors obvious. Commented Jan 21, 2023 at 21:48
• Why do you say that you would "implement the discrete system in a control software"? I am not sure to understand that. What you implement is a controller. The model is just there for design purposes.
– KBS
Commented Jan 21, 2023 at 22:24
• @KBS thank you for your reaction. The above mentioned model is intended to be used in the state observer role. The observer is necessary because it provides the information about the unmeasurable state variables of the controlled system. Commented Jan 22, 2023 at 9:49
• What is $\omega_m$?
– KBS
Commented Jan 22, 2023 at 12:13
• @KBS the $\omega_m$ is a mechanical speed in $(\mathrm{rad}\cdot \mathrm{s}^{-1})$. Commented Jan 22, 2023 at 14:39

I see two possible ways to address this problem.

The first one consists of considering the discretized system but implementing each coefficient as a triplet of integers $$(a,b,c)$$ which corresponds to the number $$\dfrac{a}{b}\cdot 10^{c}$$ and you do the calculations in this form. Check out fixed-point arithmetic.

The second consists of using the so-called delta transform. If we consider a continuous-time system

$$\dot{x}(t)=Ax(t)+Bu(t)$$

and we assume that the system is controlled with a zero order hold with sampling time $$T_s$$, we get the following discrete-time representation

$$x_{k+1}=e^{AT_s}x(kT)+\left(\int_0^{T_s}e^{As}Bds\right)u(kT_s),$$

where $$x_k:=x(kT_s)$$ and $$u_k:=u(kT_s)$$.

The issue with that is when $$T_s$$ is small, $$e^{AT_s}\approx I$$ and $$\int_0^{T_s}e^{As}Bds\approx 0$$.

One way to circument this, is to consider the operator $$\delta$$ defined as $$\delta x_k=(x_{k+1}-x_k)/T_s$$, which yields the representation

$$\delta x_{k}=\dfrac{e^{AT_s}-I}{T_s}x_k+\dfrac{1}{T_s}\left(\int_0^{T_s}e^{As}Bds\right)u_k.$$

The benefit is that when $$T_s$$ is small, one has $$\dfrac{e^{AT_s}-I}{T_s}\approx A$$ and $$\dfrac{1}{T_s}\int_0^{T_s}e^{As}Bds\approx B$$. In this respect it is much easier to consider this representation when the sampling period is small.

I have tried with your numerical values, and the matrix we get is very close to A, which contains values of reasonable magnitude.

You can check the book by Middleton and Goodwin, "Digital control and estimation: a unified approach" on the topic as well as the many published papers.

In any way, let us consider the following system in $$\delta$$-form

$$\begin{array}{rcl} \delta x_{k}&=&A_\delta x_k+A_\delta u_k\\ y_k&=&Cx_k = \displaystyle C\left(T_s\sum_{i=0}^{k-1}\delta x_i+x_0\right). \end{array}$$

One can define the following observer

$$\begin{array}{rcl} \delta \hat x_{k}&=&A_\delta \hat x_k+A_\delta u_k+L(y_k-C\hat x_k) \\ \hat x_k &=& \displaystyle T_s\sum_{i=0}^{k-1}\delta \hat x_i+\hat x_0. \end{array}$$ An important point here is that $$\hat x_k$$ can be computed recursively using fixed point arithmetic (by avoiding to multiply all the time by $$T_s$$).

The dynamics of the error $$e:=x-\hat x$$ is given by

$$\delta e_{k}=(A_\delta-LC)e_k$$

and is stable if and only if the spectrum of $$A_\delta-LC$$ lies in the disc centered about $$-1$$ and with radius $$1/T_s$$. This can be done with LMI methods.

Now, if we consider a state-feedback controller of the form

$$u_k=K\hat x_k$$

this can be expressed as

$$u_k=K\left(T_s\sum_{i=0}^{k-1}\delta\hat x_i+\hat x_0\right).$$

• Thank you very much for your answer. As far as your second approach. My understanding is that you suggest to calculate: $x_{k+1} = x_k + \delta x_k$, where the $\delta x_k$ calculation is a "trick" how to avoid the original ill-conditioned system of ode's (other interpretion of the system poles very close to the unit circle in the z-plane). Is my understanding correct? Commented Jan 23, 2023 at 21:30
• Yes, this can be roughly summarized this way (the sampling period factor is missing in your expression). You can also work directly in this formalism.
– KBS
Commented Jan 23, 2023 at 21:53
• @Steve I have updated the answer with more technical details.
– KBS
Commented Jan 24, 2023 at 14:56
• Thank you for your updated answer. I am still not able to justify for myself whether usage of more terms in the state transition matrix $\mathbf{\Psi}$ calculation can prevent the discretized set of equations against divergency. Please can you tell me your opinion? Commented Jan 26, 2023 at 16:06
• If you are talking of using a Taylor expansion for computing the exponential this a bad idea because it may converge poorly. Also, powers of $T_s$ will go to small values quickly and reach machine precision rather fast, yielding rounding errors. Your problem is not novel, people have had it in the past and has been solved.
– KBS
Commented Jan 26, 2023 at 16:36