Efficient simulation of DC-AC power converters using Multirate Partial Differential Equations

TL;DR

This paper introduces a multirate MPDE-based simulation method for DC-AC power converters that efficiently handles PWM excitation with varying duty cycles, reducing computational cost and eliminating switch event detection.

Contribution

It develops a novel MPDE approach for simulating switch-mode power converters with varying duty cycles, improving efficiency and accuracy over traditional methods.

Findings

Significantly reduces simulation time compared to conventional methods.

Eliminates the need for switch event detection in simulations.

Effectively handles PWM excitation with varying duty cycles.

Abstract

Switch-mode power converters are used in various applications to convert between different voltage (or current) levels. They use transistors to switch on and off the input voltage to generate a pulsed voltage whose arithmetic average is the desired output voltage of the converter. After smoothening by filters, the converter output is used to supply devices. The simulation of these switch-mode power converters by conventional time discretization is computationally expensive since a high number of time steps is necessary to properly resolve the unknown state variables and detect switch events of the excitation. This paper proposes a multirate method based on the concept of Multirate Partial Differential Equations (MPDEs), which splits the solution into fast varying and slowly varying parts. The method is developed to work with pulse width modulated (PWM) excitation with a constant switching cycle and varying duty cycle. The important case of varying duty cycles in the MPDE framework is addressed for the first time. Switching event detection is no longer necessary and a much smaller number of time steps for a decent resolution are required, thus leading to a highly efficient method.

Tables

Table 1. Table 1 : Accuracy of the MPDE approach and the conventional simulation using time discretization in MATLAB and PLECS.

	MPDE		MATLAB		PLECS
Simulation setting	Tol.	Error voltage/current	Tol.	Error voltage/current	Tol.	Error voltage/current
lowest order	${10}^{-3}$	$8.5\cdot {10}^{-3}$ / $2.5\cdot {10}^{-3}$	${10}^{-2}$	$1.2\cdot {10}^{-3}$ / $2.5\cdot {10}^{-3}$	${10}^{-2}$	$6.5\cdot {10}^{-3}$ / $8.1\cdot {10}^{-3}$
medium order	${10}^{-4}$	$6.9\cdot {10}^{-4}$ / $6.2\cdot {10}^{-4}$	${10}^{-4}$	$1.6\cdot {10}^{-4}$ / $1.6\cdot {10}^{-4}$	${10}^{-4}$	$2.3\cdot {10}^{-4}$ / $3.8\cdot {10}^{-4}$
high order	${10}^{-7}$	$8.1\cdot {10}^{-7}$ / $1.1\cdot {10}^{-6}$	${10}^{-6}$	$3.2\cdot {10}^{-6}$ / $2.7\cdot {10}^{-6}$	${10}^{-6}$	$1.6\cdot {10}^{-6}$ / $3.5\cdot {10}^{-6}$

Table 2. Table 2 : Speedup of the MPDE approach compared to the conventional simulation using MATLAB and PLECS in different simulation settings.

	Setting 1 (lowest order)			Setting 2 (medium order)			Setting 3 (high order)
	MPDE	MATLAB/ PLECS	Speedup (approx.)	MPDE	MATLAB/ PLECS	Speedup (approx.)	MPDE	MATLAB/ PLECS	Speedup (approx.)
time steps	235	10533/4243	44/18	474	12410/9940	26/21	4516	18476/17508	4/4
failed steps	50	36	0.7	76	805	11	111	2641	24
LU decom.	95	2535	27	157	4369	28	728	7435	10
solution lin. systems	557	21938	39	1062	27230	26	7553	43034	6

Equations

$$\mathbf{A}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{x}(t)+\mathbf{Bx}(t)=\mathbf{c}(t),\mathbf{x}(0)={\mathbf{x}}_0,$$

$$\mathbf{A}\left(\frac{\partial \hat{\mathbf{x}}(t_1,t_2)}{\partial t_1}+\frac{\partial \hat{\mathbf{x}}(t_1,t_2)}{\partial t_2}\right)+\mathbf{B}\hat{\mathbf{x}}(t_1,t_2)=\hat{\mathbf{c}}(t_1,t_2).$$

$$\hat{\mathbf{x}}(t,t)$$

$$\hat{\mathbf{c}}(t,t)$$

$${\hat{x}}_j(t_1,t_2)\doteq \sum _{k=0}^{N_{\mathrm{p}}}{w}_{j,k}(t_1)p_k(\tau (t_2),d(t_1))={\mathbf{p}}^{\top }(\tau (t_2),d(t_1)){\mathbf{w}}_j(t_1),$$

$$\mathbf{p}(\tau)=\left[\begin{array}[]{c}p_{0}\\ p_{1}(\tau)\\ p_{2}(\tau)\\ \vdots\\ p_{N_{\mathrm{p}}}(\tau)\end{array}\right]\ ,\ \ \quad\mathbf{w}_{j}(t_{1})\ =\ \left[\begin{array}[]{c}w_{j,0}(t_{1})\\ w_{j,1}(t_{1})\\ w_{j,2}(t_{1})\\ \vdots\\ w_{j,N_{\mathrm{p}}}(t_{1})\end{array}\right]\,.\vspace{1em}$$

