Nonlinear Modelling, Flatness-Based Current Control, and Torque Ripple Compensation for Interior Permanent Magnet Synchronous Machines

: A nonlinear mathematical model for the dynamics of permanent magnet synchronous machines with interior magnets is discussed. The model of the current dynamics captures saturation and dependency on the rotor angle. Based on the model, a ﬂatness-based ﬁeld-oriented closed-loop controller and a feed-forward compensation of torque ripples are derived. Effectiveness and robustness of the proposed algorithms are demonstrated by simulation results.


Introduction
As modern three-phase current-fed machines are increasingly used for industrial applications, efficient and accurate current control of such machines is of great economic interest. This is due to the fact that the current controller influences the operational behaviour of the machine to a very large extent, namely the energy efficiency and the occurrence of undesirable side-effects, e.g., torque ripples. The aforementioned electric machines comprise induction machines (IMs) and permanent magnet synchronous machines (PMSMs), where the latter are more energy-efficient and have a higher power density [1].
The present work considers a PMSM with interior magnets (IPMSM). Such IPMSMs are particularly advantageous in terms of their achievable torque density [2]. Due to increasingly involved designs of the rotor [3,4] and the stator [5], and due to higher power density [6], effects such as magnetic saturation of the iron material, cross-coupling [7], and the dependency of the system dynamics on the rotor position become more important for high-precision control design [8]. Typically, the influence of the rotor position is neglected. However, some modelling approaches exist in the literature that emphasise the importance of considering the rotor dependence in the machine model, see e.g., [8][9][10][11][12][13].
For PMSMs, the consideration of the rotor position in modelling and control design is useful for compensating ripples on the torque and the current, see e.g., [8]. In [9] and [10], stator-fixed PMSM flatness-based control schemes are proposed, which rely on a graph-theoretic approach to magnetic equivalent circuit PMSM modelling from [14]. The proposed method to derive the mathematical model is rather involved, which might contradict industrial practice. Moreover, the dependency of the model on the rotor angle is greatly simplified for control design.
In contrast to control schemes based on a stator-fixed model, other flatness-based approaches for PMSMs exist in the literature that build upon the field-oriented control method-ology. That is, the machine is described in rotor-fixed coordinates, see e.g., [15][16][17][18][19][20][21]. However, usually only fundamental wave models are considered that impede the handling of physical interdependencies of the rotor and the stator field, affecting both the current dynamics and the electromagnetic torque. In saturated electrical machines, those interactions manifest themselves as torque and current ripples.
In [8,22], the dynamics of a PMSM, including the dependence on the rotor angle, is modelled based on a rigorous investigation of the physical phenomena. The model is derived in rotor-fixed coordinates, and a feed-forward control scheme based on look-up tables is proposed to compensate ripples on the current and the torque. In [23] torque ripples of a flux-switching permanent magnet machine are reduced by means of a feedforward control scheme, where polynomial approximations for the Fourier coefficients of an inverse torque model are derived.
The contribution of the present paper is twofold. On the one hand, a nonlinear mathematical model of the dynamics of an IPMSM is presented (Section 2). The model includes saturation, cross-coupling, and the dependency on the rotor angle. In industrial practice, the model can easily be transferred to various types of machines without much additional development effort. Moreover, the frequently used rotor-fixed coordinates are preserved. Following the notion of imperfect dynamical systems of [24], the angledependent proportion of the dynamics may be considered an imperfection and, thus, a conceptually straightforward extension of the well-known fundamental wave model.
On the other hand, the proposed machine model is exploited in order to derive a flatness-based field-oriented control (FOC) scheme (Section 3.1) as well as a feed-forward compensation of torque ripples (Section 3.2). That is, the aforementioned imperfection is systematically compensated. To keep the embedded control algorithms computationally efficient and to reduce the required memory, they rely on a simplified version of the model from Section 2. Here, physical considerations are exploited to achieve an efficient reduction at reasonable losses in accuracy, and the parameters of the reduced model are calibrated by using existing data from finite element method (FEM) simulations. The proposed scheme for current control and compensation of torque ripples explicitly decouples the dq-currents. This significantly simplifies tuning of controller gains and yields an accurate tracking of even highly dynamic reference currents and a good robustness against model uncertainties, which is demonstrated in Section 4.

