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

Abstract

Numerical simulation of complex and heterogeneous electronic systems can be a very challenging issue. Circuits composed of a combination of analog, mixed-signal and digital blocks or even radio frequency (RF) blocks, integrated in the same substrate, are very difficult to simulate as a whole at the circuit level. The main reason is because they contain a lot of state variables presenting very distinct properties and evolving in very widely separated time scales. Examples of practical interest are systems-on-a-chip (SoCs), very common in mobile electronics applications, as well as in many other embedded electronic systems. This chapter is intended to briefly review some advanced circuit-level numerical simulation techniques based on circuit-block partitioning schemes, which were especially designed to address the simulation challenges brought by this kind of circuits into the computer-aided-design (CAD) field.

1. Introduction

Electronic circuit simulation has emerged in the 1970s, triggered by the necessity of engineers having a tool to help in design and analysis of integrated circuits (ICs). Since probing internal nodes of semiconductor chips is extremely difficult, or even prohibitive in almost all cases, manufacturing integrated circuits having no help of a simulation tool would lead to an unbearable set of successive physical prototypes until a final solution was achieved. A simulation package combining device modeling and numerical simulation will help engineers to verify correctness and debug circuits during their design, avoiding physical prototyping and reducing product development expenses. Over the years, the continuous scaling of semiconductor devices and the increasing complexity of electronic architectures have been making CAD tools more and more important for circuits and systems designing. Demands for continuously providing new systems’ functionalities, lower-power consumptions or higher transmission rates (e.g., RF and microwave communication systems) are typical requisites that have led to a complex scenario of highly heterogeneous electronic systems. Conventional algorithms are not capable of simulating such kind of electronic systems in an efficient way. The justification for such ineffectiveness relies on the fact that standard simulation techniques do not perform any distinction between nodes or blocks within the circuits, treating all the variables in the same manner. This causes all the blocks of the circuit (analog, mixed-signal, digital or RF blocks) to be computed with the same numerical scheme, without taking their nature into consideration. To cope with this scenario, some advanced numerical simulation algorithms based on circuit-block partition have been proposed in recent years in the scientific literature. The most important ones are briefly reviewed in this chapter.

2. Review of basic circuit simulation concepts

2.1. Mathematical modeling of electronic systems

Dynamic behavior of electronic systems is modeled as systems of differential algebraic equations (DAEs) involving voltages, currents, charges and fluxes. These systems are usually obtained via nodal analysis, or modified nodal analysis (MNA), which consists of applying the Kirchhoff’s current law (KCL) to each electrical node and writing the branch currents in terms of the circuit node voltages using the corresponding constitutive relations to each circuit element. Such systems have, in general, the following form:

$$\mathit{p}\left(\mathit{y}\left(t\right)\right)+\frac{d\mathit{q}\left(\mathit{y}\left(t\right)\right)}{\mathit{dt}}=\mathit{x}\left(t\right),$$E1

in which $\mathit{x}\left(t\right)\in {\mathbb{R}}^n$ is the vector of independent stimuli (voltage or current sources) to the circuit, $\mathit{y}\left(t\right)\in {\mathbb{R}}^n$ is the vector of unknowns (voltages and currents waveforms) and n is the total number of unknowns. $\mathit{p}:{\mathbb{R}}^{\mathit{n}}\to {\mathbb{R}}^{\mathit{n}}$ describes the memoryless elements in the circuit (linear and nonlinear elements, as resistors, nonlinear voltage-controlled current sources, etc.), while $\mathit{q}:{\mathbb{R}}^{\mathit{n}}\to {\mathbb{R}}^{\mathit{n}}$ models all linear and nonlinear reactive circuit elements, as capacitors or inductors, represented as voltage-dependent electric charges or current-dependent magnetic fluxes, respectively.

Applying the chain differentiation rule to the reactive term of the DAE system of (1), we are able to obtain

$$\frac{d\mathit{q}\left(\mathit{y}\right)}{d\mathit{y}}\frac{d\mathit{y}\left(t\right)}{\mathit{dt}}=\mathit{x}\left(t\right)-\mathit{p}\left(\mathit{y}\left(t\right)\right),$$E2

or,

$$\mathit{M}\left[\mathit{y}\left(t\right)\right]\frac{d\mathit{y}\left(t\right)}{\mathit{dt}}=\mathit{x}\left(t\right)-\mathit{p}\left(\mathit{y}\left(t\right)\right),$$E3

where $\mathit{M}\left[\mathit{y}\left(t\right)\right]$ is habitually denoted as the mass matrix. If $\mathit{M}\left[\mathit{y}\left(t\right)\right]$ is invertible then it is possible to convert (3) into the following ordinary differential equations’ (ODE) system,

$$\frac{d\mathit{y}\left(t\right)}{\mathit{dt}}=\mathit{M}{\left[\mathit{y}\left(t\right)\right]}^{-1}\left(\mathit{x}\left(t\right)-\mathit{p}\left[\mathit{y}\left(t\right)\right]\right),$$E4

which can be rewritten in the classical form

$$\frac{d\mathit{y}\left(t\right)}{\mathit{dt}}=\mathit{f}\left(t,\mathit{y}\left(t\right)\right),$$E5

ordinarily utilized in the mathematical literature. When $\mathit{M}\left[\mathit{y}\left(t\right)\right]$ is singular, the DAE system of (3) will not degenerate into a ODE system, but it is often possible to express it as a set of algebraic equations combined with a set of differential equations of the form of (5).

2.2. Classic SPICE-like simulation

The most natural way of simulating the dynamic behavior of an electronic circuit is to numerically time-step integrate, in time domain, the DAE system, or the ODE system, modeling its operation. This means that the solution of (1), or (5), has to be computed over a specified time interval $\left[t_0,t_{\mathit{\text{Final}}}\right]$ from a specific initial condition $\mathit{y}\left(t_0\right)={\mathit{y}}_0$, leading to the so-called initial value problems

$$\mathit{p}\left(\mathit{y}\left(t\right)\right)+\frac{d\mathit{q}\left(\mathit{y}\left(t\right)\right)}{\mathit{dt}}=\mathit{x}\left(t\right),\mathit{y}\left(t_0\right)={\mathit{y}}_0,t_0\le t\le t_{\mathit{\text{Final}}},\mathit{y}\left(t\right)\in {\mathbb{R}}^n,$$E6

or

