An Accurate Wide-Speed Range Control Method of IPMSM Considering Resistive Voltage Drop and Magnetic Saturation

This paper deals with the high accurate current set-points solution for interior permanent-magnet synchronous motors (IPMSM) in wide-speed range applications. Considering voltage and current constraints, the operating regions can be divided into maximum torque per ampere, maximum current, field weakening, and maximum torque per voltage regions, which requires solving different nonlinear functions in real time to obtain optimal current set-points. Traditional methods including curve-fitting methods and polynomial approximation (PA) methods are not easy to obtain these solutions, especially involving magnetic saturation problems. In this paper, Newton–Raphson algorithm for improving the control accuracy of the current set-points is proposed. Meanwhile, parameters influence including magnetic saturation and resistive voltage drop is fully investigated. Compared with PA method, the proposed method is able to converge to accurate solutions in few numbers of iterations with reduced execution time, which can be easily implemented on an off-the-shelf digital signal processor. Both simulation results and experimental results on an 8-kW IPMSM rig are conducted showing good agreement with the expected results.


I. INTRODUCTION
nterior Permanent Magnet Synchronous Motors (IPMSM) are widely used in industrial applications such as in hybrid and electric vehicles, thanks to their high efficiency, wide speed range and power density capabilities. In some applications such as electric vehicles, where a wider speed range is required, the operating regions of an IPMSM needs to be extended from constant-torque region to constant-power region [1] [2] . In general, the operating conditions can be summarized in four regions: Maximum Torque Per Ampere (MTPA), Maximum Current (MC), Field Weakening (FW), and Maximum Torque Per Voltage (MTPV).
When the motor is operating under the base speed, the maximum torque per ampere (MTPA), which aims at minimising the copper loss, becomes more attractive. The works proposing MTPA control strategies can be mainly divided into three categories: 1) Mathematical-model-based MTPA calculations [3] [4] ; 2) Signal injection-based MTPA point tracking [5]- [8] ; 3) Searching-method-based MTPA [9] [10] . First methods including curve-fitting [3] and polynomial approximation [4] , are proposed to solve nonlinear MTPA formula directly. However, inductance is regarded as linear variation with current and cross-saturation effect is ignored. Second methods utilized high frequency current injections [5] [6] or voltage injections [8] , MTPA points can be detected either from the torque or the speed response. However, the injected current signals may result in torque ripple problems. Third methods, aiming at minimising copper losses, present gradient method in [9] [10]. However, most of the time, the variation of inductances with respect to the current is neglected for computation and implementation simplicity.
When the speed of IPMSM goes up above base speed in MC region and FW region, the flux-weakening methods are needed, presented in [11]- [16]. Defined by the way to acquire more negative d-axis current, two approaches for implementing the flux-weakening control are presented as feedforward [11]- [14] and feedback [15]- [17] methods. In the feedforward approaches [11]- [14], the intersections (solutions) involving current circle, voltage ellipse and torque hyperbola requires solving complex nonlinear equations online. To simplify the solving process, resistive voltage drop and magnetic saturation are often neglected [3] [11], which loses the accuracy of the current setpoints. The resistive voltage drop, which causes the voltage ellipse shifting, can be compensated from control part in [11] or from calculation part using "Chord method" proposed in [12]- [14]. In the feedback approaches [15]- [17] , the negative d-axis current is automatically adjusted by tracking the voltage An Accurate Wide-speed Range Control Method of IPMSM Considering Resistive Voltage Drop and Magnetic Saturation constraint while the speed is increasing. These methods are considered robust because they do not need the prior knowledge of motor parameters, but this is at the expense of reduced performance during transients [17] . When the motor is operating in deep flux-weakening region, so-called MTPV region mentioned in [18]- [20], the MTPV is the intersection of the voltage ellipses with the highest torque hyperbola at saturated stator voltage. The common method proposed in [19] and [20] are mathematical model-based MTPV calculations. The same problems are parametric sensitive. All aforementioned open loop methods have an advantage on high dynamic response. However, in most of cases, the variation of inductances with respect to the current is neglected for computation simplicity. In fact, there are always some parameters uncertainties including inductance variations and permanent magnetic linkage variations [21]- [26] . Two common solutions including PA methods and parameter identification methods are investigated in related works. The former one applies PA methods for curve-fitting, such as first-order PA method in [21] [22]. The inductance is considered linear variation with current. However, cross-saturation is ignored. The second-order PA methods, to achieve inductance approximation, are introduced to obtain MTPA and FW curves, which makes the curves a fourth order complicate nonlinear equations, discussed in [23]- [25]. The latter one such as PSO (Particle Swarm Optimization) are presented in [27]- [29], the identification results are compared with experimental results obtained by current decay test. The experimental accuracy is also key issue for practical use.
Nevertheless, to the best knowledge of the authors, a complete theory, which (a) covers all operation regions (such as MTPA, MC, FW, MTPV) ,(b) allows for an analytical computation of optimal currents set-points, (c) considers stator resistance and variable inductance , is not yet available.
In this paper, a novel and accurate feedforward Newton-Raphson searching method in all operating regions is proposed, while taking into account both the resistive voltage drop and the magnetic saturation. The proposed method is characterized by the following:(a) Analysis of all the operating regions (MTPA, MC, FW, MTPV); (b) Implicit problem formulation (the nonlinear solutions involving current circle, voltage ellipse and torque hyperbola); (c) Real-time implementation on a digital microcontroller with few iterations and reduced computation burdens. Both simulation and experimental results under different conditions prove the effectiveness and rapidity of the proposed methods.
The paper is organized as follows: in Section II, the different operation modes for IPMSM are described. In Section III, the proposed N-R searching method is elaborated. The influence of resistive voltage drop and magnetic saturation is explained in Section IV. Simulation and experimental results are presented in Section V and Section VI respectively. Conclusions are given in Section VII.

