4. Multitime partitioned simulation engines

This section is devoted to briefly review some advanced circuit-block partitioning numerical simulation techniques operating in a multivariate time-domain framework. Section 4.1 introduces the multivariate formulation theory, in which the 1-D time is converted into a set of artificial time variables. Section 4.2 addresses some fundamental aspects of the numerical simulation algorithms. Finally, Section 4.3 describes some techniques for automatic circuit-block partition.

### 4.1. Multivariate formulation

The multivariate formulation is a useful stratagem that plays an important role in the EDA scientific community, especially in the RF and microwave areas. It was initially introduced in 1996 [12] as a sophisticated derivation of quasi-periodic HB, and it has been adopted by other researchers (e.g., [13–19]), which have demonstrated that it can be an efficient strategy when dealing with electronic circuits operating on widely distinct time scales. The success of multivariate formulation relies on the fact that voltages and currents containing components that evolve themselves at two, or more, widely separated rates of variation can be represented much more efficiently if we define them as functions of two, or more, time variables (artificial time scales). With this stratagem all signals (stimuli and responses) will be represented as multivariate functions, which will imply that dynamic behavior of the circuits will no longer be modeled by DAE, or ODE, systems formulated in the 1-D time. It will be described by partial differential systems.

In order to provide a simple and illustrative mathematical description of the multivariate formulation let us consider the classical example of a generic nonlinear RF circuit driven by an envelope modulated signal of the form

$$\mathbf{x}(t) = \mathbf{e}(t)\cos\left(\omega\_{\mathbb{C}}t + \phi(t)\right),\tag{26}$$

where tmodT<sup>2</sup> represents the remainder of division of t by T2.

∂t<sup>1</sup>

4.2. Partitioned algorithms for envelope following computation

<sup>p</sup>ðbyð Þ <sup>t</sup>1; <sup>t</sup>2;…; tm Þ þ <sup>∂</sup>qð Þ <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup>2;…; tm

by setting t<sup>1</sup> ¼ t<sup>2</sup> ¼ ⋯ ¼ tm ¼ t.

reviewed in the following.

The generalization of the bivariate strategy to a multidimensional problem with more than two time scales is straightforward. In fact, if the signals in the circuit present m separate rates of

þ ⋯ þ

and the univariate solution, <sup>y</sup>ð Þ<sup>t</sup> , may be recovered from its multivariate form, <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup>2;…; tm ,

The main advantage of the earlier described MPDAE approach is that it can result in significant improvements in simulation speed when compared to DAE-based alternatives [13–15]. However, by itself this approach does not perform any distinction between nodes or blocks in the circuit under analysis. In fact, in the first multivariate circuit simulation schemes initially introduced in [12], and then in [13], the same numerical algorithm was used to compute all the unknowns of the circuit. Only a few years later other advanced multivariate algorithms were proposed (e.g., [16–19]) in way to take into account possible circuit's heterogeneities. These algorithms are based on circuit-block partitioning stratagems defined within the multivariate frameworks. The most important ones regarding pure time-domain operations are briefly

As stated earlier, envelope modulated regimes are typical cases of practical interest. Computing responses to excitations of the form of (26) is a technique generally referred to as envelope following, which correspond to a combination of initial and periodic boundary conditions in the

þ

defined on the rectangular domain 0; tFinal ½ �� ½ � 0; T<sup>2</sup> . gð Þ� is a vector-valued initial-condition satisfying <sup>g</sup>ð Þ¼ <sup>0</sup> <sup>g</sup>ð Þ¼ <sup>T</sup><sup>2</sup> <sup>y</sup>ð Þ<sup>0</sup> , and <sup>b</sup>yð Þ¼ <sup>t</sup>1; <sup>0</sup> <sup>b</sup>yð Þ <sup>t</sup>1; <sup>T</sup><sup>2</sup> is the periodic boundary condition due to the periodicity of the problem in the t<sup>2</sup> fast time scale. In order to solve this initial-boundary value problem let us begin to consider the following semi-discretization of 0; tFinal ½ �� ½ � 0; T<sup>2</sup> in

0 ¼ t1,<sup>0</sup> < t1,<sup>1</sup> < ⋯ < t1,i�<sup>1</sup> < t1,i < ⋯ < t1,K<sup>1</sup> ¼ tFinal, h1,i ¼ t1,i � t1,i�<sup>1</sup>, (33)

<sup>∂</sup>qð Þ <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> ∂t<sup>2</sup>

