 # Accelerated Steady-State Torque Computation for Induction Machines using Parallel-In-Time Algorithms

This paper focuses on efficient steady-state computations of induction machines. In particular, the periodic Parareal algorithm with initial-value coarse problem (PP-IC) is considered for acceleration of classical time-stepping simulations via non-intrusive parallelization in time domain, i.e., existing implementations can be reused. Superiority of this parallel-in-time method is in its direct applicability to time-periodic problems, compared to, e.g, the standard Parareal method, which only solves an initial-value problem, starting from a prescribed initial value. PP-IC is exploited here to obtain the steady state of several operating points of an induction motor, developed by Robert Bosch GmbH. Numerical experiments show that acceleration up to several dozens of times can be obtained, depending on availability of parallel processing units. Comparison of PP-IC with existing time-periodic explicit error correction method highlights better robustness and efficiency of the considered time-parallel approach.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

Induction machines are used in a wide variety of industrial applications. Their operation covers the power ranges from hundreds of watts to several megawatts. For the design of these motors, engineers are often interested in the steady-state operating characteristics such as, e.g., mean torque, efficiency, periodically changing currents and voltages at certain revolution speeds. During initial design stages, information about these quantities can be obtained from numerical analysis of transient eddy current problems using space and time discretization , see, e.g., Figure 1. Application of an implicit time stepping leads to a nonlinear system of equations at each time step. Numerical simulation of induction motors is often computationally expensive because a lot of time steps might need to be calculated to reach the steady state. Therefore, efficient numerical methods are necessary to accelerate these calculations.

Different methodologies have been developed in order to achieve the periodic steady-state solution. For instance, the time-periodic finite element method  requires solution of a coupled periodic system, which could be prohibitively expensive for modern real-life applications. Accelereted convergence to the steady state with the simplified time-periodic explicit error correction method (TP-EEC) [3, 4] is based on correction of the sequential solution after each (half-) period. However, it performs well only for a specific class of problems, particularly when the time constant is big enough.

The Parareal algorithm [5, 6] is a powerful approach, which allows to speed up the sequential solution of the underlying evolution problem via parallelization in the time domain. The method has been recently applied to simulation of an induction machine by the authors in , where quick convergence and efficiency of Parareal is illustrated. The periodic Parareal algorithm with initial-value coarse problem (PP-IC), introduced in , is a natural extension of the original Parareal to the class of time-periodic problems. A multirate version of PP-IC was introduced and applied to a simplified electrical machine model, operated at synchronous speed, by the authors in . In this paper we investigate the industrial use case of an asynchronous (induction) machine and apply PP-IC to the underlying time-periodic problem. Performance of the method as well as its comparison with the simplified TP-EEC method is illustrated via steady-state analysis of an induction machine model, developed at the Robert Bosch GmbH, Germany. Figure 1: Finite element model of an induction motor of electric vehicle drive, provided by Robert Bosch GmbH.

## Ii Problem Setting and Discretization

The eddy current problem in A

-formulation with magnetic vector potential

is given by

 σ∂t→A+∇×(ν∇×→A)=→Js(t) (1)

on with computational domain , depicted in Figure 2, and time interval The regions and denote the conductive () and stranded conductors subdomains. Input currents are homogeneously distributed by means of the stranded-conductor winding functions  and form the current density .

Nonlinearity of the eddy current equation is given by the reluctivity function . To obtain a well-posed problem we complete (1) with gauging (e.g., the Coulomb gauge ) and a suitable boundary condition such as Dirichlet , where . We also prescribe an initial value .

### Ii-a Computation of the Torque

The torque-speed characteristic is one of the most important technical features of an electrical machine. The torque can be calculated with the eggshell method , . Within this approach the moving rigid piece is surrounded by a hull

, whose thickness does not need to be constant. Movement is then described by the deformation of the eggshell region only. Using the Maxwell stress tensor

and the velocity field associated to a displacement of the moving body by an infinitesimal distance one can write

 ˙WEM=−F⋅δ˙u=∫SσEM:∇vdS (2)

for the mechanical power in the eggshell . The velocity and its gradient are determined by and respectively, where is any smooth function that is equal to on the inner surface of and is on its outer surface. With these definitions the resultant force on the moving piece is given by

 F=−∫SσEM⋅∇γdS. (3)