II. ANALYSIS OF OPERATING MODES FOR IPMSM
In this Section, the operating loci introduced in the introduction are extended and described by mean of analytical equations and figures. Within the rotating reference frame, the equations for an IPMSM can be written as follows: Where, ud and uq are the d-q-axis stator voltages; id and iq are dq-axis current; Ld and Lq are the dq-axis inductance, Rs, ψf , p are the stator resistance, the permanent-magnet flux linkage and the pole pairs, respectively. ωe is the electrical angular frequency.
There are also two additional constraints to be taken into account. The first one is related to the current constraints of the motor or the inverter (3). The second one is related to the voltage constraints (4).
where Umax (Udc/√3) is the maximum voltage for SVPWM. As can be seen in (3), the current constraints trajectory in the dq-axis forms a circle (curve AH in Fig.1), whose center is the origin and the radius is Imax,pk. Substituting (1) and (2) into (3), and neglecting the resistance drop, (4) can be rewritten as the following: Equation (5) describes the iso-voltage constraints ellipse shrunk with the increasing speed in the dq-axis current plane, whose center Q is (-ψf/Ld, 0). The electromagnetic torque is given by the following equation: The effective operating regions are shown in Fig.1, which consists of MTPA region (curve OA), MC region (curve AH), FW region (area OAHQ) and MTPV region (curve HQ). Defining superscript * as any generic set-point, * is the torque set-point, which at steady state is equal to the load torque L . * is the speed set-point, n is the measured speed, b is the base speed, max is the maximum speed defined by the mechanic components (i.e. shaft, bearings, etc.), and n is the boundary speed defined by the iso-voltage constraint ellipse passing through the origin. n is the deep-flux-weakening speed defined by iso-voltage constraint ellipse passing through the intersection of MTPV locus and current circle (point H in Fig.  2(g)). Defining also L0 as the torque hyperbola at no load, D as the cut-off torque identified by the intersection of the MTPA locus with the iso-voltage ellipse at speed set-points (point D in Fig.1(c) and Fig.1(d)), I as the cut-off torque identified by the intersection of the MTPV locus with the iso-voltage ellipse at speed set-points (point I in Fig. 1(g) and Fig.1(h)), and MC as the maximum torque for a measured speed in MC region. The operating mode selector provides the right mode taken into account together with the proposed Newton-Raphson searching method explained in Section. III, which presents four main branches selecting one mode among the following ones: MTPA, MC, FW and MTPV. The relationship for speed and the torque is explained in the following paragraphs with aid of Fig.  1 and Fig. 2. In the following paragraphs, the four branches are described as follows: In this branch, the speed is lower than the base speed b and the IPMSM is operating on the MTPA curve OA, shown in Fig.1(a). The red star is the steady operating point C and it is defined by the intersection of the MTPA and the load torque curves. The MTPA solution is formulated as (7) aiming at copper loss minimization problem: The Lagrange multiplier is applied to obtain the solution: The partial derivative of (8) is obtained as (9): Solving the first two equations in (9) by eliminating the parameter λ , the MTPA relationship for dq-axis currents can be expressed as (10): Branch II: ( < ≤ ) If the speed is higher than base speed b and less than the boundary speed o , voltage constraint ellipse shrinks to point A. More negative d-axis current is needed to make the dq-axis current set-points move along MC circle along curve AB shown in Fig.1(b), limited by both voltage constraint and current constraint, which can be expressed in (11).
can be calculated at the same time in (11): When the measured speed is approaching the speed setpoint * and if the condition < * ＜ is satisfied, the dqaxis current set-point will move along the voltage ellipse curve BD towards point G in FW region shown in Fig.1(c). Point D is the voltage constraint ellipse at speed set-point * intersected with MTPA curve marked in Fig.1(c) and Fig.1(d).
can be calculated in (12) and the first equation in (12) comes from (10): In this case, the dq-axis current set-points are limited by both torque hyperbola and voltage constraint expressed in (13): If the condition T e * < T D is satisfied, the point F, the intersection between voltage constraints ellipse and the load torque curve, is located outside the top left quadrant, the final operating point can't go across the boundary MTPA curve OA and finally it settles at point E, shown in Fig.1(d). Curve DE is part of MTPA curve can be expressed as (10) as well.
Branch III: ( < ≤ ) When the measured speed is higher than and if the measured speed is approaching the set-point speed * , The dq-axis current set-point will move along the MC circle shown in Fig.1(e). In this case, the dq-axis current set-point is limited by both current constraint and voltage constraint can be expressed in (11). If the condition T * ≤ T MC is satisfied, the dqaxis current set-point will move along the voltage ellipse at speed * curve BD and balance at point G shown in Fig.1(f), as can be expressed as (13) as well.
Branch IV: ( ≤ ＜ ) When the measured speed is higher than (iso-voltage constraint ellipse passing through the intersection of MTPV locus and current circle), the motor is operating in MTPV region (deep flux-weakening region). If the measured speed is lower than the set-point speed * , The dq-axis current setpoint will move along the MTPV curve shown in Fig.1(g). MTPV solution is formulated as (14) aiming at output power minimization: Solving (14), the MTPV relationship for dq-axis currents can be expressed as first equation in (15). In this MTPV mode, the dq-axis current set-point is limited by both MTPV curve and voltage constraint: If the condition T * ＜T I is satisfied, the dq-axis current setpoint will move along the voltage ellipse at speed * (curve IK) and balance with the load torque (point K) shown in Fig.2(h), as can be expressed as (13) as well.
is cut-off torque, identified by the intersection of the MTPV locus with the iso-voltage ellipse at speed set-points. can be calculated in (16). The first equation of (16) comes from MTPV solution (15): The challenge in solving (10)-(16) is discussed below: 1) The solutions of (10)-(16) are nonlinear with a square root operation, which is quite complex to solve by real-time controller directly. In addition, the equations (11)-(16) contain electrical angular velocity that requires real-time measurements, which increases the difficulty of solving equations.
2) The parameters in the equations cannot remain constant. When the magnetic saturation and cross saturation is significant, the inductance Ld and Lq will decrease with the increase of the dq-axis currents. The voltage constraints ignore the effect of resistance, which also affects the control accuracy, which will be further discussed in Section IV.
3) The solutions of these equations can be listed in advance in a pre-made table; however, considering inductance variation and resistance variation, these tables tend to be huge because it is necessary to create many separate tables.

