    # Multirate PWM balance method for the efficient field-circuit coupled simulation of power converters

The field-circuit coupled simulation of switch-mode power converters with conventional time discretization is computationally expensive since very small time steps are needed to appropriately account for steep transients occurring inside the converter, not only for the degrees of freedom (DOFs) in the circuit, but also for the large number of DOFs in the field model part. An efficient simulation technique for converters with idealized switches is obtained using multirate partial differential equations, which allow for a natural separation into components of different time scales. This paper introduces a set of new PWM eigenfunctions which decouple the systems of equations and thus yield an efficient simulation of the field-circuit coupled problem. The resulting method is called the multirate PWM balance method.

## Authors

##### This week in AI

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

## 1 Introduction

Switch-mode power converters are used in various devices from small-scale applications like mobile phone chargers to industrial large-scale applications like welding devices . These converters use transistors to switch on and off the input voltage to produce an output voltage, which, in average, has the desired amplitude. A filter circuit is used to smoothen the output. The simulation of these devices is computationally expensive since, through the transistor switching, steep transients occur in the converter. Furthermore often a switch-event detection is necessary to avoid step size rejection or even solver failures . A multirate method has been developed in [3, 4] which uses the concept of Multirate Partial Differential Equations (MPDEs) [5, 6] and a combination of a Galerkin ansatz and conventional time discretization to efficiently solve problems with pulsed excitation. The method is applicable to power converters in which the switching behavior is idealized and known a-priori. It is particularly efficient in the case of linear elements. Some circuit elements may only be accurately represented by field models. For example the induced currents in the conducting materials of an inductor usually cause eddy current losses, which can easily be accounted for in a field model but not in a circuit model. In this paper the multirate method from [3, 4] is applied to a linear buck converter circuit (see Fig. 1) in which the inductor is represented by a 2D finite element model. This substantially increases the size of the strongly coupled system of equations. To still ensure an efficient simulation, a basis transformation is applied to the pulse-width modulation (PWM) basis functions  leading to decoupled systems of equations which can be solved efficiently in parallel. The resulting method is called the multirate PWM balance method in analogy with the harmonic balance method where harmonic functions take the place of the PWM basis functions. Numerical results on the buck converter show the efficiency and accuracy of the proposed method in field-circuit coupled problems.

The paper is structured as follows. Section 2 introduces the concept of MPDEs and explains the solving procedure using Galerkin approach and conventional time discretization. Subsequently section 3 presents the original PWM basis functions as described in . In section 4 the PWM eigenfunctions are developed and their advantageous properties for the solving process are highlighted. Finally section 5 summarizes numerical results and compares the three different solution approaches, i.e., conventional time discretization and the MPDE approach with PWM basis functions on the one hand and PWM eigenfunctions on the other hand. The paper is concluded by summarizing its content in section 6.

## 2 Multirate formulation

Let the field-circuit coupled model  of the converter be described by the system of ordinary differential or differential-algebraic equations

 =c(t), (1)