<sup>¼</sup> <sup>b</sup>xð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup>

(32)

bivariate framework, leading to the following initial-boundary value problem

<sup>b</sup>yð Þ¼ <sup>0</sup>; <sup>t</sup><sup>2</sup> <sup>g</sup>ð Þ <sup>t</sup><sup>2</sup> <sup>b</sup>yð Þ¼ <sup>t</sup>1; <sup>0</sup> <sup>b</sup>yð Þ <sup>t</sup>1; <sup>T</sup><sup>2</sup>

<sup>∂</sup>qð Þ <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> ∂t<sup>1</sup>

<sup>p</sup>ð Þþ <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup>

the t<sup>1</sup> slow time dimension defined by

<sup>∂</sup>qð Þ <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup>2;…; tm ∂tm

Numerical Simulation of Electronic Systems Based on Circuit-Block Partitioning Strategies

<sup>¼</sup> <sup>b</sup>xð Þ <sup>t</sup>1; <sup>t</sup>2; …; tm , (31)

33

http://dx.doi.org/10.5772/intechopen.75490

change, then m time scales will be used. In that case (29) assumes the generic form

where e tð Þ and ϕð Þt are, respectively, the amplitude and phase slowly varying baseband signals (envelope), modulating the cos ð Þ ωCt fast-varying carrier. Computing the numerical solution of such circuit using conventional 1-D time-step integrators (RK methods or LMM methods) is very expensive. This is so because the response of the circuit has to be evaluated during a long time window defined by the slowly varying envelope, wherein the time-step size is severely restricted by the high frequency carrier. Taking into consideration the disparity between the baseband and the carrier time scales, and assuming that they are also uncorrelated, which is typically true, we are able to rewrite (26) as a bivariate function

$$
\widehat{\mathfrak{X}}(t\_1, t\_2) = e(t\_1) \cos \left( \omega\_\zeta t\_2 + \phi(t\_1) \right),
\tag{27}
$$

where t<sup>1</sup> is the slow baseband artificial time scale and t<sup>2</sup> is the fast carrier artificial time scale. It must be noted that <sup>b</sup>x tð Þ <sup>1</sup>; <sup>t</sup><sup>2</sup> is a periodic function with respect to <sup>t</sup><sup>2</sup> but not to <sup>t</sup>1, that is,

$$
\widehat{\mathfrak{X}}(t\_1, t\_2) = \widehat{\mathfrak{X}}(t\_1, t\_2 + T\_2), \quad T\_2 = 2\pi/\omega\_{\mathbb{C}}.\tag{28}
$$

This means that a generic 1-D 0; tFinal ½ � time interval will be mapped into a 2-D 0; tFinal ½ �� ½ � 0; T<sup>2</sup> rectangular domain. It is easy to attest that, in general, the number of points required to represent <sup>b</sup>x tð Þ <sup>1</sup>; <sup>t</sup><sup>2</sup> in 0; tFinal ½ �� ½ � <sup>0</sup>; <sup>T</sup><sup>2</sup> is much less than the number of points needed to represent x tð Þ in 0; tFinal ½ �. This is especially evident when the t<sup>1</sup> and t<sup>2</sup> time scales are widely separated [13, 14].

Let us now consider the DAE system of (1) modeling the dynamic behavior of a generic RF circuit excited by the envelope modulated signal of (26). Taking the abovementioned into account, we will adopt the following procedure for the vector-valued functions xð Þt and yð Þt : t is replaced by t<sup>1</sup> in the slowly varying entities (envelope time scale) and t is replaced by t<sup>2</sup> in the fast-varying entities (RF carrier time scale). The application of this stratagem will turn the DAE system of (1) into the following partial differential system

$$
\hat{\mathfrak{p}}(\hat{\mathfrak{y}}(t\_1, t\_2)) + \frac{\partial \mathfrak{q}(\widehat{\mathfrak{y}}(t\_1, t\_2))}{\partial t\_1} + \frac{\partial \mathfrak{q}(\widehat{\mathfrak{y}}(t\_1, t\_2))}{\partial t\_2} = \widehat{\mathfrak{x}}(t\_1, t\_2), \tag{29}
$$

