Analytical Design of Self-Sensing Control for PMSM Using Quasi-Direct Calculation

The analytical description and parameterization of a self-sensing control (SSC) for an electrical machine is an important step toward easier commissioning of these systems. In this article, the advantages of high bandwidth position estimation via numerical optimization and the filtering characteristics of a phase-locked loop are combined in the quasi-direct (QD) calculation. The QD calculation uses two parameters for estimation. With the help of the maximum possible acceleration of the drive train, an interdependency between these two parameters is derived. The remaining degree of freedom is used to tune the dynamics of the estimation. Using the transfer function of the estimator, which is derived analytically, the parameters of the speed control are selected, and a specified phase-margin is implemented. With the help of the analytical parameterization, no empirical or numerical tuning needs to be done, which is unique for SSC. All results are experimentally validated.


INDEX TERMS
Analytical parameterization, self-sensing control (SCC), permanent magnet synchronous machine (PMSM). (αβ) x vector x in the αβ-reference frame. (dq) x vector x in the dq-reference frame. (dq) x vector x in thedq-reference frame. X (αβ) matrix X in the αβ-reference frame. X (dq) matrix X in the dq-reference frame. X (dq) matrix X in thedq-reference frame. x d , x q scalar x of the d-and q-axis. xd, xq scalar x of thed-andq-axis. x ana , x num analytical and determined variable.

NOMENCLATURE
x AS , x PS variable of the active switching state and passive switching state. x exp , x fit experimental determined and fitted variable. x full , x red variable x with a full and reduced value of acceleration. x max maximum value of x. x r rated value of x. x stab stabilitiy limit of x. (dq) T (αβ) rotation matrix from the dq to the αβ-reference frame. (dq) T (dq) rotation matrix from the dq to thedq-reference frame.x estimated quantity. x measured quantity.
factor of the symmetrical optimum. D damping factor. f PWM PWM frequency. f QD calling rate of the QD calculation. G cor corrected transfer function. G i , G ω transfer function of the current and speed controller. G mech transfer function of the mechanical system. K gain factor. k discrete sample. K M torque proportional factor. L inductance matrix. L differential inductance matrix. 0 l number of iterations. L dq mutual inductance of the d-and q-axis. L dq differential mutual inductance of the d-and qaxis. L d , L q self inductance of the d-and q-axis. L q , L d differential self inductance of the d-and q-axis. n rotational speed. n r rated rotational speed. Laplace transform of the rotor position. γ el , ω el electrical rotor position and angular speed. γ err rotor position estimation error. γ step postion difference over a PWM period. η initial step width. Ψ stator flux linkage. Ψ PM flux linkage of the magnet. Laplace transform of the rotor angular speed. ω angular frequency.

