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.


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 [1]. 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 [2]. 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 [7] 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 [7]. 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.

Multirate formulation
Let the field-circuit coupled model [8] of the converter be described by the system of ordinary differential or differential-algebraic equations where A ∈ R Ns×Ns is a possibly singular matrix, B ∈ R Ns×Ns is assumed to be a regular matrix, x(t) ∈ R Ns is the unknown solution, c(t) ∈ R Ns is the excitation, and t ∈ (0, T ] is the simulation interval. The initial state of the model is given by consistent [9] initial values is used as input of the power converter circuit. We denote by τ (t) = t Ts modulo 1 the relative time, T s is the switching cycle and D is the duty cycle. The system of Multirate Partial Differential(-Algebraic) Equations (MPDEs or MP-DAEs) with two time scales corresponding to (1) is given by [5,6,10] A where x(t 1 , t 2 ) is the unknown multivariate solution and c(t 1 , t 2 ) is the multivariate excitation. Choosing the multivariate excitation such that c(t, t) = c(t), the solution of (1) and (3) are related by x(t, t) = x(t). 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 x(0, t 2 ) = h(t 2 ), i.e., at t 1 = 0, where h with h(0) = x 0 is a function defining the initial values for all t 2 . The solution along the fast time scale t 2 is periodic, i.e., x(t 1 , t 2 + T s ) = x(t 1 , t 2 ). The multivariate right-hand side is chosen as c(t 1 , t 2 ) = c(t 2 ), 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.
To solve the MPDEs (3), a Galerkin approach and time discretization is applied [4,11]. The j-th solution component x j (t 1 , t 2 ) is approximated by expanding it into periodic basis functions p k depending on the fast time scale t 2 and coefficients w j,k depending on the slow time scale t 1 where the periodicity of the basis functions is accounted for by using the relative time τ (t 2 ) = t 2 Ts modulo 1. Applying the Galerkin approach with respect to t 2 and over one period of the excitation [0, T s ] leads to and right-hand side p denotes the complex conjugate of p and ⊗ denotes the Kronecker product.

PWM basis functions
The PWM basis functions developed in [7] are built up starting from the zero-th constant basis function p 0 (τ ) = 1 and the piecewise linear basis function which includes the duty cycle D of the excitation by construction. The higher-order basis functions p k (τ ), 2 ≤ k ≤ N p are recursively obtained by integrating the basis functions of lower order p k−1 (τ ) and orthonormalizing them using the Gram-Schmidt algorithm. The generated basis functions are depicted in Fig. 2.
For the PWM basis functions, the matrices J and Q from (6) are given by T s 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).

PWM eigenfunctions
To enable an easy parallelization of the method, the equations (5) can be decoupled, for example by diagonalizing Q, i.e., a basis transformation. We define new basis functions g k (τ ) as linear combinations of the PWM basis functions, i.e., where v k,l are unknown coefficients with k ∈ {0, . . . , N p }, and g k (τ ) are eigenfunctions of the time derivative operator We enforce this property in a weak sense by a Galerkin approach, i.e., where integration by parts and the periodicity of the basis functions is used. Inserting the expansion of the basis functions into (11) gives Since J is T s multiplied by the identity matrix (thanks to the orthonormality of the PWM basis functions), the λ k and v k are the eigenvalues and eigenvectors of the matrix Q, respectively. Furthermore since Q is real-valued and skew symmetric, and therefore a normal matrix, the eigenvectors v k are orthonormal. The new basis functions (complexvalued) are depicted in Fig. 3 for N p = 4. Note that the basis consists of pairs of conjugate complex basis functions. Inserting the transformed basis functions instead of the PWM basis functions into (6) and (7) leads, using the orthonormality of the eigenvectors, to where A as in (5), and Λ is a diagonal matrix with diagonal entries λ 0 , λ 1 , . . . , λ Np . Thus the resulting matrices in (13) are block-diagonal and the degrees of freedom can be block-wisely decoupled. This leads to N p + 1 independent systems of equations given by where w = [w 0 , w 1 , . . . , w Np ] 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".

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, σ fe = 250 S/m). The coils are modeled as stranded conductors. The simulation interval is given by Ψ = [0, 10] ms. The switching frequency is f s = 1 Ts = 1000 Hz. For the pulsed excitation (2) we use V 0 = 24 V. All calculations are performed in MATLAB. The partial differential equations governing the magnetoquasistatic problem are given by where r is the position vector, t is the time, A m is the modified magnetic vector potential [12], J s are the imposed currents, µ = 4π × 10 −7 H/m is the magnetic permeability and σ is the conductivity which is only non-zero in the ferrite core (σ fe ). The problem is considered on a 2D planar domain with homogeneous Dirichlet boundary conditions. Correspondingly, the Finite Element magnetoquasistatic [8] discretization of the magnetoquasistatic inductor model is given by the differential-algebraic system of equations [13] M σ d dt where M σ is the singular conductivity matrix, K is the stiffness matrix, a(t) gathers the degrees of freedom (DOFs) related to the magnetic vector potential, P is the discretization of the winding function [8] and i L (t) is the current through the inductor. The field-circuit coupling is expressed as follows. An additional variable is introduced for the magnetic flux linkage Φ(t) = P a(t). All equations are coupled monolithically into the index-1 differential-algebraic system of equations [14] P a − Φ = 0, which for the example in Fig. 1 contains a total of 11053 DOFs. The initial conditions are given by v C (0) = 0, i L (0) = 0 and a(0) = 0.

The initial condition for the MPDEs (3) can be written as
where h j is the j-th element of h. It only has to satisfy the condition h(0) = x(0, 0) = x 0 . Consequently there is a high degree of freedom in choosing the initial values w(0) 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., Secondly, the initial coefficients for k = 1, . . . , N p are extracted from the steady-state solution w j,k (0) = w s j,k for all j. The remaining coefficients are calculated by solving the solution expansion (4) for w j,0 (0) and using the condition x(0, 0) = x 0 . In summary the initial coefficients are given by 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.
The multivariate solution x(t 1 , t 2 ) 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 where E is the electric field strength, Ω is the spatial computation domain, the superscript H denotes the Hermitian, i.e., the complex conjugate transposed, and e(t) = − d dt a(t) 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 w(t 1 ), exemplary for the current through the inductor i L . As one can see, using the initial values (25), only the coefficient w j,0 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., abstol = reltol ∈ [10 −6 , 10 −1 ]; 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 abstol = reltol = 10 −7 and the number of basis functions N p ∈ [1, 10] is changed. The accuracy is measured for the voltage output of the converter, i.e., the voltage at the capacitor. The relative L 2 error is given by where v C,ref is the reference solution and v h C 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.
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 10 −6 in Fig. 6 for values larger than N p = 7 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 N p,pwmbal = 1, 2, 4, 6, 8, 10. 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 g k with k = 0, . . . , N p , a complexvalued 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 (N p + 1). The overall runtime is then determined only by the initial value problem that takes the longest to integrate. For this numerical example it is k = 0. 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 N p,pwmbal > 1 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 [4]. If the amplitude of the ripples is small compared to the amplitude of the envelope, the particular efficient approach described in [3] 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.

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.