Performance Analysis of Speed-Sensorless Induction Motor Drive Using Discrete Current-Error Based MRAS Estimators

In the literature on sensorless control of induction motors, many algorithms have been presented for rotor flux and speed estimation. However, all these algorithms have been developed in the continuous–time domain. The digital realization of the control systems, requires the implementation of those estimation methods in a discrete–time domain. The main goal of this article is comparison of the impact of different numerical integration methods, used in analogue emulation under the digital implementation of the control systems, to the operation of classical Model Reference Adaptive System; CC-based on two current models (MRASCC) speed estimator and its three modified versions developed for the extension of the estimator stability region. In this paper the generalized mathematical model of MRASCC estimator is proposed, which takes into account all known methods for the extension of the stability region of classical speed estimator of this type. After the short discussion of the discretization methods used for the microprocessor implementation of control algorithms the impact of different numerical integration methods on the stable operation range of the classical and modified MRASCC estimators is analyzed and validated in simulation and experimental tests. It is proved that Modified Euler discretization method is much more accurate than forward and backward Euler methods and gives almost as accurate results as Tustin method, however is much less complicated in practical realization.


Introduction
Nowadays, various emerging applications of electric drives such as: transportation systems, industrial mechatronics systems just like tool spindle drives and robotic systems, pose new challenges and problems for designers. In many scientific centers, research is underway on new topologies and techniques for controlling electric drives, especially induction motor (IM) and permanent-magnet synchronous motor (PMSM) drives. Particular attention is paid to torque and angular velocity or position control using advanced methods of control theory [1,2]. Among these issues are sensorless drives, where attention is paid to various methods of inaccessible state variables estimation [3]. It should also be noted that when designing fault-tolerant systems [4], in the event of a failure of the AC current sensors, this signal must be reproduced also, because no control structure can work without the stator current information [5].
In the case of control of IMs, which due to low cost, maintenance-free, ruggedness and reliability are willingly used in various applications, the estimation of the rotor flux and speed based on easily available electrical signals such as voltage and phase current is a significant problem. Many algorithms based on the IM mathematical model for state variables estimation have been proposed in the literature in the last decades, which are briefly discussed and evaluated among others in [6][7][8]. Among different rotor speed and flux estimation methods, such as: open and closed-loop observers, Kalman filters, Model Reference Adaptive Systems (MRAS)-based speed estimators are most popular in practice [1,3,6,7]. A recent review on the most important features of these techniques is given in [8], taking into account among others stability problems in different operating ranges of IM drive.
Currently, the vast majority of electric drive control systems are implemented in microprocessor technology. For this reason, the issue of discretization is increasingly discussed in the literature and two methodologies have been proposed: discrete design and analogue emulation [9][10][11]. The first method consists in obtaining the discrete-time domain model for the previously developed continuous-time domain model of the controlled object and next designing the control algorithm in the discrete-time domain [12][13][14]. The disadvantage of this method is the difficulty in calculating the sampled model, especially for non-linear systems [11,12]. In special cases, such as for IM, this procedure can be simplified by assuming constant angular velocity during the sampling period, and then a discrete model can be formulated in a relatively simple way, which can next be used in the design of the control or estimation algorithm [13,14]. Analogue emulation is based on the discretization of the control law (including the observer of state variables, if necessary), previously designed on the basis of a continuous-time model. The discrete control algorithm thus obtained may have worse properties, especially because of the restrictions on the lowest value of the sampling step that can be used. However, it is much simpler than the first method and is therefore very popular in practical applications. Various methods of numerical integration are used in an analogue emulation [15][16][17].
The forward Euler (FE) method is the most commonly used numerical method during discretization of the continuous model, but stability problems may occur, especially for high speed drives and relatively low sampling step used for the discretization of different state variables estimators [18]. That is why several articles have recently appeared about discrete state variable estimators for sensorless IM drives: the rotor flux [19][20][21][22] or motor speed [20][21][22][23][24][25][26].
As it results from [8], almost all MRAS-type speed estimators are unstable in some ranges of regenerating operation. Therefore, a number of attempts were made to modify their mathematical models in order to obtain stable operation in the entire range of the speed and load torque changes. In the case of the AFO estimator, originally proposed in [27], these issues were raised in [28][29][30][31][32][33], while for the MRAS CC estimator in [32][33][34][35][36], but most of these solutions have only been tested in the continuous-time domain in the regenerating mode. Only in the case of the MRAS CC estimator, originally developed in [37], the authors of the conference paper [38] analyzed the impact of two discretization methods (forward and modified Euler) on the stability of one of its modified versions, namely MRAS CC+ϕ .
Therefore, this article summarizes the influence of selected numerical integration methods used in analogue emulation on the stability issues of the classical MRAS CC speed estimator, based on the previous works of the authors [25,26], however the special attention is paid to all modified stabilized versions of MRAS CC [33][34][35][36][37] and their operation after discretization, especially in the regenerating mode. The main goal of the research is to check whether the three modified versions of the MRAS CC -type speed estimators are stable also in the discrete-time domain and to what extent, after discretization realized with the use of the four most popular numerical integration methods. So the previous research presented in conference papers of the authors [25,26,38] has been extended significantly in this paper. The comparative study of all these discretized MRAS-type speed estimators, namely: classical MRAS CC , Energies 2020, 13, 2595 3 of 23 and three modified versions: MRAS CC+G , MRAS CC+ϕ and MRAS CC+µ (all details on these estimators will be given in the next sections) is presented in motoring and regenerating modes (mechanical characteristics) and transient behaviour of these estimators is demonstrated under the slow reversal and constant speed operation in regenerating mode, as well as in simulation and experimental tests.
The paper is organized as follows: In Section 2 the mathematical model of the induction motor and classical as well as all modified MRAS CC -type speed estimators are presented. It is worth noting, that the generalized mathematical model is proposed, which takes into account all known methods for the extension of the stability region of the classical MRAS CC speed estimator. Next, in Section 3 the short overview of the discretization methods used for the microprocessor implementation of control algorithms is given. The main simulation and experimental results are presented in following sections. In Section 4 the descriptions of the direct rotor field oriented control (DRFOC) structure and experimental set-up are presented, while in the Section 5 the impact of the numerical integration methods on the stable operation range of the classical MRAS CC estimator is analyzed and validated in simulation tests. In the next section the influence of the discrete realization of three modified MRAS CC speed estimators is tested in simulation as well as validated by experimental tests, especially in the regenerating mode. Finally, some concluding remarks close out the paper.