usually denoted as multirate partial differential algebraic equations' (MPDAE) system [13, 14]. It is easy to demonstrate that, if <sup>b</sup>xð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> and <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> satisfy (29), then the univariate forms <sup>x</sup>ðÞ¼ <sup>t</sup> <sup>b</sup>xð Þ <sup>t</sup>; <sup>t</sup> and <sup>y</sup>ðÞ¼ <sup>t</sup> <sup>b</sup>yð Þ <sup>t</sup>; <sup>t</sup> satisfy (1) [13]. Thus, <sup>y</sup>ð Þ<sup>t</sup> may be retrieved from its bivariate form <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> by simply setting <sup>t</sup><sup>1</sup> <sup>¼</sup> <sup>t</sup><sup>2</sup> <sup>¼</sup> <sup>t</sup>, meaning that univariate solutions of (1) are available on diagonal lines t<sup>1</sup> ¼ t, t<sup>2</sup> ¼ t, along the bivariate solutions of (29). In truth, attending to the periodicity of the problem in the t<sup>2</sup> time dimension, the univariate version of the vector containing the circuit responses is obtained from its bivariate representation on the rectangular domain 0; tFinal ½ �� ½ � 0; T<sup>2</sup> as

$$\hat{y}(t) = \hat{\bar{y}}(t, t \text{mod} \, T\_2), \tag{30}$$

where tmodT<sup>2</sup> represents the remainder of division of t by T2.

x tðÞ¼ e tð Þ cos <sup>ω</sup>Ct <sup>þ</sup> <sup>ϕ</sup>ð Þ<sup>t</sup> � �, (26)

<sup>b</sup>x tð Þ¼ <sup>1</sup>; <sup>t</sup><sup>2</sup> <sup>b</sup>x tð Þ <sup>1</sup>; <sup>t</sup><sup>2</sup> <sup>þ</sup> <sup>T</sup><sup>2</sup> , T<sup>2</sup> <sup>¼</sup> <sup>2</sup>π=ωC: (28)

� �, (27)

<sup>¼</sup> <sup>b</sup>xð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> , (29)

where e tð Þ and ϕð Þt are, respectively, the amplitude and phase slowly varying baseband signals (envelope), modulating the cos ð Þ ωCt fast-varying carrier. Computing the numerical solution of such circuit using conventional 1-D time-step integrators (RK methods or LMM methods) is very expensive. This is so because the response of the circuit has to be evaluated during a long time window defined by the slowly varying envelope, wherein the time-step size is severely restricted by the high frequency carrier. Taking into consideration the disparity between the baseband and the carrier time scales, and assuming that they are also uncorrelated, which is

<sup>b</sup>x tð Þ¼ <sup>1</sup>; <sup>t</sup><sup>2</sup> e tð Þ<sup>1</sup> cos <sup>ω</sup>Ct<sup>2</sup> <sup>þ</sup> <sup>ϕ</sup>ð Þ <sup>t</sup><sup>1</sup>

where t<sup>1</sup> is the slow baseband artificial time scale and t<sup>2</sup> is the fast carrier artificial time scale. It

This means that a generic 1-D 0; tFinal ½ � time interval will be mapped into a 2-D 0; tFinal ½ �� ½ � 0; T<sup>2</sup> rectangular domain. It is easy to attest that, in general, the number of points required to represent <sup>b</sup>x tð Þ <sup>1</sup>; <sup>t</sup><sup>2</sup> in 0; tFinal ½ �� ½ � <sup>0</sup>; <sup>T</sup><sup>2</sup> is much less than the number of points needed to represent x tð Þ in 0; tFinal ½ �. This is especially evident when the t<sup>1</sup> and t<sup>2</sup> time scales are widely

Let us now consider the DAE system of (1) modeling the dynamic behavior of a generic RF circuit excited by the envelope modulated signal of (26). Taking the abovementioned into account, we will adopt the following procedure for the vector-valued functions xð Þt and yð Þt : t is replaced by t<sup>1</sup> in the slowly varying entities (envelope time scale) and t is replaced by t<sup>2</sup> in the fast-varying entities (RF carrier time scale). The application of this stratagem will turn the

þ

usually denoted as multirate partial differential algebraic equations' (MPDAE) system [13, 14]. It is easy to demonstrate that, if <sup>b</sup>xð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> and <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> satisfy (29), then the univariate forms <sup>x</sup>ðÞ¼ <sup>t</sup> <sup>b</sup>xð Þ <sup>t</sup>; <sup>t</sup> and <sup>y</sup>ðÞ¼ <sup>t</sup> <sup>b</sup>yð Þ <sup>t</sup>; <sup>t</sup> satisfy (1) [13]. Thus, <sup>y</sup>ð Þ<sup>t</sup> may be retrieved from its bivariate form <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> by simply setting <sup>t</sup><sup>1</sup> <sup>¼</sup> <sup>t</sup><sup>2</sup> <sup>¼</sup> <sup>t</sup>, meaning that univariate solutions of (1) are available on diagonal lines t<sup>1</sup> ¼ t, t<sup>2</sup> ¼ t, along the bivariate solutions of (29). In truth, attending to the periodicity of the problem in the t<sup>2</sup> time dimension, the univariate version of the vector containing the circuit responses is obtained from its bivariate representation on the rectangular