$$\frac{d\mathit{y}\left(t\right)}{\mathit{dt}}=\mathit{f}\left(t,\mathit{y}\left(t\right)\right),\mathit{y}\left(t_0\right)={\mathit{y}}_0,t_0\le t\le t_{\mathit{\text{Final}}},\mathit{y}\left(t\right)\in {\mathbb{R}}^n.$$E7

Computing the solution of (6), or (7), is frequently referred to as transient analysis and can be done by using initial value solvers, such as linear multistep methods (LMM) [1] or one-step methods (the popular Runge-Kutta (RK) schemes) [2, 3]. Either LMM or RK methods offer a large variety of explicit and implicit schemes, with very distinct properties in terms of order (accuracy) and numerical stability.

This time-step integration technique was used in the first digital computer programs of circuit analysis initially developed at the Electronics Research Laboratory of the University of California, Berkeley, in the early 1970s, and is still today the most widely used technique for such purpose. It is the core of all SPICE (which means Simulation Program with Integrated Circuit Emphasis) or SPICE-like computer programs, available in many commercial simulators.

2.3. Periodic steady-state simulation

As described earlier, SPICE-like simulation tools are primarily focused on transient analysis. However, in some cases electronics designers are essentially interested in obtaining circuits’ steady-state responses and not their transient regimes. The reason for that is because some specific properties of the circuits are better characterized, or simply only defined, in steady-state (e.g., impedance, voltage gain, current gain, harmonic distortion, signal to noise ratio, etc.).

SPICE tools are not adequate for computing steady-state responses of circuits presenting very different time constants, or high Q resonances, as is typically the case of RF and microwave circuits. This is so because they have to pass through the lengthy process of integrating all transients, and expecting them to vanish. Indeed, in such cases time-step integration can be extremely inefficient, since the number of discretization time steps used by the numerical integration scheme will be dramatically enormous. This is because the time interval over which the differential equations must be numerically integrated is set by the lowest frequency, or by how long the circuit takes to achieve steady-state, while the length of the time steps is constrained by the highest frequency component.

Periodic steady-state response of an electronic circuit is a regime where its unknowns are a set of generic waveforms (node voltages and branch currents) presenting a common period. Computing the periodic steady-state response, without having to first integrate all the transients, consists of finding the initial condition, $\mathit{y}\left(t_0\right)$, for the DAE, or ODE, system describing the circuit’s operation, such that the values of the unknowns at the end of one period T match their values at the beginning of that period, that is to say, $\mathit{y}\left(t_0\right)=\mathit{y}\left(t_0+T\right)$. These problems (evaluating the solution to a differential system that satisfies constraints at two or more distinct time instants) are referred to as boundary value problems. In this case, we have a periodic boundary value problem that can be formulated as:

$$\mathit{p}\left(\mathit{y}\left(t\right)\right)+\frac{d\mathit{q}\left(\mathit{y}\left(t\right)\right)}{\mathit{dt}}=\mathit{x}\left(t\right),\mathit{y}\left(t_0\right)=\mathit{y}\left(t_0+T\right),t_0\le t\le t_0+T,\mathit{y}\left(t\right)\in {\mathbb{R}}^n,$$E8

or, in the ODE form, as

$$\frac{d\mathit{y}\left(t\right)}{\mathit{dt}}=\mathit{f}\left(t,\mathit{y}\left(t\right)\right),\mathit{y}\left(t_0\right)=\mathit{y}\left(t_0+T\right),t_0\le t\le t_0+T,\mathit{y}\left(t\right)\in {\mathbb{R}}^n,$$E9

in which $\mathit{y}\left(t_0\right)=\mathit{y}\left(t_0+T\right)$ is known as the periodic boundary condition.

The most widely used techniques for computing the periodic steady-state solution of electronic circuits are briefly reviewed in the following: the shooting-Newton method [4] and the harmonic balance method [5, 6].

2.3.1. The shooting-Newton method

The shooting-Newton method [4] is the time-domain technique most commonly used by the electronic design automation (EDA) community for computing periodic steady-state solutions of electronic circuits. Shooting-Newton is an iterative technique that: (1) starts by computing the solution of the circuit for one period T using a LMM or RK integrator, considering some guessed initial condition (which is in general determined from a previous DC analysis); (2) then the computed solution at the end of the period is checked, and if it does not agree with the initial condition, the initial condition is wisely updated and (3) the circuit is then re-simulated for one period with the adjusted initial condition, and this process is repeated until the solution at the end of the period matches the initial condition.

Shooting is an iterative solver that uses an initial value technique to solve a boundary value problem. With the purpose of providing some technical details on the implementation of the shooting-Newton method, let us consider (8). Since with shooting, we perform numerical time-step integration of the DAE system from $t=t_0$ until $t=t_0+T$, the main difficulty is that we do not know a priori for which initial condition $\mathit{y}\left(t_0\right)$ has to be considered that will lead to the steady-state solution, that is, that will satisfy the periodic boundary condition $\mathit{y}\left(t_0\right)=\mathit{y}\left(t_0+T\right)$. Thus, the key aspect of shooting-Newton relies on finding the solution of

$$\mathit{y}\left(t_0\right)=\mathit{y}\left(t_0+T\right)\iff \mathit{y}\left(t_0\right)-\mathit{y}\left(t_0+T\right)=0.$$E10

Let us now define $\varphi \left(\mathit{y}\left(t_0\right),T\right)=\mathit{y}\left(t_0+T\right)$, where $\varphi$ is known as the state-transition function [4, 7], and rewrite (10) as

$$\varphi \left(\mathit{y}\left(t_0\right),T\right)-\mathit{y}\left(t_0\right)=0.$$E11

Although electronic circuits may operate in strongly nonlinear regimes, their state-transition functions are often moderately nonlinear (or even quite linear). This means that slight perturbations on the initial condition (starting state) produce almost proportional perturbations in the subsequent time states. Taking this aspect into account, it is straightforward to conclude that (11) can be iteratively solved in an efficient way with the Newton’s method, which will lead to the following iterative scheme:

$$\varphi \left({\mathit{y}}^{\left[r\right]}\left(t_0\right),T\right)-{\mathit{y}}^{\left[r\right]}\left(t_0\right)+{\left.\left[\frac{\partial \varphi \left({\mathit{y}}^{\left[r\right]}\left(t_0\right),T\right)}{\partial \mathit{y}\left(t_0\right)}-I\right]\right|}_{\mathit{y}\left(t_0\right)={\mathit{y}}^{\left[r\right]}\left(t_0\right)}\left[{\mathit{y}}^{\left[r+1\right]}\left(t_0\right)-{\mathit{y}}^{\left[r\right]}\left(t_0\right)\right]=0,$$E12