III. THE PROPOSED NEWTON-RAPHSON SEARCHING METHOD FOR DIFFERENT OPERATING MODES
The second-order N-R method is a fast convergence procedure to solve nonlinear equations f(x,y) = 0 and g(x,y)=0, in our particular application, where x is * and y is * : Equation (17) is our target function to find a value of x,y. Assuming that function (17) is continuous and there exists a continuous second-order partial derivative in the neighborhood at the point (x0, y0), (17) can be expanded into Taylor series in (18) and equation (19) is the Jacobian matrix of (17): , , Substituting ∆ = +1 − , ∆ = +1 − , 0 = and 0 = into (18), which the solutions in (17) can be expressed as (20) in iterative form: If the condition (21) is satisfied When the inverse Jacobian matrix is substituted into (20), a second order N-R iterative form can be listed as (22):   1   1   , , , ) x g x y f x y f x y g x y g x y f x y f x y g x y y y g x y f x y f x y g x y Where xk+1 and yk+1 ( +1 * and +1 * ) are the new estimate of x and y, xk and yk ( * and * ) are the previous estimate of x and y ( * and * ) respectively, f(xk, yk) and g(xk, yk) ( ( , * * ) and ( , * * )) are the target functions using xk and yk ( * and * ) ( in most cases, f(xk, yk) ≠ 0 and g(xk, yk) ≠ 0 because xk is not the correct solution). k is an integer iteration index that starts with 1. The iterative procedure starts by substituting a first guess x(0) and y(0) into (22) to get a second estimate. This second estimate is then substituted into (22) to get a third estimate. This process is repeated until the geometric distance of the two iterations results are very small (E is defined as current setting precision discussed in Section V): The N-R algorithm flow chart is shown in Fig.3. The initial value for the first guess is set at first. N represents the maximum number of iterations. The Jacobian matrix and the initial target function values f(id0 ,iq0) and g(id0 ,iq0) are used for the N-R iterative algorithm. Then id1 and iq1 are update using (22). The iterative process are repeated until the error reaches a small value (i.e. minor than 0.01 2 ), the current set-points rapidly converge to the optimal values. If the maximum number of iterations N is exceeded, equation (22) is not converging and the searching algorithm has to be re-initialized.
According to four different branches shown in Fig.3, the N-R searching method uses different target functions and Jacobian matrices in Table I to obtain the optimal current set-points. The proposed method is capable of considering the magnetic saturation and the resistive voltage drop as it will be later shown in the Section IV.

