Reduced order multirate schemes in industrial circuit simulation

In this paper the industrial application of Reduced Order Multirate (ROMR) schemes is presented. This paper contains the mathematical foundations of the ROMR schemes and elaborates on the construction of these schemes using specific Model Order Reduction (MOR) techniques. Especially the Maximum Entropy Snapshot Sampling method for generating a reduced basis and reduction by Gauß–Newton with Approximated Tensors (GNAT). This basis generation method is also used for generating the basis for the gappy hyper-reduction method used for nonlinear function evaluation. For the multirate integration part, a Backward Differentation Formula approach to integration is used in conjunction with a coupled-slowest-first multirate approach. After introducing the numerical approach to industrial circuit simulation validation experiments are performed. First a simple academic model is used, and then an industrial test case is simulated as presented by STMicroelectronics. A significant speedup in simulation time is achieved whilst accuracy and convergence is kept.


Introduction
In integrated circuit design, there are a significant number of design possibilities under which the internal components need to be guaranteed to work. This leads to a whole range of explorations to ensure sound functionality of the design. These explorations are performed by numerical simulations of the circuits mathematical model. Due to the ever increasing number of components, and thus the degrees of freedom in the model, the required simulation times may become prohibitively large.
Besides the sheer number of components inside the integrated circuit, a large contribution to the complexity of the mathematical model originates form the method of deriving these equations. As these models grow, generating a state-space model with a minimal set of unknowns cannot be generated in an automatic way. Therefore, the mathematical models have to be derived through use of algorithmic analysis. This automation comes at a cost. The resulting system of differential-algebraic equations (DAE) is numerically harder to solve, and may contain redundant network variables.
To decrease these ever increasing simulation costs a multitude of different approaches have been proposed in the past decades. For instance, the large original system can be partitioned into subsystems which each have their own characteristic rate of evolution through time. This property is capitalised upon by using multirate (MR) time integration. Another approach is to exploit redundancy in the mathematical model originating from the network analysis by using model order reduction. This technique aims to solve a model of reduced size that still approximates the solution of the original model. With the main driver being that such large coupled systems eligible for both model order reduction (MOR) and MR techniques are encountered more often in industry, for instance by STMicroelectronics. Thus is this paper specifically aimed at the combination of the two previously mentioned techniques, and provides a review of the definition, with new integration approaches.
In Sect. 2, the continuous DAE problem is introduced and is discretized to be solved numerically. The backward differentiation formula as numerical method is covered and the notion of multirate time-integration, [9], is introduced.
The model order reduction techniques considered in this paper are outlined in Sect. 3. The specific type of model order reduction we discuss is nonlinear model order reduction by reduced basis approaches. Specifically the Gauß-Newton with approximated tensors method [6]. This robust model order reduction method is a nonlinear Petrov-Galerkin projection method equipped with a function-sampling hyper-reduction scheme. It operates at the level of the nonlinear arising at each time step, after discretization. In the last part of this section the combination of multirate and model order reduction techniques is described.
In the penultimate section numerical experiments are performed on both academic and industrial test-cases. The academic test case is used to verify the implementation of the circuit parser, and check the validity of the numerical scheme. Then, the implemented circuit parser and simulator is applied to an industrial test case as provided by STMicroelectronics. Finally, the last section offers conclusions and provides an outlook.