$$\frac{\partial {\hat{x}}_j(t_1,t_2)}{\partial t_1}+\frac{\partial {\hat{x}}_j(t_1,t_2)}{\partial t_2}$$

$$+{\mathbf{p}}^{\top }(\tau (t_2),d(t_1))\frac{\mathrm{d}{\mathbf{w}}_j(t_1)}{\mathrm{d}{t}_1}+\frac{\partial {\mathbf{p}}^{\top }(\tau (t_2),d(t_1))}{\partial \tau }\frac{\mathrm{d}\tau (t_2)}{\mathrm{d}{t}_2}{\mathbf{w}}_j(t_1),$$

$$\underset{0}{\overset{T_{\mathrm{s}}}{\int }}\left(\mathbf{A}\left(\frac{\partial \hat{\mathbf{x}}(t_1,t_2)}{\partial t_1}+\frac{\partial \hat{\mathbf{x}}(t_1,t_2)}{\partial t_2}\right)+\mathbf{B}\hat{\mathbf{x}}(t_1,t_2)\right)p_l(\tau (t_2),d)\mathrm{d}{t}_2$$

$$=\underset{0}{\overset{T_{\mathrm{s}}}{\int }}\hat{\mathbf{c}}(t_1,t_2)p_l(\tau (t_2),d)\mathrm{d}{t}_2\forall l=0,\dots ,N_{\mathrm{p}}.$$

$$\underset{0}{\overset{T_{\mathrm{s}}}{\int }}\mathbf{A}\frac{\partial \hat{\mathbf{x}}(t_1,t_2)}{\partial t_1}$$

$$\displaystyle+\Big{(}\mathbf{\widehat{x}}(t_{1},t_{2})p_{l}\Big{)}\Big{|}_{t_{2}=0}^{t_{2}=T_{\mathrm{s}}}=\int\limits_{0}^{T_{\mathrm{s}}}\mathbf{\widehat{c}}(t_{1},t_{2})\,p_{l}\,{\mathrm{d}}t_{2}\qquad\forall l=0,\dots,N_{\mathrm{p}}.$$

$$\displaystyle\int\limits_{0}^{1}\left(\mathbf{A}\left[\begin{array}[]{c}\frac{\partial\mathbf{p}^{\top}}{\partial\,{d}}\frac{{\mathrm{d}}\,{d}}{{\mathrm{d}}t_{1}}\mathbf{w}_{1}\\ \vdots\\ \frac{\partial\mathbf{p}^{\top}}{\partial\,{d}}\frac{{\mathrm{d}}\,{d}}{{\mathrm{d}}t_{1}}\mathbf{w}_{N_{\mathrm{s}}}\end{array}\right]+\mathbf{A}\left[\begin{array}[]{c}\mathbf{p}^{\top}\frac{{\mathrm{d}}\mathbf{w}_{1}}{{\mathrm{d}}t_{1}}\\ \vdots\\ \mathbf{p}^{\top}\frac{{\mathrm{d}}\mathbf{w}_{N_{\mathrm{s}}}}{{\mathrm{d}}t_{1}}\end{array}\right]+\mathbf{B}\left[\begin{array}[]{c}\mathbf{p}^{\top}\mathbf{w}_{1}\\ \vdots\\ \mathbf{p}^{\top}\mathbf{w}_{N_{\mathrm{s}}}\end{array}\right]\right)p_{l}\,T_{\mathrm{s}}{\mathrm{d}}\tau$$

$$\displaystyle-\int\limits_{0}^{1}\left(\mathbf{A}\left[\begin{array}[]{c}\mathbf{p}^{\top}\mathbf{w}_{1}\\ \vdots\\ \mathbf{p}^{\top}\mathbf{w}_{N_{\mathrm{s}}}\end{array}\right]\right)\frac{\partial p_{l}}{\partial\tau}\underbrace{\frac{1}{T_{\mathrm{s}}}\,T_{\mathrm{s}}}_{=1}{\mathrm{d}}\tau=\int\limits_{0}^{T_{\mathrm{s}}}\mathbf{\widehat{c}}(t_{1},t_{2})p_{l}\,{\mathrm{d}}t_{2}\qquad\forall l=0,\dots,N_{\mathrm{p}}.$$

$$\text{cA}(t_1)\frac{\mathrm{d}\mathbf{w}}{\mathrm{d}{t}_1}+\text{cB}(t_1)\mathbf{w}(t_1)=\text{cC}(t_1),$$

$$\text{cA}(t_1)$$

$$\text{cC}(t_1)$$

$$\displaystyle\cI({d})\$$

$$\displaystyle\cQ({d})\$$

$$\displaystyle\cU({d})\$$