IV. THE INFLUENCE OF THE RESISTIVE VOLTAGE DROP AND MAGNETIC SATURATION A. Magnetic Saturation and Cross Saturation
In most torque control applications, accurate inductance information is required for torque estimation and for an optimal current set-points selection. However, due to the magnetic saturation, the inductances vary nonlinearly depending on the current. In Fig.4, five measured dq-axis inductance values acquired by static inductance experiments discussed in [30] [31] are compared with FEA results.

and k<N
Select and calculate Jacobi matrixes and target function in Table I

TABLE I TARGET FUNCTIONS AND JACOBI MATRICES FOR DIFFERENT OPERATING MODES (DQ-AXIS INDUCTANCES ARE COMING FROM ID AND IQ LUTS)
Target function for MTPA, MC, FW，MTPV Jacobi matrix It is quite difficult and complicate to give an analytic relationship to describe both magnetic saturation and crosssaturation directly, so in order to reduce measurement repetition, a small table with averaged values from experimental measurements and FEA results is used in practice. The size of original tables is 11×11 shown in Fig.6(a).  Fig.6 Expansion of the inductance data table (a)dq-axis inductance LUT without using SBIM (b)dq-axis inductance LUT using SBIM In order to reduce the execution time discussed in Section VI, the second-order bilinear interpolation method (SBIM) is utilized to expand a small inductance tables to a larger one. In Fig.5, the blue points are four known points listed as: Q11 = f (x1, y1). , Q12 = f (x1,y2), Q21 = f (x2, y1) and Q22 = f (x2, y2), and the blue points are one unknown point P. If the value is required from the unknown function f(x,y) at point P. The SBIM can be expressed in the following steps: Step 1. Using linear interpolation in the x direction Step 2. Using linear interpolation in the y direction Step 3. The final unknown point P is introduced as (26)   11  21  2  2  1  2  2  1  2  1  2  1  2  1   12  22  2  1  1  1  2  1  2  1  2  1  2 Expansion results of the inductance data table can be simulated, as shown in Fig.6(b).