Problem formulation
Since the set of equations used to describe the electrical circuits is constructed according to the topological structure of the network. This often results in a coupled system of implicit differential and nonlinear equations, or more general a system of differentialalgebraic equations (DAEs) This system may represent ill-posed problems and is in general more difficult to solve numerically than the more standard systems of ordinary differential equations (ODEs).
To solve these systems numerically we start from the consistent initial values. Then the time domain is discretised into time points t 0 , t 1 , . . . , t N , and the solution for each of these time points is approximated by an implicit linear numerical integration formula. A direct approach, as proposed in [13], is by applying backward differentiation formula (BDF) method. In this paper we restrict ourselves to using a first order BDF, but the idea can be generalised for higher order methods taking some considerations into account about the smoothness. This multistep method is applied to a DAE system by using the -embedding method.
However, we want to solve system (2.1), which can be an implicit differential algebraic system. Therefore, the multistep system for an implicit DAE system, Mẋ = f (x), is given by In general form, applying Equation (2.2) to an implicit nonlinear system of DAEs at time step t n yields This gives that the numerical solution of the system is thus reduced to the solution of the system of nonlinear Equations (2.3). This system is solved iteratively for x n by applying Newton's method.
To set up the system of DAEs defined by Equation (2.2), Modified Nodal Analysis (MNA), [8], is used to obtain network equations. The Kirchoff Current Law is applied to each node, except the node that is considered the ground node. The incidence matrix is then defined as a collection of incidence matrices related to each different type of element, with A ∈ {0, +1, -1} n u ×n , where n is the cardinality the set of each type of network element. Using these incidence matrices, we can relate the branch voltages in a loop and the currents accumulating in a node by applying both the Kirchoff Current Law and the Kirchoff Voltage Law to each node, resulting in The unknowns q, φ, u, j L and j V are the charges, fluxes, node voltages, inductor currents and voltage source currents, respectively. All these quantities are time dependent, and are combined into one state vector x(t) ∈ R m of unknowns. The dimension of which is given by the cumulative dimensions of the quantities. The network equations can now be stated in compact form where q and j are mappings from R m to R m related to the network elements.
The previously seen numerical integration method is considered a singlerate timeintegration method, as it integrates each part of the equation with the same step-size. Opposed to this classical approach there are the multirate time-integration methods, which use a different step-size, or even integration method, for parts of the equations with different dynamical behaviour.
First, it is necessary for a multirate method to partition the variables of the dynamical system into sets with different temporal characteristic. This partition can be made either by manual selection, or automatically. For now, we partition the full system into a fast and a slow partition. This can be extended to a partitioning into k subsystems, but for convenience, a two subsystem partitioning, i.e. a fast and slow subsystem, is used as example here.
Using this partitioning on the compact form, Equation Applying this partition to the network equations (2.5) results in the following systems where we have the following definitions To apply multirate time-integration methods to DAEs it is assumed that all the subsystems in itself are well posed, which means that the DAE-index should be less than or equal to that of the full DAE. As we are partitioning network equations of circuit models, it is known that they are composed of sub-circuits or natural phenomena in a hierarchical way. If we partition based on this hierarchy, we can check easily if the subsystems fulfil the index requirements and obtain viable partitions. As there are several different approaches available to implement multirate timeintegration methods, we specify in this subsection which approach is used. As has been shown in [1] and [2] a feasible method is given by the coupled-slowest-first integration approach coupled with the implicit Euler method specified by a first order BDF method. First the whole system is solved for the macro-step, t n → t n+1 = t n + H, However, the fast components, x F from the full state solution x are inaccurate, thus they are discarded. Then by using interpolated values for the slow system between t n and t n+1 , the following system is solved for the fast solutions. Here we have that H = h · f mr where f mr is known as the multirate factor.
For stability reasons, the interpolated valuesx S,n+(l+i)/m are obtained by constant interpolation based on x S,n+k , then the coupled-slowest-first Euler approach is unconditionally A-stable.

Nonlinear model order reduction
Due to the differential-algebraic structure of the generated network equations, standard techniques used for ODE model order reduction, such as a direct reduction through a Galerkin projection scheme, may yield unsolvable reduced order models.
To circumvent this and preserve the solvability of the numerical model, the Gauß-Newton with Approximated Tensors (GNAT) model order reduction approach is utilized, [6]. This is combined with a gappy data reconstruction method, [14], for hyper-reduction. These reduced basis model order reduction approaches utilise a basis obtained by the Maximum Entropy Snapshot Sampling (MESS) method, [10], as proposed in [2]. To illustrate this approach, this section contains an overview of the reduction, hyper-reduction and basis-construction methods.

Gauß-Newton with approximated tensors
As GNAT operates on a discrete level, we assume that Equation (2.1) is solved by a linear implicit time integrator, see the next section for a more detailed description. Then for n t time steps, a sequence of n t systems of nonlinear equations are to be solved, with each system defined at time step where x ∈ R M and the residual mapping R : R M → R M . For ease of notation, the time dependencies have been omitted as we consider only one time instance and one input vector. The GNAT approach for reduction is to use a projection to search the approximated solution in the incremental affine trial subspace x 0 + V ⊂ R M , with the initial condition x 0 ∈ R M . The incremental subspace is used for consistent Petrov-Galerkin projections, Appendix A [6]. This reduced state vector is then given bỹ where V x ∈ R M×r is the r-dimensional projection basis for V, which is not yet defined, and x r denotes the reduced incremental vector of the state vector. Substituting Equation This approach by residual minimisation is equivalent to performing a Petrov-Galerkin projection where te test basis is given by ∂R ∂x (V x ), [5]. This nonlinear least-squares problem is solved by the Gauß-Newton method, leading to then iterative process for k = 1, . . . , K , solving and updating the search value w k r with where K is defined through a convergence criterion, initial guess w 0 Here J k is the full order Jacobian of the residual at each iteration step k. Since the computation of this Jacobian, which is used in circuit simulation to solve for the next step, scales with the original full dimension of Equation (3.1) this is a computational bottleneck. This bottleneck can be circumvented by the application of hyper reduction methods, for which this paper utilises a gappy data reconstruction method.