where $I$ is the $n\times n$ identity matrix. The cost of the solution of the linear system of (12) is dominated by the computational effort required to evaluate the derivative of the state-transition function (usually referred to as the sensitivity matrix). This matrix is computed taking into account the chain differentiation rule, that is, taking into consideration that $\varphi \left(\mathit{y}\left(t_0\right),T\right)$ is, in fact, the numerical vector ${\mathit{y}}_K$, with $K$ being the total number of time steps in the interval $\left[t_0,t_0+T\right]$, which depends on the previous step value ${\mathit{y}}_{K-1}$, which, itself, depends on ${\mathit{y}}_{K-2}$, and so forth. The sensitivity matrix is then given by

$$\frac{\partial \varphi \left(\mathit{y}\left(t_0\right),T\right)}{\partial \mathit{y}\left(t_0\right)}=\frac{\partial {\mathit{y}}_K}{\partial {\mathit{y}}_{K-1}}\cdot \frac{\partial {\mathit{y}}_{K-1}}{\partial {\mathit{y}}_{K-2}}\cdot \cdots \cdot \frac{\partial {\mathit{y}}_1}{\partial {\mathit{y}}_0}.$$E13

Although solving (12) and computing the sensitivity matrix (13) involve substantial computational cost, shooting-Newton converges to the steady-state solution much faster than classic time-step integration. This is the reason why it is the time-domain steady-state engine most widely used in the circuit simulation field.

2.3.2. The harmonic balance method

The harmonic balance (HB) method [4, 5, 6] is a frequency domain steady-state simulation technique widely used by the EDA community. In contrast to time domain tools, which represent waveforms as a set of time samples, frequency domain techniques represent periodic signals using coefficients in a sum of complex exponentials (or sines and cosines) harmonically related. The main advantage of HB over time-domain techniques (e.g., shooting-Newton) is that it can represent steady-state solutions (voltage and current waveforms) very accurately using a small number of coefficients. This is especially evident for moderately nonlinear circuits excited by smooth waveforms, in which significant reductions in the computational cost are achieved when HB is used as the simulation tool. However, it must be noted that HB is not suitable for dealing with strongly nonlinear regimes producing waveforms with sharps transitions. In such case, the large number of terms required in the Fourier series expansions will make HB very inefficient.

With the purpose of providing a brief and intuitive explanation of the HB method, let us consider (8). For achieving simplicity in our formulation, let us consider a very small circuit driven by a single periodic source $x\left(t\right)\in \mathbb{R}$ satisfying $x\left(t\right)=x\left(t+T\right)$ and whose dynamic behavior is described by a unique unknown, $y\left(t\right)\in \mathbb{R}$ (the generalization to $x\left(t\right)\in {\mathbb{R}}^n,y\left(t\right)\in {\mathbb{R}}^n$ is straightforward). Given that the steady-state regime of the circuit will be periodic with the same period $T$, both the stimulus and the steady-state solution can be represented as the Fourier series

$$x\left(t\right)=\sum _{k=-\infty }^{+\infty }{X}_k{e}^{\mathit{jk}{\omega }_{0}t},y\left(t\right)=\sum _{k=-\infty }^{+\infty }{Y}_k{e}^{\mathit{jk}{\omega }_{0}t},$$E14

in which ${\omega }_0=2\pi /T$ is the fundamental frequency. If we substitute (14) into (8), and consider an appropriate harmonic truncation at some order $k=K$, we will obtain

$$p\left(\sum _{k=-K}^{+K}{Y}_k{e}^{\mathit{jk}{\omega }_{0}t}\right)+\frac{d}{\mathit{dt}}\left[q\left(\sum _{k=-K}^{+K}{Y}_k{e}^{\mathit{jk}{\omega }_{0}t}\right)\right]=\sum _{k=-K}^{+K}{X}_k{e}^{\mathit{jk}{\omega }_{0}t}.$$E15

The HB method converts the differential problem of (8) into the frequency domain, obtaining the $\left(2K+1\right)$ algebraic equations system

$$\mathbf{F}\left(\mathbf{Y}\right)=\mathbf{P}\left(\mathbf{Y}\right)+j\mathrm{\Omega }\mathbf{Q}\left(\mathbf{Y}\right)-\mathbf{X}=0$$E16

where $\mathbf{Y}={\left[Y_{-K},\dots ,Y_0,\dots ,Y_K\right]}^T$, $\mathbf{X}={\left[X_{-K},\dots ,X_0,\dots ,X_K\right]}^T$, $j\mathrm{\Omega }=diag\left(-\mathit{jK}{\omega }_0,\dots ,0,\dots ,\mathit{jK}{\omega }_0\right)$.

and $\mathbf{P}$, $\mathbf{Q}$ are vectors containing the Fourier coefficients of $p\left(y\left(t\right)\right)$ and $q\left(y\left(t\right)\right)$, respectively. The algebraic equations system of (16) is usually denoted as the harmonic balance system, which can be iteratively solved using the Newton method

$$\mathbf{F}\left({\mathbf{Y}}^{\left[r\right]}\right)+{\left.\frac{d\mathbf{F}\left(\mathbf{Y}\right)}{d\mathbf{Y}}\right|}_{\mathbf{Y}={\mathbf{Y}}^{\left[r\right]}}\left[{\mathbf{Y}}^{\left[r+1\right]}-{\mathbf{Y}}^{\left[r\right]}\right]=0,$$E17

in which

$$\frac{d\mathbf{F}\left(\mathbf{Y}\right)}{d\mathbf{Y}}=\mathbf{J}\left(\mathbf{Y}\right)=\mathbf{G}\left(\mathbf{Y}\right)+j\mathrm{\Omega }\mathbf{C}\left(\mathbf{Y}\right)$$E18

is the $\left(2K+1\right)\times \left(2K+1\right)$ composite conversion matrix, known as the Jacobian matrix of the error function $\mathbf{F}\left(\mathbf{Y}\right)$. $\mathbf{G}$ and $\mathbf{C}$ denote the $\left(2K+1\right)\times \left(2K+1\right)$ conversion matrices (Toeplitz) [7] corresponding to $g\left(y\left(t\right)\right)=\mathit{dp}\left(y\left(t\right)\right)/\mathit{dy}$ and $c\left(y\left(t\right)\right)=\mathit{dq}\left(y\left(t\right)\right)/\mathit{dy}$.

3. Univariate time-domain partitioned simulation engines

3.1. Time-domain latency