<sup>∂</sup>qð Þ <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> ∂t<sup>2</sup>

<sup>y</sup>ðÞ¼ <sup>t</sup> <sup>b</sup>yð Þ <sup>t</sup>; <sup>t</sup>modT<sup>2</sup> , (30)

must be noted that <sup>b</sup>x tð Þ <sup>1</sup>; <sup>t</sup><sup>2</sup> is a periodic function with respect to <sup>t</sup><sup>2</sup> but not to <sup>t</sup>1, that is,

typically true, we are able to rewrite (26) as a bivariate function

32 Numerical Simulations in Engineering and Science

DAE system of (1) into the following partial differential system

<sup>∂</sup>qð Þ <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup> ∂t<sup>1</sup>

<sup>p</sup>ð Þþ <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup><sup>2</sup>

separated [13, 14].

domain 0; tFinal ½ �� ½ � 0; T<sup>2</sup> as

The generalization of the bivariate strategy to a multidimensional problem with more than two time scales is straightforward. In fact, if the signals in the circuit present m separate rates of change, then m time scales will be used. In that case (29) assumes the generic form

$$p(\widehat{\boldsymbol{y}}(t\_1, t\_2, \ldots, t\_m)) + \frac{\partial q(\widehat{\boldsymbol{y}}(t\_1, t\_2, \ldots, t\_m))}{\partial t\_1} + \cdots + \frac{\partial q(\widehat{\boldsymbol{y}}(t\_1, t\_2, \ldots, t\_m))}{\partial t\_m} = \widehat{\mathbf{x}}(t\_1, t\_2, \ldots, t\_m), \tag{31}$$

and the univariate solution, <sup>y</sup>ð Þ<sup>t</sup> , may be recovered from its multivariate form, <sup>b</sup>yð Þ <sup>t</sup>1; <sup>t</sup>2;…; tm , by setting t<sup>1</sup> ¼ t<sup>2</sup> ¼ ⋯ ¼ tm ¼ t.

### 4.2. Partitioned algorithms for envelope following computation

The main advantage of the earlier described MPDAE approach is that it can result in significant improvements in simulation speed when compared to DAE-based alternatives [13–15]. However, by itself this approach does not perform any distinction between nodes or blocks in the circuit under analysis. In fact, in the first multivariate circuit simulation schemes initially introduced in [12], and then in [13], the same numerical algorithm was used to compute all the unknowns of the circuit. Only a few years later other advanced multivariate algorithms were proposed (e.g., [16–19]) in way to take into account possible circuit's heterogeneities. These algorithms are based on circuit-block partitioning stratagems defined within the multivariate frameworks. The most important ones regarding pure time-domain operations are briefly reviewed in the following.

As stated earlier, envelope modulated regimes are typical cases of practical interest. Computing responses to excitations of the form of (26) is a technique generally referred to as envelope following, which correspond to a combination of initial and periodic boundary conditions in the bivariate framework, leading to the following initial-boundary value problem

$$\begin{aligned} \mathfrak{p}(\widehat{\mathfrak{y}}(t\_1, t\_2)) + \frac{\partial \mathfrak{q}(\widehat{\mathfrak{y}}(t\_1, t\_2))}{\partial t\_1} + \frac{\partial \mathfrak{q}(\widehat{\mathfrak{y}}(t\_1, t\_2))}{\partial t\_2} &= \widehat{\mathfrak{x}}(t\_1, t\_2) \\ \widehat{\mathfrak{y}}(0, t\_2) = \mathfrak{g}(t\_2) \\ \widehat{\mathfrak{y}}(t\_1, 0) = \widehat{\mathfrak{y}}(t\_1, T\_2) \end{aligned} \tag{32}$$