Hyper-reduction by gappy data reconstruction
The evaluation of the nonlinear function R(x 0 +V x w k r ) has a computational complexity that is still dependent on the size of the full system. To reduce the complexity of this evaluation a gappy data reconstruction, based on [7], is applied. Like the GNAT approach gappy data reconstruction uses a reduced basis to reconstruct the data. Gappy data reconstruction starts by defining a mask vector n for a solution state x as n j = 0 if u j is missing, where j denotes the j-th element of each vector. The mask vector n is applied point-wise to a vector by (n, x) j = n j x j . This sets all the unobserved values to 0. Then, the gappy inner product can be defined as (a, b) n = ((n, a), (n, b)), which is the inner product of the each vector masked respectively. The induced norm is then ( x n ) 2 = (x, x) n . Considering the reduction base obtained by MESS V gap = {v i } r i=1 , now we can construct an intermediate "repaired" full size vectorũ from a reduced vector u with only g elements bỹ where the coefficients b i need to minimise an error E between the original and repaired vector, which is defined as This minimisation is done by solving the linear system From this solutionũ is constructed. Then the complete vector is reconstructed by mapping the reduced vectors elements to their original indices and filling the rest with the reconstructed values.

Reduced basis generation
As stated previously, both the nonlinear model order reduction method and the hyperreduction method use a reduced basis. To generate such a basis, generally snapshot backended reduced basis methods are used. A prominent method for this in literature is the proper orthogonal decomposition method. However, as shown in [3,10], the way the proper orthogonal decomposition framework extracts information from the high-fidelity snapshot matrix is inherently linear. When used for nonlinear problems, it removes high-frequency components that are present and relevant to the dynamical evolution of these systems. Therefore, the MESS method for reduced basis generation is used.

Maximum entropy snapshot sampling
Let m and n be positive integers and m n > 1. Define a finite sequence X = (x 1 , x 2 , . . . , x n ) of numerically obtained states x j ∈ R m at time instances t j ∈ R, with j ∈ {1, 2, . . . , n}, of a dynamical system governed by either ODEs or DAEs. Provided probability distribution p of the states of the system, the second-order Rényi entropy of the sample X is H (2) p (X) = -log n j=1 p(x j ) 2 = -log E p(x j ) , (3.10) with p j ≡ p(x j ) and where E(p(X)) is the expected value of the probability distribution p with respect to p itself. When n is large enough, according to the law of large numbers, the average of p 1 , p 2 , . . . , p n almost surely converges to their expected value, 1 n n j=1 p(x j ) → E p(X) as n → ∞, (3.11) thus each p(x j ) can be approximated by the sample's average sojourn time or relative frequency of occurrence. To obtain this frequency of occurrence, considering a norm · on R m . Then the notion of occurrence can be translated into a proximity condition. In particular, for each x j ∈ R m define the open ball that is centred at x j and whose radius is > 0, (3.12) and introduce the characteristic function with values Under the aforementioned considerations, the entropy of X can be estimated bŷ (3.14) Provided that the limit of the evolution ofĤ (2) p exists and measures the sensitivity of the evolution of the system itself [4, §6.6]. Then, for n large enough, a reduced sequence X r = (x j 1 ,x j 2 , . . . ,x j r ), with r ≤ n, is sampled from X. This is done by requiring that the entropy of X r is a strictly increasing function of the index k ∈ {1, 2, . . . , r} [11]. The state vectorx j k added to sampled snapshot space is the average value of all states in the selected -ball. A reduced basis V x is then generated from X r with any orthonormalization process. It has been shown [10] that any such basis guarantees that the Euclidean reconstruction error of each snapshot is bounded from above by , while a similar bound holds true for future snapshots, up to a specific time-horizon. See Theorems 3.3 and 3.4, [10].
To estimate the parameter , which determines the degree of reduction within the MESS framework, the following optimisation approach is employed [3]. As no user input is now required this optimality requirement makes the MESS parameter free reduction method.

