A Driftless Estimation of Orthogonal Magnetic Flux Linkages in Sensorless Electrical Drives

Estimations of magnetic flux linkages, either between the stationary windings of the stator for the direct torque control (DTC), or between the stationary windings and the rotor for the sensorless field-oriented control (FOC), are based on integration of corresponding voltages. Integration of voltages with offsets that come from improperly calibrated measurements as well as from transient states generally produces unwanted drifts in the resulting magnetic flux linkages, which when used within any type of control of sensorless electrical drives results in instability. This paper addresses that problem and proposes a simple self-contained solution based on orthogonal properties of waveforms of input voltages and resulting magnetic flux linkages in the frame of reference fixed to the geometry of the stator. The proposed solution requires only two periodic orthogonal input waveforms with a distinct common fundamental harmonic, which as such is independent of the type and parameters of the used machine. The idea of the proposed solution is presented analytically, its stability is proven by means of the quadratic Lyapunov theory, and its functionality is demonstrated by standalone simulations and experiments within the sensorless FOC of a permanent magnet synchronous machine (PMSM).


Introduction
Variable frequency electrical drives (VFED) are typically operated either by the direct torque control (DTC) or the field-oriented control (FOC), where the FOC requires the angular position of the rotor for its operation. The angular position of the rotor is obtained via a position sensor that requires an additional space and wiring that, besides an additional cost, might not be applicable in certain environments. Thus, the demand for more robust and affordable electrical drives has spurred the sensorless FOC that can be found in household and industrial applications. As it is described in [1], the DTC is inherently a sensorless control that for its operation requires the total magnetic flux linkage of the stationary windings of the stator (often referred to as stator flux linkage), while some versions of the sensorless FOC, as those described in [2,3], are based on estimations of the magnetic flux linkage between the stationary windings and the rotor (often referred to as rotor flux linkage).
Estimations of magnetic flux linkages are generally based on integration of corresponding voltages that typically have offsets caused by improperly calibrated measurements as well as transient states because of which the resulting waveforms of those magnetic flux linkages over time horizontally drift away. A great majority of solutions to that problem is addressed in the frame of reference fixed to the geometry of the stator ([αβ]), most of which are in applications with squirrel cage induction machines (SCIMs). The simplest solution presented in [4] is based on the subtraction of the accumulated average value of the resulting waveforms of the magnetic flux linkages from themselves cyclically every period. While this solution works in steady states and for small transients, its dynamics are severely

An Electromagnetic Model of a PMSM
Three-phase SCIMs and PMSMs are the two most commonly utilized types of electrical machines. Since a PMSM compared to a SCIM has a fixed position of the magnetic flux of the rotor and fewer parameters required for its detection without a position sensor, its electromagnetic model is used for demonstrations in this paper. For that purpose, under the assumption of Y-connected stationary windings of the stator, an electromagnetic model of a permanent magnet synchronous machine can be represented by its voltage equation in terms of the spatial phasor of the voltages across the stationary windings (v s ), the spatial phasor of the electrical currents in the stationary windings (i s ), the spatial phasor of the total magnetic flux linkage of the stationary windings (λ s ), and the electrical resistance of each stationary winding (R s ). This assumption does not limit the following model from being applied to ∆-connected stationary windings, it only requires the ∆ to Y transformation of the equivalent electromagnetic network. The general form of a spatial phasor of electrical or magnetic quantities of the stationary windings ( f s ) based on scalar functions of electrical or magnetic quantities of the stationary windings ( f sa , f sb , f sc ) can be expressed via the base of natural logarithms (e) and the imaginary unit (j) in terms of the electrical angular position of the rotor (θ e ), which is observed as a function of time (t), and an arbitrary electrical angular position (ρ e ) as The quantities of the voltage equation expressed in [αβ] are obtained from Equation (1) for ρ e = θ e , whereby v s can be expressed via R s , i s , and λ s . Thus, the projection of i s onto [αβ] (i s [αβ] ) is defined in terms of the magnitude of the fundamental harmonic of each electrical current in each of the stationary windings (i s,1 ), θ e , and the electrical angular displacement of a direct axis of the synchronous magnetic flux of the stator from the subsequent direct axis of the magnetic flux of the rotor (δ e ) as where in the motoring mode δ e is commonly referred to as the torque angle, while in the generating mode it is referred to as the load angle. To account for the magnetic saliency of the machine by assuming the cross-sectional symmetry of the stationary windings, the projection of λ s onto [αβ] (λ s[αβ] ) can be defined in terms of the stray inductance of each stationary winding (L sσ ), the magnitude of the zeroth harmonic with respect to θ e of the salient inductance of each stationary winding (L s,0 ), and the magnitude of the second harmonic with respect to θ e of the salient inductance of each stationary winding (L s,2 ) including i s[αβ] , the complex conjugate of i s[αβ] (ī s[αβ] ), and the magnetic flux linkage between the stationary windings and the permanent magnets of the rotor (λ m ) as where λ m is commonly referred to as the rotor flux constant. The derivations of L sσ , L s,0 , and L s,2 are available in [24]. Based on R s and the definitions of i s[αβ] and λ s[αβ] , provided by Equations (2) and (3)