As highlighted in the Introduction of this chapter, dynamic regimes of operation of some electronic systems may involve signals (voltages and currents) presenting widely distinct time evolution rates. Typical examples of that are coupled analog-digital systems or combined technologies of baseband analog, digital and RF blocks, in the same circuit. In these examples, very fast signals and slowly varying signals cohabit in the same framework. This property, the one of having in the same problem signals presenting rapid time rates of change, while others evolve in a very slow way (or remain approximately constant within a certain time window) is usually denoted as time-domain latency. It must be pointed out that time-domain latency is a phenomenon that is not restricted to heterogeneous circuits. For instance, regimes of operation of pure digital electronic systems typically present a set of variables remaining practically constant within a specific time interval, while others evidence quick variations (fast transitions) in that interval.

In Section 2, we have seen that SPICE-like simulation engines (which are based on time-step integration schemes) are widely used for computing the numerical solution of electronic systems. However, when dealing with circuits presenting time-domain latency, that is, containing node voltages and brunch currents evolving at very different rates, traditional SPICE simulators become inefficient because they expend unnecessary work on the computation of the slowly varying components. This is so because classic initial value solvers (as LMM or RK schemes) integrate all the DAE, or ODE, unknowns with the same step size.

3.2. Partitioned algorithms for time-step integration

To deal with the earlier described time-domain latency in an effective way, some modern partitioned algorithms for univariate time-step integration have been proposed in the recent years in the scientific literature [8, 9, 10]. These powerful techniques, denoted as multirate Runge-Kutta (MRK) schemes, split the ODE system of (5) into coupled fast and slow (latent) subsystems, obtaining

$$\begin{array}{l}\frac{d{\mathit{y}}_F\left(t\right)}{\mathit{dt}}={\mathit{f}}_F\left(t,{\mathit{y}}_F,{\mathit{y}}_L\right),{\mathit{y}}_F\left(t_0\right)={\mathit{y}}_{F,0}\\ \frac{d{\mathit{y}}_L\left(t\right)}{\mathit{dt}}={\mathit{f}}_L\left(t,{\mathit{y}}_F,{\mathit{y}}_L\right),{\mathit{y}}_L\left(t_0\right)={\mathit{y}}_{L,0},\end{array}$$E19

with

$$\mathit{y}=\left[\begin{array}{l}{\mathit{y}}_F\\ {\mathit{y}}_L\end{array}\right],{\mathit{y}}_F\in {\mathbb{R}}^{n_F},{\mathit{y}}_L\in {\mathbb{R}}^{n_L},n_F+n_L=n,$$E20

where ${\mathit{y}}_F$ is the vector containing the fast-varying signals and ${\mathit{y}}_L$ is the vector holding the slowly varying (latent) ones. The former will be integrated with a small step length $h$ (microstep), whereas the latter will be integrated with a large step size $H$ (macrostep). The number of microsteps within each macrostep is an integer that we will denote by $m$. Hence, $H=m\cdot h$. The generic formulation of a MRK scheme is given in the following [8, 9].

Consider two Runge-Kutta methods of $s$ and $\overline{s}$ stages represented by their Butcher tableaus [2] $\left(\mathit{b},\mathit{A},\mathit{c}\right)$ and $\left(\overline{\mathit{b}},\overline{\mathit{A}},\overline{\mathit{c}}\right)$,

E21

The MRK method conceived for efficiently computing the numerical solution of the partitioned differential system of (19), using a microstep $h$ for integrating the fast unknowns and a macrostep $H$ for integrating the latent unknowns, is defined as follows [8, 9].

The fast-varying vector is obtained by

$$\begin{array}{l}{\mathit{y}}_F\left(t_0+\mathit{\lambda h}+h\right)\simeq {\mathit{y}}_{F,\lambda +1}={\mathit{y}}_{F,\lambda }+h\sum _{i=1}^s{b}_i{\mathit{k}}_{F,i}^{\lambda },\\ \lambda =0,1,\dots ,m-1,\end{array}$$E22

$$\begin{array}{l}{\mathit{k}}_{F,i}^{\lambda }={\mathit{f}}_F\left(t_0+\mathit{\lambda h}+c_{i}h,{\mathit{y}}_{F,\lambda }+h\sum _{j=1}^s{a}_{\mathit{ij}}{\mathit{k}}_{F,j}^{\lambda },{\tilde{\mathit{Y}}}_{L,i}^{\lambda }\right),\\ i=1,2,\dots s,\end{array}$$E23

with ${\tilde{\mathit{Y}}}_{L,i}^{\lambda }\simeq {\mathit{y}}_L\left(t_0+\mathit{\lambda h}+c_{i}h\right).$

The slowly varying vector is given by

$${\mathit{y}}_L\left(t_0+H\right)\simeq {\mathit{y}}_{L,1}={\mathit{y}}_{L,0}+H\sum _{i=1}^{\overline{s}}{\overline{b}}_i{\mathit{k}}_{L,i},$$E24

$$\begin{array}{l}{\mathit{k}}_{L,i}={\mathit{f}}_L\left(t_0+{\overline{c}}_{i}H,{\tilde{\mathit{Y}}}_{F,i},{\mathit{y}}_{L,0}+H\sum _{j=1}^{\overline{s}}{\overline{a}}_{\mathit{ij}}{\mathit{k}}_{L,j}\right),\\ i=1,2,\dots ,\overline{s},\end{array}$$E25

with ${\tilde{\mathit{Y}}}_{F,i}\simeq {\mathit{y}}_F\left(t_0+{\overline{c}}_{i}H\right)$.

From this definition we can attest that numerical coupling between fast and slow differential subsystems is achieved by ${\tilde{\mathit{Y}}}_{L,i}^{\lambda }$ and ${\tilde{\mathit{Y}}}_{F,i}$. Effective stratagems for computing these intermediate stage values are proposed in [8, 9].

3.3. Circuit-block partitioning strategy

We must be aware that circuit-block partition into fast and slow subsystems can possibly change during the time-step integration process. Hence, it will be very helpful if the simulator is capable of automatically detect the slow and fast variables in the circuit. This automatic classification can be achieved using embedded RK methods [2] and error estimates usually evaluated for step size control and stiffness detection [3]. Fast-slow partitioning strategies, step size control tools, number of microsteps within a macrostep, stiffness detection stratagems and many other technical details of MRK code implementation are thoroughly addressed in [8, 9, 10] and also in [11].