Voltage Equations
In the following, a three-phase PMSM with star-connected stator winding is considered. The three phase voltages u abc = (u a , u b , u c ) T result from ohmic losses over the stator windings (resistance R s ) and the change in the corresponding flux linkages ψ abc = (ψ a , ψ b , ψ c ) T according to Faraday's law of induction where i abc = (i a , i b , i c ) T are the phase currents. The flux linkages are generated by the permanent magnets and the coil currents. Note that in the abc-frame the current dynamics depend naturally on the electrical rotor angle θ, which complicates the design of control algorithms, see e.g., [25]. In order to simplify the latter, the system variables are transformed into a rotor-fixed dq-coordinate frame aligned with the magnetic flux of the permanent magnet. To this end, the Park transformation (·) dqn = T(θ)(·) abc (2a) with is used, see e.g., [25]. Applying this transformation to the phase voltages, Equation (1) yields the well-known voltage equations in the dq-frame with the zero component where ω = dθ/dt. In the fault-free case, the neutral point of a PMSM with symmetric star-connected stator winding is isolated such that u n = 0 [26]. Hence, Equation (4) is omitted in the following for the sake of brevity. When saturation and the dependency on the rotor angle are neglected, the well-known linear fundamental wave model is obtained, and the voltage Equations (3) simplify to with the constant inductances L d and L q as well as the permanent magnet flux ψ m . In practice, the inductances depend on both currents i d and i q and the electrical rotor angle θ. The dependency on the currents is a result of the saturation of the iron material, which yields a varying slope of the ψ-i-characteristics over the current. Additionally, the d-current affects the q-component of the flux linkage and vice versa due to common flux paths in the stator yoke, which is known as cross-coupling, see [7].
Besides the saturation and cross-coupling of the magnetic flux linkage, the rotor angle plays a significant role that should be taken into account. In order to obtain the constant inductances L d , L q and ψ m as in the simplified model in Equation (5), the flux linkages in the abc-frame must be purely sinusoidal over the rotor angle. However, in practice, several physical effects lead to deviations from pure sinusoidal distributions of the flux linkages over the rotor angle. On the one hand, the air gap field excited by the stator has a staircase shape due to concentrated fluxes in the slots [26]. On the other hand, the rotor shape and the position of the magnets often yield a significant variation of the magnetic flux distribution from an ideal sinusoidal shape caused by the arrangement of the magnets and the resulting magnetic reluctance, see e.g., [3]. It should be emphasised that due to these reasons, even in the dq-frame, the dependency of the flux linkages on the rotor angle is significantly pronounced.
Altogether, the time derivatives of the magnetic fluxes in Equation (3) can be expanded in order to achieve with the abbreviations

Torque Characteristics
As the current dynamics, the torque generated by the motor depends on the rotor angle, too. This dependency yields undesired fluctuations in the motor torque, which can be classified as follows [27]:

1.
Torque ripple due to the interaction of the stator field distribution with the electromagnetic rotor properties, which is excited by (a) a non-sinusoidal rotor field distribution and (b) a varying reluctance distribution over the circumference of the rotor.

2.
Cogging torque due to the interaction of the rotor field distribution with the circumferential variation of the stator reluctance.

3.
Pulsating torque: sum of the torque ripple and the cogging torque.
A frequently used approach to derive the torque characteristics of an electric machine is to calculate the magnetic coenergy W * m , see e.g., [28]. Then, the torque T e is achieved as where p is the number of pole pairs. Alternatively, the power balance can be exploited as follows [28]. The electrical power P e supplied to the motor reads where u d and u q are given by the voltage Equation (6). Since the mechanical power output P m can be described as with the rotor speed Ω = ω/p, the electrical torque T e is calculated by splitting the power balance up into copper losses, change in magnetic energy, and mechanical power, see e.g., [25]. Doing so, is obtained, see also [28] or [29]. If the dependency of the rotor angle is neglected, the second summand on the righthand side of the torque expression in Equation (10) is omitted and the well-known torque equation is obtained, representing a constant torque when the currents i d and i q are constant [25].
Notice that Equation (10) yields a vanishing torque for vanishing currents and, thus, is insufficient to explain the existence of the cogging torque. Hence, Equation (10) is extended analogous to [13] as follows: where T cog is the cogging torque.

Simplified Model for Embedded Control
In order to exploit the current dynamics in Equation (6) for model-based control, partial derivatives of the flux linkages must be known. To this end, either the characteristics of the flux linkages or the magnetic coenergy can be stored on the embedded control hardware. Then, the required partial derivatives of the flux linkages are obtained by (repeated) partial differentiation, see e.g., [10,13,30]. Doing so, the achievable control accuracy is the better the more accurate the embedded numerical model of the flux linkages or coenergy is.
In principle, the magnet characteristics can be represented on the embedded hardware by data-based surrogate models such as look-up tables (in combination with a suitable interpolation algorithm), artificial neural networks, or Gaussian processes. However, such universal approximators may be inefficient in terms of required memory or computational effort since they do not exploit physical interdependencies.
In the following, an alternative approach for representing the relevant characteristics of a PMSM is described. Initially, the torque Equation (11) is approximated aŝ where the mean flux linkages over one electrical revolution replace ψ d and ψ q , respectively. Note thatψ d andψ q are independent of the rotor angle. The partial derivatives ψ dθ and ψ qθ and the cogging torque T cog are replaced by the following particular truncated Fourier series: These expressions are justified as follows: 1.
Due to the interaction of the nonsinusoidal stator field distribution with both that of the rotor as well as that of the reluctance, only integer multiples of the 6 th fundamental stator frequency f s are present in ψ dθ and ψ qθ . This is due to the interaction of the three phases, see e.g., [11,12,27,31]. For the motor studied in this research, it has turned out that the 18th harmonic and frequencies greater than the 24th harmonic almost vanish. Additionally, linear regression models for the phase shifts are sufficiently accurate for the particular type of motor.