Reduced order multirate
Combining the previously seen techniques, we apply the model order reduction techniques only in the macro-step of the multirate integration scheme. However, due to the partitioned nature, this full system needs a special type of reduced basis matrix as it needs to reduce only the slow part. To illustrate an implementation of a reduced order multirate scheme, a first order implicit Euler scheme is described. The new reduction matrix is now given by ∈ R r+n F ×M . Considering the split system defined by Equations (2.7a)-(2.7c) the reduction matrix is defined by Then in each macro time-step we only need to solve the following reduced system. This nonlinear least-squares problem is solved by the Gauß-Newton method, using , Where R is defined for the whole macro step system as Note however that now we have J k r as the full order Jacobian approximated by using the gappy-MESS function evaluations. Then once the solution for the macro time-step is ob-tained, the slow partition of this solution x S are used for the interpolated values needed in the intermediate fast steps.

Experiments
The verification of the reduced order multirate method, utilising the discussed methods, and the netlist parser to create the network equations is done through numerical experiments. For the first experiment a simple academic circuit consisting of resistors, capacitors and diodes is considered. Then a real world industrial circuit is considered as provided by STMicroelectronics. Although some simplification assumptions have been made about the parameters of the underlying network components.

Academic experiment
The academic circuit shown in Fig. 1 is a combination of a short diode chain and then a long ladder of diodes and resistors. Similar to a standard diode chain model [12], this model contains sufficient redundancy to make it eligible for model order reduction. Furthermore, increasing the resistance for each R i with R < R i < R i+1 makes that the ladder part of the circuits behaves on a slower timescale. This makes the circuit excellent for a time integration with the reduced order multirate approach. The simulation parameters are given in Table 1.
With the multirate factor f mr = 20, it is evident that the multirate integration approach is more accurate than the single rate method for the same macro step sizes, see Fig. 2. The order however is equal due to the pollution introduced by the macro step. This increase of accuracy comes at a cost. The right figure shows that the multirate technique carries a

Solar panel simulation
As final industrial application STMicroelectronics has provided a test case from their development team. The test case concerns a transient analysis of a solar panel circuit. In this case we consider a silicon photovoltaic cell that is composed of two layers of semiconducting material, e.g. silicon or gallium arsenide, with different doping. The process of converting solar radiation into electricity is based on the photovoltaic effect. To model the photovoltaic effects taking place in the solar panel, a mathematical model is used based on the construction of an equivalent circuit. When the cell is lit, the generated power can be modelled as a current source. Therefor, an equivalent circuit of the photovoltaic cell is shown in Fig. 3. The full solar panel circuit is a conglomeration of photovoltaic cells, Fig. 3, that are linked in series and parallel. This grid of photovoltaic cells is then connected to a DC-DC buck converter to stabilise the output. The full schematic of the solar panel circuit is shown in Fig. 4. As the voltages and currents generated by the solar panel vary significantly slower than the operating voltages in the buck converter there is an inherent multirate advantage. Due to the replication in the structure of the solar panel, model order reduction can be used to significantly reduce the complexity of the panel simulation.
The full test case for the silicon photovoltaic solar panel is run for a series and parallel grid that is 35 by 35, which makes the full dimensional system consist of 2422 equations. The reduction parameter selection procedure provides us with = 0.43 which leads to   Figure 5 Computational effort of the numerical schemes applied to the solar panel quite a significant reduction. As shown in Table 2 the dimension of the slow subsystem of photovoltaic cells is decreased from 2416 to only 6 equations, and the hyper-reduced dimension of the nonlinear function is reduced to only 8 equations. From the figures in Fig. 5, we see that the academically achieved results are replicated for an industrial test case, Fig. 2. Due to the model order reduction there is a substantial decrease in computation time, with the reduced order multirate scheme being roughly 6.35 times faster than the multirate integration scheme, without losing accuracy.

Conclusion and outlook
In this paper the numerical approach for implementing the reduced order multirate timeintegration method has been outlined in the context of circuit simulation with the novel BDF approach. Using Modified Nodal Analysis for the generation of the network equations has been shown to be an effective and efficient approach. The resulting equations are suitable for nonlinear model order reduction, as has been shown by the application of the Gauß-Newton with Approximated Tensors method extended with a gappy data reconstruction approach. The reduced bases for both model order reduction techniques are obtained by the Maximum Entropy Snapshot Sampling method. The reduction parameter is generated using a parameter free optimisation method. The whole mathematical framework has been applied to both an academic and an industrial test case, supplied by STMicroelectronics. A significant reduction in computation time can be seen whilst the accuracy of the multirate approach closely approximated. For actual use in the optimisation flow of STMicroelectronics a more robust framework would need to be created but the numerical results show a significant potential.