I. INTRODUCTION
The advantages of omitting the speed or position sensor of an electrical machine by means of self-sensing control (SSC), or using it as a redundant speed information, make SSC a topic of major research interest [1]. The focus of this work is to simplify the commissioning of the SSC algorithm, and thus, to make it suitable for a wider range of applications. The results are shown for the low speed of permanent magnet synchronous machine (PMSM), using the anisotropic inductance matrix of the machine for estimation of the rotor position and speed.
The early observers or estimators used in SSC estimate the position and speed of the PMSM based on a phaselocked loop (PLL) [2], [3]. A phase-locked loop (PLL) is easy to implement, but the closed-loop behavior results in an oscillatory second-order transfer function, which is quite difficult to parameterize analytically [4]. Most of times, a numerical solution for the small signal transfer function of the SSC is needed.
In recent years, an approach using numerical optimization to estimate the rotor position and speed of a PMSM has gained interest [5], [6], [7], [8]. When applying numerical optimization, the stator voltage equation of the PMSM is evaluated. The best fit of the estimated rotor position and estimated speed is sought to minimize the residual (dq) r of the stator voltage equation, which should be zero with ideal parameters and correct values (1) A system based on numerical optimization promises a higher bandwidth of position and speed estimation [8], as the best fit of estimated rotor position and estimated speed is obtained within one sampling cycle. However, the downside is that additional filtering of the estimated parameters can become necessary and high computational capacity is required. The aim of our work is to combine the advantages of numerical optimization and those of a PLL. For this purpose, the stator voltage equation is brought into a form, that leads to a simple calculation of the residual (and its partial derivative, which is also needed for the optimization process). For this purpose, current oversampling according to [9] and [10] is used to measure the current slopes within a pulsewidth modulation (PWM) period. As excitation, a square-wave injection (SWI) [11] is used. The estimator based on this specific expression of the residual has been named flux derivative estimator (FDE) [4].
In addition, the numerical optimization is performed with a very simple and easy-to-tune method, instead of classical approaches like the Gauss-Newton algorithm (GNA) or gradient descent method [12]. The method filters the speed similar to a PLL, but can be described analytically. This method can be parameterized by using two parameters and is called quasi-direct (QD) calculation. The resulting transfer function is a first-order low-pass filter and therefore not capable of oscillation by itself.
Compared to the state of the art, the simple description of the transfer function leads to advantages during commissioning. In the majority of the publications, the design of the parameters of the SSC and the speed control is not described [2], [3], [5], [6], [7], [8], [9], [10], [11]. While the publication [13] gives rules for the parametrization of the position estimation, speed control is not addressed. In [14], the achieved bandwidth of speed and position control is given, but not how the parameters of the speed control are obtained. The straightforward design of the speed controller, and thus, the simplification of commissioning is the aim of this work.
To achieve the straightforward design of the speed controller, the analytic description of the QD calculation is used. The mechanical limitation, that speed cannot change faster than maximum torque and moment of inertia allow, is used to derive an interdependency between the parameters of the estimator. The remaining degree of freedom is used to tune the dynamics of the estimation while maintaining its stability. Next, the transfer function of the estimator is used to parameterize the speed controller for a specified amount of damping in closed-loop speed control. As a result, the entire commissioning can be designed analytically. Only the desired dynamics must be selected according to the requirements of the application.
The rest of this article is organized as follows: In the next section, the underlying model of the PMSM is given, as well as a description of the machine under test and the test bench. In Section III, the position estimation via the QD method is derived, as well as the specific residual as used in the FDE. Section IV is dedicated to the analytical derivation of the transfer function and the stability limits of QD, while the results are applied to the design of the speed control loop in Section V. Section VI reports on the experimental validation. Finally, Section VII concludes this article.

II. FUNDAMENTALS
In this section, the fundamentals of the SSC with FDE-QD are described. In addition, the experimental test bench is described.

A. MODEL OF THE PMSM
The estimator is based on the simplified model of the PMSM, as described in [4], which takes only the most important effects into account. Here, only the results of the simplification and the assumptions made are presented. The stator voltage differential equation of the PMSM in the rotor-fixed dq-reference frame results in with The stator terminal voltage is (dq) v, R is the stator resistance, (dq) i is the current in the dq-reference frame, and γ el is the rotor angle. The flux linkage of the permanent magnet is Ψ PM . For the inductances, the abbreviation with the differential inductance (4) and the secant inductance are used. The secant inductance is introduced with using the flux linkage of the stator in the dq-reference frame (dq) Ψ . To transform the differential equation from the statorfixed αβ-reference frame into the dq-reference frame is used. In this model, spatial harmonics of the PMSM have been neglected. For the model of the PMSM in estimated coordinates, the position estimation error is defined. In the estimateddq-reference frame, the differential equation leads to The inductance matrix L (dq) accounting for rotor position estimation errors results in + L dq -sin (2γ err ) cos (2γ err ) cos (2γ err ) sin (2γ err ) with The back electromotive force (EMF) in estimateddqcoordinates is To give an easy-to-implement residual of the estimator, the inductance matrix is linearized for small estimation errors, which yields The differential inductance matrix L (dq) is of the same shape, and is marked with .