defined on the rectangular domain 0; tFinal ½ �� ½ � 0; T<sup>2</sup> . gð Þ� is a vector-valued initial-condition satisfying <sup>g</sup>ð Þ¼ <sup>0</sup> <sup>g</sup>ð Þ¼ <sup>T</sup><sup>2</sup> <sup>y</sup>ð Þ<sup>0</sup> , and <sup>b</sup>yð Þ¼ <sup>t</sup>1; <sup>0</sup> <sup>b</sup>yð Þ <sup>t</sup>1; <sup>T</sup><sup>2</sup> is the periodic boundary condition due to the periodicity of the problem in the t<sup>2</sup> fast time scale. In order to solve this initial-boundary value problem let us begin to consider the following semi-discretization of 0; tFinal ½ �� ½ � 0; T<sup>2</sup> in the t<sup>1</sup> slow time dimension defined by

$$0 = t\_{1,0} < t\_{1,1} < \dots < t\_{1,i-1} < t\_{1,i} < \dots < t\_{1,K\_1} = t\_{\bar{\text{Final}}} \qquad h\_{1,i} = t\_{1,i} - t\_{1,i-1} \tag{33}$$

in which K<sup>1</sup> is the total number of steps in t1. Now, using a backward differentiation formula (BDF) [1] to approximate the t<sup>1</sup> derivatives of the MPDAE (for simplicity let us consider here the Gear-2 rule [1]), we obtain for each slow time instant t1,i, from i ¼ 1 to i ¼ K1, the periodic boundary value problem defined by

$$p(\hat{y}\_i(t\_2)) + \frac{3q(\hat{y}\_i(t\_2)) - 4q(\hat{y}\_{i-1}(t\_2)) + q(\hat{y}\_{i-2}(t\_2))}{2h\_{1,i}} + \frac{dq(\hat{y}\_i(t\_2))}{dt\_2} = \hat{\mathfrak{x}}\ (t\_{1,i}, t\_2),\tag{34}$$

$$\hat{y}\_i(0) = \hat{y}\_i(T\_2),$$

detecting the fast-varying and the slowly varying signals. The approach adopted in [16, 17] for

As mentioned earlier, each of the periodic boundary value problems defined by (34) is solved using shooting-Newton based on multirate time-step integrators (MRK schemes). Now, taking into account that shooting is an iterative technique, the key idea of this partitioning strategy is to use a uni-rate scheme on the first shooting iterations needed for each time line t1,i � ½ � 0; T<sup>2</sup> of the 0; tFinal ½ �� ½ � 0; T<sup>2</sup> rectangular domain. For that, it starts by considering m ¼ 1 in the adopted MRK scheme, so that it degenerates into a standard uni-rate RK scheme. This means that all the unknowns (voltages or currents) will be integrated in t<sup>2</sup> with the same microstep h2. After that, the differential system of (34) is partitioned into (35) according to the variations in the t<sup>2</sup> derivatives of the unknowns. Each unknown that practically evidences no variations in its t<sup>2</sup> time derivatives for the entire time line t1,i � ½ � 0; T<sup>2</sup> , that is, each unknown that satisfies the

> � � � min slopej � � � �

where and Tol is a small specified tolerance, will be treated as slow on the next shooting iterations. The remaining unknowns will be treated as fast. In a cheap scheme within a uniform

All the subsequent shooting iterations on the time line t1,i � ½ � 0; T<sup>2</sup> will be conducted in a multirate way using different step sizes. In a nonuniform grid, step sizes h<sup>2</sup> and H<sup>2</sup> ¼ m � h<sup>2</sup> may be chosen using any step-size control tool. In a uniform grid, they can be predefined or

The robustness of this partitioning strategy can be improved if more than one shooting iteration with the uni-rate scheme is considered for each time line t1,i � ½ � 0; T<sup>2</sup> . However, such tactic will obviously conduct to some extra computational cost, turning the overall algorithm

Finally, to conclude this section, it must be highlighted that although efficiency gains are essentially dependent on the ratio between the number of slow and fast signals, significant reductions on the computational effort were reported in the scientific literature (e.g., [16, 17]) for the simulation of several illustrative examples of practical interest using the envelope

This section is intended to provide a brief description of powerful circuit-block partitioning numerical simulation techniques operating in multivariate hybrid time-frequency frameworks. Section 5.1 introduces the fundaments of multivariate hybrid envelope following. Section 5.2

� �

Numerical Simulation of Electronic Systems Based on Circuit-Block Partitioning Strategies

� � � � <sup>=</sup>h2, t2,j <sup>∈</sup>f g <sup>0</sup>; <sup>h</sup>2; <sup>2</sup>h2; …; <sup>T</sup><sup>2</sup> � <sup>h</sup><sup>2</sup> : (38)