where is a possibly singular matrix, is assumed to be a regular matrix, is the unknown solution, is the excitation, and is the simulation interval. The initial state of the model is given by consistent  initial values . The ideal pulsed excitation

 vi(t)={V0for all τ(t)≤D0otherwise, (2)

is used as input of the power converter circuit. We denote by the relative time, is the switching cycle and is the duty cycle.

The system of Multirate Partial Differential(-Algebraic) Equations (MPDEs or MPDAEs) with two time scales corresponding to (1) is given by [5, 6, 10]

 A(∂ˆx∂t1+∂ˆx∂t2)+Bˆx(t1,t2)=ˆc(t1,t2), (3)

where is the unknown multivariate solution and is the multivariate excitation. Choosing the multivariate excitation such that , the solution of (1) and (3) are related by . Thus, the solution of (1) can be calculated solving the MPDEs and extracting the solution along a diagonal through the computation domain. To solve the MPDEs, additional conditions need to be specified. For the present application, a combination of initial and boundary values is applied. Initial values are supplied by , i.e., at , where with is a function defining the initial values for all . The solution along the fast time scale is periodic, i.e., . The multivariate right-hand side is chosen as , i.e, the pulses of the excitation occur along the fast time scale. It is possible to use MPDEs with more than two time scales. However, in the applications of this paper, it is not necessary and furthermore often not feasible since the dimension of the computation domain increases and thus also the computational effort to calculate the solution. Figure 1: Simplified circuit of the buck converter in continuous conduction mode with C=10μF and R=30Ω. The field model of the pot inductor (axisymmetric around z-axis) is designed to have an inductance of L=65mH and a series resistance of RL=800mΩat DC. The figure shows the equipotential lines of the magnetic vector potential.

To solve the MPDEs (3), a Galerkin approach and time discretization is applied [4, 11]. The -th solution component is approximated by expanding it into periodic basis functions depending on the fast time scale and coefficients depending on the slow time scale

 ˆxhj(t1,t2):=Np∑k=0wj,k(t1)pk(τ(t2)), (4)

where the periodicity of the basis functions is accounted for by using the relative time Applying the Galerkin approach with respect to and over one period of the excitation leads to

with block matrices , where

 \cI=Ts1∫0¯p(τ)p⊤(τ)dτ,\cQ=−1∫0d¯p(τ)dτp⊤(τ)dτ, (6)

and right-hand side

 \cC(t1)=Ts∫0¯p(τ(t2))⊗ˆc(t1,t2)dt2. (7)

denotes the complex conjugate of and denotes the Kronecker product.

## 3 PWM basis functions

The PWM basis functions developed in  are built up starting from the zero-th constant basis function and the piecewise linear basis function

 p1(τ) = ⎧⎨⎩√3 2τ−DDif   0≤τ≤D√3 1+D−2τ1−Dif   D≤τ≤1, (8)

which includes the duty cycle of the excitation by construction. The higher-order basis functions , are recursively obtained by integrating the basis functions of lower order and orthonormalizing them using the Gram-Schmidt algorithm. The generated basis functions are depicted in Fig. 2.

For the PWM basis functions, the matrices and from (6) are given by

multiplied by the identity matrix (due to the orthonormality of the basis functions) and a square matrix with around 25 % of non-zero entries, respectively. Solving the problem requires time stepping of the entire system (

5).

## 4 PWM eigenfunctions

To enable an easy parallelization of the method, the equations (5) can be decoupled, for example by diagonalizing , i.e., a basis transformation. We define new basis functions as linear combinations of the PWM basis functions, i.e.,

 gk(τ):=Np∑l=0vk,lpl(τ), (9)

where are unknown coefficients with , and are eigenfunctions of the time derivative operator

 ddτgk(τ)=λkgk(τ). (10)

We enforce this property in a weak sense by a Galerkin approach, i.e.,

 −1∫0gk(τ)dpm(τ)dτdτ=λk1∫0gk(τ)pm(τ)dτ, (11)

where integration by parts and the periodicity of the basis functions is used. Inserting the expansion of the basis functions into (11) gives

 Ts\cQvk=λk\cIvk. (12)

Since is multiplied by the identity matrix (thanks to the orthonormality of the PWM basis functions), the and

are the eigenvalues and eigenvectors of the matrix

, respectively. Furthermore since

is real-valued and skew symmetric, and therefore a normal matrix, the eigenvectors

are orthonormal. The new basis functions (complex-valued) are depicted in Fig. 3 for . Note that the basis consists of pairs of conjugate complex basis functions. Figure 3: PWM eigenfunctions gk(τ), k∈{0,1,2,3,4}, i.e., Np=4. (top) real part. (bottom) imaginary part.

Inserting the transformed basis functions instead of the PWM basis functions into (6) and (7) leads, using the orthonormality of the eigenvectors, to

where as in (5),

 ˜\cB =\cI⊗B+Λ⊗A, (14) ˜\cC(t1) =Ts∫0¯g(τ(t2))⊗ˆc(t1,t2)dt2, (15)

and is a diagonal matrix with diagonal entries . Thus the resulting matrices in (13) are block-diagonal and the degrees of freedom can be block-wisely decoupled. This leads to independent systems of equations given by

where Note that if a diagonal entry in is complex, there is also a complex conjugate counterpart. The solutions of the decoupled system of equations corresponding to this complex eigenvalue and its conjugate complex counterpart, are, as a result, complex conjugate to each other. Therefore it is sufficient to solve one of them. This is similar to harmonic balance methods in which the harmonic basis functions are given by pairs of complex conjugates leading to similar systems of equations. In analogy to “harmonic balance method”, we call the developed method the “multirate PWM balance method”.

## 5 Test case and numerical results

The method is applied to the buck converter from Fig. 1, where the pot inductor is represented by a 2D field model with conducting core material (ferrite, S/m). The coils are modeled as stranded conductors. The simulation interval is given by ms. The switching frequency is Hz. For the pulsed excitation (2) we use V. All calculations are performed in MATLAB. The partial differential equations governing the magnetoquasistatic problem are given by

 σ(r)∂Am(r,t)∂t+∇×((μ(r)−1)∇×Am(r,t))=Js(r,t), (17)

where is the position vector, is the time, is the modified magnetic vector potential , are the imposed currents, H/m is the magnetic permeability and is the conductivity which is only non-zero in the ferrite core (). The problem is considered on a 2D planar domain with homogeneous Dirichlet boundary conditions.

Correspondingly, the Finite Element magnetoquasistatic  discretization of the magnetoquasistatic inductor model is given by the differential-algebraic system of equations 

 Mσddta(t)+Ka(t)=PiL(t), (18)

where is the singular conductivity matrix, is the stiffness matrix, gathers the degrees of freedom (DOFs) related to the magnetic vector potential, is the discretization of the winding function  and is the current through the inductor. The field-circuit coupling is expressed as follows. An additional variable is introduced for the magnetic flux linkage . All equations are coupled monolithically into the index-1 differential-algebraic system of equations 

 P⊤a−Φ =0, (19) CddtvC−iL+1RvC =0, (20) Mσddta−PiL+Ka =0, (21) ddtΦ+RLiL+vC =vi(t), (22)

which for the example in Fig. 1 contains a total of DOFs. The initial conditions are given by , and .

The initial condition for the MPDEs (3) can be written as

 hj(t2)≈Np∑k=0wj,k(0)pk(τ(t2)) (23)

where is the -th element of . It only has to satisfy the condition . Consequently there is a high degree of freedom in choosing the initial values for the system of equations (5). However not all choices lead to an efficient simulation, i.e., low dynamics in the slow time scale. The following choice of initial values has proven advantageous. First, the steady-state solution is calculated, i.e.,

 ws=\cB−1\cC(0). (24)

Secondly, the initial coefficients for are extracted from the steady-state solution for all . The remaining coefficients are calculated by solving the solution expansion (4) for and using the condition . In summary the initial coefficients are given by

 wj,k(0)=⎧⎪ ⎪⎨⎪ ⎪⎩wsj,kfor k=1,…,Np and for all j=1,…,Nsxj(0)−Np∑l=1wsj,lpl(0)for k=0 and for all j=1,…,Ns. (25)

The initial conditions for the system of equations (13) are computed similarly using the PWM eigenfunctions. Other choices of initial values may still lead to the correct solution but might impair the efficiency of the method.

To calculate the reference solution with a conventional adaptive time discretization, the MATLAB solver ode15s

is used. It is modified to restart the simulation at the known switching instances. Consistent initial values for the restart of the solver are calculated by using a Newton-Raphson algorithm to solve the set of algebraic equations. The required differential variables are taken from the solution at the end of the prior solution interval. After finding the new set of initial values, the initial slopes of the differential variables are calculated by solving the subsystem of ordinary differential equations for the slope. Figure 4: Multivariate voltage at the capacitor calculated using the multirate PWM balance method. The solution component corresponding to the original system of equations (1) is extracted along a diagonal and marked as a black curve. Figure 5: (top) Reference solution calculated using conventional adaptive time discretization compared to the solution obtained by the MPDE approach with Np,pwmbal=4 PWM eigenfunctions. The relative L2 error of the current through the inductor similar to (27) is approximately 3×10−5. (bottom) Joule losses in the core material due to eddy currents.

The multivariate solution calculated using the multirate PWM balance method, i.e., solving (13) with ode15s, is reconstruced using (4) and the multivariate voltage at the capacitor is depicted in Fig. 4. The corresponding solution component of the original system of equations (1) is extracted along a diagonal through the computation domain. Fig. 5 shows the current through the inductor along with the reference solution. The agreement between the multirate PWM balance method solution and the reference solution is excellent. The Joule losses in the core material due to eddy currents are calculated by

 Peddy(t)=∫ΩE(r,t)⋅σ(r)E(r,t)dr=(e(t)H)Mσe(t), (26)

where is the electric field strength, is the spatial computation domain, the superscript denotes the Hermitian, i.e., the complex conjugate transposed, and is the line-integrated discrete electric field. The Joule losses are plotted as well in Fig. 5. Fig. 7 depicts the solution of (13), i.e., the coefficients , exemplary for the current through the inductor . As one can see, using the initial values (25), only the coefficient corresponding to the zero-th basis function varies and the others stay constant. To quantify the accuracy and efficiency of the multirate PWM balance method, it is compared to conventional time discretization and to the MPDE approach with the original PWM basis functions. Different settings are considered: To analyze the performance of the conventional time discretization, the relative and absolute tolerance setting of the solver is changed, i.e., ; For the case of the multirate PWM balance method and the MPDE approach with the original PWM basis functions, relative and absolute tolerances are fixed at and the number of basis functions is changed. The accuracy is measured for the voltage output of the converter, i.e., the voltage at the capacitor. The relative error is given by

 ϵ(tol,n)=||vC,ref(t)−vhC(tol,n,t)||L2(Ψ)||vC,ref(t)||L2(Ψ), (27)

where is the reference solution and is the solution using the multirate PWM balance method and the MPDE approach with the original PWM basis functions. The norm is approximated using mid-point quadrature. Fig. 6 shows the error plotted as a function of the solution time, i.e., the time that ode15s needs. For conventional time discretization the time for solving consists of the time that is needed to calculate consistent initial values and slopes after switching, and the actual time that ode15s needs. The time to calculate consistent initial values and slopes depends on the number of switching instants and is thus constant if the switching frequency or simulation interval does not change. It is given by approximately 16 seconds. The total time displayed in Fig. 6 is the sum of both contributions. Figure 6: Error ϵ as defined in (27) over time for solving the systems of equations. The MPDE approach with PWM eigenfunctions (multirate PWM balance method) is considerably faster than the MPDE approach with the original PWM basis functions and the conventional time discretization. Figure 7: Coefficients w1,k for the inductor current calculated by solving (13) with Np,pwmbal=4. (top) real part. The coefficients w1,1,…,w1,4 are approximately the same therefore they are hard to distinguish visually. (bottom) imaginary part.

As one can see the MPDE approach with the original PWM basis functions is considerably slower than conventional time stepping. This is due to the fact that the already large systems of equations (1) (due to field-circuit coupling) are even further increased in size through the Galerkin approach. The stagnation of the error at in Fig. 6 for values larger than is caused by the chosen accuracy of ode15s. Furthermore one can see that when adding another basis function the error does decrease with every second basis function. This was already observed in [4, 7]. For this reason the error for the PWM eigenfunctions is only plotted for . Since the systems of equations resulting from the multirate PWM balance method are decoupled, they can be solved efficiently in parallel. For each basis function with , a complex-valued initial value problem of the form (16) has to be solved. The size of these systems of equations is the same as that of the original system of equations (1). However, the time for solving is considerably smaller since less time steps are necessary for the same solution accuracy. Note that due to the choice of the initial values (25) most coefficients in (13) for this numerical example do not change and only those corresponding to the zero-th basis function vary. This means that only the decoupled system of equations which corresponds to the zero-th basis function takes considerable computational effort to solve. In a parallel computing environment one would choose as many processor cores as basis functions (). The overall runtime is then determined only by the initial value problem that takes the longest to integrate. For this numerical example it is . The communication overhead between processors is not taken into account since it is highly implementation and machine dependent. The slightly decreasing time to solution when is owed to the fact that initial values according to (25

) take more a-priori information into account which leads to smaller number of time steps and faster simulation. The overall accuracy of the method is problem-specific and always depends on both the tolerance for the solver and the number of basis functions. An a-priori determination of the number of basis functions and the solver tolerance is not yet available. An a-posteriori estimator can be constructed by increasing the number of basis functions and comparing the solutions. The resulting error is also related to the time stepping error.

The MPDE approach works also for nonlinear problems. However, similarly as for the harmonic balance case, the decoupling is not straightforward anymore. Furthermore the PWM basis functions and thus also the PWM eigenfunctions might not be able to represent the solution of problems with nonlinear elements . If the amplitude of the ripples is small compared to the amplitude of the envelope, the particular efficient approach described in  can be applied. It uses only the slowly varying envelope to evaluate the nonlinearities. Although the assembly of the field model matrices for a new envelope can not be parallelized, the matrices in (13) can still be decoupled and calculations to obtain the following time step can be run in parallel.

## 6 Conclusion

A new efficient technique was presented for field-circuit coupled models of DC-DC power converters, in which the switches are idealized and the filtering circuit is linear. The already existing MPDE technique with PWM basis functions splits the solution into fast varying and slowly varying parts. In this paper this method has been improved by introducing a new set of PMW basis functions which decouple the systems of equations similar as in the harmonic balance method. The new method, now called multirate PWM balance method, enables a parallel solution of all PWM modes resulting in a speed-up amounting to a factor 4 for the test example.

## 7 Acknowledgements

This work is supported by the “Excellence Initiative” of German Federal and State Governments and the Graduate School CE at TU Darmstadt. The authors thank Johan Gyselinck for fruitful discussions. Further thanks go to Jonas Bundschuh and Erik Skär for their contribution to the first implementation of the PWM eigenfunctions.

## References

•  N. Mohan, T. M. Undeland, and W. P. Robbins, Power electronics: converters, applications and design. Wiley, 3 ed., 2003.
•  J. Tant and J. Driesen, “On the numerical accuracy of electromagnetic transient simulation with power electronics,” IEEE Trans. Power Deliv., vol. 33, pp. 2492–2501, Oct. 2018.
•  A. Pels, J. Gyselinck, R. V. Sabariego, and S. Schöps, “Solving nonlinear circuits with pulsed excitation by multirate partial differential equations,” IEEE Trans. Magn., vol. 54, pp. 1–4, Mar. 2018.
•  A. Pels, J. Gyselinck, R. V. Sabariego, and S. Schöps, “Efficient simulation of DC-DC switch-mode power converters by multirate partial differential equations,” IEEE J. Multiscale Multiphys. Comput. Tech., vol. 4, no. 1, pp. 64–75, 2019.
•  H. G. Brachtendorf, G. Welsch, R. Laur, and A. Bunse-Gerstner, “Numerical steady state analysis of electronic circuits driven by multi-tone signals,” Electr. Eng., vol. 79, no. 2, pp. 103–112, 1996.
•  J. Roychowdhury, “Analyzing circuits with widely separated time scales using numerical PDE methods,” IEEE Trans. Circ. Syst. Fund. Theor. Appl., vol. 48, pp. 578–594, May 2001.
•  J. Gyselinck, C. Martis, and R. V. Sabariego, “Using dedicated time-domain basis functions for the simulation of pulse-width-modulation controlled devices – application to the steady-state regime of a buck converter,” in Electromotion 2013, (Cluj-Napoca, Romania), Oct. 2013.
•  S. Schöps, H. De Gersem, and T. Weiland, “Winding functions in transient magnetoquasistatic field-circuit coupled simulations,” COMPEL, vol. 32, pp. 2063–2083, Sept. 2013.
•  R. Lamour, R. März, and C. Tischendorf, Differential-Algebraic Equations: A Projector Based Analysis. Differential-Algebraic Equations Forum, Heidelberg: Springer, 2013.
•  R. Pulch, M. Günther, and S. Knorr, “Multirate partial differential algebraic equations for simulating radio frequency signals,” Eur. J. Appl. Math., vol. 18, no. 6, pp. 709–743, 2007.
•  K. Bittner and H. G. Brachtendorf, “Adaptive multi-rate wavelet method for circuit simulation,” Radioengineering, vol. 23, Apr. 2014.
•  C. R. I. Emson and C. W. Trowbridge, “Transient 3d eddy currents using modified magnetic vector potentials and magnetic scalar potentials,” IEEE Trans. Magn., vol. 24, pp. 86–89, Jan. 1988.
•  A. Nicolet and F. Delincé, “Implicit Runge-Kutta methods for transient magnetic field computation,” IEEE Trans. Magn., vol. 32, pp. 1405–1408, May 1996.
•  A. Bartel, S. Baumanns, and S. Schöps, “Structural analysis of electrical circuits including magnetoquasistatic devices,” APNUM, vol. 61, pp. 1257–1270, Sept. 2011.