Finally, it must be pointed out that significant efficiency gains in computation speed for the simulation of several illustrative examples have been reported in the scientific literature, which demonstrate an effective reduction on the MRK computational cost in comparison to traditional SPICE simulation engines.

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, 14, 15, 16, 17, 18, 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

$$x\left(t\right)=e\left(t\right)cos\left({\omega }_{C}t+\varphi \left(t\right)\right),$$E26

where $e\left(t\right)$ and $\varphi \left(t\right)$ are, respectively, the amplitude and phase slowly varying baseband signals (envelope), modulating the $cos\left({\omega }_{C}t\right)$ 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{x}\left(t_1,t_2\right)=e\left(t_1\right)cos\left({\omega }_C{t}_2+\varphi \left(t_1\right)\right),$$E27

where $t_1$ is the slow baseband artificial time scale and $t_2$ is the fast carrier artificial time scale. It must be noted that $\widehat{x}\left(t_1,t_2\right)$ is a periodic function with respect to $t_2$ but not to $t_1$, that is,

$$\widehat{x}\left(t_1,t_2\right)=\widehat{x}\left(t_1,t_2+T_2\right),T_2=2\pi /{\omega }_C.$$E28

This means that a generic 1-D $\left[0,t_{\mathit{\text{Final}}}\right]$ time interval will be mapped into a 2-D $\left[0,t_{\mathit{\text{Final}}}\right]\times \left[0,T_2\right]$ rectangular domain. It is easy to attest that, in general, the number of points required to represent $\widehat{x}\left(t_1,t_2\right)$ in $\left[0,t_{\mathit{\text{Final}}}\right]\times \left[0,T_2\right]$ is much less than the number of points needed to represent $x\left(t\right)$ in $\left[0,t_{\mathit{\text{Final}}}\right]$. This is especially evident when the $t_1$ and $t_2$ 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 $\mathit{x}\left(t\right)$ and $\mathit{y}\left(t\right)$: $t$ is replaced by $t_1$ in the slowly varying entities (envelope time scale) and $t$ is replaced by $t_2$ 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

$$\mathit{p}\left(\widehat{\mathit{y}}\left(t_1,t_2\right)\right)+\frac{\partial \mathit{q}\left(\widehat{\mathit{y}}\left(t_1,t_2\right)\right)}{\partial t_1}+\frac{\partial \mathit{q}\left(\widehat{\mathit{y}}\left(t_1,t_2\right)\right)}{\partial t_2}=\widehat{\mathit{x}}\left(t_1,t_2\right),$$E29

usually denoted as multirate partial differential algebraic equations’ (MPDAE) system [13, 14]. It is easy to demonstrate that, if $\widehat{\mathit{x}}\left(t_1,t_2\right)$ and $\widehat{\mathit{y}}\left(t_1,t_2\right)$ satisfy (29), then the univariate forms $\mathit{x}\left(t\right)=\widehat{\mathit{x}}\left(t,t\right)$ and $\mathit{y}\left(t\right)=\widehat{\mathit{y}}\left(t,t\right)$ satisfy (1) [13]. Thus, $\mathit{y}\left(t\right)$ may be retrieved from its bivariate form $\widehat{\mathit{y}}\left(t_1,t_2\right)$ by simply setting $t_1=t_2=t$, meaning that univariate solutions of (1) are available on diagonal lines $t_1=t,t_2=t,$ along the bivariate solutions of (29). In truth, attending to the periodicity of the problem in the $t_2$ time dimension, the univariate version of the vector containing the circuit responses is obtained from its bivariate representation on the rectangular domain $\left[0,t_{\mathit{\text{Final}}}\right]\times \left[0,T_2\right]$ as

$$\mathit{y}\left(t\right)=\widehat{\mathit{y}}\left(t,tmod{T}_2\right),$$E30

where $tmod{T}_2$ represents the remainder of division of $t$ by $T_2$.

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

$$\mathit{p}\left(\widehat{\mathit{y}}\left(t_1,t_2,\dots ,t_m\right)\right)+\frac{\partial \mathit{q}\left(\widehat{\mathit{y}}\left(t_1,t_2,\dots ,t_m\right)\right)}{\partial t_1}+\cdots +\frac{\partial \mathit{q}\left(\widehat{\mathit{y}}\left(t_1,t_2,\dots ,t_m\right)\right)}{\partial t_m}=\widehat{\mathit{x}}\left(t_1,t_2,\dots ,t_m\right),$$E31

and the univariate solution, $\mathit{y}\left(t\right)$, may be recovered from its multivariate form, $\widehat{\mathit{y}}\left(t_1,t_2,\dots ,t_m\right)$, by setting $t_1=t_2=\cdots =t_m=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, 14, 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, 17, 18, 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{array}{l}\mathit{p}\left(\widehat{\mathit{y}}\left(t_1,t_2\right)\right)+\frac{\partial \mathit{q}\left(\widehat{\mathit{y}}\left(t_1,t_2\right)\right)}{\partial t_1}+\frac{\partial \mathit{q}\left(\widehat{\mathit{y}}\left(t_1,t_2\right)\right)}{\partial t_2}=\widehat{\mathit{x}}\left(t_1,t_2\right)\\ \widehat{\mathit{y}}\left(0,t_2\right)=\mathit{g}\left(t_2\right)\\ \widehat{\mathit{y}}\left(t_1,0\right)=\widehat{\mathit{y}}\left(t_1,T_2\right)\end{array}$$E32

defined on the rectangular domain $\left[0,t_{\mathit{\text{Final}}}\right]\times \left[0,T_2\right]$. $\mathit{g}\left(\cdot \right)$ is a vector-valued initial-condition satisfying $\mathit{g}\left(0\right)=\mathit{g}\left(T_2\right)=\mathit{y}\left(0\right)$, and $\widehat{\mathit{y}}\left(t_1,0\right)=\widehat{\mathit{y}}\left(t_1,T_2\right)$ is the periodic boundary condition due to the periodicity of the problem in the $t_2$ fast time scale. In order to solve this initial-boundary value problem let us begin to consider the following semi-discretization of $\left[0,t_{\mathit{\text{Final}}}\right]\times \left[0,T_2\right]$ in the $t_1$ slow time dimension defined by

$$0=t_{1,0}<t_{1,1}<\cdots <t_{1,i-1}<t_{1,i}<\cdots <t_{1,K_1}=t_{\mathit{\text{Final}}},h_{1,i}=t_{1,i}-t_{1,i-1},$$E33

