528 IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN OF INTEGRATED CIRCUITS AND SYSTEMS, VOL. 29, NO. 4, APRIL 2010 FSSA: Fast Steady-State Algorithm for the Analysis of Mixed Analog/Digital Circuits Angelo Brambilla, Giambattista Gruosso, Member, IEEE, and Giancarlo Storti Gajani Abstract —The shooting method is largely employed to deter- mine the steady-state working condition of both autonomous and nonautonomous circuits. In general, the conventional shooting method employs the Newton algorithm to estimate a better approximation of the steady-state working condition. The Newton algorithm requires the computation of the Jacobian matrix and this seriously limits the use of the conventional shooting method to solve medium/large scale circuits. In this paper, an approach to efficiently determine the shooting matrix is presented. It is shown that the approach is also adequate to deal with mixed analog/digital circuits. Index Terms—Leading dynamics, mixed analog/digital circuits, shooting method, steady-state analysis. I. Introduction I N THIS PAPER, an improved version of the shooting method is presented. The shooting method is a powerful and well-known tool to compute the steady-state behavior of periodic circuits. Its basic implementation for circuit applica- tions has been presented, among the others, in the classical works by Aprille and Trick [1], Kundert et al. [2], Ogrodzki [3], and Vlach and Singhal [4] and it is still a fundamental method for the analysis and design of oscillators, see, e.g., [5], [6] and [7]. Basically, the shooting method solves a boundary value problem where the unknowns are the values of the state variables at the beginning and at the end of one working period. In the case of autonomous circuits, the working period is added among the electrical unknowns. In general, the nonlinear equation implementing the boundary problem is solved with the Newton iterative method. The Newton method requires the computation of the Jacobian matrix and, at each iteration, of the residue vector, i.e., the difference of the state variables at the beginning with respect to the end of the working period. The residue is computed in a straightforward way by performing a time domain analysis along the working period. The main component of the solving matrix in the conventional shooting method (that we will here on call “shooting matrix” J S ) is the sensitivity of the state variables at the end of the working period with respect to variations of the state variables at the beginning of the period. Manuscript received March 29, 2009; revised September 6, 2009 and November 13, 2009. Current version published March 19, 2010. This paper was recommended by Associate Editor, G. Gielen. The authors are with the Dipartimento di Elettronica e Informazione, Politecnico di Milano, Milano I20133, Italy (e-mail: brambill@elet.polimi.it; gruosso@elet.polimi.it; storti@elet.polimi.it). Digital Object Identifier 10.1109/TCAD.2010.2042886 In a formal way, this sensitivity is represented by a square matrix, which, in general, is full. This matrix comes “for free” from the transient analysis. It is formed by a left matrix product of partial matrices, each constituted by the inverse of the Jacobian matrices that are computed by the Newton method within the time domain analysis at each time step. The left product of full matrices is numerically expensive since its cost is proportional to N 3 , where N is the dimension of the square matrix. Furthermore, the matrix product is composed of S terms where S is the number of time samples employed by the time domain analysis. This aspect impacts on the efficiency of the conventional shooting method and prevents its usage to compute the steady-state solution of medium/large circuits. In this context, the time spent to LU-factorize matrix J S is only a small fraction of the total time spent to perform the matrix products needed to build it. This drawback can be mitigated by using the generalized minimal residual (GMRES) method [8]. In this case, each partial matrix forming J S is multiplied by a vector so that the total cost per iteration is reduced to SN 2 . The GMRES method is based on the formation of a Krylov base spanning the leading dynamics subspace and incremental QR factorizations to ensure orthogonality. Another drawback is given by the fact that some of the Jacobians composing the matrix product exist if and only if the solution satisfies the Lipschitz condition. In the standard IEEE floating point arithmetic, the finite relative precision is about  = 10 −16 , therefore, there is no way to accurately “follow” relative variations greater than 1/. This means for example that “ideal” behavioral models of digital ports can lead to accuracy problems or converge failures of the conventional shooting method. In general, as shown in this paper, circuits with switching elements, e.g., when some electrical variables go above or below some given threshold, yield to cusps in the state variable waveforms, that is to discontinuities of the first kind in the derivatives of state variables [9]. If we think of mixed mode simulation of analog/digital circuits it is immediate to see that the digital portion of the circuit leads to discontinuities in some electrical variables or to cusps in the state variables of the analog portion of the circuit. This means that the conventional shooting method, based on the Newton method (and possibly adopting the GMRES iterative method), cannot be employed to compute the steady-state behavior for the class of circuits mentioned above. The problem of steady-state analysis of circuits with switching elements is well known, in particular in the Power 0278-0070/$26.00 c  2010 IEEE