$${P}_{i,0}(\xi)=\left\{\begin{array}[]{cc}1&\text{for }\xi_{i}\leq\xi<\xi_{i+1}\\ 0&\text{otherwise}\end{array}\right.$$

$$P_{i,\bar{p}}(\xi )=\frac{\xi -{\xi }_i}{{\xi }_{i+\bar{p}}-{\xi }_i}{P}_{i,\bar{p}-1}(\xi )+\frac{{\xi }_{i+\bar{p}+1}-\xi }{{\xi }_{i+\bar{p}+1}-{\xi }_{i+1}}{P}_{i+1,\bar{p}-1}(\xi ).$$

$$\Xi =\{\underset{\bar{p}+1}{\underbrace{a,\dots ,a}},{\xi }_{\bar{p}+1},\dots ,{\xi }_{m-\bar{p}-1},\underset{\bar{p}+1}{\underbrace{b,\dots ,b}}\},$$

$${\Xi }_{\bar{p}}=\{\underset{\bar{p}+1}{\underbrace{0,\dots ,0}},\underset{\bar{p}}{\underbrace{d,\dots ,d}},\underset{\bar{p}+1}{\underbrace{1,\dots ,1}}\}.$$

$${\Xi }_{\bar{p},K}=\{\underset{\bar{p}+1}{\underbrace{0,\dots ,0}},{\alpha }_{1}d,\dots ,{\alpha }_{K}d,\underset{\bar{p}}{\underbrace{d,\dots ,d}},{\beta }_1(1-d)+d,\dots ,{\beta }_K(1-d)+d,\underset{\bar{p}+1}{\underbrace{1,\dots ,1}}\},$$

$$\{P_{0,\bar{p}}(\xi ,d),\dots ,P_{2\bar{p}+2K,\bar{p}}(\xi ,d)\}$$

$$p_k(\tau ,d)=P_{k,\bar{p}}(\tau ,d)$$

$$P_{i,\bar{p}}(\xi ,d)=\frac{\xi -{\gamma }_{1,i}d}{{\gamma }_{2,i}d}{P}_{i,\bar{p}-1}(\xi ,d)+\frac{{\gamma }_{3,i}d-\xi }{{\gamma }_{4,i}d}{P}_{i+1,\bar{p}-1}(\xi ,d),$$

$$P_{i,\bar{p}}(\xi ,d)={\bar{\gamma }}_{\bar{p},i}\frac{{\xi }^{\bar{p}}}{d^{\bar{p}}}+{\bar{\gamma }}_{\bar{p}-1,i}\frac{{\xi }^{\bar{p}-1}}{d^{\bar{p}-1}}+\cdots +{\bar{\gamma }}_{0,i}\frac{{\xi }^0}{d^0}$$

$${\hat{\Xi }}_{\bar{p},K}={\Xi }_{\bar{p},K}-d=\{-d,\dots ,d,{\alpha }_{1}d-d,\dots ,{\alpha }_{K}d-d,0,\dots ,0,{\beta }_1(1-d),\dots ,{\beta }_K(1-d),d,\dots ,d\}.$$

$${\hat{P}}_{i,\bar{p}}(\xi ,d)=\frac{\xi -{\chi }_{1,i}(1-d)}{{\chi }_{2,i}(1-d)}{\hat{P}}_{i,\bar{p}-1}(\xi ,d)+\frac{{\chi }_{3,i}(1-d)-\xi }{{\chi }_{4,i}(1-d)}{\hat{P}}_{i+1,\bar{p}-1}(\xi ,d),$$

$${\hat{P}}_{i,\bar{p}}(\xi ,d)={\bar{\chi }}_{\bar{p},i}\frac{{\xi }^{\bar{p}}}{(1-d{)}^{\bar{p}}}+{\bar{\chi }}_{\bar{p}-1,i}\frac{{\xi }^{\bar{p}-1}}{(1-d{)}^{\bar{p}-1}}+\cdots +{\bar{\chi }}_{0,i}\frac{{\xi }^0}{(1-d{)}^0}$$

$${\text{cI}}_{i,j}(d)=\underset{{\text{cI}}_{i,j}^{0}d}{\underbrace{T_{\mathrm{s}}\underset{0}{\overset{d}{\int }}{P}_{i,\bar{p}}(\xi ,d)P_{j,\bar{p}}(\xi ,d)\mathrm{d}\xi }}+\underset{{\text{cI}}_{i,j}^1(1-d)}{\underbrace{T_{\mathrm{s}}\underset{0}{\overset{1-d}{\int }}{\hat{P}}_{i,\bar{p}}(\xi ,d){\hat{P}}_{j,\bar{p}}(\xi ,d)\mathrm{d}\xi }}$$

$${\text{cQ}}_{i,j}=\underset{{\text{cQ}}_{i,j}^0}{\underbrace{\underset{0}{\overset{d}{\int }}\frac{\partial P_{i,\bar{p}}(\xi ,d)}{\partial \xi }{P}_{j,\bar{p}}(\xi ,d)\mathrm{d}\xi }}+\underset{{\text{cQ}}_{i,j}^1}{\underbrace{\underset{0}{\overset{1-d}{\int }}\frac{\partial {\hat{P}}_{i,\bar{p}}(\xi ,d)}{\partial \xi }{\hat{P}}_{j,\bar{p}}(\xi ,d)\mathrm{d}\xi }}$$

Full text

1 Introduction

Switch-mode power converters play nowadays a vital role in various applications [1, 2]. They are used in domestic devices, e.g., in mobile phone chargers and computer power supply, as well as in industrial applications, e.g., in high-voltage DC (HVDC) power transfer and speed control of electrical motors, to convert voltage (or current) between different levels. Switch-mode power converters use transistors to periodically switch on and off the input voltage to generate a pulsed voltage, whose arithmetic average is the desired output of the converter. The pulsed voltage is smoothened by filter circuits, which also act as energy buffer before it is used to supply the appliance.

An exemplary circuit of such a switch-mode power converter is depicted in Fig. 1 along with its solution in DC-DC conversion mode. As can be seen the converter output voltage, even after filtering, still comprises of high-frequency components generated by the pulsed excitation, which are usually referred to as ripples [2]. Since the circuit is uncharged at the beginning, the converter needs some time to reach the steady state. As a result, the solution consists of a slowly varying envelope which is modulated with the fast periodically varying ripples. The method used to generate the control signals for the transistors is called pulse width modulation (PWM). There are different kinds of PWM [2], e.g., natural sampling with different carrier signals like sawtooth and triangle, and regular sampling. Since in this paper, we focus on natural sampling with a sawtooth carrier, we refer the interested reader to Vasca et al. [2] for more details on other modulations. Using natural sampling with a sawtooth carrier (trailing edge) leads to a pulsed signal as depicted in Fig. 2, where the switching (pulse) period $T_{\mathrm{s}}=1/f_{\mathrm{s}}$ and the duty cycle ${d}$ are the quantities defining the pulses. The switching frequency is constant while the duty cycle varies with time.

Large-scale simulations of power converters require simplifications, e.g., ideal switching behaviour of the transistors and ideal diode [1]. Nonetheless, the simulation of power converters with conventional time discretization is still computationally expensive since a high number of time steps is necessary to properly resolve the ripples in the solution. Special techniques are commonly applied to detect switching events and up to now most software uses techniques which reduces the order of the applied time discretization to lowest order, see Tant et al. [3]. In this paper we develop a method which alleviates the need of high number of time steps by splitting the solution into fast and slowly varying parts. The resulting method is a (high-order) multirate approach, which is based on a concept called Multirate Partial Differential Equations (MPDEs) [4, 5]. For DC-DC power converters the method has already been proposed in earlier papers by Pels et al. [6, 7]. We extend the concept to DC-AC power converters (also called inverters), in the following. In the process we restrict ourselves to power converter circuits which can be described by linear electrical elements and thus by a system of linear ordinary differential equations (ODEs). An extension to nonlinear problems can be achieved similarly as described in Pels et al. [6]. The method is verified using the example of the converter in Fig. 1, which is referred to as inverter in the following since it is operated in DC-AC mode. The main contribution is the efficient treatment of (sinusoidally) varying duty cycle.

The paper is structured as follows. Section 2 introduces the concept of MPDEs and their relation to the original ODEs describing the application. The focus of Section 3 lies on the modeling of the DC-AC converters using MPDEs. It describes how an efficient simulation can be achieved. The methods used for solving the MPDEs, i.e., a Galerkin approach for fast periodic parts of the solution and a conventional time discretization for the slowly varying parts are presented. Section 4 introduces B-spline basis functions used for the solution expansion and discusses their properties. It is shown that the basis functions are well suited for use with the proposed method since they allow to exploit smoothness almost everywhere but still lead to a cheap assembly of the arising matrices. Finally in Section 5, the method is applied to the example of the inverter and the computational efficiency is discussed. Section 6 concludes the paper.

2 Introduction to Multirate Partial Differential Equations

Let the mathematical model of the converter be described by linear ordinary differential equations in the form

[TABLE]

where $\mathbf{A},\mathbf{B}\in\mathbb{R}^{N_{\mathrm{s}}\times N_{\mathrm{s}}}$ are matrices, which depend on the topology of the circuit and the electrical elements, $\mathbf{x}(t)\in\mathbb{R}^{N_{\mathrm{s}}}$ is the unknown solution consisting of voltages and currents, $\mathbf{x}_{0}$ is the vector of initial conditions and $\mathbf{c}(t)\in\mathbb{R}^{N_{\mathrm{s}}}$ is the excitation of the circuit.

To obtain the multirate formulation, two artifical time scales $t_{1}$ and $t_{2}$ are introduced and (1) is rewritten in terms of the corresponding multivariate solution $\mathbf{\widehat{x}}(t_{1},t_{2})$ and excitation $\mathbf{\widehat{c}}(t_{1},t_{2})$ which yields the MPDEs [4, 5]

[TABLE]

The relation between the ODEs (1) and MPDEs (2) is given by [4, 5]

[TABLE]

Therefore, if any right-hand side $\mathbf{\widehat{c}}(t_{1},t_{2})$ can be found which fulfills the relation $\mathbf{\widehat{c}}(t,t)=\mathbf{c}(t)$, then the solution of the original system of ODEs can be extracted from the solution of the MPDEs by using $\mathbf{x}(t)=\mathbf{\widehat{x}}(t,t)$, i.e., evaluating the multivariate solution along a diagonal line through the computational domain. It is important to note that the simulation efficiency depends on the choice of $\mathbf{\widehat{c}}(t_{1},t_{2})$ and of course the methods which are used for solving the MPDEs. In our case, we require fast changes to occur along the “fast time scale” $t_{2}$ and slow changes along the “slow time scale” $t_{1}$. The multivariate right-hand side $\mathbf{\widehat{c}}(t_{1},t_{2})$ has then to be chosen appropriately, which depends on the application and is detailed in Section 5.

3 Multirate Modeling of DC-AC converters

The fast varying components, i.e., the ripples in the solution of the power converters, are assumed to be representable by basis functions $p_{k}(\tau(t_{2}),{d}(t_{1}))$, which depend on the fast time scale $t_{2}$ and on the duty cycle ${d}(t_{1})$ assumed to be slowly varying. The first assumption is necessary for the convergence of the method but obvious for most classical bases. The second assumption on the dynamics of the duty cycle is not crucial but will be important to obtain an efficient method. The basis functions are periodic (with period $T_{\mathrm{s}}$), which is realized using the function $\tau(t_{2})=\frac{t_{2}}{T_{\mathrm{s}}}\text{ modulo }1$, also called the relative time. The multivariate solution is expanded into the basis functions and slowly varying coefficients $w_{j,k}(t_{1})$ yielding [7, 6, 8]

[TABLE]

where

[TABLE]

Inserting the solution expansion into the partial derivatives from (2) leads to

[TABLE]

To solve the MPDEs, two methods are applied, one method for each time scale. To resolve the fast changes represented by the basis functions $p_{k}$, a Galerkin approach is applied. The slowly varying coefficients are solved afterwards with a conventional time discretization. Applying the Ritz-Galerkin approach to the MPDEs (2) with respect to $t_{2}$ and on the interval of periodicity $[0,T_{\mathrm{s}})$ yields

[TABLE]

Using integration by parts gives (function arguments are omitted in the following for the sake of readability)

[TABLE]

Since the solution is periodic with respect to $t_{2}$ in the interval $[0,T_{\mathrm{s}})$, the boundary term vanishes. Inserting (6) and substituting $t_{2}$ with $\tau$ leads to

[TABLE]

Assembling all equations finally yields a linear system of ODEs with time-varying coefficients

[TABLE]

where

[TABLE]

and

[TABLE]

The equation system (10) can be solved using conventional time discretization with much larger time steps than for the original problem, since fast varying changes are already taken into account by the Galerkin approach and only the envelope is resolved. A disadvantage of (10) is the larger equation systems, which are $N_{\mathrm{p}}+1$ times larger than the original ones (1).

To solve the system (10) numerically, initial values have to be specified. To achieve an efficient simulation the following choice has proven advantageous:

Calculate the steady-state of the system (10), i.e., \left[\begin{array}[]{ccc}\mathbf{w}^{s}_{1}(0)&\dots&\mathbf{w}^{s}_{N_{\mathrm{s}}}(0)\end{array}\right]^{\top}=\cB^{-1}(0)\cC(0) and use the solution expansion ${\widehat{x}}^{s}_{j}(0,t_{2})=\mathbf{p}^{\top}\!(\tau(t_{2}),{d}(0))\,\mathbf{w}^{s}_{j}(0)$, $j=1,\dots,N_{\mathrm{s}}$ to reconstruct the solution (ripple). 2. 2.

Since the reconstructed solution in steady-state does (usually) not fulfill the initial condition of the original ODEs (1), i.e., $\mathbf{\widehat{x}}^{s}(0,0)\neq\mathbf{x}_{0}$, it has to be shifted by a constant $\mathbf{c}^{\mathrm{s}}$ as such that $\mathbf{\widehat{x}}(0,0)=\mathbf{\widehat{x}}^{s}(0,0)-\mathbf{c}^{\mathrm{s}}=\mathbf{x}_{0}$. The shift is accomplished by modifying the coefficients $\mathbf{w}^{s}(0)$. In case of B-splines, due to partition of unity, this is achieved by choosing the coefficients $\mathbf{w}_{j}(0)=\mathbf{w}^{s}_{j}(0)-c^{\mathrm{s}}_{j}$ as initial values for (10), where $c^{\mathrm{s}}_{j}$ is the $j$-th component of $\mathbf{c}^{\mathrm{s}}$. This corresponds to the zero initial conditions as we use them for the original set of ODEs (1).

Other choices may still lead to the correct solution but may require higher effort from the time discretization algortihm.

4 Choice of basis functions

For the solution expansion (4) B-spline basis functions are employed. They allow a high-order basis while still capturing the $C^{0}$ continuity as it appears in the current ripples by construction. In the literature splines have also been used in the context of MPDE methods but for high-frequency problems. Brachtendorf et al. [9] for instant use cubic and exponential splines, which lead to a better approximation than Fourier basis functions for steep transients, while still enabling a simple extraction of the frequency spectrum and Bittner et al. [10] use an adaptive spline-wavelet to approximate solutions with steep transients. Both do not deal with $C^{0}$ solutions by construction. In contrast we take advantage of the a-priori knowledge of the excitation switching instant to contruct appropriate basis functions for the solution representation.

4.1 Introduction to B-splines

In this subsection we briefly introduce the most important properties and definitions concerning B-splines for this work. The information is taken from Piegl et al. [11], to which the interested reader is referred for more details.

To define B-splines basis functions we first setup a knot vector $\Xi=\{\xi_{0},\dots,\xi_{m}\}$, which is sorted in ascending order, i.e., $\xi_{i}\leq\xi_{i+1}$, $i=0,\dots,m-1$. The basis functions of degree ${\bar{p}}$ are now build up recursively from the piecewise constant basis function ${\bar{p}}=0$

[TABLE]

by the Cox-DeBoor recurrence formula

[TABLE]

For our purposes we assume the knot vector to be open (also called clamped), i.e., it is of the form

[TABLE]

where the first and the last knots appear ${\bar{p}}+1$ times. The regularity $r_{j}$ of the B-spline basis functions across the knots, i.e., their continuity $C^{r_{j}}$ across the knot $\xi_{j}$, $j=0,\dots,m$, can be controlled using knot repetitions. The maximum regularity for degree ${\bar{p}}$ B-splines is given by $r_{j,\mathrm{max}}={\bar{p}}-1$, which corresponds to knot repetitions of one for all knots except the knots at the boundary. To obtain less regularity, a knot repetition is introduced. Continuity $C^{r_{j}}$ is achieved by a knot repetition of $s_{j}={\bar{p}}-r_{j}$. For instant, to represent a $C^{0}$ continuity across knot $j$, the knot repetition is given by $s_{j}={\bar{p}}-0={\bar{p}}$.

4.2 Choice of knot vector

For use in this work we define a knot vector as such, that the basis functions have a $C^{0}$ continuity at the point where the excitation changes its state. Defining the basis in terms of the relative time, i.e., $\tau\in[0,1]$, and the (fixed) duty cycle ${d}\in(0,1)$ the simplest knot vector depending on the degree ${\bar{p}}$ is given by

[TABLE]

Additional knot refinement (corresponding to $h$-refinement in Finite Element Methods) leads to

[TABLE]

where $K$ is the number of additional knots inserted before and after the $C^{0}$ continuity, and $\alpha_{k}\in[0,1]$, $\beta_{k}\in[0,1]$, $\forall k=1,\dots,K$ in ascending order.

Using the refined knot vector (20) leads to the set of basis functions

[TABLE]

as depicted exemplary for $K=1$ and ${\bar{p}}=2$ in Fig. 3, i.e., we have in total $2({\bar{p}}+K)+1$ basis functions. The basis functions for the solution expansion (4) are finally given by

[TABLE]

The periodicity of the set of basis functions is ensured using periodic boundary conditions in the implementation.

Using the set of B-spline basis functions as introduced above leads to a cheap calculation of the matrices arising in the MPDE approach since their elements depend only linearly on the duty cycle. For a low-order basis like classical Finite Element hat functions (i.e., ${\bar{p}}=1$) this is rather obvious whilst for higher order B-splines it is more involved due to their non-local support. Hence this property is analyzed in the following theorem.

Theorem 1** (Dependency of the matrices $\cI,\cQ,\cU$ on the duty cycle).**

Using the B-spline basis functions (21), only the matrix $\cI$ from (13) depends linearly on the duty cycle, i.e., $\cI({d})=\frac{\cI({d}_{0})-\cI({d}_{1})}{{d}_{0}-{d}_{1}}({d}-{d}_{0})+\cI({d}_{0})$, with ${d}_{0},{d}_{1}\in(0,1)$ and ${d}_{0}\neq{d}_{1}$, and the matrices $\cQ,\cU$ from (14), (15), respectively, are independent of the duty cycle.

Proof.

The set of basis functions (21) can be split into three parts: The basis functions which are supported on the knots in the interval $[0,{d}]$, i.e., the basis functions $\{{P}_{0,{\bar{p}}},\dots,{P}_{{\bar{p}}+K-1,{\bar{p}}}\}$; the basis functions on the interval $[d,1]$, i.e., $\{{P}_{{\bar{p}}+K+1,{\bar{p}}},\dots,{P}_{2{\bar{p}}+2K,{\bar{p}}}\}$; and the remaining $C^{0}$ basis function which depends on knots from the entire interval $[0,1]$ , i.e. $\{{P}_{{\bar{p}}+K,{\bar{p}}}\}$.

Calculating the basis functions $\{{P}_{0,{\bar{p}}},\dots,{P}_{{\bar{p}}+K-1,{\bar{p}}}\}$ using the Cox-DeBoor formula leads to a recursion of the form

[TABLE]

where $\gamma_{1,i},\dots,\gamma_{4,i}$ are constants depending on $\alpha_{1},\dots,\alpha_{K}$, but independent of the duty cycle ${d}$. Therefore the resulting polynomials ${P}_{i,{\bar{p}}}(\xi)$ can be written as

[TABLE] 2. 2.

To calculate the B-splines $\{{P}_{{\bar{p}}+K+1,{\bar{p}}},\dots,{P}_{2{\bar{p}}+2K,{\bar{p}}}\}$, we first redefine the knot vector for ease of notation by shifting it

[TABLE]

The resulting B-splines are the same as with the original knot vector $\Xi_{{\bar{p}},K}$ but shifted by $-{d}$. Using the Cox-DeBoor formula leads to

[TABLE]

The final polynomials $\widehat{{P}}_{i,{\bar{p}}}(\xi,{d})$ can be written as

[TABLE] 3. 3.

The single basis function ${P}_{{\bar{p}}+K,{\bar{p}}}$ consists of two terms due to the knot repetition. In the Cox-DeBoor formula, the first term in the sum stems from basis function in the interval $[0,d]$, the second term stems from basis functions in the interval $[1,d]$. Therefore as well as for the other basis functions, the left term can be written as a polynomial of the form (24), the right term can be written as a polynomial of the form (27).

Using the above knowlegde, the matrix (13) is of the form

[TABLE]

where $\cI_{i,j}^{0},\cI_{i,j}^{1}$ are constants. Therefore the matrix $\cI$ depends only linearly on the duty cycle.

The matrices (14) and (15) are given by

[TABLE]

and

[TABLE]

where $\cQ_{i,j}^{0}$, $\cQ_{i,j}^{1}$, $\cU_{i,j}^{0}$ and $\cU_{i,j}^{1}$ are constants. The matrices are thus independent of the duty cycle ${d}$. ∎

This is a helpful result since it allows a cheap calculation of the matrices $\cI$, $\cQ$ and $\cU$ for different duty cycles and time steps.

5 Numerical results

To verify the proposed method and measure its efficiency, we test the method on the inverter example as depicted in Fig. 1. To measure the accuracy of solutions they are compared to a reference solution calculated in Simulink using PLECS. In the following we use three simulation setups: 1.) The MPDE approach; 2.) A conventional time discretization with switch detection implemented in MATLAB; 3.) A simulation in Simulink using PLECS. The conventional approaches 2.) and 3.) exploit a-priori knowledge on discontinuities, e.g. by a pre-defined event function.