2.
The fundamental frequency of the cogging torque can be derived from the number of pole pairs of the rotor and the stator slots, as presented e.g., in [32]. For the present machine, 36 f s is achieved. It has turned out that subharmonics also need to be considered. Here, constant approximations of the phase shifts are sufficient.
Note that no coefficient of the problem-tailored surrogate models of Equation (13) depends on more than two variables. This allows efficient implementation of these models on an embedded hardware by simple look-up tables.
In order to approximate the current dynamics of Equation (6) as well, surrogate models for the partial derivatives ψ dd , ψ dq , ψ qd , and ψ qq of the flux linkages are required in addition to Equation (13). To this end, averaged flux linkages in Equation (13a) are numerically differentiated on the embedded hardware.

Model Calibration
In order to calibrate the torque model in Equation (13), training datasets can be generated either based on laboratory measurements or on a numerical simulation model such as a finite element model. The first approach would require a physical prototype of the motor and an accurate test bench for dynamic torque measurement. In contrast to this, the second approach can be pursued in an earlier project phase, when no physical prototype of the motor is available. Since the utilisation of high-fidelity FEM computations has become common practice in industry, we propose to simply reuse the available simulation models as they are for data generation. This also avoids the necessity for involved torque measurement equipment.
For the present machine with 2p = 6 poles, representative normalised magnetic flux linkages and their normalised partial derivatives, respectively, are shown in Figures 1 and 2, respectively. It can be seen that the magnetic flux linkage and its partial derivatives vary over the electrical rotor angle θ. This variation is mostly pronounced for small currents and gets smaller with increasing magnetic saturation, see Figure 1. Moreover, the saturation influences the q-direction more than the d-direction, see Figure 2.  Given the training data on a sufficiently fine grid, the standard particle swarm optimisation algorithm [33] has been used for least-squares fitting of the model parameters.

Flatness-Based Current Control
Given a nonlinear system where u is the input and x is the state, whether the system is differentially flat or not is of great interest from a control perspective. This is due to the fact that a differentially flat system is controllable and, moreover, a so-called flat output y ∈ R m can be found such that the design of multi-variable decoupling feed-forward and feedback control laws for y becomes particular simple [34].
Considering the nonlinear machine model in Equation (6), that may be written in the form of Equation (14) with the input u = (u d , u q ) T , the state x = (i d , i q ) T , and the output y = (i d , i q ) T , it can easily be seen that for the current dynamics in Equation (6), y is a flat output. In order to design a feedback law, the new input v = (v d , v q ) T with v d = di d /dt, v q = di q /dt is introduced, which allows us to rewrite the current dynamics given in Equation (6) as Choosing the components of v as PI control laws and where t → y * j (t) denotes the (at least once differentiable) reference trajectory for y j , the decoupled closed-loop dynamics of the tracking errors e j = y * j − y j is achieved, which is asymptotically stable for K Pj , K I j > 0. Note that the closed-loop error dynamics in the dq-coordinates are linear and decoupled, which simplifies the tuning of the controller gains. The actual controller gains can, e.g., be chosen by pole placement.