Mathematical Description of the Induction Motor
The IM mathematical model can be written in the form of linear nonstationary state equation, using well-known simplifications (constant and lumped parameters of the windings, symmetrical windings and the air gap, neglected hysteresis and eddy currents), using spatial vectors, in per unit [p.u.] system [2]. In this article this model is presented in a reference frame d-q, rotating with an arbitrary angular velocity ω k and all spatial vectors are marked with a bold font as k = k d + jk q . This mathematical model consists of state equation for electromagnetic state variables: with: and motion equation together with the electromagnetic torque expression: where: x, y, u-state, output and control vectors, respectively, A, B, C-state, control and output matrices, respectively, u s , i s , ψ s , ψ r -stator voltage, stator current, stator and rotor flux spatial vectors, ω m -electrical rotor speed, l s , l r , l m ,-stator, rotor and mutual inductance, r s , r r -stator and rotor winding resistance, k r = l m /l r , l σ = l s σ, r 1 = r s + r r k r 2 , τ r = l r /r r , σ = 1-l m 2 /l s /l r -total leakage factor, and rotor flux-current-based estimator: where: i s ,î s ,ψ r -measured stator current, estimated stator current and estimated rotor flux vectors, respectively.
Both these models are adapted by the estimated rotor speed,ω m , obtained from a special adaptive rule designed using Lyapunov or Popov stability theory [27,37]: where: K pω , K iω -are the parameters of the PI-type speed adaptation algorithm,e i = i s −î s -is stator current estimation error, i.e. the difference between the estimated and measured stator current signals. The schematic diagram of the MRAS CC estimator is shown in Figure 1. stator current estimator: and rotor flux-current-based estimator: where: ˆ, , s s r i i ψ -measured stator current, estimated stator current and estimated rotor flux vectors, respectively. Both these models are adapted by the estimated rotor speed, ˆm ω , obtained from a special adaptive rule designed using Lyapunov or Popov stability theory [27,37]: where: K pω , K iω -are the parameters of the PI-type speed adaptation algorithm,ŝ s = − i e i i -is stator current estimation error, i.e. the difference between the estimated and measured stator current signals. The schematic diagram of the MRAS CC estimator is shown in Figure 1. This estimator is similar to well-known adaptive full-order observer (AFO), described in [27] and widely tested in the literature [28][29][30][31]. In the case of the original AFO, the feedback from the This estimator is similar to well-known adaptive full-order observer (AFO), described in [27] and widely tested in the literature [28][29][30][31]. In the case of the original AFO, the feedback from the stator current estimation error with a suitable gain matrix G has been proposed [37]. However, in most research works till early 2000, G = 0 was assumed. In [28], it was mentioned, that for such assumption the AFO is unstable in some part of regenerating mode. Moreover, new expression for G matrix have been proposed and thus the first stabilization method of this speed estimator has been developed. This issue was extensively addressed in next research works and new stabilization methods for AFO have been proposed in [29][30][31].
Both these speed estimators, AFO and MRAS CC , use the IM as a reference model and their speed adaptation mechanism is the same. However, in the MRAS CC the rotor flux estimator uses the measured stator current as the input, instead of the estimated current as in the case of the full-order observer (AFO) with G = 0. Although the difference is small, their performances are quite different, especially in the case of regenerating mode.
The original scheme of the MRAS CC estimator and its basic properties can be found in [37], while in this article its generalized version, which includes classical and modified MRAS CC -type estimators is presented. This speed estimator in its classic version has many advantages, e.g. design simplicity, very good operation in motoring mode and relatively good robustness to motor parameter changes, however its main drawback is some region of instable operation in the regenerating mode [32][33][34][35]. Detailed analyses on its instability issues can be found in [34,35]. In a few papers [32][33][34]36] some improvement methods of stability range for this speed estimator have been developed, like: proper design of the gain matrix, G [35], -appropriate modification of the speed adaptation algorithm through the shift angle, ϕ between the stator current error and the rotor flux vectors for the MRAS CC [32,33,35], -introduction of the auxiliary variable adapted on-line [36].
Thus, the generalized mathematical model of the MRAS CC estimator, covering all of the above modifications (previously proposed in the cited literature) can be created: where now the generalized state matrix of the estimator is: and: K = 0 r r k r , G = g s g r = g sd + jg sq g rd + jg rq , e i = i s −î s (12) where: G-the gain matrix of the estimator, e i -stator current estimation error,ω m -estimated angular speed of the IM,μ-auxiliary variable. The generalized speed adaptation algorithm can be formulated as: where ϕ is the additional shift angle between the stator current error and the rotor flux vectors, and the auxiliary variable is given by the following algorithm: Energies 2020, 13, 2595 6 of 23 where: K pω , K iω , K pµ , K iµ are the parameters of the adaptation algorithms for speed and auxiliary variable, while ε ω and ε µ are described by error Equations (14) and (16), derived from the Lyapunov theory [32,36]. The general schematic diagram of the mathematical model for the MRAS CC speed estimator with all three stability improvement methods is presented in Figure 2.
{ } * μ ε = ℜ r i ψ e (16) where: K pω , K iω , K pμ , K iμ are the parameters of the adaptation algorithms for speed and auxiliary variable, while ε ω and εμ are described by error Equations (14) and (16), derived from the Lyapunov theory [32,36]. The general schematic diagram of the mathematical model for the MRAS CC speed estimator with all three stability improvement methods is presented in Figure 2.
where: k-parameter to be selected appropriately, -MRAS CC+φ -when G = 0, μ = 0 and the shift angle, φ in the speed adaptation mechanism (13)- (14) is calculated as follows [32,33]: Substituting suitable parameters G, ϕ, µ of the generalized mathematical model of statorcurrent-error based speed estimator, the following classical and modified MRAS CC systems can be created: classical MRAS CC -when G = 0, ϕ = 0, µ = 0; -MRAS CC+G -when ϕ = 0, µ = 0 and elements of the gain matrix G, are calculated according to [35,36]: g sd = k r r l r , g sq = −kω r , g rd = − r s k r , g rq = −l r k rωr (17) where: k-parameter to be selected appropriately,ω r = i sy ψ 2 ref k r r r -is the value of the rotor slip angular velocity calculated in the field-oriented control structure of the IM drive system; i sy -the torque component of the stator current vector in the synchronously rotating reference frames (x-y); -MRAS CC+ϕ -when G = 0, µ = 0 and the shift angle, ϕ in the speed adaptation mechanism (13) and (14) is calculated as follows [32,33]: -MRAS CC+µ -when G = 0, ϕ = 0 and the auxiliary variable µ is on-line adapted according to (15) and (16) [36].
All modified estimators MRAS CC+G , MRAS CC+ϕ and MRAS CC+µ developed in continuous-time domain are stable in the motoring and regenerating modes, as it was proved in [36]. However, in a drive system implemented in digital technology, the observer algorithm designed in the continuous-time domain must be discretized. The commonly used FE method can lead to instability of the observer at a higher speed range, especially for larger sampling steps. Therefore, the discretization method and sampling step should be carefully selected [23][24][25][26].
In the next sections all discussed MRAS CC -type speed estimators will be discretized, their performance and accuracy will be compared depending of the discretization method used and next tested in the discrete implementation of the whole DRFOC drive system.