B. The Influence of the Resistive Voltage Drop
If the stator resistive voltage drop is considered, the voltage equation (1) and (2) (27) Noting that the temperature variation was negligible during the experiments presented in Section VI, the resistance is assumed to be constant. Substituting (27) into (4), the equation can be rewritten as (28), which indicates the voltage constraints trajectory accounting for the stator resistive voltage drop.
The corrected curve for (28) forms a series of slant ellipses with speed, whose centers are moving with the increase of the motor speed. If the general equation of the ellipse is defined as 2 As shown in Fig.7, the non-standard ellipse trajectory is composed of four elements: Center coordinates (xc,yc), Semimajor axis a, Semi-minor axis b, and the inclination angle θc, which can be solved with the coefficients (31) : In conclusion, if the resistive voltage drop is considered, the intersection of the current constraint circle and the voltage constraint ellipse, marked as stars in Fig.7, should be moved from point B1 to point B2. Compared with the state-of-art, the N-R method combined with variable inductance LUT proposed in this paper has some advantages listed as follows: 1) The proposed N-R method is based on iterative operations, the intersections of complex nonlinear equations (10)- (16) involving current circle, voltage ellipse and torque hyperbola can be solved in real-time digital controllers.
2) According to different operating region, such as MTPA, MC, FW and MTPV region, the N-R searching method can switch different target functions and Jacobian matrices in Table  I to obtain the optimal current set-points.
3) Thanks to the LUTs, the proposed N-R method can cope with the actual motor parameter variation, such as dq-axis inductance, resistance and permanent magnetic flux linkage. The controllability, accuracy improvement, and reduced computation burden achieved by using the LUTs will be demonstrated in Section VI.