B. EXPERIMENTAL TEST BENCH
The parameters of the PMSM used are given in Table 1, and the current-dependent differential inductances are shown in Fig. 2. The measured quantities of parameters are marked with˘.
For the experimental investigation, the estimator is implemented using a rapid prototyping system. It consists of a Xilinx Zynq 7000 SoC combining an ARM Dual Core Cortex A9 with an Artix-7 field programmable gate array (FPGA). The shunt-based current measurement is sampled with a frequency of 1 MSPS. The estimator is executed once per PWM period ( f QD = f PWM = 10 kHz). The speed controller is executed with f speed = 5 kHz. The PWM generator controls the MOSFET-based inverter with a dc-link voltage of 48 V and uses FPGA-based high dynamic dead time compensation [15]. For current control, a field-oriented control (FOC) with maximum torque per current (MTPC) is implemented, as explained in [4]. The control scheme used with SSC is shown in Fig. 1.

III. POSITION ESTIMATION WITH THE FDE USING A QD CALCULATION
In this section, the FDE-QD is presented. The numerical optimization of the QD calculation is presented and the residual function is derived.

A. PRINCIPLE OF NUMERICAL OPTIMIZATION
This chapter presents the algorithm used to estimate the rotor position γ el . In this algorithm, the GNA [12] and the bisection method (BM) [16] are combined in order to minimize the residual of a nonlinear equation of the PMSM. The optimization criterion is to minimize the norm of the residuals min x∈R n with r = (r 1 , . . ., r m ) : R n → R m and the vector of the estimated parameters x. Starting with an initial guess x(0), the iterative process of the GNA searches for the best fit using the formula In the case of the position estimation, the estimated parameter is the estimated rotor position errorγ err and a scalar x =γ err (k). The residual of the FDE-QD is derived in Section III-B [cf.
(26)]. The variable k is used as an index to count the number of calls. Instead of using a variable step width, as described in (16), a fixed step width is used to optimize the estimated parameter. Only the sign of the adaption law is evaluated, leading tô The fixed step width η is analytically derived to represent the maximum value of the acceleration, as described in Section IV. Due to the mechanical boundary, the interval [γ el (k) + η,γ el (k) − η] defines the region in which the rotor position can be located in this call of the QD calculation. A search for the rotor position is performed within this interval, by an iterative process. During each iteration, the residual is updated using the currently estimated rotor position errorγ err and after each iteration the step size is halved for the next iteration, which forms the basis for the analogy to the BM.
An example function of the QD calculation with two iterations is shown in Fig. 3. Here, the rotor position of the PMSM is depicted in the dq-reference frame, which rotates in the stator-fixed αβ-reference frame. In addition, the rotor position of the next call of the QD calculation is shown in the d q -reference frame. The difference between the two positions is the step with the calling rate of the QD calculation f QD . Starting from the initial estimated rotor position of the previous optimization 0 , the search direction [cf. (19)] is evaluated for the first time, leading to a step toward the rotor position γ el , marked by 1 . After reevaluating the residual using the updated estimated parameter, another step toward the rotor position is taken. In this case, the search direction is inverted. By halving the step width, the estimation error is reduced, which results in 2 . Here, the iterative process is stopped by reaching a fixed number of iterations l = 2, as the termination criterion. After the iterative process has  been completed, the estimated position is fed-forward by a speed-dependent stepγ step which estimates γ step , leading to 0 '. The stepγ step is estimated by integration ofγ err with the time constant T QD,I , yieldinĝ The values of the initial step width η and the time constant T QD,I are set according to the analytical equations in Section IV. In Fig. 4, the structure of the QD calculation is shown. The estimated speed is calculated usinĝ

B. DERIVING THE RESIDUAL
The residual of the FDE is derived using a separate evaluation of the voltage differential equation in the different switching states of the converter. As injected signal the SWI of [11] is used. An example PWM period of the converter is shown in Fig. 5. In the upper plot, the voltage of the d-axis with a superimposed SWI is shown. In the bottom plot, the example current is depicted. For the passive switching state (PS), the differential equation of the PMSM results in The differential equation of the active switching state (AS) results in In the FPGA, the current derivative during the PS (dq) i PS is measured with a least mean squares (LMS) algorithm. The voltages (dq) v fm,PS and (dq) v fm,AS make up the fundamental components of the differential equation. The current derivative during the AS can be calculated using the results of two consecutive PSs.
To receive the equation of the residual, the FDE uses a consecutive PS and AS [4]. The assumption is made. This is valid as long as the switching frequency is higher than the time constants of the mechanical and electrical system. Substituting the voltage (dq) (26) The estimated inductance matrixL (dq) uses the estimated rotor position errorγ erȓ leading to an updated residual within each iteration of the QD calculation. In the next step, the Jacobian matrix is derived analytically, resulting in The adaption law of the parameter optimization giveŝ In (29), the denominator of the equation within the signum function is neglected, due to it always being positive, leading toγ err (i + 1) =γ err (i) − (30)

