## I Introduction

The ongoing miniaturization of modern IC devices has led to extremely complex circuits. This results in the increase of the problems associated with the analysis and simulation of their physical models. In particular, the performance and reliable operation of ICs are largely determined by several critical subsystems such as the power distribution network, multi-conductor interconnections, and the semiconductor substrate. The electrical models of the above subsystems are very large, consisting of hundreds of millions or billions of electrical elements (mostly resistors R, capacitors C, and inductors L

), and their simulation is becoming a challenging numerical problem. Although their individual simulation is feasible, it is completely impossible to combine them and simulate the entire IC in many time-steps or frequencies. However, for the above subsystems it is often not necessary to fully simulate all internal state variables (node voltages and branch currents), as we only need to calculate the responses in the time or frequency domain for a small subset of output terminals (ports) and given excitations at some input ports. In these cases, the very large electrical model can be replaced by a much smaller model whose behavior at the input/output ports is similar to the behavior of the original model. This process is called Model Order Reduction (MOR).

MOR methods are divided into two main categories. System theoretic techniques, such as Balanced Truncation (BT) [1], provide very satisfactory and reliable bounds for the approximation error. However, BT techniques require the solution of Lyapunov matrix equations which are very computationally expensive, and also involve storage of dense matrices, even if the system matrices are sparse. On the other hand, moment-matching (MM) techniques [2] are well established due to their computational efficiency in producing reduced-order models. Their drawback is that the reduced-order model depends only on the quality of the Krylov subspace.

The majority of MM methods exploit the standard or the rational Krylov subspace in order to approximate the original model. Authors in [3, 4]

employ rational Krylov MM methods to reduce power delivery networks. Using this projection subspace requires a heuristic and expensive parameter selection procedure, while the approximation quality is usually very sensitive to an inaccurate selection of these parameters. Moreover, in

[2, 5] a standard Krylov subspace is employed for the reduction of regular and singular systems, respectively. Generally, established MM methods construct the subspace only for positive directions, usually leading to a large approximated subspace to obtain a satisfactory error. Recent developments in a wide range of applications have shown that the approximation quality of the Extended Krylov Subspace (EKS) outperforms the one of the standard Krylov subspace [6]. However, the application of EKS in the context of circuit simulation is not trivial. In several problems, EKS computation involves singular circuit models and dense matrix manipulations, which can hinder the applicability of this subspace.In this paper, we introduce an EKS Moment-Matching (EKS-MM) method that greatly decreases the error induced by MM methods by approximating both ends of the spectrum. To enable the simulation of many-port models, the proposed method exploits the superposition property. More specifically, we develop a procedure for applying EKS-MM to large-scale regular and singular models, by implementing computationally efficient transformations in order to preserve the original form of the sparse input matrices. Finally, we evaluate our methodology on industrial IBM power grids. The rest of the paper is organized as follows. Section II presents the theoretical background of MM methods for the reduction of regular and singular circuit models. Section III presents our main contributions on the application of EKS to MM methods, as well as its efficient implementation by sparse matrix manipulations for both regular and singular circuit models. Section IV presents the experimental results, while conclusions are drawn in Section V.

## Ii Background

### Ii-a MOR by Moment-Matching

Consider the Modified Nodal Analysis (MNA) description of an n-node, m-branch (inductive), p-input, and q-output RLC circuit in the time domain:

(1) | |||

where (node conductance matrix), (node capacitance matrix), (branch inductance matrix), (node-to-branch incidence matrix),

(vector of node voltages),

(vector of inductive branch currents), (vector of input excitations from current sources), (input-to-node connectivity matrix), (vector of output measurements), (node-to-output connectivity matrix), (input-to-output connectivity matrix). Without loss of generality, in the above we assume that any voltage sources have been transformed to Norton-equivalent current sources, and that all outputs are obtained at the nodes as node voltages. Furthermore, and .If we now denote the model order as , the state vector as , and also:

then expression (1) can be written in the following generalized state-space form, or so-called descriptor form:

(2) |

The objective of MOR is to produce a reduced-order model:

(3) |

where , , . The reduced model has order , and the output error is bounded as for given input and given small . The bound in the output error can be equivalently written in the frequency domain as via the Plancherel’s theorem [7]. If

are the transfer functions of the original and the reduced-order model, respectively, then the output error in frequency domain is:

(4) | |||

where is the induced matrix norm, or norm, of a rational transfer function. Therefore, in order to bound the output error, we need to bound the distance between the transfer functions .

The most important and successful MOR methods for linear systems are based on MM. They are very efficient in circuit simulation problems and are formulated in a way that has a direct application to the linear model of (2).

By applying the Laplace transform to (2), we obtain the domain equations as:

(5) | |||

Assuming that and that an impulse response is applied to (i.e. ), then the above system of equations can be written as follows:

(6) |

and by expanding the Taylor series of around zero, we derive the below equation:

(7) |

The transfer function of (2) is a function of , and can be expanded into a moment expansion around as follows:

(8) |

where , , , , are the moments of the transfer function. Specifically, in circuit simulation problems, is the DC solution of the linear system. This means that the inductors of the circuit are considered as short circuits, and the capacitors as open circuits. Moreover, is the Elmore delay of the linear model, which is defined as the time required for a signal at the input port to reach the output port. Finally, is related to the system matrices as:

(9) |

The goal of MM reduction techniques is the derivation of a reduced-order model where some moments of the reduced-order transfer function match some moments of the original transfer function .

Let us now denote the two projection matrices onto a lower dimensional subspace as and , respectively. These matrices can be derived from the associated moment vectors using one or more expansion points. As a result, if we assume that , then the matrices and are defined as follows:

(10) | |||

The computed reduced-order model matches the first moments and is obtained by the following matrices:

(11) |

This reduced model provides a good approximation around the DC point. Finally, in case we employ an one-sided Krylov method, which is usually the case, the matrix can be set equal to , an equality that also holds for symmetric systems.

### Ii-B Handling of Singular Descriptor Models

In certain circuit simulation problems, the matrix

might be singular. A method for dealing with such models is to compute spectral projections onto the left and right deflating subspaces corresponding to the finite eigenvalues of the model, which is computationally prohibitive for large-scale systems. However, singular descriptor models typically result when there are some nodes, say

, where no capacitance is connected, leading to corresponding all-zero rows and columns in the submatrix . Note that in case the circuit contains no voltage sources, the submatrix of inductive branches is always nonsingular. If the nodes with no capacitance connection are enumerated last, and the remaining nodes first, then (1) can be partitioned as follows:(12) | |||

where , , , , , , , , , , , and .

Assuming now that the submatrix is nonsingular
(a sufficient condition for this is at least one resistive connection from any of the non-capacitive nodes to ground), the second row of (12) can be solved for as follows:

(13) |

The above can be substituted to the first and third row of (12), as well as the output part of (12), to give:

This can be put together in the following descriptor form:

(14) | |||

The above is a nonsingular (i.e. regular) state-space model which can be reduced normally.

## Iii Extended Krylov Subspace for MOR

### Iii-a EKS Moment-Matching (EKS-MM)

The essence of MM methods is to iteratively compute a projection subspace, and then project the original system into this subspace in order to obtain the reduced-order model of (3). The dimension of the projection subspace is increased in every iteration, until an a-priory selection of the moments is matched. More specifically, if is the desired order for the reduced system and is the number of moments, then () is a projection matrix whose columns span the -dimensional Krylov subspace:

where

Then, the reduced-order model is obtained through the following matrix transformations:

(15) |

with , , .

The projection process is independent of the subspace selection, but its effectiveness is critically dependent on the chosen subspace. As a result, one choice is to consider the rational Krylov subspace [3, 4]. However, this projection subspace requires the input of a number of shift parameters, whose choice greatly affects the produced reduced-order model. The reason for this is that it relies on unclear heuristics and is highly problem-dependent. In order to address this issue, the standard Krylov subspace [2, 5] must be enriched with information from the subspace , which corresponds to the inverse matrix , leading to EKS:

(16) |

The Arnoldi procedure [10] that computes EKS begins with the pair , and then generates a sequence of extended subspaces in order to compute the matrix and produce the reduced-order model as described in (15). EKS can be considered a special case of the rational Krylov subspace with two expansion points, one expansion point at zero and one at infinity. The complete procedure is given in Algorithm 1.

At this point, we can elaborate on some aspects regarding the efficient implementation of the proposed EKS procedure:

#### Iii-A1 Sparse matrix inputs

It is worth mentioning that Algorithm 1 does not require matrices , as inputs, but only the sparse system matrices , are necessary. This is due to the fact that the generally dense inverse matrices are only needed in products with vectors (initially in step 3) and vectors (in step 7 at every iteration, where the iteration count is typically very small and thus ). These products can be implemented as sparse linear solves ( and ) by employing any sparse direct [8] or iterative [9] algorithm.

#### Iii-A2 Orthogonalization in steps 3 and 9

A modified Gram-Schmidt procedure [10] is employed to implement the corresponding qr() procedures.

#### Iii-A3 Orthogonalization in step 8

### Iii-B Sparse Implementation for Singular Descriptor Models