The eggshell formula (3) can be then used for calculation of the resultant electromagnetic torque as

 T=−∫Sr×(∇γ⋅σEM)dS. (4)

### Ii-B Spatial Discretization

The Ritz-Galerkin approach leads to the following weak formulation for

 ∫Ω→w⋅σ∂t→A+∇×→w⋅(ν∇×→A)dΩ =∫Ω→w⋅→JsdΩ

for all . For a rotating machine, approximation by edge elements 

 →A(→x,t)≈n∑i=1→wi(→x)ai(t)

yields the following system of differential algebraic equations (DAEs):

 Mσdta(t)+kν(a(t),θ(t)) =js(t) (5)

for the unknown (line-integrated) magnetic vector potentials Here denotes the (singular) mass matrix, function is given by with the curl-curl matrix , which depends on the rotor angle and is the discretized source current density.

Rotation is modeled with the moving band approach [14, 15] and is described by the mechanical equations for the angle and angular velocity

 dtθ(t)=ω(t)andIdtω(t)+Cω(t)=T(a(t)), (6)

where

is the moment of inertia,

denotes the friction coefficient, and is the mechanical excitation determined by the magnetic field.

Assigning initial conditions and together with and combining (5) and (6) we obtain the coupled problem

 Mdtu(t)+K(u(t))u(t) =f(t) (7)

with respect to the -dimensional unknown . The matrices and are given by

 M=⎛⎜⎝Mσ0001000I⎞⎟⎠,K=⎛⎜⎝kν(⋅,θ)0000−1T(⋅)0C⎞⎟⎠

and the right-hand side

### Ii-C Time-Integration Scheme

Since equation (7) is an index-1 DAE 

, by means of the implicit Euler method it can essentially be treated as an ordinary differential equation

. In this case the time-stepping scheme propagating the solution from to is written as

 (8)

This time-integration formula defines a numerical solution operator

 F:I×I×Rm→Rm (9)

such that . Similarly, we introduce a coarse propagator

 G:I×I×Rm→Rm, (10)

which, as solves the initial-value problem (IVP) for (7), however uses a lower precision, i.e., a time step

## Iii Parallel-in-Time Methods

### Iii-a Parareal Algorithm

The main idea of the Parareal method is to parallelize sequential time stepping by distributing the calculations among available CPUs. First, the time interval is divided into subintervals with . Equation (7) has to be solved on each subinterval with initial value to obtain the final value Applying the propagator from (9) to each IVP the matching conditions at the synchronization points are imposed as

 H(U):=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩U0−u0=0,U1−F(T1,T0,U0)=0,⋮UN−F(TN,TN−1,UN−1)=0. (11)

In fact, equation (11) is a root-finding problem for the nonlinear operator with respect to its argument Application of the Newton method and a finite difference approximation of the Jacobian [17, 7] with the coarse solver defined in (10) gives the Parareal update formula:

 U(k)0:= u0, (12) U(k)j:= F(Tj,Tj−1,U(k−1)j−1) +G(Tj,Tj−1,U(k)j−1)−G(Tj,Tj−1,U(k−1)j−1) (13)

for and The pseudo code for (12)-(13) and a detailed explanation of the iterative procedure is presented in .

### Iii-B PP-IC Iteration Figure 3: Comparison of PP-IC and time stepping when visualizing the resulting waveform of iteration k on time interval [t0+(k−1)T,t0+kT] for k=1,…,10 with t0=0 and T=0.02 s. Illustrated data is based on the solution of a linear RL-circuit model excited with a 50 Hz sinusoidal current source.

Aiming at the steady state of an induction machine we would now like to exploit a version of the Parareal method, adapted to time-periodic problems. Such an algorithm was introduced in  and was called PP-IC. A time-periodic formulation for (7) on period requires initial value to be equal to the final one . Within the Parareal setting it means substitution of the first equation in (11) for the periodicity condition Analogous derivations as those performed for (11) together with an additional relaxation on the coarse grid give the PP-IC iteration:

 U(k)0:= U(k−1)N, (14) U(k)j:= F(Tj,Tj−1,U(k−1)j−1) +G(Tj,Tj−1,U(k)j−1)−G(Tj,Tj−1,U(k−1)j−1) (15)