General Remarks
As it is mentioned in the introduction, there are two ways to implement a continuous-time models and control methods in microprocessor-based control systems: (1) Discrete design method, which consists in obtaining a discrete-time model of the object based on the continuous-time domain model and designing next the control law in the discrete-time domain [12][13][14]; (2) Analogue emulation, where the entire control system is first designed in a continuous time-domain and next the control law (regulator or/and state estimator if necessary) is discretized using selected numerical integration method [15][16][17].
The disadvantage of the first method consists in difficulties of the derivation a digitized model of an object, especially for nonlinear systems, as mentioned in the introduction. In contrast, the model obtained using the second method may not have high performance due to restrictions on the lower value of sampling time and the numerical method that can be used. However, it is much simpler to use in various applications.
In analogue emulation, various methods of numerical integration are used, such as simple Euler methods (Forward Euler-FE, Backward Euler-BE) or bilinear method-Tustin (TU). The TU method gives much more accurate results than both Euler methods, but is much more complicated in terms of computational complexity. Therefore other methods have been proposed, as e.g., modified Euler method [9,10,24].

Numerical Integration Methods
As mentioned above, in order to implement the continuous control law in a discrete system, the selected method of approximation of differential equations should be used. For the numerical methods listed above, they take the following form [10,26]: for the forward Euler (FE) method: for the backward Euler (BE) method: for the Tustin (TU) method: for the modified Euler (ME): with: where: T p -a sampling time.
If the mathematical model of the object is written in the state equation form ( . x = Ax + Bu), its discretized version is as follows: Energies 2020, 13, 2595 8 of 23 -for the FE method: for the BE method: for the TU method: for the ME method: All approximations (24)-(27) of the state equation of the controlled object can be written in a form: where S-is an object state matrix in a discrete form, which depends on the discretization method selected. This matrix can then be used to analyze the position of the eigenvalues (poles) of the object based on the characteristic equation of the discrete model (in z-domain) [9,10]: The presented discretization methods will be used in this paper for the performance analysis of the MRAS CC speed estimators in discrete-time domain, including their stability in a wide speed range.