Algorithm 1 is computationally inefficient for the reduction of the model given in (14), which results from the regularization of a singular descriptor model, since the inversion of renders the matrices dense and hinders the solution procedure. In this subsection, we present efficient ways to implement the EKS algorithm by preserving the original sparse form of the system matrices.

#### Iii-B1 Construction of RHS

#### Iii-B2 Sparse linear system solutions

The system matrix

(18) |

of the model given in (14) is rendered dense due to the inversion of . The linear system solutions with in steps 3, 7 of Algorithm 1 can be handled by partitioning the RHS of these systems conformally to , i.e. with , , and implementing their solution efficiently by keeping all the sub-blocks in their original sparse form as follows:

(19) |

where is a temporary sub-matrix.

#### Iii-B3 Sparse matrix-vector products

The matrix-vector products with in step 7 of Algorithm 1 can be implemented efficiently by observing that:

(20) | |||

Therefore, the product with vectors can be carried out by a sparse solve , followed by a sum of products .

#### Iii-B4 Construction of system matrix

In order to construct and then reduce the dense system matrix of (18), we need to employ sparse solves with the submatrix . Since usually , it is better to first compute the left-solves and , followed by products with and . The left-solves can be performed as and , where contains the rows of each left-solve.

### Iii-C Superposition Property of LTI Models

While in the previous subsections we emphasized on the efficient execution of the proposed methodology, it still can not handle many-terminal models. To this end, we consider the superposition principal of LTI models. Using the superposition property, the output response of the initial multi-input multi-output (MIMO) descriptor model of (2) can be computed as the sum of the output responses of the following single-input multi-output (SIMO) subsystems as:

(21) |

where is a matrix with only one nonzero column of the input-to-node-connectivity matrix , and . From these relations, it can be derived that and .

This property can be employed for the parallel computation of the reduced-order model. Each partitioned model of (21) can be reduced by a projection matrix whose columns span the k-dimensional EKS:

(22) |

with , and similarly the reduced-order model is obtained by:

(23) |

Moreover, each reduced-order transfer function is computed as:

(24) |

Finally, the approximate transfer function of the reduced-order model is computed as:

(25) |

It must be noted that there is no guarantee that the passivity of the reduced-order models obtained using the superposition property is preserved. In recent years, however, the focus of MOR has been shifted from provably passive models to passivity enforcement after efficient reduction. A wealth of passivity enforcement techniques, such as [11], have been developed to assure passivity of the reduced-order models obtained using the superposition property.

Ckt | Dimension | #ports | ROM Order | Moment-Matching (MM) | EKS Moment-Matching (EKS-MM) | |||||

#moments | Max Error | Runtime(s) | #moments | Max Error | Error Reduction | Runtime(s) | ||||

percentage | ||||||||||

ibmpg1 | 44946 | 600 | 1200 | 2 | 0.037 | 0.146 | 1 | 0.014 | 62.16% | 0.146 |

ibmpg2 | 127568 | 500 | 2000 | 4 | 0.233 | 1.206 | 2 | 0.131 | 43.78% | 1.277 |

ibmpg3 | 852539 | 800 | 1600 | 2 | 0.253 | 11.029 | 1 | 0.146 | 42.29% | 11.060 |

ibmpg4 | 954545 | 600 | 2400 | 4 | 0.233 | 16.642 | 2 | 0.038 | 83.69% | 17.981 |

ibmpg5 | 1618397 | 600 | 1200 | 2 | 0.242 | 10.228 | 1 | 0.063 | 73.97% | 10.998 |

ibmpg6 | 2506733 | 1000 | 6000 | 6 | 0.161 | 19.155 | 3 | 0.130 | 19.25% | 21.780 |

ibmpg1t | 54265 | 400 | 800 | 2 | 4.767 | 0.259 | 1 | 1.814 | 61.95% | 0.273 |

ibmpg2t | 164897 | 800 | 3600 | 4 | 0.785 | 0.250 | 2 | 0.411 | 47.64% | 0.268 |

## Iv Experimental Results

For the experimental evaluation of the proposed methodology, we used the available IBM power grid benchmarks [12]. Their characteristics are shown in the first three columns of Table I. Note that for the transient analysis benchmarks, ibmpg1t and ibmpg2t, a matrix of energy storage elements (capacitances and inductances) is provided. However, in order to perform transient analysis for the DC analysis bechmarks, ibmpg1 to ibmpg6, we had to add a (typical for power grids) diagonal capacitance matrix with random values on the order of picofarad. In order to evaluate our methodology on singular benchmarks, we enforced the capacitance matrix of ibmpg2 and ibmpg4 to have at least one node that was missing a capacitance connection. These benchmarks along with ibmpg1t and ibmpg2t were represented as singular descriptor models of (12), thus we applied the techniques described in Section III-B for their efficient sparse handling.