Fig.8 Wide-speed range control method based on N-R method
In order to better verify the proposed control strategy in different operation modes shown in Fig.2, the control loop is constructed like the block diagram shown in Fig.8 and 8kW IPMSM parameters are presented in Table II. Noting that the MTPV locus does not intersect the MC circle, it will not be discussed in this Section V and Section VI. A. The Whole Convergence Process for MTPA Mode Taking initial value id0 iq0 (-30A, 20A) and (-4A, 12A) as the first guess, the iterative convergence process of the N-R algorithm in MTPA mode is shown in Table III. Comparing the operation condition I and II, due to the magnetic saturation effect of the inductance, the inductance is going down with a larger given torque. Comparing the operation condition I and III, symbol E is defined as current setting precision in the equation (23). When E is set to be 0.01 2 , three steps of iteration are required to reach the MTPA set-points (-16.20A, 75.78A). When E is set to be 0.001 2 , it takes four steps of iteration to reach the MTPA set-points (-16.208A, 75.768A). In this case, the accuracy is increased as the cost of more resources occupied. In this particular work, according to current sensor 0.2A/mV, the current setting precision E is set to be 0.01 2 . Fig.9(a) corresponds to Fig.1(a) for the algorithm branch I, the speed set-points is set to be 500rpm less than the base speed 2400rpm. The IPMSM is ramping from 0rpm to 500rpm at max torque 32N.m, which corresponds to point A in dq-axis current plane and then go to steady state at no load, which corresponds to point Lo. The speed command is then stepped to 1000rpm, the IPMSM is ramping to point A and again settles down at point Lo. If the load torque is added at 1s, the operating point is moving from point Lo to point C along MTPA curve and balance at 10N.m finally. Fig.9(b) corresponds to Fig.1(a)(b)(c)(d) for the algorithm branch II, the speed command * is set to be 2800rpm (between base speed 2400rpm and 2953rpm). IPMSM is ramping with max torque 32N.m at point A. When the measured speed is higher than 2400rpm, the operating mode selector switch to the MC region and the flux is weakened by increasing the daxis current set-points in the negative directions. As the measured speed reaches 2800rpm, the dq-axis current starts to settle down to point Lo passing by point D along the voltage ellipses * . If the load torque is added at 1.15s, the operating point is moving from point Lo to point E along MTPA curve. Fig.9(c) corresponds to Fig.1(a)(e)(f) for the algorithm branch III, the speed set-points is firstly set to be 1000rpm below the base speed, the IPMSM is ramping with max torque 32N.m at point A and balance to steady state. When the speed set-point * steps to 3600rpm (higher than the boundary speed   n 0 2953rpm), the d-axis current will decrease towards point B along the MC circle. When the measured speed reaches steady state at 3600rpm, the operating will go along the voltage ellipse towards point Lo and, as soon as the load torque is added at 1.2s, it will balance at point G with the final torque.

VI. EXPERIMENTAL RESULTS
Several experimental tests are carried out to verify the effectiveness of the proposed method. The test bench is made up of a rated 8kW target IPMSM (peak power 20kW) with resolver as position sensor, a 100Nm torque sensor, a 5:1 reducer and a magnetic powder brake as load, as shown in Fig.10. A DSP (Type:TMS320F28335) is used as core device for the controller. The inverter is a three-phase IGBT power module (Type: Infineon FS400R07A1E3), two current sensors, DC-link voltage sensor and some protection circuits. The switching frequency is 10kHz (Ts=100us).

A. Experimental Results For Different Operation
Regions The experimental results for testing MTPA region (branch I, Fig.1(a)) dealing with magnetic saturation are shown in Fig.11(a)(d). The initial speed set-points n * in orange and the measured speed n in pink in Fig.11(a) are 500rpm. When the speed set-points step to 1000rpm (n ≤ n b ), the dq-axis current set-points are -16A and 76.8A (at max torque 32N.m). The dqaxis inductance goes saturation at 0.326mH (decline 2.7%) and 0.522mH (decline 4.2%), respectively, due to the high starting current. After reaching 1000rpm, the dq-axis inductance will be 0.334mH and 0.545mH, quite close to the rated inductance value (Ld=0.335mH and Lq=0.545mH).When the load torque is added to 15N.m, the dq-axis current set-points are -4.8A and 38A, while the dq-axis inductance will be 0.332mH (decline 0.9%) and 0.537mH (decline 1.5%).
The experiment results for testing branch II (corresponding to Fig.1(a)(b)(c)(d)) are shown in Fig.11(b)(e). The initial speed is 1500rpm at steady state. After the speed set-points step to 2700rpm (n b < n ≤ n o ), IPMSM began to accelerate and the given torque (max torque 32N.m) is assigned by N-R searching method to obtain dq-axis current set-points at -16A and 76.8A (point A). When the speed increases higher than the base speed 2400rpm(switch to MC region), deeper flux-weakening d-axis current is required to shift the point from A to B along the current circle. Later on, when the speed reaches the speed setpoints, the reference torques finally balance with a load torque 5N.m along the voltage ellipse in FW region and finally balance on MTPA curve (point E).
The experiment results for testing branch III (corresponding to Fig.1(a)(e)(f)) are shown in in Fig.11(c)(f). The measured speed n is 2000rpm at steady state. After the speed set-points steps to 3600rpm (n > n o ), IPMSM began to enter the speed transition. The given torque (at max torque 32N.m) is assigned by N-R searching method to -16A and 76.8A at point A. When the speed increases higher than the base speed 2400rpm, deeper flux-weakening is required to limit voltage, shifting the point L0 D D D from A to B along the current constraints circle. Later on, the speed is approaching the speed set-points 3600rpm, the reference torques finally balance with a load torque 5 N.m along the voltage constraints ellipse towards point G. (c) branch III (f) branch III Fig.11 The experimental results for different branches: waveform of torque , speed set-points , speed measurement and phase B current(a) branch I (b) branch II (c) branch III; waveform of dq-axis current set-points and measurement (d) branch I (e) branch II (f) branch III