A Brief Description of the Drive Control Structure and Test Bench
All versions of the MRAS CC -type speed estimators, classical and those with the extended stability range (in regenerating mode) have been tested in the discrete implementation in the DRFOC structure, presented in Figure 3. Four PI (proportional-integral) controllers are used to control motor speed, rotor flux amplitude and stator current components (expressed in the x-y reference coordinate system, rotating synchronously with the rotor flux vector, with ω s electrical angular velocity). In this control structure the amplitudeψ r and the positionγ ψr of the rotor flux vector are obtained using well-known current model of the rotor flux [2]. The decoupling block is used to obtain the independent control of the rotor flux and motor torque, and its inputs are the i sx and i sy components, rotor flux angle and its amplitude. The subordinated PI stator current i sx and i sy components generate the control signals f x and f y , which supported by the decoupling signals e x and e y define the reference values of the stator voltage u sx and u sy components. Next, they are transformed to the stationary frame, to be useful for the space voltage modulator (SVM). The voltage vector components are calculated using the DC-bus voltage u DC and the PWM signals S A , S B , S C .
Space Vector Modulator (SVM) was used for controlling the switches of VSI, both in simulation and experimental tests. The nominal data of the tested IM are given in the Appendix A.
Energies 2020, 13, 2595 9 of 23 vector are obtained using well-known current model of the rotor flux [2]. The decoupling block is used to obtain the independent control of the rotor flux and motor torque, and its inputs are the i sx and i sy components, rotor flux angle and its amplitude. The subordinated PI stator current i sx and i sy components generate the control signals f x and f y , which supported by the decoupling signals e x and e y define the reference values of the stator voltage u sx and u sy components. Next, they are transformed to the stationary frame, to be useful for the space voltage modulator (SVM). The voltage vector components are calculated using the DC-bus voltage u DC and the PWM signals S A , S B , S C . In simulation studies, obtained using in the MATLAB/Simulink environment, the entire control structure was modelled in a discrete-time domain with a switching frequency equal to the sampling frequency. The IM and voltage source inverter (VSI) are modelled in simulation studies in "continuous-time" domain (with a sampling time of 1e-6 s), while in the experimental test both are real objects. The laboratory drive system consists of two IMs (1.1 kW-working as a driving motor and 1.5 kW-as a load torque generator), powered from separate VSIs. The control structure was digitally implemented using the Rapid Prototyping System dSpace1103. Speed measurement was carried out using an incremental encoder, and current and voltage measurements using current and voltage sensors of the LEM type. The voltage inverters were controlled with a switching frequency equal to the sampling frequency, as well. Due to the limitations of the load torque generation on the laboratory set-up, tests were performed up to 150% of nominal torque. A general view of the laboratory test bench is presented in Figure 4. In simulation and experimental studies, the algorithms of the rotor flux and speed estimators were implemented in stationary reference frame αβ, (ω k = 0 in mathematical models of the respective MRASs). In simulation studies, obtained using in the MATLAB/Simulink environment, the entire control structure was modelled in a discrete-time domain with a switching frequency equal to the sampling frequency. The IM and voltage source inverter (VSI) are modelled in simulation studies in "continuous-time" domain (with a sampling time of 1e-6 s), while in the experimental test both are real objects. The laboratory drive system consists of two IMs (1.1 kW-working as a driving motor and 1.5 kW-as a load torque generator), powered from separate VSIs. The control structure was digitally implemented using the Rapid Prototyping System dSpace1103. Speed measurement was carried out using an incremental encoder, and current and voltage measurements using current and voltage sensors of the LEM type. The voltage inverters were controlled with a switching frequency equal to the sampling frequency, as well. Due to the limitations of the load torque generation on the laboratory set-up, tests were performed up to 150% of nominal torque. A general view of the laboratory test bench is presented in Figure 4. In simulation and experimental studies, the algorithms of the rotor flux and speed estimators were implemented in stationary reference frame α-β, (ω k = 0 in mathematical models of the respective MRASs).

Impact of the Numerical Integration Methods on the Stable Operation Range of the Classical MRAS CC Estimator
In order to simplify the stability analysis of the mathematical model of the MRAS CC -type estimators after the discretization, the state Equations (10)- (12) can be rewritten in the following

Impact of the Numerical Integration Methods on the Stable Operation Range of the Classical MRAS CC Estimator
In order to simplify the stability analysis of the mathematical model of the MRAS CC -type estimators after the discretization, the state Equations (10)- (12) can be rewritten in the following modified form: where the general modified state matrix of the estimator is: and: The influence of the numerical integration methods on the stable operation range of the classical MRAS CC estimator has been tested. First the poles placement of this estimator discretized using discussed methods was analyzed for different sampling times (T p = 0.125 ms, T p = 0.25 ms and T p = 0.5 ms) in a wide speed range ω m = (0 ÷ 10)ω mN , for the estimator model realized in stationary (ω k = 0), α-β and field-oriented (ω k = ω s ), x-y reference frames. The sample results for α-β model are presented in Figure 5.