for and It can be seen that PP-IC (14)-(15) is based on the Parareal iteration (12)-(13) with the only difference in the update of the initial value at iteration We note that in contrast to the classical Parareal method, which convergences superlinearly in iterations as it is proved in , convergence of PP-IC is only linear .

The iterative procedure of PP-IC is summarized with the pseudo code in Algorithm 1. Step updates the solution at the beginning of the period by assigning the end value from the previous iteration. Starting from the corrected sequential solutions are performed by means of the coarse propagator at step and are expected to be computationally cheap due to a large step size Calculations on the fine grid, obtained via discretization with a small time step can be performed for each subinterval in parallel (step ), starting from initial values already given from steps and . We choose both coarse and fine solvers to be the implicit Euler method, using low and high fidelity, respectively. Figure 4: Torque-speed characteristic and some representative operating points of the induction motor from Figure 1.

PP-IC is an iterative approach applied to a fixed time interval, on which the periodicity constraint is to be satisfied. However, PP-IC can be reinterpreted as applying the classical Parareal method on period with initial guess at and using only a single iteration. PP-IC then proceeds to subsequent time intervals with a single iteration each until the periodicity constraint together with the matching conditions at every for are fulfilled at some We illustrate this interpretation of PP-IC as a forward-in-time Parareal iteration in Figure 3. It is visible that the solution after the first iteration of Algorithm 1, when only initial fine solves have been executed on subintervals in parallel, contains jumps at the synchronization points. However, these discontinuities quickly decrease already after the next PP-IC update at iteration and the obtained solution nearly replicates the classical time stepping solution on the second period ( s here). Further PP-IC iterations, in the same way as the standard time-stepping, eventually reduce the periodicity error up to a prescribed tolerance and deliver the periodic steady-state solution.

### Iii-C Estimation of Computational Costs

For a clear estimation of computational costs of the Parareal algorithms considered in this paper we calculate the number of effective time steps instead of the direct wall clock time measurement. By effective time steps we mean the following. First, assume the number of available processing units is equal to

Then we consider the splitting of the time domain into subintervals with exactly one coarse step per subinterval. Second, parallelization among the CPUs allows us to take into account fine computations on one subinterval only. We denote the number of fine steps per parallel process by . Effective time steps after iterations could then be counted using the formula

 Ne=It⋅(N+Np). (16)

At every iteration time steps are calculated sequentially. Multiplication with the number of iterations results in the overall time steps which are not executed in parallel.

We note that the described approach to measure the computational costs is legitimate only if the solution at a time step is always performed with a comparable effort and communication costs are negligible. In general, this may not hold, especially since both fine and coarse time steps are considered. Within the current implementation for the examples presented in the next section, the solver always requires - Newton iterations to solve a time step. Estimate with formula (16) is therefore justifiable.

## Iv Numerical Results

The torque-speed characteristic of an induction motor is determined by evaluation of some tens to hundreds of operating points, i.e., transient eddy current simulations at a constant revolution speed and mostly sinusoidal coil currents. In electric vehicles such operating points correspond to certain driving conditions. In the following we apply the PP-IC algorithm to the induction motor whose two-dimensional mesh view with degrees of freedom is depicted in Figure 1. Due to symmetry only half of the geometry is modeled and discretized. The torque-speed characteristic as well as some representative operating points of the motor are shown in Figure 4.

To demonstrate the capability of PP-IC we consider two representative operating points of the machine. OP1 describes the operation limit at maximum speed of rpm while OP2 corresponds to the maximum acceleration at low speed. The torque evolutions, obtained with the classical time stepping, are depicted in Figure 5. We note that within the performed simulations the rotational speed was prescribed by the operating points, which omits solution of mechanical equations (6).

Within the time-marching procedure we consider the steady-state solution to be attained up to a prescribed tolerance at instant which belongs to time interval Period is defined here based on the (relative) periodicity error given in terms of torque by

 err(k)=|T(t0+(k−1)T)−T(t0+kT)||T(t0+kT)|. (17)

We denote the periodicity error of the calculated torque in the steady state by Performance of the considered acceleration methods will be estimated in comparison to the sequential time stepping until

### Iv-a Operating point 1