� <sup>&</sup>lt; Tol, (37)

http://dx.doi.org/10.5772/intechopen.75490

35

max slopej

� � � <sup>b</sup>yi <sup>t</sup>2,j

�

grid each slope can be simply given by

to be a little less efficient.

slopej <sup>¼</sup> <sup>b</sup>yi <sup>t</sup>2,j <sup>þ</sup> <sup>h</sup><sup>2</sup>

successively refined to achieve a desired accuracy.

following bivariate partitioned techniques.

5. Hybrid time-frequency partitioned techniques

that purpose is briefly described in the following.

condition

where <sup>b</sup>y<sup>i</sup> ð Þ <sup>t</sup><sup>2</sup> <sup>≃</sup> <sup>b</sup><sup>y</sup> <sup>t</sup>1,i ð Þ ; <sup>t</sup><sup>2</sup> . This means that, once <sup>b</sup>y<sup>i</sup>�<sup>1</sup>ð Þ <sup>t</sup><sup>2</sup> is known, the solution on the next slow time instant, <sup>b</sup>y<sup>i</sup> ð Þ t<sup>2</sup> , is obtained by solving (34). Hence, a set of K<sup>1</sup> boundary value problems have to be solved for obtaining the whole bivariate solution in the entire domain 0; tFinal ½ �� ½ � 0; T<sup>2</sup> . A very powerful technique has been proposed in the literature for solving each of the periodic boundary problems defined by (34) in an efficient way [16, 17]. This technique uses shooting-Newton based on MRK schemes. It splits the circuits into two distinct subsets according to the time evolution rates of their voltages and currents and performs time-step integration with different step lengths in each of the consecutive shooting iterations needed to solve (34). For that (34) is firstly divided into the following coupled fast-slow subsystems

$$\begin{split} &p\_{\mathbf{F}}(\widehat{\boldsymbol{y}}\_{F,i}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i}(t\_{2})) \\ &+\frac{3q\_{F}(\widehat{\boldsymbol{g}}\_{F,i}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i}(t\_{2}))-4q\_{F}(\widehat{\boldsymbol{g}}\_{A,i-1}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i-1}(t\_{2}))+q\_{F}(\widehat{\boldsymbol{g}}\_{F,i-2}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i-2}(t\_{2}))}{2h\_{1,i}} \\ &+\frac{dq\_{A}(\widehat{\boldsymbol{g}}\_{F,i}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i}(t\_{2}))}{dt\_{2}} = \widehat{\mathbf{x}}(t\_{1,i},t\_{2}) \\ &p\_{\mathbf{L}}(\widehat{\boldsymbol{\mathcal{Y}}}\_{F,i}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i}(t\_{2})) \\ &+\frac{3q\_{L}(\widehat{\boldsymbol{g}}\_{F,i}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i}(t\_{2}))-4q\_{L}(\widehat{\boldsymbol{g}}\_{F,i-1}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i-1}(t\_{2}))+q\_{L}(\widehat{\boldsymbol{g}}\_{F,i-2}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i-2}(t\_{2}))}{2h\_{1,i}} \\ &+\frac{dq\_{L}(\widehat{\boldsymbol{g}}\_{F,i}(t\_{2}),\widehat{\boldsymbol{y}}\_{L,i}(t\_{2}))}{dt\_{2}} = \widehat{\mathbf{x}}(t\_{1,i},t\_{2}). \end{split} \tag{35}$$

with

$$
\widehat{\boldsymbol{y}}\_{i}(t\_{2}) = \begin{bmatrix}
\widehat{\boldsymbol{y}}\_{\boldsymbol{f},i}(t\_{2}) \\
\widehat{\boldsymbol{y}}\_{\boldsymbol{L},i}(t\_{2})
\end{bmatrix}, \quad \widehat{\boldsymbol{y}}\_{\boldsymbol{F},i}(t\_{2}) \in \mathbb{R}^{n\_{\boldsymbol{f}}}, \quad \widehat{\boldsymbol{y}}\_{\boldsymbol{L},i}(t\_{2}) \in \mathbb{R}^{n\_{\boldsymbol{L}}}, \quad n\_{\boldsymbol{F}} + n\_{\boldsymbol{L}} = n. \tag{36}
$$