Impact of the Numerical Integration Methods on the Stable Operation Range of the Classical MRAS CC Estimator
In order to simplify the stability analysis of the mathematical model of the MRAS CC -type estimators after the discretization, the state Equations (10)- (12) can be rewritten in the following modified form: (30) where the general modified state matrix of the estimator is: and: The influence of the numerical integration methods on the stable operation range of the classical MRAS CC estimator has been tested. First the poles placement of this estimator discretized using discussed methods was analyzed for different sampling times (Tp = 0.125 ms, Tp = 0.25 ms and Tp = 0.5 ms) in a wide speed range m = (0  10)mN, for the estimator model realized in stationary (k = 0), αβ and field-oriented (k = s), x-y reference frames. The sample results for α-β model are presented in Figure 5. In the first row of this figure the poles placement of this estimator is presented in the surrounding of the unit circle and BE method circle, for the considered speed range and for different sampling times. In the second row the enlarged parts of these drawings are presented to show the exact location of the observer eigenvalues in a surrounding of the stability border. As expected from the literature analysis, it can be seen that FE discretization leads to unstable estimator operation In the first row of this figure the poles placement of this estimator is presented in the surrounding of the unit circle and BE method circle, for the considered speed range and for different sampling times. In the second row the enlarged parts of these drawings are presented to show the exact location of the observer eigenvalues in a surrounding of the stability border. As expected from the literature analysis, it can be seen that FE discretization leads to unstable estimator operation (eigenvalues location outside the unit circle) at ever lower drive speeds as the sampling step increases. The BE and TU methods ensure completely stable operation in the whole tested range; observer eigenvalues are located along the unit circle-TU or inside-BE. However, it should be emphasized that these methods require a large amount of calculations related to additional matrix products and the inverse matrices Equations (25) and (26). The last method-ME, can also cause unstable operation for certain speed values, however, it should be noted that these values are large enough to be considered unreachable under normal operating conditions of the drive. In addition, this method is much simpler to implement than the BE and TU methods, as it does not require calculating the inverse of the system state matrix Equation (27). This makes the implementation much easier.
It is also worth noting that the estimator stability is also affected by the coordinate system in which it was modelled; it is illustrated in the Table 1. From the presented results it can be concluded that modelling the estimator in the synchronously rotating x-y coordinate system provides much more stable operation than in the case of using the stationary α-β coordinate system, which is applied in most cases. Thus, if low sampling frequencies are taken into account, the MRAS CC speed estimators should be modelled in the synchronous frame. As it results from the Table 1, for each of the sampling steps the influence of the coordinate system is significant, i.e., in the case of a stationary α-β coordinate system the MRAS CC estimator loses stability already at a much lower speed value than in the case of a reference system x-y rotating synchronously with the machine field. However, the implementation of the discrete algorithm of the MRAS CC estimator in the x-y coordinate system requires the calculation of the actual value of the slip pulsation, which is usually burdened with an error resulting from the variability or incorrect identification of IM parameters. However, it can be seen that when using the ME method to discretize the estimator's mathematical model, the range of stable operation of this algorithm in the α-β system is much larger, even for a relatively big sampling step value. The next Figure 6 presents the results of the estimator tests for each of the discretization methods, operating in an open-loop, in the motoring mode, for reference speeds ω m = {0.2, 0.5, 0.8, 1.1, 1.4, 1.7} p.u., with no load. These results confirm the previous theoretical analysis of the unstable work of the estimator digitized using the FE method. As can be seen from the graphs in Figure 6, this estimator loses stability earlier than it results from Table 1, because its work is also influenced by the speed adaptation mechanism, which was not taken into account when analyzing the location of the estimator poles. The next Figure 6 presents the results of the estimator tests for each of the discretization methods, operating in an open-loop, in the motoring mode, for reference speeds ωm = {0.2, 0.5, 0.8, 1.1, 1.4, 1.7} p.u., with no load. These results confirm the previous theoretical analysis of the unstable work of the estimator digitized using the FE method. As can be seen from the graphs in Figure 6, this estimator loses stability earlier than it results from Table 1, because its work is also influenced by the speed adaptation mechanism, which was not taken into account when analyzing the location of the estimator poles. Additionally, as the sampling step increases, the estimator discretized by the FE method loses stability the faster the larger the sampling step. Even in the conditions of stable operation for the FE and BE methods, a fixed speed estimation error appears, which increases with the increase of the sampling step and motor speed value. Depending on the discretization method used, the estimated speed value is higher or lower than the actual value. On the contrary, for TU and ME methods, the accuracy is almost always perfect (only for higher speeds at a fairly large sampling step, a very small error of tenths of percent appears- Figure 6c, row III).
In the following Figure 7 the estimation of the rotor flux magnitude is presented. It shows the same features that can be seen when estimating the rotor speed. However, in this case, for the stable operating range of the estimator discretized using FE and BE methods, the flux estimation error is Additionally, as the sampling step increases, the estimator discretized by the FE method loses stability the faster the larger the sampling step. Even in the conditions of stable operation for the FE and BE methods, a fixed speed estimation error appears, which increases with the increase of the sampling step and motor speed value. Depending on the discretization method used, the estimated speed value is higher or lower than the actual value. On the contrary, for TU and ME methods, the accuracy is almost always perfect (only for higher speeds at a fairly large sampling step, a very small error of tenths of percent appears- Figure 6c, row III).
In the following Figure 7 the estimation of the rotor flux magnitude is presented. It shows the same features that can be seen when estimating the rotor speed. However, in this case, for the stable operating range of the estimator discretized using FE and BE methods, the flux estimation error is very large, which indicates a very poor accuracy of these methods (except for a small sampling step and low speed values- Figure 7, row I). In the case of the TU and ME methods, estimation of the rotor flux magnitude is almost always ideal-at high speed values and with increasing sampling steps, the value of the steady-state estimation error increases, but very slightly.
Energies 2020, 13, x FOR PEER REVIEW 13 of 22 very large, which indicates a very poor accuracy of these methods (except for a small sampling step and low speed values- Figure 7, row I). In the case of the TU and ME methods, estimation of the rotor flux magnitude is almost always ideal-at high speed values and with increasing sampling steps, the value of the steady-state estimation error increases, but very slightly. Analyzing the waveforms presenting the estimation of the rotor flux component in α axis in Figure 8, it can be stated that the reason for the poor accuracy of the estimation of the rotor flux magnitude is the incorrect estimation of the rotor flux vector components. Figure 8 shows that the estimated value of the rotor flux in the α axis is incorrect both in amplitude and in phase when FE and BE methods are used to discretize the MRAS CC estimator. Of course, as in the case of speed estimation, the accuracy is the higher the smaller the sampling step and the smaller the estimated speed value. For both TU and ME methods, the quality of the rotor flux estimation in the α axis is excellent.
All the properties of four analyzed discretization methods have been collected in Table 2. As the table shows, the FE method is not computationally complex and is simple to implement, however it may cause unstable operation of the digitized system. If the estimator works stably, despite the fact that the quality of the speed estimation is good, the accuracy of the estimation of the rotor flux is very poor. Analyzing the waveforms presenting the estimation of the rotor flux component in α axis in Figure 8, it can be stated that the reason for the poor accuracy of the estimation of the rotor flux magnitude is the incorrect estimation of the rotor flux vector components. Figure 8 shows that the estimated value of the rotor flux in the α axis is incorrect both in amplitude and in phase when FE and BE methods are used to discretize the MRAS CC estimator. Of course, as in the case of speed estimation, the accuracy is the higher the smaller the sampling step and the smaller the estimated speed value. For both TU and ME methods, the quality of the rotor flux estimation in the α axis is excellent.
All the properties of four analyzed discretization methods have been collected in Table 2. As the table shows, the FE method is not computationally complex and is simple to implement, however it may cause unstable operation of the digitized system. If the estimator works stably, despite the fact that the quality of the speed estimation is good, the accuracy of the estimation of the rotor flux is very poor.  Compared to the FE method, the BE method presents almost the same quality of the estimation of state variables, however ensures the estimator's stability over the full range of operation. However, due to the need to calculate a large number of products and the inverse matrix, this method is much more computationally complex and difficult to implement in hardware.
The TU method, similarly to the BE method, ensures the stability of the estimator in the full range of work, but it is very complex computationally and difficult in microprocessor implementation. The advantage of this method is a significant improvement in the quality of the estimation of all state variables.
The last method-a modified Euler method-just like TU, provides excellent quality of state variables estimation. On contrary, this method does not require as much calculation as the BE and TU methods (it is not necessary to calculate the inverse matrices), and thus it is not so difficult to implement in hardware. The disadvantage of this method is the possibility of losing the stability of the discrete speed estimator for high speed and relatively big sampling time. It should be noted, however, that the range of stability of the estimator discretized using this method is very wide,  Compared to the FE method, the BE method presents almost the same quality of the estimation of state variables, however ensures the estimator's stability over the full range of operation. However, due to the need to calculate a large number of products and the inverse matrix, this method is much more computationally complex and difficult to implement in hardware.
The TU method, similarly to the BE method, ensures the stability of the estimator in the full range of work, but it is very complex computationally and difficult in microprocessor implementation. The advantage of this method is a significant improvement in the quality of the estimation of all state variables.
The last method-a modified Euler method-just like TU, provides excellent quality of state variables estimation. On contrary, this method does not require as much calculation as the BE and TU methods (it is not necessary to calculate the inverse matrices), and thus it is not so difficult to implement in hardware. The disadvantage of this method is the possibility of losing the stability of the discrete speed estimator for high speed and relatively big sampling time. It should be noted, however, that the range of stability of the estimator discretized using this method is very wide, because theoretically determined speed limit for which the MRAS CC estimator discretized by ME method may work unstable is very large (see Table 1).
To sum up, it can be stated that the ME method can successfully replace not only the simple FE method, which causes problems with the stability of the estimator for relatively low speed values, but also more computationally complex BE and TU methods.
In this section the four most popular numerical integration methods are analysed. Although some of them are considered to be highly complex regarding the calculation effort and difficult in microprocessor implementation, as shown in Table 2, it must be mentioned here that they are still much simpler in terms of the complexity and implementation than the approximate sample data models of induction motor.