The Angular Position and Speed of the Rotor
Since θ e points to the direct axis of the magnetic flux of the rotor that coincides with the angular position of the last term in Equation (3), this term needs to be expressed preferably together with any other coincident phasorial components. For that purpose, the inductive terms next to i s [αβ] in Equation (3) can be denoted by the total amount of L s,0 with L sσ in [αβ] (L s[αβ],0 ) defined as Since the conjugate of the product of two complex numbers is equal to the product of the conjugates of those two complex numbers,ī s [αβ] in Equation (3)  (t) e jθ e (t) + λ m e jθ e (t) .
The real parts of Equation (9) (t) sin(θ e (t)) + λ m cos(θ e (t)), (10) while the imaginary parts define the quadrature component of λ (t) cos(θ e (t)) + λ m sin(θ e (t)). (11) According to the relationship defined by Equation (7) as well as in terms of i s[β] and i s[d] as By substituting Equation (12) into Equation (10) and Equation (13) into Equation (11), λ s[αβ] can be expressed as Moreover, L s[αβ],0 and L s[αβ],2 can be expressed in terms of the direct synchronous inductance of the stationary windings (L s [d] ) and the quadrature synchronous inductance of the stationary windings (L s[q] ), which are typically used because they are easily measurable. Thus, L s[αβ],0 can be expressed as and L s[αβ],2 as Based on Equations (15) and (16), Equation (14) can be rewritten as where the projection of the component of λ s whose argument equals θ e onto [αβ] (λ sθ [αβ] ) can be denoted as and the electrical angular speed of the rotor (ω e ) as According to Equations (17) and (18), λ s[αβ] can be expressed as from where can λ sθ[αβ] based on Equation (4) be expressed as This way of expressing θ e and ω e eliminates the need for L s [d] , which simplifies the observer that can be constructed based on Equations (19), (20) and (22) as it is shown by the diagram in Figure 1. The simplicity and filtering properties of integration make this observer highly desirable and frequently used in selsorless applications of PMSMs whose examples can be found in [25][26][27][28], however, the drift caused by the integration represents the major obstacle in its practical implementations. In the case of significant changes in R s that might be caused by extreme changes in temperature as well as in L s[q] that, depending on the design of the machine, might saturate by high currents, actual values of R s and L s[q] can be estimated using the method described in [29]. The overall diagram of the sensorless FOC for PMSMs with a voltage source inverter (VSI) controlled by means of the space-vector modulation (SVM) that is regulated by proportional-integral controllers (PIs) is shown in Figure 2.

The Problem of the Drift
The drift can be defined as the accumulated offset caused by integration of an offset in the voltage represented by the integrand in Equation (22). Such an offset is primarily caused by initial conditions at transient states and additionally often by improperly calibrated shunts whose influence should be and typically is negligible. To present the problem in [αβ], it is sufficient to consider v s[αβ] and i s[αβ] in a steady state. Therefore, v s[αβ] with an offset can be defined in a steady state according to Equation (20) in terms of the magnitude of the zeroth harmonic of v s[αβ] (v s[αβ],0 ), which represents an offset in v s[αβ] , the magnitude of the fundamental harmonic of each voltage across each of the stationary windings (v s,1 ), ω e , δ e , and the angular shift in the phase of Similarly, i s[αβ] with an offset can be defined in a steady state based on Equation (2) by adding the magnitude of the zeroth harmonic of i s[αβ] (i s[αβ],0 ), which represents an offset in i s [αβ] , and replacing θ e according to Equation (20) by ω e t as Based on Equations (23) and (24), the derivative of λ s[αβ] with respect to t in Equation (4) can be expressed as From Equation (26) (25)

