I Introduction
Multirate behaviour can be observed in a number of technical applications. In highfrequency electrical circuit simulation, e.g., [1, 2, 3], the solution often consists of widely separated frequencies with slow and fast varying components. Furthermore, a division of the circuit into subcircuits whose state variables are either latent or active is often possible, especially in highly integrated circuits with many electrical elements [4]. The same holds for fieldcircuit coupled simulations describing the same physical phenomenon, e.g. an electrical circuit coupled to a magnetoquasistatic field model of an electrical machine. In coupled multiphysical simulations, different physical phenomena exhibit different characteristic time constants (e.g. electrothermal problems) and thus also lead to different rates of variation in the unknowns [5].
The solution of the abovementioned problems by conventional time discretization is inefficient as it enforces a step size to resolve the dynamics of the most active components of the system. Thus, the latent parts of the system are resolved with a much smaller time step size than necessary. This results in long simulation intervals and high computational effort. To efficiently solve these problems, various multirate methods have been developed.
Problems described by ordinary differential and differential algebraic equations (ODEs and DAEs) can be split into subsystems [4, 6, 7, 8]
, i.e., into several systems of equations describing the latent and active components, respectively. The subsystems are coupled, e.g., by extrapolation and interpolation of the state variables and resolved by different time step sizes and/or methods. A method of this kind for the simulation of power electronics has been proposed by Pekarek et al.
[9]. The coupling variables are synchronized in certain intervals. Between the synchronization time instants, the solution of the slow subsystem within the fast subsystem is calculated using a predictor and interpolation. Within the slow subsystem the solution of the fast subsystem is calculated by averaging.Another recent concept to deal with multirate phenomena is the reformulation of the ODEs or DAEs describing the problem into multirate partial differential equations (MPDEs) [1, 2]. The concept of the MPDEs allows to split the solution into components associated with different explicitly stated time scales . MPDEs have already been successfully applied in highfrequency circuit simulation using different solution approaches, e.g., multitone harmonic balance [1], multivariate finite difference time domain, hierarchical shooting [2] or combinations of different time stepping methods [3].
In this paper we focus on the efficient simulation of problems which are excited by periodic 2level pulsed (control) signals, as is often the case in power electronics, e.g., in DCDC switchmode power converters [10]. The system of differential equations describing the application is reformulated into a system of MPDEs. The MPDEs are solved by a combination of a RitzGalerkin approach and conventional adaptive time discretization. For the solution expansion of the unknowns, two types of basis functions are proposed: 1. the standard nodal FE basis functions; 2. the excitationspecific pulse width modulation (PWM) basis functions, introduced in [11]. The latter are designed to represent the ripple component in the output of switchmode power converters. These are, for the first time, interpreted in the framework of MPDEs and their approximation properties are mathematically analyzed. Since the concept of MDPEs allows to explicitly split the solution into components of different rates and solve them with different methods, it is possible to efficiently simulate not only the steadystate (as done in [11]), but also the transient behaviour of an application. The MPDE approach is validated in the example of a buck converter [11] in continuous conduction mode, as depicted in Fig. 4. Its solution, shown in Fig. 5, consists of a fast periodic ripple component and a slowly varying envelope.
The paper is structured as follows. Section II introduces the concept of MPDEs as described in the literature and establishes a link between the original system and the MPDEs. Section III is devoted to the solution of the MPDEs using a combination of a RitzGalerkin approach and conventional time discretization. In Section IV the two different types of basis functions for the solution expansion of the unknowns are presented. Finally in Section V the method is numerically validated on the simplified buck converter and convergence, accuracy and computational efficiency are analyzed. Section VI concludes the work by briefly summarizing the proposed approach and the main results.
Ii Introduction to Multirate Partial Differential Equations
In the following the proposed method is developed starting from a general linear circuit model. Let the vector of
unknown state variables consisting of node voltages and branch currents be given as(1) 
Modified nodal analysis [12] can be used to determine the system of firstorder linear DAEs or ODEs governing the circuit. This leads to the initial value problem (IVP)
(2)  
where are matrices, is the vector of excitations, is the vector of initial values and determines the simulation interval.
For , i.e., if is continuously differentiable, (2) can be written equivalently as system of MPDEs [1, 2], introducing different time scales
(3)  
where , are the multivariate forms of , , respectively, and determine the simulation domains. The vector of state variables is given by
(4) 
In the following, a relation between the solution and excitation of (2) and (3) is established, which was first introduced by Brachtendorf et al. [1]. They developed a socalled multitone harmonic balance method using MPDEs to efficiently simulate highfrequency circuits with more than one fundamental frequency.
Let be a solution of the MPDEs (3) and the corresponding excitation. Then the solution and excitation of the DAEs or ODEs (2) and MPDEs (3) are related by and , respectively, for any fixed [1, 13].
To proof this statement, the chain rule of differentiation is applied to (
2), which yields [1, 13](5)  
Thus if a solution of the MPDEs (3) can be found for a multivariate righthand side fulfilling , the solution of the DAEs or ODEs (2) can be extracted from the multivariate solution using . To solve the MPDEs (3), initial and boundary conditions have to be imposed. As only IVPs are considered, whose solution can be separated into periodic and nonperiodic parts, the setting of envelopemodulated solutions [2] is appropriate. We therefore define initial and boundary conditions to the MPDEs (3) as
(6)  
where is a function specifying the initial conditions and are time intervals of periodicity.
For the sake of simplicity, hereafter, we restrict to two time scales () leading to the mixed initial boundary value problem
(7)  
In the following section, we will apply the MPDE framework to problems with discontinuous righthand sides. Existence and uniqueness can still be assured by the Carathéodory conditions. The interested reader is referred to, e.g., [14]. However, in the case of the pulsed excitations introduced in the next section, a piecewise analysis is possible and a detailed discussion is not needed.
Iii Solution of the Multirate Partial Differential Equations
In this section we focus on the solution of the MPDEs for applications with PWM (pulsed) excitation. The switching cycle and duty cycle are assumed to be constant. We propose the following procedure for solving: 1.) a RitzGalerkin approach is applied to one dimension of the MPDEs (7); 2.) the remaining linear system of ODEs or DAEs is solved with conventional time discretization.
Iiia Solution expansion by basis functions
In a first step the multivariate solution is expanded into a finite set of basis functions and coefficients. It reads
(8) 
where , is the th approximated state variable, , are periodic basis functions and are coefficients. The superscript in denotes that it is an approximation to . By defining the expansion as above we associate the slowly varying envelope with the time scale , which will be therefore referred to as the slow time scale, and the fast periodically varying ripples with the time scale , which will be referred to as fast time scale. The basis functions are periodic with switching cycle , which can be accounted for by introducing the relative time
(9) 
The switching cycle is related to the switching frequency by . For simplicity the basis functions will in the following be expressed as functions of the relative time .
IiiB Galerkin approach
The RitzGalerkin approach is applied to the MPDEs with respect to the fast time scale in the interval
(16)  
i.e., the MPDEs are weighted by the same basis functions used for the solution expansion.
Let the matrices and be given as
(17) 
Inserting the relation (15) into (16) leads to
(18) 
where
(19) 
is the unknown vector of coefficients and and are, using the Kronecker product, given by
(20)  
(21)  
(22) 
Note that the sparsity pattern of the matrices and depends on the choice of the basis functions.
IiiC Time discretization
The equations (18) are now formulated only in . According to (22) their righthand side naturally depends on the righthand side of the MPDEs, which only needs to satisfy the relation according to Section II. As a result infinitely many choices for are possible. However to minimize the dynamic of the system (18) and thus maximizing the efficiency of the approach, it is reasonable to head for a constant righthand side. As is periodic with switching cycle for the considered problems, we choose . Inserting this into (22), the time scales and vanish which leads to
(23) 
Note that the system (18) is times larger than the original one (2).
Iv Choice of basis functions for solution expansion
We propose the use of standard FE nodal functions as in classical finite element methods (FEM) or the PWM basis functions introduced in [11].
Iva Finite element nodal basis functions
To start with the FE nodal functions of first order, let us introduce a division of the relative time interval into elements, such that , where the are the nodes defining the elements. The nodal basis functions are piecewise linear functions defined as
(24) 
To enforce periodicity on the interval we set the basis functions at the boundary, i.e., and , to zero
(25)  
(26) 
To resolve the envelope, we introduce an additional constant basis function
(27) 
The FE nodal basis as defined above is depicted in Fig. 1. As the FE nodal functions offer local support, except , the matrices are sparsely populated matrices. Due to the constant basis function supporting the entire relative time interval , the matrices are not purely banded matrices as in classical FE methods.
Instead of setting the boundary functions to zero and defining an additional constant basis function, it is also possible to enforce periodic boundary conditions on the set of standard FE nodal functions in the final system of equations.
IvB PWM basis functions
A problemspecific choice of basis functions in case of apriori known duty cycle is the PWM basis functions, which were developed in [11]. The apriori knowledge enables us to build the basis functions such that they mimic the shape of the ripple components in the solution by construction. The zeroth basis function is which resolves the envelope as in the case of nodal basis functions. The PWM basis is iteratively built starting from the normalized, zero average, piecewise linear basis function defined as [11]
(28) 
The higherorder basis functions , are obtained recursively by integrating the basis functions of lower order ensuring continuity
(29) 
This extended set of basis functions is successively orthonormalized, starting from , by orthogonalizing
(30) 
and normalizing
(31) 
which corresponds to a GramSchmidt orthonormalization [15]. Note that it is possible to calculate the PWM basis functions analytically. The basis functions of order up to 3 are depicted in Fig. 2.
Opposed to the FE nodal functions, the PWM basis functions are global polynomials on the relative time interval
as in spectral methods and offer the same accuracy with less degrees of freedom compared to the nodal basis functions
[16]. Due to the orthonormality of the PWM basis functions the matrixis the identity matrix. The matrix
is dense, however only 25% are nonzero elements.The approximation properties of these specific basis functions have been studied up to now only numerically [11]. The basis functions are by construction restricted to represent piecewise exponential solutions. Their properties are studied analytically in the Appendix for a duty cycle of .
V Numerical results
In this section the method is numerically validated. Computational efficiency, accuracy and convergence results are presented. All calculations have been performed in GNU Octave [17]. For solving equation (18
), an implicit RungeKutta method of order 5 with 6 stages is used. For step size prediction the estimated error is measured in the infinity norm instead of the 2norm as originally proposed in
[18], p. 124, i.e.,(32) 
where is the dimension of the equation system and is the estimated error of the th solution component in each step. The quantity depends on the relative and absolute tolerance. For more information the reader is referred to [18]. The absolute tolerance is fixed at so that the error estimation is controlled by the relative tolerance . The solver supports dense output which is used in reconstructing the MPDE solution.
Va Test case
The test case is a buck converter circuit [11] as depicted in Fig. 3. The buck converter consists of a DC voltage source , a switch (e.g. an IGBT), a diode, an inductor (consisting of inductance and resistance ) and a capacitor (capacitance ). At the output a load resistance is connected. The switch is controlled by a 2level pulsed signal, which closes and opens the switch at switching frequency and with a duty cycle .
Assuming continuous conduction mode (), an ideal switch, and an ideal diode, the buck converter can be simplified as depicted in Fig. 4[10]. The switch and diode have been removed and the voltage source has been replaced by a pulsed voltage source , which output voltage alternates between and , i.e.
(33) 
The circuit can be described by two state variables, namely the current through the coil and the voltage across the capacitor , which is also the output voltage of the buck converter. Using Kirchhoff’s circuit laws leads to the firstorder linear ODEs
(34) 
where the following parameter values are chosen:

V;

Hz;

;

mH, m;

F;

.
The initial conditions are set to and .
As reference solution a closedform analytic solution of the buck converter ODE (34) is used. Figs. 5 and 6 show the voltage at the capacitor and current through the coil of the buck converter for Hz and Hz, respectively. The solution consists of a slowly varying envelope and ripple components which are periodic. Increasing the switching frequency, the magnitude of the ripples decreases.
VB Multirate solution
To obtain the multirate solution of the buck converter circuit the MPDE approach as described before is applied. For the FE nodal functions equidistant spacing between the nodes dividing the relative time interval into elements is used. The number of basis functions is always chosen such that the jump of the excitation as defined in (33) occurs at a time instant which coincides with a node. Thereby the continuity in the solution coincides exactly with a node and is properly represented. For a duty cycle of this corresponds to . For the PWM basis functions no special care is needed to choose as they take the duty cycle into account by construction.
The equation system (18) is solved for the vector of coefficients . To find the initial values , the steadystate solution of the system is calculated
(35) 
The coefficients corresponding to the constant basis function are set as such that the solution satisfies the initial condition and .
The coefficients for 12 basis functions (11 FE nodal functions + 1 constant function), i.e., , are exemplary depicted in Fig. 7 for the capacitor voltage after solving. All coefficients except stay constant during the simulation time, i.e. the coefficients controlling the shape of the ripples do not change.
The multivariate solution is reconstructed using the solution expansion (8). As the solver often uses less time steps than for the original equations (34), it is taken advantage of dense output to extract a reasonably fine sampled solution. Fig. 8 shows the result in a 3D plot. Along the time axis the slow dynamic resolved by time discretization can be observed while along the time axis the high dynamic resolved by the Galerkin approach is visible. The solution of the original equations (34) is marked as black line and can be extracted using according to Section II.
VC Convergence
To compare the two types of basis functions used for the solution expansion, the convergence of the solution with respect to and the tolerance of the solver is examined. We consider the simulation time interval ms. As reference solution for the buck converter, a closedform analytic solution is calculated. The following analysis is restricted to the output voltage of the buck converter, i.e., the voltage at the capacitor. The convergence behaviour of the current through the inductor is similar. Let and define the relative error of the solution by
(36) 
where is the voltage at the capacitor calculated by the MPDE approach for different relative tolerance, number of basis functions and time instants and is the respective reference solution. The norm is approximated by numerical quadrature using the midpoint rule. In the following will simply be referred to as error.
The error is evaluated at a fixed number of 500 samples per period . For this, again, the dense output feature of the solver is used.
Fig. 9 shows the convergence of the error using nodal basis functions (nodal BFs) with refinement and the PWM basis functions (PWM BFs) with refinement for a fixed relative tolerance of for the time stepper. It’s tolerance determines a limit for the accuracy of the solution. To ensure that the employed tolerance is small enough, the error is compared for and . The absolute difference between the errors for a maximum number of basis functions, for the PWM basis, and for the FE basis, is several orders of magnitude smaller than the obtained error . A relative tolerance of is therefore adequate for all calculations.
According to Fig. 9, the solution using PWM basis functions converges significantly faster than the solution using nodal basis functions with respect to the number of basis functions .
VD Computational efficiency
To validate the efficiency of the method, the MPDE approach with nodal and PWM basis functions is compared to classical time discretization of the original ODEs (34). For the time discretization the accuracy is controlled by varying the relative tolerance of the solver while for the MPDE approach we fix the relative tolerance (at ), and vary the number of basis functions to achieve a certain accuracy. Fig. 10 shows that the efficiency in terms of time for solving the differential equation systems of the MPDE approach depends on the choice and number of basis functions. While the PWM basis functions yield excellent efficiency, the MPDE approach using nodal functions becomes inferior than time discretization for about . To better understand this effect, two additional quantities are examined.
Fig. 11 shows the error versus number of function evaluations. For the time discretization this number increases to reach higher accuracy as more time steps are necessary. For the MPDE approach, due to the slow dynamics of the equation system (18), much less time steps and thus less function evaluations are needed. For higher accuracy (i.e. increasing ) the number of function evaluations even decreases. This effect results from adding additional basis functions by which there is more apriori information on the solution already taken into account. The envelope stored in the zeroth coefficient including its initial value to ensure the initial conditions of the buck converter therefore varies with different and the ODE solver needs less time steps and thus less function evaluations. In Fig. 12 the error versus the average time per function evaluation is shown. In this plot the effect of larger equation system in the MPDE approach becomes visible. The average time increases dramatically for the MPDE approach with nodal basis functions as the number of basis functions is large while for the PWM basis functions the effect is much smaller due to smaller . For conventional time discretization the average time per function evaluation is constant as the size of the equation system does not change and therefore the computational effort per step is constant. The effects visible in Figs. 11 and 12 determine the overall efficiency depicted in Fig. 10. In conclusion, for the FE nodal functions this means that the effect of increasing size of equation systems and therefore more effort per step begins to outweigh the advantage of less required time steps for and larger.
The reconstruction of the solution using the solution expansion (8) is not taken into account in the above efficiency measurements. The time for evaluation depends mainly on the number of samples at which the solution is reconstructed. If the number of samples per period is known, the evaluation of the basis functions can be done apriori. As a result, the reconstruction of the solution is cheap. In the case of 500 samples per period it takes considerably less than ms and can therefore be neglected.
Note that the speedup of the MPDE approach compared to time discretization can be expected to increase if larger time intervals are considered or higher switching frequencies are used. The higher the frequency, the more ripples have to be resolved. Time discretization therefore needs more and more time steps in the same time interval while for the MPDE approach the number of time steps do not change as the periodically varying ripples are resolved by the Galerkin approach. The same happens for increasing time intervals and fixed switching frequency.
Vi Conclusion
An efficient and accurate approach to simulate PWM driven applications with constant switching and duty cycle has been presented. The linear circuit model of the power converter is first reformulated as MPDEs. A Galerkin approach and time discretization are used to solve the MPDEs. By this the fast periodically varying components of the solution are taken into account by the basis functions and the time discretization only resolves the dynamics of the envelope. This leads to a reduced number of time steps. For the solution expansion two types of basis functions have been proposed, namely FE nodal functions and PWM basis functions. The MPDE approach has been validated on the example of a simplified buck converter. The convergence of the solution in terms of solver tolerance and number of basis functions has been examined. The solution using PWM basis functions converges much faster than when using FE nodal functions. The computational efficiency of the method strongly depends on the choice and number of basis functions. By using the Galerkin approach, the size of the resulting equation system is determined by how many basis functions are used. To solve the final equation system a much smaller number of time steps is necessary however with the drawback of more time spent in each step due to the larger equation systems. A tradeoff between accuracy and speedup is therefore necessary. This becomes particularly visible for the FE nodal functions. When a large number of basis functions are used, the drawback of the approach begins to outweigh the advantage which leads to inefficient simulation. For small number of basis functions the MPDE approach is highly efficient on the presented example of the buck converter.
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 Herbert Egger from the department of mathematics and Felix Wolf from the Graduate School of Computational Engineering at TU Darmstadt for fruitful discussions.
References
 [1] H. G. Brachtendorf, G. Welsch, R. Laur, and A. BunseGerstner, “Numerical steady state analysis of electronic circuits driven by multitone signals,” Electrical Engineering (Archiv für Elektrotechnik), vol. 79, no. 2, pp. 103–112, 1996.
 [2] J. Roychowdhury, “Analyzing circuits with widely separated time scales using numerical PDE methods,” IEEE Transactions on Circuits and Systems Part I: Fundamental Theory and Applications, vol. 48, pp. 578–594, May 2001.
 [3] T. Mei, J. Roychowdhury, T. Coffey, S. Hutchinson, and D. Day, “Robust, stable timedomain methods for solving MPDEs of fast/slow systems,” IEEE Transactions on ComputerAided Design of Integrated Circuits and Systems, vol. 24, pp. 226–239, Feb. 2005.
 [4] M. Günther and P. Rentrop, “Partitioning and multirate strategies in latent electric circuits,” report tum m9301 (93), Technische Universität, München, Jan. 1993.
 [5] S. Schöps, H. De Gersem, and A. Bartel, “A cosimulation framework for multirate timeintegration of field/circuit coupled problems,” IEEE Transactions on Magnetics, vol. 46, pp. 3233–3236, July 2010.
 [6] M. Striebel, A. Bartel, and M. Günther, “A multirate ROWscheme for index1 network equations,” Applied Numerical Mathematics, vol. 59, no. 3, pp. 800–814, 2009.
 [7] C. W. Gear and D. R. Wells, “Multirate linear multistep methods,” BIT Numerical Mathematics, vol. 24, no. 4, pp. 484–502, 1984.
 [8] A. Kværnø and P. Rentrop, “Low order multirate rungekutta methods in electric circuit simulation,” 1999.
 [9] S. D. Pekarek, O. Wasynczuk, E. A. Walters, J. V. Jatskevich, C. E. Lucas, N. Wu, and P. T. Lamm, “An efficient multirate simulation technique for powerelectronicbased systems,” IEEE Transactions on Power Systems, vol. 19, pp. 399–409, Feb. 2004.
 [10] N. Mohan, T. M. Undeland, and W. P. Robbins, Power electronics: converters, applications and design. Wiley, 3 ed., 2003.
 [11] J. Gyselinck, C. Martis, and R. V. Sabariego, “Using dedicated timedomain basis functions for the simulation of pulsewidthmodulation controlled devices – application to the steadystate regime of a buck converter,” in Electromotion 2013, (ClujNapoca, Romania), Oct. 2013.
 [12] C.W. Ho, A. E. Ruehli, and P. A. Brennan, “The modified nodal approach to network analysis,” IEEE Transactions on Circuits and Systems, vol. 22, pp. 504–509, June 1975.
 [13] S. Knorr, WaveletBased Simulation of Multirate Partial DifferentialAlgebraic Systems in Radio Frequency Applications. Dissertation, Düsseldorf, 2007. VDI Verlag. FortschrittBerichte VDI, Reihe 20.
 [14] A. F. Flippov, Differential Equations with Discontinuous Righthand Sides. Mathematics and Its Applications, Springer, 1988.
 [15] L. N. Trefethen and D. Bau, Numerical Linear Algebra. Society for Industrial and Applied Mathematics, 1997.
 [16] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods: Fundamentals in Single Domains. Scientific Computation, Springer Berlin Heidelberg, 2006.
 [17] J. W. Eaton, D. Bateman, S. Hauberg, and R. Wehbring, The GNU Octave 4.0 Reference Manual 1/2: Free Your Numbers. Samurai Media Limited, Oct. 2015.