Simulation Tests of the Discrete Realization of all Tested Speed Estimators
Simulation tests included the analysis of the mechanical characteristics of the drive with classic and all modified MRAS CC -type estimators in the motoring and regenerating operation as well as the properties of these estimators during a slow reverse. This second test was aimed at analyzing the behavior of estimators in the range of near zero speed and during the transition of the drive from motoring to regenerating mode. All estimators were discretized using FE and ME methods, to validate the results and conclusion presented in the previous section.
The mechanical characteristics of the drive system with the speed estimated by a given MRAS were obtained with the methodology illustrated in Figure 9, where an example for low speed of the drive system is shown. When the motor reaches the set speed, for example 0.1 p.u. Figure 9b, the load torque after t = 5 s decreases gradually from zero value to −1.5 of its nominal value Figure 9a within 15 s. Due to the constant speed, motor and load torques can be assumed to be equal and the drive operates in the regenerating mode. The process ends at t = 20 s. Then, based on obtained speed and torque, appropriate dynamic mechanical characteristic is created Figure 9c. It should be noted that in the case of stable operation, the estimated speed coincides with the actual motor speed, and the mechanical characteristic is a straight line. As can be seen in Figure 9b,c, the classic MRAS CC loses the stability, while MRAS CC+ϕ maintains stable operation. According to the methodology described above, the torque-speed characteristics were determined for all discrete MRAS CC -type estimators and they are demonstrated in Figure 10 for various set-values of the rotor speed. This figure also compares the effects of two discretization methods, FE and ME, and two sampling times (T p = 0.125 ms or T p = 0.250 ms). The D 1 line in the Figure 10 is the unobservability line of induction motor (the synchronous angular speed ω s is equal to zero across the D 1 line [30,32]). The D 2 line appears only in the case of the classical MRAS CC estimator and shows the theoretical stability border line of this estimator [32,33,35]. The characteristics in Figure 10 are light green when the system is stable and they are dark green when the system becomes unstable.
Energies 2020, 13, x FOR PEER REVIEW 16 of 22 According to the methodology described above, the torque-speed characteristics were determined for all discrete MRAS CC -type estimators and they are demonstrated in Figure 10 for various set-values of the rotor speed. This figure also compares the effects of two discretization methods, FE and ME, and two sampling times (Tp = 0.125 ms or Tp = 0.250 ms). The D1 line in the Figure 10 is the unobservability line of induction motor (the synchronous angular speed ωs is equal to zero across the D1 line [30,32]). The D2 line appears only in the case of the classical MRAS CC estimator and shows the theoretical stability border line of this estimator [32,33,35]. The characteristics in Figure 10 are light green when the system is stable and they are dark green when the system becomes unstable. The mechanical characteristics presented in Figure 10 confirm the analysis presented in the previous section. For the FE method, in the case of a smaller sampling step, the classic MRAS CC estimator works unstable before as well as just after switching on the load torque in the motoring operating range for a speed of 0.8 p.u. Figure 10a row II. This is related to the impact of the discretization method on the estimator stability. In the case of discrete version of the classic MRAS CC estimator Figure 10a-rows I and II, an unstable area in regenerating mode can be seen in all four graphs, which confirms earlier analyzes of this estimator in the continuous-time domain [35][36][37]. However, it should be noted that with the same sampling step and using the ME method, more accurate results are obtained over the entire working range tested than with the FE method. First of The mechanical characteristics presented in Figure 10 confirm the analysis presented in the previous section. For the FE method, in the case of a smaller sampling step, the classic MRAS CC estimator works unstable before as well as just after switching on the load torque in the motoring operating range for a speed of 0.8 p.u. Figure 10a row II. This is related to the impact of the discretization method on the estimator stability. In the case of discrete version of the classic MRAS CC estimator Figure 10a-rows I and II, an unstable area in regenerating mode can be seen in all four graphs, which confirms earlier analyzes of this estimator in the continuous-time domain [35][36][37]. However, it should be noted that with the same sampling step and using the ME method, more accurate results are obtained over the entire working range tested than with the FE method. First of all, the estimated speed calculated using the ME method does not oscillate as much as using the FE method. Moreover, with the FE method and a lower switching frequency (4 kHz) there is unstable operation of the estimator even in the motoring mode, when the load torque is turned on Figure 10a-row I and II, but the estimator stabilizes when the load torque begins to increase.
Despite the analysis of the waveforms shown in Figure 2, the MRAS CC+G and MRAS CC+µ estimators discretized using the FE method do not lose stability at 0.8 p.u. Figure 10b,d, row II, in contrast to MRAS CC+ϕ Figure 10c, row II. Additionally, it should be noted that for MRAS CC+G and MRAS CC+ϕ discretized using the FE method, along with the increasing speed value and decreasing sampling step (in the case of MRAS CC+G also with increasing load torque)-the steady-state error value of the estimated speed increases. For the MRAS CC+µ estimator, both the FE method and the sampling step do not affect the accuracy of the estimation, which is very good. For all modified MRAS CC estimators discretized using the ME method, stable operation is ensured in the full operating range, regardless of the sampling step used.
In Figure 11 the behavior of the MRAS CC estimators under slow speed reverse of the drive system was shown, with nominal load torque switched on in 2.75 s. All properties resulting from the analysis of mechanical characteristics Figure 10 are also revealed during reverse operation.
Energies 2020, 13, x FOR PEER REVIEW 17 of 22 all, the estimated speed calculated using the ME method does not oscillate as much as using the FE method. Moreover, with the FE method and a lower switching frequency (4 kHz) there is unstable operation of the estimator even in the motoring mode, when the load torque is turned on Figure  10a-row I and II, but the estimator stabilizes when the load torque begins to increase. Despite the analysis of the waveforms shown in Figure 2, the MRAS CC+G and MRAS CC+µ estimators discretized using the FE method do not lose stability at 0.8 p.u. Figure 10b,d, row II, in contrast to MRAS CC+ϕ Figure 10c, row II. Additionally, it should be noted that for MRAS CC+G and MRAS CC+ϕ discretized using the FE method, along with the increasing speed value and decreasing sampling step (in the case of MRAS CC+G also with increasing load torque)-the steady-state error value of the estimated speed increases. For the MRAS CC+µ estimator, both the FE method and the sampling step do not affect the accuracy of the estimation, which is very good. For all modified MRAS CC estimators discretized using the ME method, stable operation is ensured in the full operating range, regardless of the sampling step used.
In Figure 11 the behavior of the MRAS CC estimators under slow speed reverse of the drive system was shown, with nominal load torque switched on in 2.75 s. All properties resulting from the analysis of mechanical characteristics Figure 10 are also revealed during reverse operation. The MRAS CC estimator behaves unstable after switching to the regenerating mode and crossing the unobservability line ωs = 0 (see for details on stability issues of these estimator in [32,35])-regardless of the discretization method used. In the case of modified MRAS CC -type The MRAS CC estimator behaves unstable after switching to the regenerating mode and crossing the unobservability line ω s = 0 (see for details on stability issues of these estimator in [32,35])-regardless Energies 2020, 13, 2595 18 of 23 of the discretization method used. In the case of modified MRAS CC -type estimators, special attention should be paid to the fact that for the FE method, MRAS CC+G and MRAS CC+ϕ introduce a fixed error Figure 11, columns (b) and (c). In addition, the MRAS CC+ϕ behaves unstably for a smaller sampling step after establishing the speed at the rated value in the regenerating mode Figure 11c, rows III and IV. In the case of the MRAS CC+µ estimator, the work is almost perfect, regardless of the discretization method and sampling step used. The estimated speed value oscillates only around crossing the observability line. These oscillations disappear very quickly.
To assess the accuracy of speed estimation by individual MRAS CC estimator, discretized using the numerical methods analyzed, the Integral of Time multiplied by Absolute Error (ITAE) criterion was used: where: ω m ,ω m -real and estimated (by each MRAS CC ) rotor speeds, e ω -speed estimation error, T c -calculation time.
The ITAE criterion for all analyzed discrete estimators was calculated for the studied reversal operation. These values have been collected and presented in the Table 3. Table 3. ITAE criterion values for slow speed reversal of Figure 11. All the above analyses are confirmed by the results obtained for the ITAE criterion. It should be emphasized that in the most cases the ME method provides more accurate operation of the estimator. In addition, it is also worth noting that the MRAS CC modification regarding the introduction of an auxiliary variable µ allows the stability and quality of speed estimation to be improved even for an unstable range of work resulting from the use of the FE method. The MRAS CC+µ estimator provides the most accurate work of all discrete speed estimators tested.