byF,i ð Þ t<sup>2</sup> is the vector containing the fast-varying circuit variables at the slow time instant t1,i, and <sup>b</sup>yL,i ð Þ t<sup>2</sup> is the vector holding the slowly varying ones at the same slow time instant. The former will be integrated in t<sup>2</sup> with a small step size h2, while the later will be integrated with a much larger step size H2.

### 4.3. Circuit-block partitioning strategy

Since the partition of the electronic system in fast and slow subsystems may dynamically vary with time, it will be of great usefulness if the simulation algorithm is capable of automatically detecting the fast-varying and the slowly varying signals. The approach adopted in [16, 17] for that purpose is briefly described in the following.

in which K<sup>1</sup> is the total number of steps in t1. Now, using a backward differentiation formula (BDF) [1] to approximate the t<sup>1</sup> derivatives of the MPDAE (for simplicity let us consider here the Gear-2 rule [1]), we obtain for each slow time instant t1,i, from i ¼ 1 to i ¼ K1, the periodic

ð Þ <sup>t</sup><sup>2</sup> <sup>≃</sup> <sup>b</sup><sup>y</sup> <sup>t</sup>1,i ð Þ ; <sup>t</sup><sup>2</sup> . This means that, once <sup>b</sup>y<sup>i</sup>�<sup>1</sup>ð Þ <sup>t</sup><sup>2</sup> is known, the solution on the next slow

to be solved for obtaining the whole bivariate solution in the entire domain 0; tFinal ½ �� ½ � 0; T<sup>2</sup> . A very powerful technique has been proposed in the literature for solving each of the periodic boundary problems defined by (34) in an efficient way [16, 17]. This technique uses shooting-Newton based on MRK schemes. It splits the circuits into two distinct subsets according to the time evolution rates of their voltages and currents and performs time-step integration with different step lengths in each of the consecutive shooting iterations needed to solve (34). For that

> <sup>3</sup>qFðbyF,ið Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,ið Þ <sup>t</sup><sup>2</sup> Þ � <sup>4</sup>qFðbyA,i�<sup>1</sup>ð Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,i�<sup>1</sup>ð Þ <sup>t</sup><sup>2</sup> Þ þ qFð Þ <sup>b</sup>yF,i�<sup>2</sup>ð Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,i�<sup>2</sup>ð Þ <sup>t</sup><sup>2</sup> 2h1,i

> <sup>3</sup>qLðbyF,ið Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,ið Þ <sup>t</sup><sup>2</sup> Þ � <sup>4</sup>qLðbyF,i�<sup>1</sup>ð Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,i�<sup>1</sup>ð Þ <sup>t</sup><sup>2</sup> Þ þ qLð Þ <sup>b</sup>yF,i�<sup>2</sup>ð Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,i�<sup>2</sup>ð Þ <sup>t</sup><sup>2</sup> 2h1,i

> > ð Þ <sup>t</sup><sup>2</sup> <sup>∈</sup> <sup>R</sup>nF , <sup>b</sup>yL,i

ð Þ t<sup>2</sup> is the vector containing the fast-varying circuit variables at the slow time instant t1,i,

former will be integrated in t<sup>2</sup> with a small step size h2, while the later will be integrated with a

Since the partition of the electronic system in fast and slow subsystems may dynamically vary with time, it will be of great usefulness if the simulation algorithm is capable of automatically

ð Þ t<sup>2</sup> is the vector holding the slowly varying ones at the same slow time instant. The

ð Þ t<sup>2</sup> , is obtained by solving (34). Hence, a set of K<sup>1</sup> boundary value problems have

þ

<sup>d</sup>qð Þ <sup>b</sup>yið Þ <sup>t</sup><sup>2</sup> dt<sup>2</sup>

<sup>¼</sup> <sup>b</sup><sup>x</sup> <sup>t</sup>1,i ð Þ ; <sup>t</sup><sup>2</sup> ,

ð Þ <sup>t</sup><sup>2</sup> <sup>∈</sup> <sup>R</sup>nL , nF <sup>þ</sup> nL <sup>¼</sup> <sup>n</sup>: (36)

(34)

(35)

<sup>3</sup>qð Þ� <sup>b</sup>yið Þ <sup>t</sup><sup>2</sup> <sup>4</sup>qð Þþ <sup>b</sup>y<sup>i</sup>�<sup>1</sup>ð Þ <sup>t</sup><sup>2</sup> <sup>q</sup>ð Þ <sup>b</sup>y<sup>i</sup>�<sup>2</sup>ð Þ <sup>t</sup><sup>2</sup> 2h1,i