The buck converter is described by the ordinary differential equation

[TABLE]

where $L=4\,$mH (inductance of the coil), $C=10$ $\mu$F (capacitance of the capacitor), $R_{\mathrm{L}}=10$ m$\Omega$ (coil resistance), and $R=20\,\Omega$ (load resistance) are fixed parameters, and $v_{\mathrm{C}}$, $i_{\mathrm{L}}$ and $v_{\mathrm{i}}$ are the voltage at the capacitor, the current through the coil and the PWM excitation, respectively. The excitation for the inverter is generated using natural sampling PWM with a sawtooth carrier. It can be written as

[TABLE]

where $s(t)=\frac{t}{T_{\mathrm{s}}}\text{ mod }1$ is the sawtooth carrier, $\bar{v}_{\mathrm{i}}$ is the peak excitation voltage, and $\mathrm{sgn}(t)$ is the sign function. The switching frequency is fixed at $f_{\mathrm{s}}=1/T_{\mathrm{s}}=5\,$kHz. An excerpt of the excitation is shown in Fig. 4.

For the MPDE formulation we choose the multivariate excitation as $\widehat{v}_{\mathrm{i}}(t_{1},t_{2})=v_{\mathrm{i}}(t_{2},{d}(t_{1}))$, i.e., we force the duty cycle, which is slowly varying compared to the switching of the excitation, to evolve along the time scale $t_{1}$ and the switching to occur along the time scale $t_{2}$. As it turns out the right-hand side of the ODEs after the MPDE approach, i.e., (12) is also linearly dependent on the duty cycle and thus allows a cheap evaluation. The proof is analog to the one detailled in Section 4. The duty cycle is assumed to be sinusoidal and given by