A Solution to the Problem
According to Equation (26), λ s[α] is defined as while λ s[β] is defined as where λ s[α] (0) and λ s[β] (0) are initial conditions. Since the trigonometric terms in Equations (28) and (29) Similarly, the trigonometric terms in Equations (27) and (30) Hence, the compensated derivative of λ s[α] with respect to t is then obtained as while the compensated derivative of λ s[β] with respect to t is Due to causality, ω e needed for the proposed compensation cannot be obtained by Equation (20), instead, it needs to be calculated prior to the compensation. Thus, ω e is calculated based on v sλ[α] and v sλ [β] as and as such is used instead of Equation (20) for the sensorless FOC. Since differentiation amplifies noise, it is practically always implemented in a cascade with an LPF. Thus, by expressing Equation (35) in the form of a transfer function in the Laplace domain in terms of the cut-off frequency of the low-pass filter for the filtering of ω e (ω c ) as which can be rearranged to the form and observed as a closed-loop system, the observer of ω e can be constructed as it is shown in Figure 3. Moreover, based on the diagram in Figure 3, the observer of ω e can be expressed in the time domain as ω e (t) = ω c The error.
The wrapping of the error between −π and π. , where the resulting error needs to be wrapped between −π and π to match the range of the four-quadrant inverse tangent function. The proposed compensation of the drift, presented by Equations (31)-(34) including Equation (38), represents a multiple input-multiple output (MIMO) phase-locked loop (PLL), which for an easier understanding is graphically presented by the diagram in Figure 4.

The Stability of the Proposed Compensation of the Drift
From Equations (31)-(34) it is obvious that the proposed compensation of the drift represents a nonlinear time-varying system that can be described by a pair of differential equations, namely and where the nonlinearity defined by Equation (38) is represented by ω e , which varies with t. Due to their specific structure, Equations (39) and (40) can be observed in a state-space as a linear parameter-varying system. Thus, by defining the state vector (x) as and the input vector (u) as as well as the state matrix (A) as and the input matrix (B) as the observed system can be presented in the state-space form d dt According to [30], the quadratic stability of the observed system for all uncertainties in ω e , where ω e is bounded, can be proven in terms of a square matrix (P) by constructing a quadratic Lyapunov function candidate (V) in the form so that a quadratic Lyapunov equation (Q) fulfills the condition −Q(ω e (t)) = A T (ω e (t)) P + PA(ω e (t)) < 0, where which is only true if all leading principal minors of P are positive. Accordingly, P can be selected in the diagonal form, in which case the element in the first row and the first column of P (P 1,1 ) and the element in the second row and the second column of P (P 2,2 ) must be positive definite, thus By substituting Equations (43) and (49) into Equation (47), the following inequality is obtained which is fulfilled if det(−Q 1,1 (ω e (t))) = − 2k|ω e (t)| k 2 + 1 P 1,1 < 0 (51) and det(−Q(ω e (t))) = 2k|ω e (t)| The first condition defined by Equation (51) is fulfilled if k > 0 and ω e = 0, while the second condition defined by Equation (52) is always fulfilled if k = 0, ω e = 0, and P 1,1 = P 2,2 . Therefore, by choosing P 1,1 = P 2,2 , the system is quadratically stable for arbitrarily fast time-varying uncertainties in ω e if k > 0 and ω e = 0. The validity of these conditions can be checked by Equations (31) and (32), where if ω e = 0, then v sλ[α]-corr = 0 and v sλ[α]-corr = 0 because sgn(0) = 0, which according to Equations (33) and (34) cancels the proposed compensation. Hence, the condition ω e = 0 ensures the existence of v sλ[α]-corr and v sλ[α]-corr , while sgn(ω e ) ensures the negative feedback in Equations (33) and (34). Similarly, k = 0 cancels the proposed compensation, while k < 0 turns the negative feedback into a positive feedback and therefore makes the system unstable. Thus, only k > 0 preserves the negative feedback and the stability.