ð Þ T<sup>2</sup> ,

byi ð Þ¼ <sup>0</sup> <sup>b</sup>y<sup>i</sup>

(34) is firstly divided into the following coupled fast-slow subsystems

<sup>¼</sup> <sup>b</sup><sup>x</sup> <sup>t</sup>1,i ð Þ ; <sup>t</sup><sup>2</sup> ,

<sup>¼</sup> <sup>b</sup><sup>x</sup> <sup>t</sup>1,i ð Þ ; <sup>t</sup><sup>2</sup> ,

, <sup>b</sup>yF,i

boundary value problem defined by

34 Numerical Simulations in Engineering and Science

<sup>p</sup> <sup>b</sup>y<sup>i</sup> ð Þ t<sup>2</sup> � � <sup>þ</sup>

pF <sup>b</sup>yF,i

pL <sup>b</sup>yF,i

byi ð Þ¼ t<sup>2</sup>

much larger step size H2.

4.3. Circuit-block partitioning strategy

þ

þ

þ

þ

with

byF,i

and <sup>b</sup>yL,i

ð Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,i

<sup>d</sup>qAð Þ <sup>b</sup>yF,ið Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,ið Þ <sup>t</sup><sup>2</sup> dt<sup>2</sup>

ð Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,i

<sup>d</sup>qLð Þ <sup>b</sup>yF,ið Þ <sup>t</sup><sup>2</sup> ; <sup>b</sup>yL,ið Þ <sup>t</sup><sup>2</sup> dt<sup>2</sup>

� �

� �

ð Þ t<sup>2</sup>

ð Þ t<sup>2</sup>

byF,i ð Þ t<sup>2</sup>

" #

byL,i ð Þ t<sup>2</sup>

where <sup>b</sup>y<sup>i</sup>

time instant, <sup>b</sup>y<sup>i</sup>

As mentioned earlier, each of the periodic boundary value problems defined by (34) is solved using shooting-Newton based on multirate time-step integrators (MRK schemes). Now, taking into account that shooting is an iterative technique, the key idea of this partitioning strategy is to use a uni-rate scheme on the first shooting iterations needed for each time line t1,i � ½ � 0; T<sup>2</sup> of the 0; tFinal ½ �� ½ � 0; T<sup>2</sup> rectangular domain. For that, it starts by considering m ¼ 1 in the adopted MRK scheme, so that it degenerates into a standard uni-rate RK scheme. This means that all the unknowns (voltages or currents) will be integrated in t<sup>2</sup> with the same microstep h2. After that, the differential system of (34) is partitioned into (35) according to the variations in the t<sup>2</sup> derivatives of the unknowns. Each unknown that practically evidences no variations in its t<sup>2</sup> time derivatives for the entire time line t1,i � ½ � 0; T<sup>2</sup> , that is, each unknown that satisfies the condition

$$\left| \max \left( slope\_j \right) - \min \left( slope\_j \right) \right| < \left. \text{Tol} \right| \tag{37}$$

where and Tol is a small specified tolerance, will be treated as slow on the next shooting iterations. The remaining unknowns will be treated as fast. In a cheap scheme within a uniform grid each slope can be simply given by

$$\text{slope}\_{j} = \left[ \hat{y}\_{i} (t\_{2,j} + h\_{2}) - \hat{y}\_{i} (t\_{2,j}) \right] / h\_{2}, \quad t\_{2,j} \in \{0, h\_{2}, 2h\_{2}, \dots, T\_{2} - h\_{2}\}. \tag{38}$$

All the subsequent shooting iterations on the time line t1,i � ½ � 0; T<sup>2</sup> will be conducted in a multirate way using different step sizes. In a nonuniform grid, step sizes h<sup>2</sup> and H<sup>2</sup> ¼ m � h<sup>2</sup> may be chosen using any step-size control tool. In a uniform grid, they can be predefined or successively refined to achieve a desired accuracy.

The robustness of this partitioning strategy can be improved if more than one shooting iteration with the uni-rate scheme is considered for each time line t1,i � ½ � 0; T<sup>2</sup> . However, such tactic will obviously conduct to some extra computational cost, turning the overall algorithm to be a little less efficient.

Finally, to conclude this section, it must be highlighted that although efficiency gains are essentially dependent on the ratio between the number of slow and fast signals, significant reductions on the computational effort were reported in the scientific literature (e.g., [16, 17]) for the simulation of several illustrative examples of practical interest using the envelope following bivariate partitioned techniques.