in which $K_1$ is the total number of steps in $t_1$. Now, using a backward differentiation formula (BDF) [1] to approximate the $t_1$ derivatives of the MPDAE (for simplicity let us consider here the Gear-2 rule [1]), we obtain for each slow time instant $t_{1,i}$, from $i=1$ to $i=K_1,$ the periodic boundary value problem defined by

$$\begin{array}{l}\mathit{p}\left({\widehat{\mathit{y}}}_i\left(t_2\right)\right)+\frac{3\mathit{q}\left({\widehat{\mathit{y}}}_i\left(t_2\right)\right)-4\mathit{q}\left({\widehat{\mathit{y}}}_{i-1}\left(t_2\right)\right)+\mathit{q}\left({\widehat{\mathit{y}}}_{i-2}\left(t_2\right)\right)}{2h_{1,i}}+\frac{d\mathit{q}\left({\widehat{\mathit{y}}}_i\left(t_2\right)\right)}{{\mathit{dt}}_2}=\widehat{\mathit{x}}\left(t_{1,i},t_2\right),\\ {\widehat{\mathit{y}}}_i\left(0\right)={\widehat{\mathit{y}}}_i\left(T_2\right),\end{array}$$E34

where ${\widehat{\mathit{y}}}_i\left(t_2\right)\simeq \widehat{\mathit{y}}\left(t_{1,i},t_2\right)$. This means that, once ${\widehat{\mathit{y}}}_{i-1}\left(t_2\right)$ is known, the solution on the next slow time instant, ${\widehat{\mathit{y}}}_i\left(t_2\right)$, is obtained by solving (34). Hence, a set of $K_1$ boundary value problems have to be solved for obtaining the whole bivariate solution in the entire domain $\left[0,t_{\mathit{\text{Final}}}\right]\times \left[0,T_2\right]$. 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{array}{l}{\mathit{p}}_{\mathit{F}}\left({\widehat{\mathit{y}}}_{F,i}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i}\left(t_2\right)\right)\hfill \\ +\frac{3{\mathit{q}}_{\mathit{F}}\left({\widehat{\mathit{y}}}_{F,i}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i}\left(t_2\right)\right)-4{\mathit{q}}_{\mathit{F}}\left({\widehat{\mathit{y}}}_{A,i-1}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i-1}\left(t_2\right)\right)+{\mathit{q}}_{\mathit{F}}\left({\widehat{\mathit{y}}}_{F,i-2}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i-2}\left(t_2\right)\right)}{2h_{1,i}}\hfill \\ +\frac{d{\mathit{q}}_{\mathit{A}}\left({\widehat{\mathit{y}}}_{F,i}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i}\left(t_2\right)\right)}{{\mathit{dt}}_2}=\widehat{\mathit{x}}\left(t_{1,i},t_2\right),\hfill \\ {\mathit{p}}_{\mathit{L}}\left({\widehat{\mathit{y}}}_{F,i}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i}\left(t_2\right)\right)\hfill \\ +\frac{3{\mathit{q}}_{\mathit{L}}\left({\widehat{\mathit{y}}}_{F,i}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i}\left(t_2\right)\right)-4{\mathit{q}}_{\mathit{L}}\left({\widehat{\mathit{y}}}_{F,i-1}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i-1}\left(t_2\right)\right)+{\mathit{q}}_{\mathit{L}}\left({\widehat{\mathit{y}}}_{F,i-2}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i-2}\left(t_2\right)\right)}{2h_{1,i}}\hfill \\ +\frac{d{\mathit{q}}_{\mathit{L}}\left({\widehat{\mathit{y}}}_{F,i}\left(t_2\right),{\widehat{\mathit{y}}}_{L,i}\left(t_2\right)\right)}{{\mathit{dt}}_2}=\widehat{\mathit{x}}\left(t_{1,i},t_2\right),\hfill \end{array}$$E35

with

$${\widehat{\mathit{y}}}_i\left(t_2\right)\mathit{=}\left[\begin{array}{l}{\widehat{\mathit{y}}}_{F,i}\left(t_2\right)\\ {\widehat{\mathit{y}}}_{L,i}\left(t_2\right)\end{array}\right],{\widehat{\mathit{y}}}_{F,i}\left(t_2\right)\in {\mathbb{R}}^{n_F},{\widehat{\mathit{y}}}_{L,i}\left(t_2\right)\in {\mathbb{R}}^{n_L},n_F+n_L=n.$$E36

${\widehat{\mathit{y}}}_{F,i}\left(t_2\right)$ is the vector containing the fast-varying circuit variables at the slow time instant $t_{1,i}$, and ${\widehat{\mathit{y}}}_{L,i}\left(t_2\right)$ is the vector holding the slowly varying ones at the same slow time instant. The former will be integrated in $t_2$ with a small step size h₂, while the later will be integrated with a much larger step size H₂.

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.

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 $t_{1,i}\times \left[0,T_2\right]$ of the $\left[0,t_{\mathit{\text{Final}}}\right]\times \left[0,T_2\right]$ 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_2$ with the same microstep $h_2$. After that, the differential system of (34) is partitioned into (35) according to the variations in the $t_2$ derivatives of the unknowns. Each unknown that practically evidences no variations in its $t_2$ time derivatives for the entire time line $t_{1,i}\times \left[0,T_2\right]$, that is, each unknown that satisfies the condition

$$\left|max\left({\mathit{\text{slope}}}_j\right)-min\left({\mathit{\text{slope}}}_j\right)\right|<\mathit{Tol},$$E37

where and $\mathit{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

$${\mathit{\text{slope}}}_j=\left[{\widehat{y}}_i\left(t_{2,j}+h_2\right)-{\widehat{y}}_i\left(t_{2,j}\right)\right]/h_2,t_{2,j}\in \left\{0,h_2,2h_2,\dots ,T_2-h_2\right\}.$$E38

All the subsequent shooting iterations on the time line $t_{1,i}\times \left[0,T_2\right]$ will be conducted in a multirate way using different step sizes. In a nonuniform grid, step sizes $h_2$ and $H_2=m\cdot h_2$ 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 $t_{1,i}\times \left[0,T_2\right]$. 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.

5. Hybrid time-frequency partitioned techniques

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 addresses some basic details of the simulation algorithms and finally, Section 5.3 describes the approach for automatic circuit-block partition.

5.1. Multivariate hybrid envelope following

Let us once again consider the initial-boundary value problem of (32) characterizing the bivariate nature of an electronic circuit operating at two widely separated time scales. Let us also consider the semi-discretization of the $\left[0,t_{\mathit{\text{Final}}}\right]\times \left[0,T_2\right]$ rectangular domain, which has led to the set of periodic boundary value problems defined by (34). When waveforms are not excessively demanding on the quantity of harmonic components needed for an accurate frequency domain representation, we are able to use the HB method for efficiently computing the solution of (34). In such case we will obtain for each time line $t_{1,i}\times \left[0,T_2\right]$ the following HB system

$$\mathbf{P}\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)+\frac{3\mathbf{Q}\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)-4\mathbf{Q}\left(\widehat{\mathbf{Y}}\left(t_{1,i-1}\right)\right)+\mathbf{Q}\left(\widehat{\mathbf{Y}}\left(t_{1,i-2}\right)\right)}{2h_{1,i}}+j\mathbf{\Omega }\mathbf{Q}\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)=\widehat{\mathbf{X}}\left(t_{1,i}\right),$$E39