We consider the operating point OP1 at the revolution speed of rpm and the sinusoidal three-phase input current with the amplitude of A. The resulting mean torque in the steady state has to be determined by the transient eddy current simulation.

First, we apply PP-IC on the time interval s, i.e., on one electric period in rotor, using parallel processors. iterations of the algorithm with the fine step s produced a periodic solution up to tolerance with respect to torque. The simulation required effective time steps to calculate the torque, periodic up to tolerance . This same value is also a bound for the periodicity error of the steady-state solution, obtained from the classical time stepping on period i.e., on s. The effort of the sequential computations is evaluated by time steps of size , calculated on s. Comparing this to the number of effective time steps performed within PP-IC, we get the speedup of factor due to time-parallelization. We would like to note that period s has been empirically determined, more sophisticated approaches could include as an additional variable to be determined by Parareal.

The periodic torque, obtained with the PP-IC iteration is depicted in Figure 5(a), together with the reference time-stepping behavior. We note that the calculated PP-IC solution still contains a residual oscillation remains, although being almost perfectly periodic. This is in contrast to the sequential steady-state solution, which eventually flattens out. The deviation (by Nm), however, does not exceed of the mean torque of Nm. Therefore, we consider this oscillation to be acceptable, since it is in a good agreement with the chosen periodicity tolerance .

On the other hand, refining the fine step size by the factor of i.e., using s, let the residual oscillation vanish and delivered the steady-state solution, compatible with the sequential calculation of the same time-stepping precision. We conclude that PP-IC performs well but is more sensitive to the step size of the fine propagator than an ordinary (sequential) time domain simulation.

We now compare the performance of PP-IC with the simplified TP-EEC method . The main idea of the approach is a successive reduction of the time-stepping solution at every half a period by the average of the initial value and the final value until the steady state is attained. Application of this method to OP1 shows that time stepping with the TP-EEC correction at every s converges to the steady state about times faster than the ordinary time-marching procedure of steps. Time-stepping calculation for OP1, corrected with the simplified TP-EEC, is depicted in Figure 5(b), where the steady-state solution, periodic up to tolerance , is obtained already after the second correction.

Please note, TP-EEC is not a parallel method and requires less overall computing power than Parareal. Furthermore, the speedup provided by TP-EEC could be increased by its joint application with the Parareal algorithm, e.g., via application of Parareal on half-period followed by the TP-EEC correction at

### Iv-B Operating point 2

At OP2 the rotor speed is rpm and peak value of the phase currents is A. In contrast to the previously considered example, OP2 has a different transient behavior: there is no significant overshoot of the mean torque (Figure 4(b)). PP-IC, applied on one rotor period s with cores and fine step size s converged in iterations, reaching the steady-state solution, periodic up to tolerance in terms of torque, and thereby requiring calculation of only effective time steps. On the other hand, classical time stepping led to the torque with relative periodicity error on period . Compared to sequential steps of the same (fine step) size performed on s, a speedup of factor is achieved with PP-IC.

An attempt to apply the simplified TP-EEC to OP2 unfortunately did not converge. Correction at every half of the period s set the transient solution even further apart from reaching the steady state. We summarize the computational costs of PP-IC and simplified TP-EEC for the two considered operating points, in contrast to the standard sequential time stepping until the steady state, in Table I.

## V Conclusion

A periodic Parareal algorithm with initial-value coarse problem (PP-IC) is a suitable approach for accelerated attainment of the steady state, since it is applied to a time-periodic problem directly, in contrast to the standard Parareal method, which solves an IVP without taking into account the periodicity condition. Application of PP-IC to an industry-relevant problem illustrated a significant speedup in terms of effective time steps. Numerical experiments illustrated that convergence of PP-IC is more sensitive to the fine step size, compared to the classical time stepping. Nonetheless, the algorithm could be successfully applied to all tested operating points, whereas the simplified TP-EEC method delivered feasible results only in some cases. In the applicable case (OP1) the speedup obtained with the simplified TP-EEC is smaller than that of PP-IC, however the Parareal-based approach needs more (parallel) computing power. Note that for steady-state calculation of asynchronous machines with methods, based on the periodicity condition, additional research should be dedicated to a suitable choice of the common period . Further development of efficient numerical methods for accelerated steady-state analysis could be based, e.g., on another parallel-in-time algorithm, tailored to periodic problems, called periodic Parareal algorithm with periodic coarse problem (PP-PC) , where the periodicity constraint on the coarse grid is imposed explicitly.