IV. ANALYTICAL DESCRIPTION OF THE QD CALCULATION
In this section, the analytical description of the QD calculation is used to derive its transfer function. The transfer function is then used to derive a parameter dependency of the initial step width η and time constant K QD,I of the QD calculation and so to achieve a stability criterion for a given maximum acceleration. The analytical description was first published in [17], but no experimental validation of the transfer function and stability criteria has been carried out.

A. TRANSFER FUNCTION OF THE QD CALCULATION
In the QD calculation, the rotor position is estimated directly. The speed is calculated via the integrally determined feedforwardγ step . Based on the difference equation of the previous section, the transfer function is derived next. For this purposê is inserted into (21) to substitute the change of the speeddependent fed-forwardγ step .
Further substituting for γ step with (20) andγ step with (22) yieldsω Looking at the equation above, the input value γ el in the numerator of the fraction can be seen in time-differentiated form. The rest of (33) corresponds to a first-order low-pass filter, and thus, the transfer function leads to The transfer function is validated on the test bench. For this purpose, an oscillating q-current setpoint with various frequencies is applied to the PMSM and the response of the estimator is evaluated. Fig. 6 shows the result for three different time constants K QD,I for the experimentally determined values and the analytical transfer function. Fig. 6 shows the Bode plot of the measured position (s) against the estimated speedˆ (s). Deviations only occur at high angular frequencies, confirming that the analytically derived transfer function is valid.

B. STABILITY LIMIT
To further validate the QD calculation, an experimental investigation of the stability criterion is presented in this section. Using the equations introduced in [17], stable operation for a given maximum acceleration a max of a drive train is guaranteed.
The initial step width η depends on the mechanical boundary. It takes the maximum acceleration of the drive train in the electrical reference frame into account. For the initial step width η, the position change between two samples of the estimator is used. This leads to γ er,max = a max f 2 QD and η ≥ γ er,max .
Using the tracking error of the low-passed speed, a relationship between the initial step width η and the time constant K QD,I is derived. Using a maximum acceleration and a linear increase in speed, the actual speed can be expressed by Next, the difference between the transfer functions of the real and estimated speeds obtained from the QD calculation is multiplied by the Heaviside step function. The Laplace function results in (38) This continuous-time function is transformed into a discretetime function using the bilinear transform (39) T QD = 1 /f QD is the time elapsed between two cycles of the estimator. After deriving the discrete-time function, the final value theorem is used to determine the difference between the estimated and actual electrical speeds, which results in A relationship between the initial step width and the time constant of the integrator is established by using leading to Due to the iterative solution, the maximal step width per call of the estimator results in On the test bench, l = 4 is implemented as a termination criterion. This results in a maximum step size of To investigate the stability criterion, the behavior of the QD calculation is compared using two different variants for the stability limit. One variant uses the maximal acceleration a full = a max of the drive train for parameterization, while the other variant uses a reduced acceleration a red = a max/4. To consider the maximum step size η max for the maximum representable acceleration the two variants a stab,full = a full η max η and a stab,red = a red η max η (45) are used. The parameters for the validation of the two variants [according to (36) and (42)] are shown in Table 2. For validation, the PMSM is current controlled and the measured position of the encoder is used for the FOC. As a sequence, a q-current i q is set until a speed of n = 1000 1 /min is reached. After reaching this speed, the q-current is reversed and maintained until the speed n = −1000 1 /min is achieved. Finally, a positive q-current is applied until standstill is reached. This experiment is performed for the current i q in the interval of [1, 2, . . ., 10]A. Fig. 7 shows the results for i q = 7 A and i q = 9 A. Both test scenarios are carried out with the full and reduced stability variants, which are denoted by the corresponding indexes.  For the experimental test with the current i q = 7 A, stability is achieved for both settings of the QD calculation. This is evident from the low estimation errors γ err,full and γ err,red . The speedsn full andn red are the same. At a current of i q = 9 A, stable operation can be observed for the variant with full stability. The estimation error of the variant with reduced stability γ err,red shows an increasing error with acceleration in the range with positive q-current and positive speed. For negative q-current and positive speed, the estimation error γ err,red increases until the rotor position γ red is estimated incorrectly. In this case, the frictional torque and machine torque act in the same direction and the estimation fails because the maximum representable acceleration a stab is exceeded.
For better illustration, Fig. 8 shows both the root mean square (rms) value and maximum value of the estimation error during constant acceleration. It is confirmed by both plots when a red,stab is exceeded, the estimation error increases and the estimation for the variant with reduced stability fails. The variant with full stability shows a low estimation error for all points considered, and thus, demonstrates stable behavior.