Conclusions
In this article the influence of four numerical integration methods: forward and backward Euler, Tustin and Modified Euler methods, on the discretization of classical and three modified As can be seen from the speed waveforms ( Figure 12, row I and III), each of the methods for improving MRAS CC stability derived in a continuous-time domain can be used in the discrete model of the MRAS CC estimator, because all methods ensure stable operation of the modified speed estimators. However, it is worth noting, in particular for ω m = 0.6 p.u., that the method with the use of an auxiliary variable µ (Figure 12, row II and IV) seems to be the most accurate, and is characterized by the smallest speed estimation error. Table 4 below presents the ITAE criterion values for ω m = 0.2 p.u. and ω m = 0.6 p.u., respectively. It should be noted that the increase in speed and increase of the sampling step increase the estimation inaccuracy using discrete speed observers. However, it is clearly seen that the use of the ME method allows a much more accurate speed estimation than when using the FE method to discretize the mathematical model of the estimator, with slightly bigger computational effort.

Conclusions
In this article the influence of four numerical integration methods: forward and backward Euler, Tustin and Modified Euler methods, on the discretization of classical and three modified MRAS CC -type speed estimators, with improved stability in the regenerating operation mode is discussed and compared. As it results from the presented simulation and experimental research, the ME discretization method is much more accurate than both FE methods and gives almost as accurate results as the bilinear TU method, however it is much less complicated in practical implementation. ITAE criteria values calculated for various tests of four discrete MRAS CC -type estimators confirm these properties. It has been shown that the ME method can successfully replace the inaccurate and unstable FE method in a wide range of IM drive speeds at a relatively low sampling frequency, which is especially important for higher drive system speeds.
The most important result of the presented research is the conclusion that the methods of improving stability developed for mathematical models in the continuous-time domain of the MRAS CC -type speed estimators can be effectively used in their discrete models if analogue emulation methodology with a properly selected discretization method and sampling time is used.
It should be highlighted that all modified estimators are characterized by the stable operation in the whole range of IM speed changes, and in particular in the regenerating mode, also after discretization using the TU or ME method. However, it should be also noticed that the discretized MRAS CC+µ estimator provides the most accurate operation of all discrete speed estimators tested.
Considering that the analogue emulation is one of the methods of systems' discretization, the authors plan in future works to use a discrete design method, based on an approximate sampled-data model, to create a discrete-time version of the MRAS CC estimator. It will be verified whether this method gives different or better results and allows a stable operation of the drive system in a wide range of drive speed changes at a relatively low sampling rate.