in which $\widehat{\mathbf{X}}\left(t_{1,i}\right)$ is a $\left(2K+1\right)\times n$ vector containing the Fourier coefficients of the stimuli (independent sources) and $\widehat{\mathbf{Y}}\left(t_{1,i}\right)$ a $\left(2K+1\right)\times n$ vector containing the Fourier coefficients of the unknowns (node voltages and brunch currents’ waveforms), at $t_1=t_{1,i}$. $K$ is the maximum harmonic order and $n$ is the number of unknowns. In order to obtain the solution in the entire $\left[0,t_{\mathit{\text{Final}}}\right]\times \left[0,T_2\right]$ rectangular domain a total of $K_1$ HB systems have to be solved.

With this hybrid time-frequency approach we have an envelope following technique that handles the slow variations of the unknowns in the time domain and their fast variations in the frequency domain. This technique, usually denoted as multivariate envelope transient harmonic balance [20, 21], is able to exploit the existence of time-rate disparities in circuits operating under moderately nonlinear regimes of operation.

5.2. Partitioned algorithms for hybrid envelope following computation

The hybrid time-frequency envelope following technique presented earlier does not perform any distinction between nodes or blocks in the circuit. Although it has been widely used, especially by the RF and microwave community, only a few years ago other versions regarding circuit-block partition were proposed [18, 19] in way to take into account possible circuit’s heterogeneities. The most important aspects of these partitioning techniques are briefly reviewed in the following.

Let us rewrite (39) as

$$\begin{array}{l}\mathbf{F}\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)=\mathbf{P}\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)+\frac{3\mathbf{Q}\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)-4\mathbf{Q}\left(\widehat{\mathbf{Y}}\left(t_{1,i-1}\right)\right)+\mathbf{Q}\left(\widehat{\mathbf{Y}}\left(t_{1,i-2}\right)\right)}{2h_{1,i}}\\ +j\mathbf{\Omega }\mathbf{Q}\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)-\widehat{\mathbf{X}}\left(t_{1,i}\right)=0.\end{array}$$E40

The Newton iterative solver is usually utilized to solve (40), leading to

$$\mathbf{F}\left({\widehat{\mathbf{Y}}}^{\left[r\right]}\left(t_{1,i}\right)\right)+{\left.\frac{d\mathbf{F}\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)}{d\widehat{\mathbf{Y}}\left(t_{1,i}\right)}\right|}_{\widehat{\mathbf{Y}}\left(t_{1,i}\right)={\widehat{\mathbf{Y}}}^{\left[r\right]}\left(t_{1,i}\right)}\left[{\widehat{\mathbf{Y}}}^{\left[r+1\right]}\left(t_{1,i}\right)-{\widehat{\mathbf{Y}}}^{\left[r\right]}\left(t_{1,i}\right)\right]=0,$$E41

which implies that for computing ${\widehat{\mathbf{Y}}}^{\left[r+1\right]}\left(t_{1,i}\right)$ from the previous estimate ${\widehat{\mathbf{Y}}}^{\left[r\right]}\left(t_{1,i}\right)$ we have to solve a linear system composed of $n\times \left(2K+1\right)$ equations. The system of (41) involves the so-called Jacobian matrix of $\mathbf{F}\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)$, which has a block structure, consisting of an $n\times n$ matrix of square submatrices (blocks), each one with dimensions $\left(2K+1\right)\times \left(2K+1\right)$. The general block of row $m$ and column $l$ can be expressed as

$$\begin{array}{l}\frac{d{\mathbf{F}}_m\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)}{d{\widehat{\mathbf{Y}}}_l\left(t_{1,i}\right)}=\frac{d{\mathbf{P}}_m\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)}{d{\widehat{\mathbf{Y}}}_l\left(t_{1,i}\right)}+\frac{1}{h_{1.i}}\frac{d{\mathbf{Q}}_m\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)}{d{\widehat{\mathbf{Y}}}_l\left(t_{1,i}\right)}\\ +j\mathrm{\Omega }\frac{d{\mathbf{Q}}_m\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)}{d{\widehat{\mathbf{Y}}}_l\left(t_{1,i}\right)},l,m=1,2,\dots ,n,\end{array}$$E42

where $d{\mathbf{P}}_m\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)/d{\widehat{\mathbf{Y}}}_l\left(t_{1,i}\right)$ and $d{\mathbf{Q}}_m\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)/d{\widehat{\mathbf{Y}}}_l\left(t_{1,i}\right)$ denote the Toeplitz matrices [7] of the vectors containing the Fourier coefficients of ${\mathit{dp}}_m\left(\widehat{y}\left(t_{1,i},t_2\right)\right)/d{\widehat{y}}_l\left(t_{1,i},t_2\right)$ and ${\mathit{dq}}_m\left(\widehat{y}\left(t_{1,i},t_2\right)\right)/d{\widehat{y}}_l\left(t_{1,i},t_2\right)$, respectively.

A very powerful technique has been proposed in the literature [18, 19] for solving (41) in an efficient way. This technique takes into account that, in some cases (e.g., RF heterogeneous systems), there are parts of the circuits in which there are no fluctuations dictated by the fast carrier. As a consequence, bidimensional forms $\widehat{y}\left(t_1,t_2\right)$ of voltages and currents in those parts have no dependence on $t_2$. This means that for each time line $t_{1,i}\times \left[0,T_2\right]$ each signal $\widehat{y}\left(t_{1,i},t_2\right)$ has a constant value that can be represented as a Fourier series with a unique $k=0$ coefficient. Thus, while quickly varying signals are represented in the frequency domain by a set of $\left(2K+1\right)$ Fourier coefficients, signals in latent (slowly varying) blocks are represented by a single coefficient. Having this in mind it is straightforward to conclude that the size of vector $\widehat{\mathbf{Y}}\left(t_{1,i}\right)$ will be significantly reduced, as well as the size (number of equations) of the HB system. Significant Jacobian matrix size reductions will also be achieved, since some of the submatrices will no longer be $\left(2K+1\right)\times \left(2K+1\right)$ matrices, but, instead, simple $1\times 1$ scalar elements. For instance, let us suppose that the $m\mathrm{th}$ component of (40) is a slowly varying signal exclusively dependent on other slowly varying signals. Let us also suppose that the $l\mathrm{th}$ component is also a slowly varying signal. Then, the Jacobian matrix block of row $m$ and column $l$ will be a $1\times 1$ block given by