[18]
E. Hairer, S. P. Nørsett, and G. Wanner,
Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems
. Springer Series in Computational Mathematics, Berlin: Springer, 2 ed., 2002.
Appendix
Proof.
The basis functions with duty cycle are defined as follows: The zeroth and first basis function are given piecewisely as
(37) 
and
(38) 
which essentially corresponds to a scaled and translated hat function. The subscript letter refers to the interval in which the basis function is defined, i.e., if “”, the polynomial in is considered, if “” the polynomial in is considered. If the subscript comprises only a number, the entire basis function is addressed.
The symmetry of the basis function can be expressed as follows
(39) 
The basis functions of higher order, i.e., are calculated by integrating the basis functions of lower order
(40)  
(41) 
and orthogonalizing the integrated basis functions against the constant basis function . The basis function therefore is
(42)  
Similarly is given as
(43)  
Note that orthonormalization of all basis functions against each other is possible and has been originally proposed in [11]. However, it spans the same space as the basis functions without full orthonormalization (42),(43), and is therefore neglected.
Via Symmetry properties of PWM basis functions
The symmetry properties of the basis functions as defined by (38)(43) with duty cycle are examined in the following.
ViA1 Induction hypothesis
The symmetry of the basis functions is given by
(44) 
i.e., for all basis functions with even index, and
(45) 
i.e., for all basis functions with odd index.
ViA2 Induction base
ViA3 Induction step
We calculate the basis functions and , where and their symmetry properties.
Basis function .
The basis function is given by integration and orthogonalization against the constant basis function , i.e.,
(50) 
and
(51) 
The orthogonalization term yields using the symmetry properties (45) and substitution
(52)  
Therefore the expressions for the basis function (50) and (51) simplify to
(53) 
and
(54) 
The symmetries are given by, using (45),
(55) 
and
(56) 
and thus fulfill the hypothesis (44).
Basis function .
The basis function is calculated by
(57) 
and
Comments
There are no comments yet.