Results of Simulations of the Proposed Compensation
To demonstrate the standalone functionality of the proposed compensation of the drift, without the sensorless FOC and considerations of possible states of the drive, several simulations were made in Matlab/Simulink (R2016a, MathWorks Inc., Natick, MA, USA) using the diagram shown in Figure 5. As test signals, the waveforms of v sλ [α] and v sλ[β] were set to demonstrate the performance of the proposed compensation for changes in their magnitude and frequency. Thus, the initial magnitude from 0 to 0.5 s was set to 0 V, from 0.5 to 3 s was set to 1 V, and from 3 s onward was set to 2 V, while the frequency from 0 to 6 s was set to 10 rad·s −1 and from 6 s onward was set to 20 rad·s −1 , as it is shown in Figure 6

Results of Comparative Simulations of the Proposed Compensation and Referred Methods
For qualitative comparisons of the proposed compensation with the state of the art, the standard solution with LPFs instead of the integrators presented in [5][6][7][8][9][10][11] was chosen as the first reference method, while the previous compensation introduced in [23] and presented by the diagram in Figure 9 was chosen as the second reference method.
In the simulations, the cut-off frequencies of the LPFs in the first reference method were set to 1 rad·s −1 to asymptotically match the magnitude characteristic of an integrator above 1 rad·s −1 in the frequency domain. For the most unbiased comparison with the second reference method, k in the proposed compensation was set to 1 and ω c in both models was set to be equal, since the calculation of ω e in both models is the same. The resulting waveforms of λ s[α] and λ s[β] to the same reference waveforms of v sλ [α] and v sλ [β] in Figure 6 are shown in Figure 10.  [5][6][7][8][9][10][11], the previous compensation presented in [23], and the proposed compensation.
The magnitude of λ s (|λ s |) obtained as is shown in Figure 11, while the angular shift in the phase of λ s with respect to the spatial phasor of the voltage equivalent to the uncompensated derivative of λ s with respect to t (φ λ ) is acquired as and shown in Figure 12.  Figure 11. The results of |λ s | obtained by a comparative simulation of the standard solution with LPFs presented in [5][6][7][8][9][10][11], the previous compensation presented in [23], and the proposed compensation. Figure 12. The results of φ λ obtained by a comparative simulation of the standard solution with LPFs presented in [5][6][7][8][9][10][11], the previous compensation presented in [23], and the proposed compensation.
From the results presented in Figures 10-12, it can be seen that the proposed compensation and the previous compensation presented in [23] have very similar settling times that are inversely proportional to ω e as well as significantly faster than the standard solution with LPFs. In Figure 12 it can be noticed that φ λ of both the proposed compensation and the previous compensation presented in [23] in the steady states is −90 • , as it is expected for integration, while for the standard solution with LPFs that is not the case. The proposed compensation has better filtering properties than the previous compensation in [23] due to its feedback-based structure, no discontinuities, and no divisions, which in fixed-point implementations on a DSP for small divisors can cause overflows and are generally more expensive than multiplications. To demonstrate the filtering properties, a random noise with the magnitude of 5% of the magnitude of v sλ [α] and v sλ[β] was added to them, which resulted in the noisy test signals shown in Figure 13.  Figure 13 are shown in Figure 14, from which it can be noticed that the proposed compensation preserves the dynamic of the previous compensation in [23] while having the filtering properties of LPFs, but without any undesirable influence on |λ s | and φ λ in steady states.  [5][6][7][8][9][10][11], the previous compensation presented in [23], and the proposed compensation, all with the noisy test signals.