$$\frac{d{\mathbf{F}}_m\left(\widehat{\mathbf{Y}}\left(t_{1,i}\right)\right)}{d{\widehat{\mathbf{Y}}}_l\left(t_{1,i}\right)}=\frac{{\mathit{dp}}_m\left(\widehat{y}\left(t_{1,i}\right)\right)}{d{\widehat{y}}_l\left(t_{1,i}\right)}+\frac{1}{h_{1,i}}\frac{{\mathit{dq}}_m\left(\widehat{y}\left(t_{1,i}\right)\right)}{d{\widehat{y}}_l\left(t_{1,i}\right)}.$$E43

It may be noted that this $1\times 1$ scalar block can be seen as a particular case of the general $\left(2K+1\right)\times \left(2K+1\right)$ block of (42), if $K=0$ is assumed as the maximum harmonic order. Actually, given that ${\mathit{df}}_m\left(\widehat{y}\left(t_{1,i},t_2\right)\right)/d{\widehat{y}}_l\left(t_{1,i},t_2\right)$ is a constant function evidencing no fluctuations in the $t_2$ fast time scale, the last term of (42) will vanish and there will be no more necessity of converting the right-hand side terms of (43) into the frequency domain.

Since only fast-varying signals are converted into the frequency domain, this partitioned technique can be seen as a combination of multivariate envelope transient harmonic balance (by treating the fast-varying signals in a hybrid time-frequency framework) with a pure time marching simulation engine (by treating some of the signals in a pure time domain scheme). Thus, beyond the notorious significant vector and matrix size reductions above mentioned, there is another important advantage brought by this partitioned technique. For example, in complex heterogeneous RF systems strongly nonlinear regimes of operation are in general associated to digital or baseband blocks, whereas moderately nonlinear regimes are typical of RF blocks. With this partitioned technique signals in digital and baseband blocks are appropriately computed in the time domain, while signals in RF blocks are treated in the hybrid time-frequency framework.

5.3. Circuit-block partitioning strategy

Similar to what we have mentioned for the methods discussed in the previous sections, it will be of great utility if the simulator is able to automatically distinguish the fast-varying variables from the slowly varying ones. We now briefly review the approach addressed in [19] for that purpose, which splits the circuit into distinct blocks according to the time rates of change of their voltages and currents.

Let us consider the HB system of (39). As stated earlier, this system has to be solved with the iterative scheme of (41) for each artificial time line $t_{1,i}\times \left[0,T_2\right]$. The automatic partitioning strategy proposed in [19] consists in considering all the circuit’s variables as fast on the first iteration of (41), that is, it consists in initially representing all the unknowns in the frequency domain as a set of $\left(2K+1\right)$ Fourier coefficients. This way, the algorithm starts by computing a single iteration in (41) to evaluate ${\widehat{\mathbf{Y}}}^{\left[1\right]}\left(t_{1,i}\right)$

$${\widehat{\mathbf{Y}}}^{\left[1\right]}\left(t_{1,i}\right)={\left[{{\widehat{\mathbf{Y}}}_1}^{\left[1\right]}{\left(t_{1,i}\right)}^T,{{\widehat{\mathbf{Y}}}_2}^{\left[1\right]}{\left(t_{1,i}\right)}^T,\dots ,{{\widehat{\mathbf{Y}}}_n}^{\left[1\right]}{\left(t_{1,i}\right)}^T\right]}^T,$$E44

where each ${\widehat{\mathbf{Y}}}_v\left(t_{1,i}\right)$, $v=1,\dots ,n$, is a $\left(2K+1\right)\times 1$ vector defined as

$${{\widehat{\mathbf{Y}}}_v}^{\left[1\right]}\left(t_{1,i}\right)={\left[{Y_{v,-K}}^{\left[1\right]}\left(t_{1,i}\right),\dots ,{Y_{v,0}}^{\left[1\right]}\left(t_{1,i}\right),\dots ,{Y_{v,K}}^{\left[1\right]}\left(t_{1,i}\right)\right]}^T.$$E45

Each of the ${\widehat{\mathbf{Y}}}_v^{\left[1\right]}\left(t_{1,i}\right)$ is then inspected. If its Fourier coefficients of order $k\ne 0$ are practically null (their absolute values stay under a very small prescribed tolerance) it will be classified as slow. Otherwise it will be classified as fast. After this classification, which will temporarily split the system into fast and slow subsystems, the simulator considers that signals in slow subsystems can be represented by a single Fourier coefficient for the remaining iterations needed to compute the solution of (41).

Similar to what we have discussed for the method described in Section 4, the robustness of this partitioning strategy may be improved if more than one iteration in (41) is computed before the simulator decides which signals are slow and which signals are fast. The main drawback of such approach is the loss of some efficiency due to the extra computational effort required.

As a final point of this section, we would like to point out that significant gains in computation speed have been reported in the scientific literature (e.g., [18, 19]) for the simulation of several illustrative examples of practical interest using these hybrid time-frequency envelope following partitioned algorithms.

6. Conclusions

In this chapter, we have briefly reviewed some powerful numerical simulation techniques based on partitioned stratagems. Such techniques were especially designed to cope with the simulation challenges brought by emerging electronic technologies to the EDA community, as is the case of complex heterogeneous electronic systems composed of a combination of different kinds of circuit blocks (analog, mixed-signal and digital blocks, or even radio frequency blocks) containing node voltages and brunch currents of very distinct formats and running on widely separated time scales. With these partitioned techniques signals within different blocks of the circuits are computed with distinct algorithms and/or step sizes. Considerable reductions on the computational cost have been registered in several experiments published in the scientific literature (in comparison to previously recognized techniques) without compromising the accuracy of the results.

Acknowledgments

This work is funded by National Funds through FCT—Fundação para a Ciência e Tecnologia, under the project UID/EEA/50008/2013.