## References

•  A. Arkkio, “Finite element analysis of cage induction motors fed by static frequency converters,” IEEE Trans. Magn., vol. 26, no. 2, pp. 551–554, Mar. 1990.
•  T. Nakata, N. Takahashi, K. Fujiwara, K. Muramatsu, H. Ohashi, and H. L. Zhu, “Practical analysis of 3-d dynamic nonlinear magnetic field using time-periodic finite element method,” IEEE Trans. Magn., vol. 31, no. 3, pp. 1416–1419, May 1995.
•  Y. Takahashi, T. Tokumasu, A. Kameari, H. Kaimori, M. Fujita, T. Iwashita, and S. Wakao, “Convergence acceleration of time-periodic electromagnetic field analysis by singularity decomposition-explicit error correction method,” IEEE Trans. Magn., vol. 46, no. 8, pp. 2947–2950, Jul. 2010.
•  H. Katagiri, Y. Kawase, T. Yamaguchi, T. Tsuji, and Y. Shibayama, “Improvement of convergence characteristics for steady-state analysis of motors with simplified singularity decomposition-explicit error correction method,” IEEE Trans. Magn., vol. 47, no. 5, pp. 1458–1461, May 2011.
•  J.-L. Lions, Y. Maday, and G. Turinici, “A parareal in time discretization of PDEs,” Comptes Rendus de l’Académie des Sciences – Series I – Mathematics, vol. 332, no. 7, pp. 661–668, 2001.
•  M. J. Gander and E. Hairer, Nonlinear Convergence Analysis for the Parareal Algorithm.   Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, pp. 45–56.
•  S. Schöps, I. Niyonzima, and M. Clemens, “Parallel-in-time simulation of eddy current problems using parareal,” IEEE Trans. Magn., vol. 54, no. 3, pp. 1–4, Mar. 2018.
•  M. J. Gander, Y.-L. Jiang, B. Song, and H. Zhang, “Analysis of two parareal algorithms for time-periodic problems,” SIAM J. Sci. Comput., vol. 35, no. 5, pp. A2393–A2415, 2013.
•  M. J. Gander, I. Kulchytska-Ruchka, and S. Schöps, “A new parareal algorithm for time-periodic problems with discontinuous inputs,” 2018, arXiv:1810.12372.
•  S. Schöps, H. De Gersem, and T. Weiland, “Winding functions in transient magnetoquasistatic field-circuit coupled simulations,” COMPEL, vol. 32, no. 6, pp. 2063–2083, Sep. 2013.
•  P. Monk, Finite Element Methods for Maxwell’s Equations.   Oxford: Oxford University Press, 2003.
•  F. Henrotte, G. Deliége, and K. Hameyer, “The eggshell approach for the computation of electromagnetic forces in 2D and 3D,” COMPEL, vol. 23, no. 4, pp. 996–1005, 2004.
•  F. Henrotte, M. Felden, M. van der Giet, and K. Hameyer, “Electromagnetic force computation with the eggshell method,” in 14th International IGTE Symposium on Numerical Field Calculation in Electrical Engineering, Graz, 2010.
•  B. Davat, Z. Ren, and M. Lajoie-Mazenc, “The movement in field modeling,” IEEE Trans. Magn., vol. 21, no. 6, pp. 2296–2298, Nov. 1985.
•  M. V. Ferreira da Luz, P. Dular, N. Sadowski, C. Geuzaine, and J. P. A. Bastos, “Analysis of a permanent magnet generator with dual formulations using periodicity conditions and moving band,” IEEE Trans. Magn., vol. 38, no. 2, pp. 961–964, Mar. 2002.
•  A. Nicolet and F. Delincé, “Implicit Runge-Kutta methods for transient magnetic field computation,” IEEE Trans. Magn., vol. 32, no. 3, pp. 1405–1408, May 1996.
•  M. J. Gander and S. Vandewalle, “On the superlinear and linear convergence of the parareal algorithm,” in Domain decomposition methods in science and engineering XVI, ser. Lecture Notes in Computational Science and Engineering.   Berlin: Springer, 2007, vol. 55, pp. 291–298.