B. Experimental Results For Control Accuracy And
Computation Burden Fig.12 shows the proposed method dealing with parameter mismatch in MC region as discussed in Section IV. At the initial steady state I, the IPMSM is spinning at the speed of 3000rpm. Fig. 12 shows the large q-axis current error when the fixed values of Ld = 0.335 mH, Lq = 0.545 mH and Rs=0Ω are used. The q-axis current measurement differs significantly from the current set-points, due to large mismatch of the motor parameters. Therefore, the performance of the IPMSM is degrading. When the variable parameters have been compensated by LUT later on in steady state II, the q-axis current set-points control errors reduce from about 20A to about 0A. The main reason is that the corrected voltage ellipse considering the resistive voltage drop in the dq-axis current plane forms a series of slant ellipses. The resistive voltage drop forces the current set-points to move from point B1 to point B2 corresponding to Fig.7, discussed in section IV. This proves the necessity of updating Ld and Lq and considering the influence for resistive voltage drop.
In order to test the execution time of different methods, the IO pin in the DSP is set to be high before the proposed method and reset it to low after the function has executed, as shown in Fig.13(a). The duration of the pulse can be measured in the oscilloscope, shown as histogram in Fig.13(b). In the method 1, the N-R method is applied with dq-axis inductance matrices (size:11×11) and expanded by SBIM took least computation time 45.85us. In method 2, the second-order PA method is used for the inductance curve-fitting mentioned in [23] and [24], the solution in Table I occupies 56.2us. In method 3, the N-R method needs 65.8us combined with 61 × 61 inductance matrices. In method 4, the square root operation is saved in advance using premade tables and 61×61 inductance matrices are applied for parameter compensation. Comparing method 3 and method 4, the proposed N-R method occupied 18.1us less resource than square root operation, since the proposed N-R method is based only on iteration and selection operations. Comparing method 1 and method 3, the SBIM can reduce the size of LUTs and 20us execution time margin is saved.

VII. CONCLUSION
In this paper, the operating modes for IPMSM are divided into MTPA region, MC region, FW region and MTPV region. The target non-linear functions for optimal current set-points are presented. In order to get more accurate solution, N-R feedforward searching method is proposed to find the optimal current set-points in iterative forms. It has been highlighted that neglecting the resistive voltage drop and magnetic saturation leads the actual current trajectory to deviate from the optimal one, resulting in reduced accuracy. Therefore, LUT combined with second-order bilinear interpolation method is used for variable parameters compensation.
The simulations and experimental results are demonstrating the validity and accuracy of the proposed searching method. This has been verified in various operation cases and also the case dealing with parameter mismatch. The applied method has been proved to improve the control accuracy of IPM machines in wide-speed range, while reducing the computation burden at the same time.