[TABLE]

where $\bar{v}_{\mathrm{C},\text{desired}}$ is the desired peak output voltage of the converter and ${f_{\mathrm{ac}}}$ is the desired frequency of the AC output voltage. The constants are fixed to $\bar{v}_{\mathrm{i}}=350\,$V, $\bar{v}_{\mathrm{C},\text{desired}}=325\,$V (corresponding to $230\,$V effective voltage) and ${f_{\mathrm{ac}}}=50\,$Hz. The resulting inverter output, i.e., voltage at the capacitor and current through the coil, are depicted in Fig. 4. The multivariate solution of the MPDEs is depicted in Fig. 5. The solution of the original equations (1) is marked as black line.

The MPDE approach is implemented in MATLAB and the solver ode15s is used for the Simulink, the time discretization and the MPDE approach simulation. For the MPDE approach, the maximum order MaxOrder for the solver is set to 2 whilst for the other two methods the original setting of 5 is used. To compare the solutions, the relative ${\mathrm{L}^{2}}$ error of the capacitor voltage on the simulation interval ${\Omega=[0,4{T_{\mathrm{ac}}}]}$

[TABLE]

is approximated using the mid-point quadrature rule. The error of the current through the coil $\epsilon_{\mathrm{i}}$ is defined analogously. The reference solution $v_{\mathrm{C,ref}}(t)$ is calculated using a very fine tolerance of $\mathtt{abstol}=\mathtt{reltol}=10^{-12}$ and a maximum step size of $T_{\mathrm{s}}/1000$ in Simulink.