V. PARAMETERIZATION OF THE SPEED CONTROL LOOP
In this section, the open-loop speed control using an SSC with a QD calculation is derived and used to parameterize the speed controller with a specified amount of damping. For the design of the speed control loop, the inner current control loop is approximated by a first-order low-pass element. According to [18], the current control loop, which is parameterized using the amplitude optimum with a damping For the transfer function of the open-loop speed control loop, dead times are approximated by first-order low-pass elements.
Combining the transfer function of the discrete speed control, the current control, and the transfer function of the QD calculation results in a third-order low-pass filter with an integrator To simplify the formula, the three first-order low-pass elements are approximated by one first-order low-pass element. For this purpose, the accumulated time constant is used. The results of the approximation and numerical transfer function, using dead-time elements within the control loop, are verified in Fig. 9 for various K QD,I . In the Bode plot, the amplitude and phase agree in the range up to the cutoff frequency. Here, the dominant time constant is that of the QD calculation 1 /K QD,I , confirming that the approximation is valid. For a higher K QD,I , a numerical solution or a higher order approximation is advisable. For the analytical design of the speed controller, the transfer function is assumed. According to the design based on the symmetrical optimum, the parameters are calculated via The derivation can be found in [18]. In this work, b = 2 is used. Due to the design featuring the symmetrical optimum, setpoint filtering for the speed setpoint using a first-order lowpass filter with the time constant T n,I is necessary. In Fig. 10, the result of the controller parameters obtained with a factor b = 2 of the analytical equations are compared with the numerical solution. The numerical solution is based on a linearization at the operating point of the model which considers the dead times of the system. Instead of the integration time T n,I , the reciprocal K n,I = 1 /T n,I is shown.
Up to a time constant of the QD calculation K QD,I = 400 1 /s, a high correlation is achieved for the integral part K n,I . For the proportional part V n,P , a high correspondence is obtained across the whole range shown. Thus, the analytical approximation is assumed to be valid.
In the next step, the open-loop speed control is validated on the test bench. The machine is operated in open-loop control. In this configuration, the speed control setpoint can be used to manually set the speed to n ≈ 200 1 /min by having the speed control integrator assume a q-current that overcomes the test bench friction. By superimposing an oscillating setpoint at a fixed frequency, the phase shift and open-loop gain for one frequency can be determined. The results for different time constants K QD,I are presented as a Bode plot of the speed setpoint * (s) against the estimated speedˆ (s) in Fig. 11.
For each time constant, one numerical and two experimentally determined solutions are shown. The transfer function G exp is derived from the experimental investigation. At low frequencies, a deviation in the phase of G exp and G num is evident. This deviation is corrected in the transfer function G cor . In order to achieve the correction, the transfer function of the mechanical system is evaluated. The transfer function of the mechanical system G mech is shown in Fig. 12 and measured via the mounted encoder.
To determine the transfer function for the mechanical system, the input value is defined as the q-current I q (s), which is a torque-proportional value, and the output is the measured  speed (s). The correction of the phase for angular frequencies ω ≤ 3 rad /s is necessary since a system without friction, i.e., a pure integrator (cf. G mech,num ), is assumed by the numerical solution. In reality, the friction leads to a mechanical system with a first-order low-pass behavior. In the design, the mechanical system is still assumed to be an integrator, since the time constant of the friction G mech,fit with T mech = 3.3 s turns out to be very large and is sufficiently well represented by an integrator. Using the difference in phase of the numerical solution G mech,num and also the measured phase angles, the phase of the experimental results of Fig. 11 can be corrected, yielding the transfer function G cor . For the correction of the phase is used. Thus, up to an angular frequency of ω < 6 rad /s, the results of the experimental solution G cor achieve a high correlation with the numerical solution G num . Above the angular frequency of ω ≥ 6 rad /s, the deviation between the experimental solution and the numerical solution increases. Here, the numerical transfer function lies between the solutions G exp and G cor . Looking at the transfer function shown in Fig. 12, another phase shift occurs from an angular frequency