Results of Experiments of the Proposed Compensation within the Sensorless FOC of a PMSM
The practical functionality and the robustness of the proposed compensation was tested in the fixed-point 16-bit Q15 implementation within the Sensorless FOC of a PMSM on a MOSFET-based power electronics converter controlled by the Texas Instruments (Dallas, TX, USA) TM4C123BE microcontroller at the switching frequency of 10 kHz. Parameters of the PMSM used in the experiments are listed in Table 1. To test the speed of the convergence of the proposed compensation, a step in n from 0 to 4000 min −1 for k = 0.5 and ω c matching the nominal electrical angular speed of the machine (ω n ) was performed after the alignment of the rotor without any load and the initialization described at the end of Section 2.1. The results of that experiment are shown in Figure 15. Since each settling time in |λ s | is equal to the corresponding one in φ λ , as it is demonstrated in Figures 11 and 12, φ λ can be taken as the measure of the speed of the convergence of the proposed compensation. Thus, from the plot of φ λ in Figure 15, it can be seen that the proposed compensation converges under 20 ms, while the longer settling time of the magnitude of λ sθ[α] and λ sθ [β] , which corresponds to λ m in Table 1 since i s[d]-ref was kept zero, is imposed by the rest of the control system. It should also be noticed that |λ s | is proportional to the quotient of the magnitude of v sλ [α] and v sλ[β] (|v sλ |) and ω e , where the relationship between ω e and n is defined in terms of the number of pole pairs (p p ) as ω e (t) = 2π p p 60 n(t).
Although a startup from a standstill is possible, it is advisable to perform the initialization described at the end of Section 2.1 to avoid initial uncertainties and ensure a smooth start. To demonstrate the influence of variations in ω c on the robustness of the proposed compensation and its influence on the sensorless FOC, three comparative experiments of a step in n from 1000 to 4000 min −1 for k = 0.5 and varying ω c were performed. The results of those experiments are shown in Figure 16, where the error in θ e (θ e-err ) is expressed in terms of θ e and the actual electrical angular position of the rotor (θ e-act ) obtained by a position sensor as θ e-err (t) = θ e (t) − θ e-act (t). (56) From the results presented in Figure 16 it can be noticed that lower values of ω c result in higher overshoots and longer settling times for the same requirements of the sensorless FOC regarding the dynamics. However, overshoots can be reduced by reducing the dynamics of the sensorless FOC, which can be achieved by proportionally reducing the gains of the PIs in Figure 2.  For the demonstration of the influence of variations in k on the robustness of the proposed compensation and its influence on the sensorless FOC, three additional comparative measurements of a step in n from 1000 to 4000 min −1 for ω c = ω n and varying k were done. From the results of those experiments presented in Figure 17, it can be noticed that k mostly influences the amount of oscillations in the system. Those oscillations become more noticeable for smaller values of k when v sλ[α]-corr and v sλ[β]-corr are not sufficiently high enough to adequately compensate the drift. In addition, from the results of φ λ for k = 0.2 in Figure 17, it can be noticed that at lower values of ω e , or equivalently n, the magnitude of oscillations is greater because there is more time per an electric period for the accumulation of an offset that causes the drift and, consequently, oscillations. The robustness of the proposed compensation as well as its influence on the sensorless FOC during loading of the machine was tested by loading the machine with a hysteresis brake. Four comparative experiments for k = 0.5 and ω c = 0.2ω n were carried out in the range of n from 1000 to 4000 min −1 in steps of 1000 min −1 . From the resulting measurements shown in Figure 18, it can be seen that the machine was overloaded to about 125% of the nominal torque (T n ) listed in Table 1, whereby the peaks of the error in φ λ as well as in θ e-err are more prominent at lower values of n because, as previously mentioned, there is more time per an electric period for the accumulation of an offset that causes the drift. Based on the presented results, it can be concluded that the error in φ λ and θ e-err can be reduced by increasing k.
It is important to mention that the measured values of θ e-act were used only as reference values for the evaluation of the performance of the proposed compensation within the sensorless FOC.

Conclusions
This paper has presented a simple self-contained solution that stems from the idea previously introduced in [23], which as such introduces improvements in the performance regarding fixed-point implementations on DSPs as well as filtering capabilities. The results of simulations presented in Section 4 have shown that the proposed compensation combines the dynamics of the previous solution proposed in [23] with the filtering properties of the standard solution with LPFs presented in [5][6][7][8][9][10][11]. Additionally, the results of experiments presented in Section 5 have confirmed the robustness of the proposed compensation and its functionality within the sensorless FOC.
In summary, the proposed compensation: • is efficiently designed to be computationally suitable for inexpensive applications; • for its operation requires only two periodic orthogonal input waveforms with a distinct common fundamental harmonic; • is completely independent of the type and parameters of the used machine; • provides correct values of |λ s | and φ λ in steady states; • has dynamics adjustable by ω c and k, which provide a fine-tuning without affecting |λ s | and φ λ in steady states; • can be used within the DTC and the sensorless FOC, or any other system that satisfies the requirement of orthogonality stated in the second point. the direct synchronous inductance of the stationary windings L s [q] the quadrature synchronous inductance of the stationary windings λ s the spatial phasor of the total magnetic flux linkage of the stationary windings λ s [αβ] the projection of λ s onto [αβ]