EKS-MM was implemented with the procedures described in Section III, and was compared with a standard MM method also implemented with the superposition property. The reduced-order models (ROMs) were evaluated in the frequency range with respect to their accuracy for given ROM order. For our experiments, an appropriate number of matching moments was selected such that the ROM order for both EKS-MM and MM is the same. All experiments were executed on a Linux workstation with a 3.6GHz Intel Core i7 CPU and 32GB memory using MATLAB R2015a.

The results are reported in the remaining columns of Table I, where #moments refers to the number of moments that matched in order to produce the ROMs, Max Error refers to the error between the infinity norms of the transfer functions, i.e. , Runtime refers to the computational time (in seconds) needed to generate each submatrix of (25), while Error Reduction percentage refers to the error reduction percentage achieved by EKS-MM over MM. It can be clearly verified that, compared to MM for similar ROM order, EKS-MM produces ROMs with significantly smaller error. As depicted in Table I, the Error Reduction percentage ranges from 19.25% to 83.69%. The execution time of EKS-MM is negligibly larger than standard MM for each moment computation, due to the expansion in two points, however the efficient implementation can effectively mask this overhead to a substantial extent and make the procedure applicable to very large circuit models.

To demonstrate the accuracy of our method, we compare the transfer functions of the original model and the ROMs generated by EKS-MM and MM. The corresponding transfer functions for one regular (ibmpg1) and one singular (ibmpg2t) benchmark, in the band , are shown in Fig. 1. Fig. 2 presents the transfer functions of ROMs produced by EKS-MM and MM along with the absolute errors induced over the original model for a selected benchmark in the same band. As can be seen, the response of EKS-MM ROM is performing very close to the original model, while the response of MM ROM exhibits a clear deviation. In particular, responses of ROMs produced by MM do not capture effectively the dips and overshoots that arise in some frequencies.

## V Conclusions

In this paper, we proposed the use of EKS to enhance the accuracy of MM methods for descriptor circuit models. Our method provides clear improvements in reduced-order model accuracy compared to a standard Krylov subspace MM technique. For the implementation, we made efficient computational choices, as well as adaptations and modifications for large-scale singular models. As a result, the proposed method still remains computationally efficient, introducing only a small overhead in the reduction process.

## References

- [1] J. Phillips, L. Daniel and L. Miguel Silveira, “Guaranteed passive balancing transformations for model order reduction,” in Design Automation Conference, pp 52–57, 2002.
- [2] A. Odabasioglu, M. Celik and L. T. Pileggi, “PRIMA: passive reduced-order interconnect macromodeling algorithm,” in IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 17, no. 8, pp. 645–654, 1998.
- [3] S. Mei and Y. I. Ismail, “Stable Parallelizable Model Order Reduction for Circuits With Frequency-Dependent Elements,” in IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 56, no. 6, pp. 1214–1220, 2009.
- [4] W. Zhao et al., “Automatic adaptive multi-point moment matching for descriptor system model order reduction,” in International Symposium on VLSI Design, Automation, and Test, pp. 1-4, 2013.
- [5] N. Banagaaya et al., “Implicit index-aware model order reduction for RLC/RC networks,” in Design, Automation & Test in Europe Conference & Exhibition, pp. 1–6, 2014.
- [6] L. Knizhnerman and V. Simoncini, “Convergence analysis of the extended Krylov subspace method for the Lyapunov equation, ” in Numerische Mathematik, vol. 118, no. 3, pp. 567–586, 2011.
- [7] K. Gröchenig, “Foundations of time-frequency analysis,” in Applied and Numerical Harmonic Analysis, 2001.
- [8] T. A. Davis and E. P. Natarajan, “Algorithm 907: KLU, A Direct Sparse Solver for Circuit Simulation Problems,” in ACM Trans. Math. Softw., vol. 37, no. 3, 2010.
- [9] D. Garyfallou et al., “A Combinatorial Multigrid Preconditioned Iterative Method for Large Scale Circuit Simulation on GPUs,” in International Conference on Synthesis, Modeling, Analysis and Simulation Methods and Applications to Circuit Design, pp. 209–212, 2018.
- [10] G. Golub and C. F. Van Loan, Matrix Computations, Johns Hopkins University Press, 1996.
- [11] S. G. Talocia and A. Ubolli, “A comparative study of passivity enforcement schemes for linear lumped macromodels,” in IEEE Transactions on Advanced Packaging, vol. 31, no. 4, pp. 673–683, 2008.
- [12] S. R. Nassif, “Power grid analysis benchmarks,” in Asia and South Pacific Design Automation Conference, pp. 376–381, 2008.

Comments

There are no comments yet.