VI. EXPERIMENTAL RESULTS WITH THE FDE-QD
In order to validate the system in closed loop speed control, the dynamic response of the system is considered. In Table 3, the parameters of the controllers and estimator are shown. Two scenarios are considered to evaluate the dynamic behavior of the speed control with SSC. In both scenarios, the PMSM is speed controlled, using only the estimated position and speed for control. The mounted encoder is solely used to determine the quality of the SSC. The first scenario tests the disturbance response at zero speed. The PMSM is speed controlled with a setpoint of n = 0 1 /min. The load machine applies positive rated torque, then reverses the torque and finally the torque is set to zero. This test has been proven in the past to be a reliable indicator of the dynamics of the SSC, as the limitation of the jolt is only limited by the dynamics of the current control of the load machine. Since there is no position control of the device under test, the load torque is compensated solely by the speed control, using the estimated speed as feedback. Accordingly, after the load torque has been applied, there is a change in position until the speed control has reduced the actual speed back to the setpoint, i.e., standstill. The more dynamic the speed control is configured, the smaller the change in position during the jolt. The result is shown in Fig. 13.
The results show that the electrical rotor position changes by approximately 120 • when the rated torque is applied, while standstill is already reached after approximately 40 ms. On the basis of the speed graph, the damping D < 1 can be seen, as there is a slight overshoot. When changing from positive to negative rated torque, the rotor position changes by approximately 260 • . The value is more than doubled. The cause lies in the current limitation for the q-current of 15 A, which corresponds to the maximum current of the machine, and thus, creates a nonlinearity. When the load torque is switched OFF, again a position change of approximately 120 • occurs. During  the entire test sequence, the rotor position estimation error remains less than 5 • . Thus, the stability of the rotor position estimation is proven. If the estimated speed is compared with the measured speed, a slight phase shift can be seen. Overall, the deviations remain small.
The control response is shown in Fig. 14. Here, the speed setpoint is set to n = 200 1 /min (≈ 5% of the rated speed), reversed to n = −200 1 /min, and set back to n = 200 1 /min, while rated torque is applied by the load machine. The damping of the control behavior can be determined by the transient response. In addition, the applied load torque tests both the generator and motor operation, as the speed is reversed. In combination with the test according to Fig. 13, the important operating points of standstill, generative and motor operation of the PMSM in the range of low speeds are covered.
As expected, the dynamic response shows a slight overshoot of the speed, while the q-current of the PMSM is set dynamically. The low rotor position estimation error indicates that the rotor position estimation is stable. Again, a small phase shift between the estimated and measured speed becomes apparent. All in all, the expected control behavior could be validated.

VII. CONCLUSION
In this article, an analytical description of the transfer function of the QD calculation is presented. A stable operation can be achieved for a given maximal acceleration of the drive train by deriving an interdependence between the two parameters of the QD calculation. The remaining degree of freedom is used to tune the dynamics of the estimator. By using the transfer function of the QD calculation, the speed controller for SSC can be parameterized analytically. The resulting phase margin of the open-loop speed control is validated experimentally. Based on these results, the commissioning of the SSC can be performed purely analytically. While this article has shown the application of QD calculation to the specific FDE residual, which has been designed for the low-speed operation of PMSM, the results can be applied to other estimators based on the QD calculation method as well.