MPDE approach simulation results are obtained for three different settings:

Lowest order approximation: ${\bar{p}}=1$ and $K=1$ (corresponds to Finite Element hat functions), 2. 2.

Medium order approximation: ${\bar{p}}=2$ and $K=1$, 3. 3.

High order approximation: ${\bar{p}}=3$ and $K=3$.

The error $\epsilon$ of the MPDE approach using these settings is calculated for different tolerances of the time discretization and is depicted in Fig. 6 for both the current through the coil and the voltage at the capacitor. The error is evaluated using 100 points per cycle. Reference solution values which are not available at the corresponding points are interpolated linearly. Fig. 6 shows that the error decreases with smaller tolerances for the solver until it stagnates, which is due to the fact that the approximation of the Galerkin approach bounds the accuracy. It is also visible that with better discretization settings for the Galerkin approach, the stagnation is shifted to smaller tolerances. For the three settings introduced above we now fix the tolerances for the solver as such, that we achieve highest accuracy without wasting computational effort in the stagnation region. The corresponding tolerances are marked as dots in Fig. 6 and summarized in Table 1. To be able to compare statistical data like the number of time steps, number of LU decompositions and number of function evaluations used by the solvers, the accuracy of conventional time discretization in MATLAB and PLECS are controlled by the tolerance $\mathtt{abstol}=\mathtt{reltol}$. The error of both with respect to the tolerance is shown in Fig. 7 for both voltage and current. The errors corresponding approximately to the errors of the MPDE approach are listed in Table 1 with the associated solver tolerance. Statistical data of all approaches, i.e., number of time steps, number of failed steps, number of LU decompositions, number of function evaluations and number of forward/backward substitution (solution of linear systems) are compared in Table 2 for the settings in Table 1. For the PLECS simulation only the number of time steps is supplied. The other figures of merit can to our knowledge not be extracted from the Simulink solvers. For high solver tolerance, i.e., $\mathtt{abstol}=\mathtt{reltol}=10^{-2}$, the conventional time discretization in MATLAB is slightly more accurate than the PLECS results (see Fig. 7) since the MaxStep option is used to guarantee that no switching events are missed. This explains the considerably higher number of time steps compared to PLECS in Table 2, setting 1 (lowest order). For the other settings, the number of time steps between MATLAB and PLECS is almost similar.