Torque Control and Compensation of Torque Ripples
Given a desired torque T * e , various strategies exist in the literature to choose advantageous desired combinations i * d and i * q of the dq-currents, such as the so-called maximum torque per ampere (MTPA) algorithm [35]. However, these algorithms typically do not take into account the dependency of the torque characteristics on the rotor angle. This results in torque ripples, which excite noise and vibrations during operation of the motor.
Besides a ripple-reducing design of the motor [32], torque ripples can be compensated by control algorithms. Here, it is common practice in the literature to inject a ripplecompensating current in the q-direction, see e.g., [23,36,37].
Such an algorithm based on the efficient torque model from Section 2.3 is described in the following. That is, given desired values i * d and i * q for the dq-currents, which have been computed without consideration of the torque ripples, an injection current i qc is sought such that when instead of i * q is used as the reference current for the underlying current control law, torque ripples are reduced.
Based on Equation (11), the desired torque T * e is computed as whereψ d andψ q are given by Equation (13a). By comparison of the desired torque of Equation (19) with the approximated actual torque of Equation (11), the injection current i qc is achieved as the solution of the implicit nonlinear equation Note that no analytical solution of Equation (20) for i qc can be found. Hence, an iterative solution is performed in each real-time step, where the bisection method is proposed in the present work due to its robustness and simplicity. To keep the required number of iterations small and, thus, to reduce the computation time, an initial guess i 0 qc is computed in each real-time step as whereL dq andL qq are the arithmetic means ofψ dq andψ qq over the whole range of operation of the motor. Given the initial guess i 0 qc from Equation (21), the search interval of the bisection method is reduced to i 0 qc − δ/2, i 0 qc + δ/2 , where the interval width δ is chosen to capture relevant errors of the initial guess. Starting with i 0 qc , the approximate solution i N qc of i qc is achieved after N iterations of the bisection method.
For the flatness-based control law in Equation (16), not only the desired currents i * d and i * * q are required, but also their time derivatives di * d /dt and di * * q /dt = di * q /dt + di qc /dt. If the overlying algorithm, which generates i * d and i * q , does not also provide di * d /dt and di * q /dt, it is often reasonable to assume di * d /dt ≈ di * q /dt ≈ 0 or to approximate di * d /dt and di * q /dt numerically from i * d and i * q . In contrast to this, the injection current i qc is relatively quickly fluctuating. Hence, greater care should be taken on the computation of di qc /dt to achieve a sufficient bandwidth of the ripple compensation. In the present work, a simple linear time-invariant differentiating filter with the transfer function s/(Ts + 1), where s is the differentiation operator and T is the time constant of the filter, is used to approximate di qc /dt from i N qc . If a higher accuracy is required, di qc /dt can be computed by implicit differentiation of Equation (20) at the cost of increased computational complexity.
The proposed control scheme, consisting of the flatness-based field-oriented current controller from Section 3.1 and the feed-forward compensation of torque ripples, is illustrated in Figure 4.

Simulation Results
For numerical validation, a co-simulation of the proposed algorithms with a detailed model of the physical system has been implemented in MATLAB/SIMULINK. The physical system is simulated with the ode45 variable-step solver using a relative tolerance of 10 −8 . The control algorithms have been discretised using the trapezoidal integration rule and run at a fixed sampling period of τ el /20, where the electrical time constant τ el is defined as Note that for the present motor, the selected time constant based on the d-direction is faster than its counterpart for the q-direction. The controller gains in Equation (16) To demonstrate the robustness of the control algorithms, 10 % relative parameter errors on both R s and the magnetic flux characteristics as well as a bias of 0.1 • on the measured rotor angle are simulated. Figure 5 demonstrates that the flatness-based feedback controller from Section 3.1 enables a very accurate tracking of large simultaneous set-point changes on both i d and i q within a time significantly shorter than the time constant defined by Equation (22). Furthermore, the controller shows a good robustness against the simulated parameter errors and the measurement error on the rotor angle. Figure 6 shows the effectivity of the torque ripple compensation from Section 3.2 over the whole range of achievable torques for the given motor. It can be seen that already by the initial guess of Equation (21) for the required compensation current i qc , an average relative reduction of the torque ripples by approximately 73% is achieved. By iterative refinement of the initial guess based on Equation (20), the average relative ripple reduction increases to 84%. For operating points where the torque ripples are particularly pronounced, even higher relative reductions are reached. This can be seen in more detail in Figure  7. The remaining relatively small ripple after the proposed compensation is due to the simplifications of Equation (13) in comparison to the high-fidelity FEM-based model. T e (i * d , i * q , θ) T e (i * d , i * q + i 0 qc , θ) T e (i * d , i * q + i N qc , θ) Figure 6. Averaged torque ripple relative to the respective mean torque without compensation, with compensation by i 0 qc and by i N qc .

Conclusions
Mathematical models for the current dynamics and the torque of PMSMs have been introduced. The proposed models capture saturation, cross-coupling, and dependency on the rotor angle. Based on the models, a control scheme containing a flatness-based fieldoriented closed-loop current controller and a feed-forward compensation of torque ripples has been derived. In order to improve the computational efficiency and to reduce required memory on embedded hardware, the motor characteristics used in the control algorithms have been approximated by truncated Fourier series. These application-specific surrogate models have been calibrated using data from FEM simulations, which are typically present in an industrial design process for electric motors. High accuracy of the current control and effective compensation of torque ripples has been simulatively demonstrated even under significant model uncertainties and a bias on the measurement of the rotor angle.

Conflicts of Interest:
The authors declare no conflicts of interest.