As can be seen in Table 2 the MPDE approach is several times faster than the conventional methods with respect to the considered quantities. Especially when low or medium accuracy (setting 1 and 2) is necessary, the speedup is very high. Furthermore the efficiency of the method increases with higher switching frequency, see [7]. If a switching frequency is assumed to be twice as high, for instance, the conventional time discretization will likely take twice as much time since twice as many ripples have to be resolved. The MPDE approach on the other hand will need approximately the same number of time steps since the envelope does not change. As a result the speedup compared with the conventional methods would be twice as high.

It should be noted, that the actual efficiency in terms of computing time depends on the efficiency of the solver used for the solution of the equation systems. Since the equation systems in the MPDE approach are larger than the ones of the original problem, the solver will take more time. This has to be taken into account when choosing the number of basis functions for the solution expansion. The MPDE approach will be especially efficient, if the computational effort for function evaluation is higher than the effort for solving the equation systems.

6 Conclusions

A multirate method for the efficient simulation of DC-AC switch-mode inverters has been presented. The system of equations describing the converter circuit is first formulated as a system of Multirate Partial Differential Equations (MPDEs), which allow to split the solution into components of different explicitly stated time scales. The MPDEs are efficiently solved using a combination of a Galerkin approach with B-spline basis functions for the solution expansion, and a conventional time discretization. The functionality of the proposed approach is verified on a simple inverter circuit with sinusoidal AC output. The computational efficiency is analyzed among others in terms of the number of time steps, number of LU decompositions, number of functions evaluations, compared to a conventional time discretization of the problem. It shows the high potential efficiency of the method. The method is particularly efficient in applications in which the function evaluations, i.e., the evaluations of the matrices/functions in the differential equation describing the application are computationally expensive.

Acknowledgments

This work is supported by the “Excellence Initiative” of German Federal and State Governments and the Graduate School CE at TU Darmstadt.