Fang Liu

Jacob J Flomenberg

Devaka V Yasaratne

Sule Ozev

Hierarchical Variance Analysis for Analog Circuits Based on Graph Modelling and Correlation Loop Tracing

Process variations play an increasingly important role on the success of analog circuits. State-of-the-art analog circuits are based on complex architectures and contain many hierarchical layers and parameters. Knowledge of the parameter variances and their contribution patterns is crucial for a successful design process. This information is valuable to find solutions for many problems in design, design automation, testing, and fault tolerance. In this paper, we present a hierarchical variance analysis methodology for analog circuits. In the proposed method, we make use of previously computed values whenever possible so as to reduce computational time. Experimental results indicate that the proposed method provides both accuracy and computational efficiency when compared with prior approaches.

Introduction

As semiconductor manufacturers push the boundaries of the process technology, a steady trend has been shrinking the feature sizes. One direct result from the reduction in device sizes is the increasing importance of the role of process variations. As devices are scaled down to their physical limits, precise control over process parameters becomes nearly impossible, resulting in larger percentage variations in the parameters at each level of the design hierarchy.

Process variations affect the performance of analog and digital circuits alike, as performance parameters, such as gain, frequency of operation, or path delay, are directly influenced by process parameters. Thus, increased variations in process parameters will also increase variability in performance parameters at the circuit and the system level. An analysis of such parametric variations is needed for many problems in VLSI design, design automation, testing, and fault tolerance. There is an increasing need to incorporate process variations into the design flow from early design stages and understand which parts of the design are the most important contributors.

Since statistical tolerance analysis has been widely used for many applications in industry for several decades, a plethora of approaches have been developed for various application domains over the years. Evans [START_REF] Evans | Statistical tolerancing: The state of the art, part i: Background[END_REF] provides a detailed application background on state-of-the-art tolerance analysis approaches.

A popular variance analysis approach has been the collection of statistical data through simulating many instances of the circuit. The most well known and widely used sample-and-simulate technique is the Monte Carlo analysis. Taguchi's method [START_REF] Taguchi | Performance analysis design[END_REF] and modified Taguchi's method [START_REF] Errico | Statistical tolerancing using a modification of taguchi's method[END_REF], which are based on deterministic sampling, have been proposed as more reliable alternatives to the Monte-Carlo analysis. However, the large number of required simulations (3 n ) prohibits their use beyond a handful of parameters.

Due to its tractability, worst-case min-max analysis based on a sensitivity-weighted addition of tolerance windows has been intensively used in many fields [START_REF] Tian | Worst case tolerance analysis of linear analog circuits using sensitivity bands[END_REF]. Hierarchical analysis based on Monte-Carlo and sensitivity analyses has also been introduced as a compromise [START_REF] Koskinen | Hierarchical tolerance analysis using statistical behavioral models[END_REF].

Recently, variance analysis of path delays has attracted much attention in the digital domain [START_REF] Agarwal | Statistical timing analysis for intra-die process variations with spatial correlations[END_REF][START_REF] Choi | Timing analysis of combinational circuits including capacitive coupling and statistical process variation[END_REF][START_REF] Orshansky | A general probabilistic framework for worst case timing analysis[END_REF], where the problem typically involves three levels of hierarchy (transistor level, gate level, and path level) and a single parameter (delay). While computing the delay variance of a path, correlations among the gate delays have been taken into account through their covariances [START_REF] Choi | Timing analysis of combinational circuits including capacitive coupling and statistical process variation[END_REF][START_REF] Orshansky | A general probabilistic framework for worst case timing analysis[END_REF]. In [START_REF] Orshansky | A general probabilistic framework for worst case timing analysis[END_REF], sensitivity analysis is used to derive the relations between gate delays and process variables. However, instead of min-max analysis, the covariance of each pair of gate delay is computed based on the sensitivity information and this information is used to update the tolerance bounds computed under the independent assumption. A similar approach is taken in [START_REF] Choi | Timing analysis of combinational circuits including capacitive coupling and statistical process variation[END_REF] while the covariances are assumed to be fixed between the gates and specified by the system designer. In [START_REF] Liu | Hierarchical analysis of process variations for mixed-signal systems[END_REF], we have developed a theoretical basis for variance analysis of analog circuits. The approach in [START_REF] Liu | Hierarchical analysis of process variations for mixed-signal systems[END_REF] is based on computation of correlation coefficients along a hierarchical path.

In this paper, we propose a variance analysis method that is specifically geared towards information re-use. We model the circuit representation as a tree and re-use the variance information at each hierarchical layer. The method takes advantage of identification of correlation loops in the hierarchical structure so as to avoid re-calculation of previously computed values. Such information re-use is extremely helpful during design iterations where only minor changes are made to the circuit.

Hierarchical Variance Analysis

In analog circuits, a small number of transistors constitute building blocks (such as current mirrors, diff-amp pairs, etc.), several building blocks constitute modules (such as an OPAMP), and modules are connected in a network to build the circuits (such as a filter or an ADC). The relations among the circuit parameters at various levels can be defined either through approximate behavioral models, or through simulator-based models, such as the sensitivity analysis. Thus, there is inherently a hierarchical construct followed during the design process. The behavior of the circuit can be generalized as follows:

P i (r) = f i (r) (P 1 (i-1) , P 2 (i-1) , ..., P NP (i-1) i-1
)

where P i (r) denotes the ith parameter at the hierarchy Level-r, f i (r) denotes its functional relation in terms of parameters one level down in the hierarchy, and N P r denotes the number of parameters at the hierarchy Level-(r).

In general, the covariance analysis for the linear system with normally distributed inputs can be conducted hierarchically through chainwise matrix multiplications [START_REF] Bilodeau | Theory of Multivariate Statistics[END_REF]:

Let A : R n → R m , linear, and

x ∼ N n (µ, Σ). T hen, Ax ∼ N m (Aµ, AΣA ).
The coefficients of the matrix A can be formed by the sensitivities between adjacent levels since first order Taylor approximation always results in a linear relation between variances in adjacent levels. As a result, the covariance matrix can be propagated hierarchically.

Although the matrix multiplication approach can provide the variance and covariance information for each parameter, the following drawbacks limit its usage in hierarchical variance analysis for analog circuits.

• The covariance matrices only provide lumped information about the correlation between variables. It is hard to determine the contribution in variance for each lower level variable with respect to the top level variances. However, variation analysis in analog circuits requires the auxiliary information to help determine which critical sub-circuits to focus on during design iterations.

• The covariance matrix approach lacks the flexibility to trade off the computational complexity and the accuracy. As the number of variables in each level increases, time complexity becomes the major concern in variance analysis.

• In the covariance matrix based approach, each time the hierarchical construct is changed, all the affected values need to be recomputed. This need for the complete recomputation arises from the inability to destruct the components of the variance of each parameter. During the design iterations, an analog circuit is modified multiple times and the resulting hierarchical construct can be slightly different from the previous one. It is desired to save some computational effort on the unchanged part of the hierarchical construct.

In order to overcome these challenges, we propose a hierarchical approach where structural information on the circuit is tracked to enable information re-use during design iterations. In our earlier work [START_REF] Liu | Hierarchical analysis of process variations for mixed-signal systems[END_REF], we have derived expressions for computing the variances of parameters at each hierarchical level of the circuit:

σ 2 P j (i) = k (S P j (i) P k (i-1) ) 2 σ 2 P k (i-1) + r (0) i (0) r (1) i (1) ,j (1)
...

r (n-2) i (n-2) ,j (n-2) • r (n-1) i (n-1) =j (n-1) S P i (n) P i (n-1) ...S P i (1) P i (0) S P i (n) P j (n-1)
...S P j (1) P i (0)

σ 2 P i (0) (1)
where,

S P j (k) P j (k-1) = ∂P j (k) ∂P j (k-1) µ P 1 (k-1) ,µ P 2 (k-1) ,...,µ P r (k-1)
represents the first order sensitivity of a Level-k parameter

P i (k) to a Level-(k -1) parameter P i (k-1)
. Section 4 presents a detailed discussion of this mathematical expression. In this work, we concentrate on a graph based interpretation and an efficient implementation of this processing scheme. We use a tree model for the hierarchical structure of the circuit. We focus on a layered variance analysis in which tracing to lower levels is avoided if the necessary information has been processed previously.

Modelling the Circuit Behavior and Parameter Correlations

Top-down design flow in analog circuits allows for a well-structured hierarchical representation of the circuit behavior. In our circuit model (illustrated in Figure 1), the lowest level of the hierarchy consists of process (e.g. the oxide thickness, the dopant concentration) and layout (e.g. the transistor width, the transistor length) parameters. The next level of hierarchy consists of DC parameters such as the bias current, followed by small signal parameters, such as the transconductance and the output resistance. The nodes in the circuit graph in Figure 1 correspond to the parameters. An edge between two nodes indicates a functional dependency. The weight of that edge is assigned as the sensitivity between the two parameters. We call these weights correlation coefficients.

For the relations between process and layout parameters (Level-0 parameters) and DC parameters (Level-1 parameters), we use the SPICE circuit simulator and perturbationbased sensitivity computation. For the relations among the parameters at higher levels of hierarchy, we use analytical models. Such analytical models have been derived and used by circuit and system designers [START_REF] Razavi | Design of Analog CMOS Integrated Circuits[END_REF][START_REF] Johns | Analog Integrated Circuit Design[END_REF], and extensively in analog design automation [START_REF] Antao | Archgen: Automated synthesis of analog systems[END_REF][START_REF] Alpaydin | Multi-level optimisation approach to switched capacitor filter synthesis[END_REF].

Once the hierarchical circuit model is constructed, the variances of parameters at levels above Level-0 need to be computed. At a first glance, this computation step involves nothing more than computing the variance of a response function using the variances of its input parameters. However, as the parameters at intermediate hierarchical nodes are correlated, simple statistical analysis is not sufficient for accurate computation.

Adjacent level correlations

In this work, we process the correlation between two same-level parameters through their dependency on lowerlevel parameters. This dependency is determined using the hierarchical structure of the circuit representation and the correlation coefficients. A correlation between two Level-i parameters exists if they have functional dependency on at least one common lower level parameter.

As an example, consider the hierarchy in Figure 2. The correlation between F and G is due to their dependence on the Level-0 parameter, Y. Thus, for two parameters to be correlated, there has to be a path from both of them to at least one common Level-0 parameter. In prior work, correlations among intermediate parameters are handled using covariances that are either obtained through simulations [START_REF] Choi | Timing analysis of combinational circuits including capacitive coupling and statistical process variation[END_REF] or assumed given [START_REF] Orshansky | A general probabilistic framework for worst case timing analysis[END_REF]. (Simulation-based derivation of covariance is primarily limited to 3 levels of hierarchy.) In this work, we aim at efficient processing of these correlations using the circuit representation.

Layered Processing of Variances

The variance of Level-1 parameters can easily be computed using the sensitivity information since there is no correlation among Level-0 parameters. However, for the parameters at Level-2 and above, correlations need to be taken into account. Without loss of generality, consider a partial circuit graph given in Figure 3. Parameter P i (j) at Level-i can only be processed after processing of all parameters at Level-(i -1) that P i (j) has functional dependence on. The processing of parameter P i (j) involves the following steps:

• Calculation of the variance based only on the variances of parameters at Level-(i -1) on which P i (j) has functional dependency. This step corresponds to the first term in Equation 1 and only involves the parameters one level down in the hierarchy. For parameters at Level-2 and above, correlations need to be taken into account.

• Calculation of the correlation correction term (CCT)

for the variance of P i (j) . This step corresponds to the remaining term in Equation 1and involves all parameters at Level-(i -1) and all the parameters at Level-0 that cause the correlation. Although straightforward as a mathematical expression, direct implementation of this step can be time-consuming since the circuit graph needs to be traversed through Level-0 multiple times for each parameter being processed.

Loop tracing for CCT calculation

In order to avoid processing of the same information multiple times, we aim at using the connectivity information in the circuit graph in a more efficient manner. Let us first make several definitions: Definition: A path between two parameters at distinct levels of hierarchy constitutes a correlation path. Definition: For a parameter at Level-i, all parameters at lower levels that can reached through two distinct paths constitute its correlation parameters. Definition: Two distinct paths between a parameter at Level-i and one of its correlation parameter constitute a correlation loop.

Equation 1 indicates that for a parameter, P i (j) , and each of its Level-0 correlation parameters, the correlation paths need to be traced down to Level-0 in pairs. Equation 1 also requires that parameters at (Level-(i -1)) in these path pairs be distinct. The reason for this condition is illustrated in a partial circuit graph shown in Figure 4. There exist more than one correlation paths between the parameter under processing, H, and the Level-0 parameter, X 2 , making X 2 a Level-0 correlation parameter for H. As Equation 1 suggests, the correlation paths need to be processed in pairs. The pair of correlation paths shown in Figure 4a have already been processed while processing the parameter, G 1 , thus should be omitted. However, the correlation paths shown in Figure 4b have never been processed, thus should be included in the analysis. The correlation paths shown in Figure 4a clearly include a loop excluding the parameter under processing (H). We call such loops convergent loops. The paths in Figure 4b demonstrate a non-convergent loop. Thus, the implementation of Equation 1 on the circuit graph suggests identification and processing of non-convergent loops (correlation loops) in the circuit graph.

F 3 H G 1 G 2 X 2 X 3 X 1 F 2 F 1 X 1 H G1 G2 X 2 X 3

Correlation through an arbitrary level

Correlation loops need not always go through Level-0 parameters. While one can still use Equation 1 and the aforementioned loop-based technique, it is advantageous in terms of computational efficiency to use the information at higher levels of the hierarchy. In order to achieve this goal, one can re-write the CCT term in Equation 1 by introducing the variances of the intermediate level parameters:

CCT = r (n-2) i (n-2) =1 W P i (n) P i (n-2) σ 2 P i (n-2) + • • • + r (0) i (0) =1 W P i (n) P i (0) σ 2 P i (0) (2)
where

W P i (n) P i (k) = r (k) i (k) =1 • • • r (n-1) i (n-1) =j (n-1) =1 S P i (n) P i (n-1) • • •S P i (k+1) P i (k) S P i (n) P j (n-1) • • •S P j (k+1) P i (k)
is defined as remainder correlation weight between a Leveli parameter P i (n) and a Level-k (k ≤ i -2) parameter P i (k) .

It is called remainder since the complementary part of the correlation coefficients has already been taken into account. Figure 5: A four-layer hierarchical structure Equation 2 suggests that a correction loop may exist between the parameter under processing and another parameter at an arbitrary level of hierarchy. Such loops need to be handled in the same way loops through Level-0 are handled. For example, in a five-layer hierarchical structure as shown in Figure 5a, the remainder correlation weight r H1 X2 corresponds to the correlation loop through H 1 and X 2 , as show in Figure 5b. In this case, the remainder correlation weight relating H 1 and X 2 will be r

H 1 G 2 G 3 F 3 F 2 F 4 F 1 F 5 G 1 X 2 X
H1 X2 = C 1 • C 2 • C 3 • C 4 • C 5 • C 6 .
To calculate a remainder correlation weight, one needs to identify the loops in the partial tree corresponding to the remainder correlation weight and sum up the corresponding products of correlation coefficients along the loops.

Data structure and implementation issues

For VLSI systems with many hierarchical layers, it is likely that the connectivity information will be sparse, since independent variables at one corner of the die are not likely to affect parameters on another corner of the die. Moreover, in the computation of the correction terms, the connectivity information is processed many times. Thus, an efficient data structure for the representation of the connectivity information is a bi-directional tree structure where each node may have multiple children and multiple parents.

The partial circuit graph between a parameter under processing and one of its correlation parameters resembles a directed control flow graph in the sense that it has a single start and a single end point. Identification of loops in control flow graphs is a well investigated topic in compiler research. In [START_REF] Ramalingam | Identifying loops in almost linear time[END_REF], the loop identification is conducted in almost O(k) time, where k is the number of edges in the control graph (the worst case time complexity is O(k 2 )).

The correlation path information provides valuable insight during design iterations. First, the correlation paths provide a means for determining the relative contribution of each sub-circuit. Second, unmodified correlation paths need not be processed during design iterations. Third, for large circuits, we can compare the relative importance of the sensitivity paths and only retain a set of important sensitivity paths at each hierarchical level.

Although the storage of these sensitivity paths may lead to space overhead, the overhead is actually quite small in practice since, for each path, we only need to record the index of the variables the path visits at each level. Also, dynamic programming techniques can be utilized when propagating these paths in a bottom-up manner along the hierarchical construct.

Comparison of computation time and accuracy with prior approaches

Two experiments have been conducted to evaluate the computational time and accuracy of the proposed hierarchical analysis approach. In this section, we provide comparisons with the Monte Carlo analysis, Taguchi's method, the min-max method, and the covariance based method given in [START_REF] Orshansky | A general probabilistic framework for worst case timing analysis[END_REF]. Example 1: A four-layer non-linear hierarchy Figure 6 shows an illustrative four-layer hierarchical construct with non-linear functional relationships between two adjacent layers. To evaluate the computational time and accuracy of the proposed method, means and variances of all the variables in the hierarchical construct are calculated using both the proposed method and various widely used statistical analysis approaches: the Monte Carlo analysis with both a large sample size and a small sample size, Taguchi's method, sensitivity-based min-max analysis, and covariance based analysis. (Covariance is assumed to be user given in [START_REF] Choi | Timing analysis of combinational circuits including capacitive coupling and statistical process variation[END_REF], and is assumed to be based on a constant correlation of 0.3 between each pair of intermediate layer variables in [START_REF] Orshansky | A general probabilistic framework for worst case timing analysis[END_REF]. To enable comparison, the same 0.3 correlation based analysis is used.) In the proposed method, the correlation coefficients are modelled with the first order Taylor series approximations around the nominal. The standard deviation of Level-0 parameters is given as 3.3%.

The computational time, the accuracy for the Level-3 parameter, H, and the accuracy for the Level-2 parameters (G 1 -G 3 ) have been summarized in Table 1. The 50Ksample Monte Carlo (MC) method is chosen as the base case for all accuracy comparisons.

The poor accuracy of the min-max and covariance-based approaches in conducting the hierarchical analysis for analog circuits limits their usage even though they can be done with relatively small computational efforts. These results

G 1 G 2 G 3 F 1 F 2 F 3 F 4 F 5 X 2 X 3 X 4 X 5 X 6 X 7 X 1 X 8 H F1=X7X2+X 2 2 +X2X3 X 6 X 4 F2=X 2 1 +X3X6+X 2 4 G1=F1F2+F 2 3 F3=X2X6+X3X4+X 2 8 G2=F 2 2 +F 2 4 F4=X1X8+X 2 4 +X5X7 G3=F4F5+F2F3+X 2 1 F5=X 2 1 +X 2 7 +X 2 4 /X5 H=G 2 1 +G2G3
Figure 6: A system of hierarchy with non-linear functional relations between two adjacent layers confirm that the covariance-based approaches developed for digital timing analysis [START_REF] Orshansky | A general probabilistic framework for worst case timing analysis[END_REF][START_REF] Choi | Timing analysis of combinational circuits including capacitive coupling and statistical process variation[END_REF] may not be easily extended to analog circuits where correlations may substantially vary between each pair of intermediate level parameters. Assuming a fixed correlation for all intermediate parameter pairs is effective in timing analysis since propagation delay is the only parameter and delays through various gates will have similar correlations as they consist of the same (or similar) transistor combinations. In the analog and mixed-signal circuits, however, this condition does not hold any more, thus a constant correlation modelling approach results in much higher errors. Example 2: A differential amplifier circuit Figure 7 shows a current mirror load MOS differential amplifier circuit that consists of 3 NMOS and 2 PMOS transistors. The corresponding hierarchical construct is shown in Figure 8, where Av corresponds to the low frequency gain. In this example, six process variables (W, L, V to , U o , λ, T ox ) for each of the five transistors are chosen as independent Level-0 parameters with given variances. Level-1 parameters consist of DC parameters (I D ,V DS ,V GS ). The small-signal model parameters (g m , g o ) constitute the set of Level-2 parameters, since they are usually expressed as functions of DC parameters. The top level (Level-3) of the hierarchy consists of the performance parameters of the circuit, such as gain (A v ) and cut-off frequency (f c ). In this experiment, only the variance of the gain is investigated. The variances of other performance parameters can be obtained in the same manner. As in the previous example, correlation coefficients between Level-1 and higher level parameters are derived using first order Taylor Series approximations around the nominal. Unlike the previous example, the correlation coefficients between Level-0 and Level-1 are obtained through a perturbation based approach. The standard deviation of Level-0 parameters is assumed to be 0.67%. The goals of this experiment are (a) to evaluate the impact of approximate modelling (1st order approximation in this case) on more complicated functional relations, (b) to compare the proposed method to other state-of-the-art approaches for a practical circuit, (c) to illustrate the possibil- Table 2 shows the accuracy and computation time comparison for the parameters at Level-1 and above for various approaches. A 12K-sample Monte Carlo simulation is taken as the base case for all comparisons. Taguchi's method is excluded from the experiment as the required sample size (3 30 ) prohibits its use in this circuit with 30 Level-0 parameters. While tracking the results of the Monte Carlo simulation within 10% error, the proposed method only takes a fraction of computation time required by it. On the other end of the spectrum, the proposed analysis is highly accurate compared to the min-max or covariance approaches with only a small increase in computation time. Thus, even though the functions governing the relations between variables in adjacent layers are much more complicated compared to the first example, a similar conclusion can be drawn from this example. The proposed method provides both computational efficiency and accuracy.

Parameter G1 G2 G3 H Time Method % % % % s MC-50K - - - - 2.5×10

Conclusion

In this paper, we present a hierarchical method to analyze the impact of process variations on circuit parameters. We model the hierarchical structure of the circuit as a tree wherein relations among circuit parameters are modelled as weighted edges. We start from Level-1 parameters and process each parameter at higher levels. Through identification of correlation loops, we avoid reprocessing of the same information to save time. The information reuse is particularly important during design iterations. In order to enable such reuse, we keep track of the correlation paths in the system and process only the updated information.

Two experimental cases have been used to evaluate the efficiency of the proposed approach: a highly non-linear hi- samples is taken as the base case.

erarchy of random functions, and a 5-transistor differential amplifier circuit. The results indicate that the proposed analysis provides almost the same computational efficiency as the simple min-max approach, and fixed-covariance based approach while providing much higher accuracy (9% for the proposed approach for the top-level parameter compared to 1585% for the min-max analysis and 414% for the fixed covariance approach [START_REF] Orshansky | A general probabilistic framework for worst case timing analysis[END_REF].).

Figure 1 :Figure 2 :

 12 Figure 1: Hierarchical graph model of the circuit behavior

3 P 1 P 2 P 3 P 4 P 5 P 2 P 3 P 4 P 5 P 6 Figure 3 :

 31235234563 Figure 3: Determination of correlations and processing of variance information in a hierarchical manner

H to X that have not yet been processed 2 Figure 4 :

 24 Figure 4: A partial circuit graph to process the CCT in the variance of H due to the correlation parameter X2

Figure 7 :

 7 Figure 7: A differential amplifier circuit ity of combining the analytical derivation and the perturbation based computation of correlation coefficients.Table2shows the accuracy and computation time comparison for the parameters at Level-1 and above for various approaches. A 12K-sample Monte Carlo simulation is taken as the base case for all comparisons. Taguchi's method is excluded from the experiment as the required sample size (3 30 ) prohibits its use in this circuit with 30 Level-0 parameters. While tracking the results of the Monte Carlo simulation within 10% error, the proposed method only takes a fraction of computation time required by it. On the other end of the spectrum, the proposed analysis is highly accurate compared to the min-max or covariance approaches with only a small increase in computation time. Thus, even though the functions governing the relations between variables in adjacent layers are much more complicated compared to the first example, a similar conclusion can be drawn from this example. The proposed method provides both computational efficiency and accuracy.

Figure 8 :

 8 Figure 8: Hierarchical graph for the differential amplifier

Table 1 :

 1 Comparison for the accuracy of variance computation.

	3

Monte Carlo analysis with 50K samples is taken as the base. Taguchi's method requires 3 8 (6561) samples.

Table 2 :

 2 Results for Example 2. Monte Carlo analysis with 12K

	σ	MC-12K Min-max Covariance Proposed
	% ID 1	-	681.1	-0.3	2.0
	% ID 4	-	698.9	-0.3	2.0
	% gm 1	-	1641.3	470.0	3.1
	% gm 2	-	1639.1	468.5	3.6
	% gm 3	-	1623.8	465.5	2.6
	% gm 4	-	1610.8	459.8	1.5
	% go 1	-	695.8	6.9	2.3
	% go 2	-	651.9	1.0	-3.4
	% go 3	-	695.0	6.8	2.1
	% go 4	-	625.2	-2.6	-6.7
	% go 5	-	574.6	-7.1	3.0
	% Av	-	1585.4	414.3	9.4
	Time(s) 5.0×10 4	1.7	1.7	2.4

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)


Hratch Mangassarian

Mohab Anis
email: manis@vlsi.uwaterloo.ca

On Statistical Timing Analysis with Inter-and Intra-die Variations

In this paper, we highlight a fast, effective and practical statistical approach that deals with inter and intra-die variations in VLSI chips. Our methodology is applied to a number of random variables while accounting for spatial correlations. Our methodology sorts the Probability Density Functions (PDFs) of the critical paths of a circuit based on a confidence-point. We show the mathematical accuracy of our method as well as implement a typical program to test it on various benchmarks. We find that worst-case analysis overestimates path delays by more than 50% and that a path's probabilistic rank with respect to delay is very different from its deterministic rank.

Introduction

Variability is posing an increasing challenge to timing analysis of VLSI designs implemented in the nanometer regime. Variations can be either environmental or physical. Environmental variations are due to unpredictable operating conditions such as the power supply voltage V dd and the chip's temperature, whereas physical variations arise during the fabrication of the chip and include permanent discrepancies in the oxide thickness tox, effective channel length L ef f and interconnect dimensions. These variations appear at different levels of the manufacturing process: inter-lot, inter-wafer, inter-die and intra-die. Traditionally, discrepancies in VLSI chip parameters have been accounted for using several cases (best-case, nominal and worst-case) in Static Timing Analysis. Worst-case analysis is a deterministic timing analysis that was a fair assumption to make until a few years ago, when intra-die parameter mismatches were considered to be negligible in the presence of the more significant interdie variations. Today, intra-die variations must be accounted for, as worst-case timing analysis is highly potential-cutting [START_REF] Boning | Design of High-Performance Microprocessor Circuits[END_REF].

Physical parameters are susceptible to random variations from their nominal values, and therefore should be modeled as Random Variables RVs. If the circuit to be analyzed contains thousands of gates, and because of the existing spatial correlations between RVs of the same kind, thousand-variable Joint PDFs (JPDFs) should be used for each RV kind. Getting all these joint distributions is computationally impossible for some reasonable discretization accuracy, let alone the exponential run-time needed for computing any path PDF using such JPDFs. If the number of discretization points of the JPDFs per dimension is equal to QUALITY, then an N-variable JPDF will have QUALITY N discretization points, making the computational complexity for computing the delay PDF of that path (QUALITY N ). Hence, the exact numerical solution for finding the PDF of the critical path is unrealistically complex.

Recently, a number of techniques have been developed to perform statistical timing analysis. They can be broadly classified into full-chip analysis and path-based analysis. Full-chip analysis maps the whole block to an acyclic graph and strives to propagate and merge the PDFs of the gate delays, in order to get the PDF of the critical path of the circuit [START_REF]Statistical Timing Analysis Using Bounds and Selective Enumeration[END_REF][START_REF] Devadas | Statistical Timing Analysis of Combinational Circuits[END_REF][START_REF] Jyu | Statistical Timing Optimization of Combinational Logic Circuits[END_REF][START_REF] Najm | Statistical Timing Analysis Based on a Timing Yield Model[END_REF][START_REF] Tongsima | Optimizing Circuits with Confidence Probability Using Probabilistic Timing[END_REF][START_REF] Jing-Jia | Fast Statistical Timing Analysis by Probabilistic Event Propagation[END_REF][START_REF]Computation and Refinement of Statistical Bounds on Circuit Delay[END_REF][START_REF]Statistical Timing Analysis for Intra-die Process Variations with Spacial Correlations[END_REF]. This method has an inherently exponential run-time with circuit size. Therefore, in the pursuit of reducing the run-time, advocators of full-chip analysis have to use primitive delay models, assume that gate delay PDFs are given [START_REF] Devadas | Statistical Timing Analysis of Combinational Circuits[END_REF][START_REF] Jyu | Statistical Timing Optimization of Combinational Logic Circuits[END_REF], neglect parameter correlations [START_REF]Statistical Timing Analysis Using Bounds and Selective Enumeration[END_REF][START_REF] Devadas | Statistical Timing Analysis of Combinational Circuits[END_REF][START_REF]Computation and Refinement of Statistical Bounds on Circuit Delay[END_REF] and draw on inaccurate approximations [START_REF] Jing-Jia | Fast Statistical Timing Analysis by Probabilistic Event Propagation[END_REF]. And even with this, they sometimes give bounds for the delay PDF and not the PDF itself [START_REF]Statistical Timing Analysis Using Bounds and Selective Enumeration[END_REF][START_REF]Computation and Refinement of Statistical Bounds on Circuit Delay[END_REF][START_REF]Statistical Timing Analysis for Intra-die Process Variations with Spacial Correlations[END_REF]. In the path-based approach, the PDF of one path is calculated at a time [START_REF]Statistical Delay Computation Considering Spacial Correlations[END_REF][START_REF] Gattiker | Timing Yield Estimation from Static Timing Analysis[END_REF]. First, Deterministic Timing Analysis is performed on the circuit in order to find some upper percentage of nominally critical paths. Then, each candidate critical path is analyzed statistically, with the aim of computing its delay PDF. The probabilistic critical path is decided by comparing some confidence-point on the PDFs, such as the 3σ points. This method allows for more complex delay and interconnect models [START_REF] Gattiker | Timing Yield Estimation from Static Timing Analysis[END_REF]. Some common assumptions and drawbacks, both in full-chip and path-based statistical analysis, include the consideration of only 1 RV (usually effective channel length) [START_REF]Statistical Timing Analysis for Intra-die Process Variations with Spacial Correlations[END_REF][START_REF]Statistical Delay Computation Considering Spacial Correlations[END_REF], the use of 1 correlation-space at most (usually spatial correlations) [START_REF]Statistical Timing Analysis for Intra-die Process Variations with Spacial Correlations[END_REF][START_REF]Statistical Delay Computation Considering Spacial Correlations[END_REF], the restriction to a certain kind of input PDF (usually Gaussian) [START_REF]Statistical Delay Computation Considering Spacial Correlations[END_REF][START_REF] Orshansky | A General Probabilistic Framework for Worst-Case Timing Analysis[END_REF], and the exponential run-time needed to get accurate results [START_REF]Statistical Timing Analysis Using Bounds and Selective Enumeration[END_REF][START_REF]Statistical Timing Analysis for Intra-die Process Variations with Spacial Correlations[END_REF][START_REF] Orshansky | A General Probabilistic Framework for Worst-Case Timing Analysis[END_REF].

In this work, we propose a path-based methodology which relaxes some these assumptions while maintaining a run-time of polynomial complexity. Our methodology includes all the steps that lead to finding the probabilistic critical path, from the initial deterministic analysis algorithms to the ranking of all probabilistic critical paths of the circuit. The methodology accounts for several RVs which affect gate delay variations. Spatial correlations for these RVs are taken into consideration. Our findings emphasize the critical path delay overestimation that arises due to worst-case analysis.

Modeling and Assumptions

Gate Delay Using Elmore's Model

In order to perform probabilistic timing analysis later in the paper, we first derive the delay of any path which is composed of gates by using Elmore's model. Based on Elmore's delay model, [START_REF] Wei | Design and Optimization of Dual-Threshold Circuits for Low-Voltage Low-Power Applications[END_REF] showed that the propagation delay of an n-input NAND gate tp is formalized as

tp = 0.345[R N C dN F I (F I -1) + F I R N Cn + R P Cn] ( 1
)
where FI is the number of fan-ins, C dN is the drain capacitance, Cn is the sum of the drain capacitances of this gate at the output node, plus the gate capacitances of all fan-out gates, as well as the wire capacitance, while RN and RP are the NMOS and PMOS onresistances, respectively. For short channel models for RN and RP , tp can be expressed as

tp = 0.345 toxL ef f ox [α( V dd (V dd -V T n ) 1.3 + 1 1.5V dd -2V T n ) + β( V dd (V dd -|V T p |) 1.3 + 1 1.5V dd -2|V T p | )] (2)
where

α = C dN F I (F I -1) + F I Cn µnWn (3) β = Cn µpWp ( 4
)
where ox is the oxide permittivity, Wn/Wp are the NMOS/PMOS channel widths, µn/µp are the NMOS/PMOS mobilities and VTn/VTp are the NMOS/PMOS threshold voltages. Similar delay equations can be found for an inverter, an n-input NOR gate and a 2-input XNOR gate. They all have the same form as (2) with different values for α and β. These are the gates under investigation which constitute all ISCAS85 benchmarks. Thus, a path consisting of N such gates will have a delay equation of the form

t P AT H = N X i=1 {0.345 tox i L ef f i ox [αi( V dd i (V dd i -V T n i ) 1.3 + 1 1.5V dd i -2V T n i ) + βi( V dd i (V dd i -|V T p i |) 1.3 + 1 1.5V dd i -2|V T p i | )]} (5)

Sensitivity Analysis

In order to investigate how much each parameter's variability affects the delay, we performed a first-degree sensitivity analysis for the delay of a 2-input NAND gate, an inverter, a 2-input NOR gate and a 2-input XNOR gate, all with fan-outs of 2, at the nominal values of all parameters. To avoid complicated calculations, and to focus the reader's attention to the approach, we assumed the parameters to be independent and all capacitances to be constant. We compared the values of ∂tp

∂x i |X nom • σx i
for each parameter xi, where X is the vector of parameters and σx i is the standard deviation of xi. 130nm CMOS technology values were used, and the typical variances were obtained from [START_REF] Nassif | Delay Variability: Sources, Impacts and Trends[END_REF]. The parameters that had the most impact on the delay were tox, L ef f , V dd , while VTn and VTp had less effect. Table 1 shows the linear approximations of gate delay variations as a result of varying those 5 RVs by one σx i .

Layering of Correlation-Spaces

In this section, we show a simple yet efficient method developed in [START_REF]Statistical Timing Analysis for Intra-die Process Variations with Spacial Correlations[END_REF], [START_REF]Statistical Delay Computation Considering Spacial Correlations[END_REF] that divides the spatial correlation-space into several levels. We use this technique to get rid of the correlated RVs while keeping the correlation information. The model first replicates the die on several layers, then divides each layer i of the die into 4 i rectangle regions. Using variabilities in L ef f for each layer's partitionsize, the model replaces the L ef f of each gate by a sum of RVs, such

x i ) (σ tox =0.15nm,σ L ef f =15nm,σ V dd =40mV,σ V T n =13mV,σ V T p =14mV)
that the correlation information between two L ef f 's is in number of common RVs they have and their variances. This method was used only for spatial correlations in L ef f . We include other RVs to the technique proposed in [START_REF]Statistical Timing Analysis for Intra-die Process Variations with Spacial Correlations[END_REF], to account for any other chip variations.

Since we are accounting for more than 1 RV, we will first provide a general formulation to the spatial correlation problem. Let χ refer to a parameter (RV) with a known marginal probability distribution, χi refers to a RV of type χ of any partition in a certain level i, and χi,j refers to the RV of type χ and of partition j in level i. χ0 is assigned a PDF whose mean is χ nominal . Each subsequent χi is assigned another PDF, with a mean of zero, such that the PDF of

P L-1
i=0 χi is equal to the PDF of χ, which is given. If χ has a Gaussian PDF, this simplifies to assigning a Gaussian PDF to each χi, such that

σ 2 χ = L-1 X i=0 σ 2 χ i (6)
where L is the number of hierarchical layers, and layer 0 represents the whole correlation space. The RV of a certain gate, χgate for instance, will be set to the sum over all layers of the partition RVs that it belongs to, yielding

χgate = L-1 X i=0 4 i -1 X j=0 ξ(i, j, gate)χi,j (7)
where ξ = 1 only if j is the partition the gate belongs to on level i, and ξ = 0 otherwise. It is the inter-die variations in some χ that decide the chip-mean of χ. The remaining layers correspond to different levels of intra-die variations. Intra-die variations are only deviations around the chip-mean. This is why all layers except layer 0 were assigned zero-mean PDFs. As for the mean of the inter-PDF of some χ (the mean of the PDF of χ0), it is equal to χ nominal , because the average of the means of χ, in all chips of all wafers of all lots, is χ nominal . Using the above spatial correlation description, the path delay ( 5) can be re-written as

t P AT H = N X i=1 {0.345 P i,u,w tox u,w P i,u,w L ef f u,w ox × [αi( P i,u,w V ddu,w ( P i,u,w V ddu,w -P i,u,w V T nu,w ) 1.3 + 1 1.5 P i,u,w V ddu,w -2 P i,u,w V T nu,w ) + βi( P i,u,w V ddu,w ( P i,u,w V ddu,w -P i,u,w |V T p u,w |) 1.3 + 1 1.5 P i,u,w V ddu,w -2 P i,u,w |V T p u,w | )]} (8)
where u and w represent the layer number and the partition number, respectively. Even though the RVs in the equation of a path are considered as independent, we still cannot use this to our advantage in order to find the PDF of any tPATH. This is because the same RV, for instance, say, L ef f 1,2 , can belong to the tp of many gates in the same path. Therefore, we cannot separately calculate the delay PDF of each tp i along a path, since the same RVs may re-occur in the path, causing correlation between tp i 's. In the following section, we will use the Taylor series first-order approximation for the gate delay, in order to linearize part of tp i , and therefore transform (8) into a suitable form for PDF calculation in low run-time complexity.

Taylor Series First-Order Approximation for Path Delay

Let Xi represent the vector of RVs of gate i, Xinter the vector of inter-RVs, ∆Xintra i the vector of intra-RVs of gate i, and for any gate i: Xi = Xinter + ∆Xintra i . Using the Taylor-series first order approximation of tp i taken at the bias point Xinter, the gate delay equation is formulated in vector and scalar forms as follows

tp i (Xinter + ∆Xintra i ) ≈ tp i (Xinter) + ∇tp i |X inter • ∆Xintra i (9)
It should be noted that the mean of Xinter is X nominal , the mean of ∆Xintra i is zero, and that this approximation is accurate if ∆Xintra i Xinter. Since the intra-die RV standard deviations are much smaller than the nominal values of those RVs, the approximation is valid. Replacing the intra and inter RVs of ( 9) by the layer RVs described in the previous subsection, we will get

tp i (Xinter + ∆Xintra i ) ≈ tp i (X0,0) + ∇tp i |X 0,0 • X u,w Xu,w (10)
where X0,0 is Xinter.

All summation RVs are assumed to be independent, both inside each sum and across sums. This is very helpful because, if their coefficients were constants, we would have an intra-part that is simply a linear combination of independent RVs. Finding the intra-delay PDF would then be very easy. However, their coefficients are the partial derivatives of delay evaluated at Xinter, which is a random vector. Hence, we will take a final zeroth-order approximation that evaluates the partial derivatives at nominal value [START_REF]Statistical Delay Computation Considering Spacial Correlations[END_REF], and thus making them constants:

∂tp i ∂χ |X inter ≈ ∂tp i ∂χ |X nominal ( 11
)
yielding

tp i (Xi) ≈ tp i (X0,0) + ai X u,w tox u,w + bi X u,w L ef f u,w + ci X u,w V ddu,w + di X u,w V T nu,w + ei X u,w |V T p u,w | (12)
where ai, bi, ci, di and ei are constants equal to the derivatives of the delay of gate i at nominal value. Thus, the path-delay equation in ( 8) changes when we use this approximation. Using [START_REF] Orshansky | A General Probabilistic Framework for Worst-Case Timing Analysis[END_REF], for a path of N gates, (8) yields

t P AT H = N X i=1 tp i (tox 0,0 , L ef f 0,0 , V dd 0,0 , V T n 0,0 , |V T p 0,0 |) + X u,w au,wtox u,w + X u,w bu,wL ef f u,w + X u,w cu,wV ddu,w + X u,w du,wV T nu,w + X u,w eu,w|V T p u,w | (13)
where au,w, bu,w, cu,w, du,w and eu,w are constants corresponding to layer u and partition w. The first summation is the tPATH inter , while the following terms are the tPATH intra . Therefore, the intradelay of a path is simply a linear combination of RVs. If the input RVs were Gaussian, then finding the intra-PDF would be equivalent to finding its variance

σ 2 t P AT H intra = X u,w a 2 u,w σ 2 tox u,w + X u,w b 2 u,w σ 2 L ef f u,w + X u,w c 2 u,w σ 2 V ddu,w + X u,w d 2 u,w σ 2 V T nu,w + X u,w e 2 u,w σ 2 |V T p u,w | (14)

Convexity Analysis

For the delay equation for the gates used, we calculated the values of

∂ 2 tp ∂x i 2 |X nom • σx i
for each parameter xi. This represents the change of the derivative of delay with respect to that parameter for a one standard deviation change of the parameter. The first derivative of gate delay with respect to those RVs is in the order of tens of pico-seconds per Volts, which is significantly higher that the values of

∂ 2 tp ∂x i 2 |X nom • σx i
for all parameter xi. Therefore, we can conclude that even a worst case change in the derivative, for a 3σ change in those RVs, would still be an order of magnitude less than the actual value of the derivative. Therefore, the convexity is small enough to ensure an acceptable accuracy for this approximation.

Doing the numerical computation of the inter-PDF at once would require a complexity of (QUALITYinter R ), where QUALITYinter is the discretization of the PDFs of the inter-RVs and does not have to be equal to QUALITYintra (the discretization of the PDFs of the intra-RVs), and R is the number of different parameters being varied. One should try to separate as many variables, in order to reduce the complexity of the PDF-computation. It should be noted however, that with more complex models, inter-delay equation will tend to become inseparable in almost all its parameters. In such a case, one should identify the parameters with the least inter-variabilities, and when using the Taylor-series approximation, instead of using their inter-RVs, one should use their nominal values. Another possible compromise of accuracy for faster run-time could be to reduce QUALITYinter.

In the next section, we explain our full methodology, which makes use of the two assumptions analyzed in this section, in order to find and rank some upper percentage of probabilistic critical paths.

Methodology

A high-level description of our methodology flowchart is the following. First, we set the number of RVs (R), the number of layers (L) for the spatial correlation-space, as well as the variabilities of each layer in each RV. The circuit is then mapped to a timing graph, and we evaluate all gate deterministic delays as well as derivatives with respect to all RVs that are being considered, at their nominal values. The inputs to the methodology are a description of the circuit connections, gates, inputs, outputs as well as any necessary information related to correlations, like (x, y) coordinates to compute the spatial correlations. These are one time calculations. Next, the deterministic critical path is found, using Bellman-Ford. We perform the probabilistic timing analysis of the deterministic critical path delay: We find its intra-delay PDF and its inter-PDF, using the techniques described in the previous section, and finally its total delay PDF (accounting for inter and intra variations). From its delay PDF, its standard deviation σC is extracted in order to have an idea about the variability of our circuit. Using that standard deviation and some arbitrary confidence constant C that the user inputs, we find all the next deterministic critical paths that are within C • σC of the first deterministic critical path delay. We perform the same probabilistic timing analysis to each of those paths sequentially. In the end, we can compare any confidence point, such as the 3σ's, on the critical path PDFs, in order to rank them and obtain the probabilistic critical path. Fig. 1 shows the graphical flowchart of our methodology.

Deterministic Delay Computation

The Bellman-Ford algorithm was used to find the deterministic critical path of a circuit. The weight of each edge in the algorithm is to be the nominal delay of the node before the edge, since the graph is simple, directed and acyclic. The worst-case complexity of Bellman-Ford is (|N | × |E|), where |N | is the number of nodes in the graph and |E| is the number of edges. This can be quite big if the graph is very dense, however, the fact that our graph is simple and acyclic makes it highly unlikely to reach that worst-case.

Probabilistic Timing Analysis to Deterministic Critical Path

Probabilistic timing analysis is applied to the deterministic critical path to find its standard deviation. The probabilistic analysis is separated into intra-and inter-delay calculations. The variance of the intra-die variations are formulated using [START_REF] Rabaey | Digital Integrated Circuits[END_REF], and the PDF of the intra-delay is computed (assuming it's Gaussian). The complexity of such a PDF computation is simply (QUALITYintra). For the inter-delay PDF calculation, the first term in ( 13) is used. Finally, the intra-and inter-PDFs are convolved to evaluate the total delay PDF. Supposing that both PDFs have a discretization of QUALITY, the complexity of the convolution is obviously (QUALITY 2 ).

We use the standard deviation (σC ) of the deterministic critical path total delay PDF as an indicator of variability in the circuit. We will choose and analyze all the paths whose delays are larger than D -C • σC , where D is the deterministic delay of the critical path and C is a constant the user specifies. The larger this constant, the more confident we get in that there is no other path in the circuit that is probabilistically longer than the one we got, but the more run-time is needed. In order to find all the next critical paths within CσC of the deterministic critical path, we used the recursive algorithm shown in Fig. 2, where Wi is the weight of edge i and LABELn i is the delay label of node ni, equal to the maximum arrival time to ni from ns, that was calculated using Bellman-Ford. In broad terms, the algorithm starts from a root node, which is n f , and checks for fan-in nodes that have labels larger than the label of the root minus the weight of the edge between them. If there is such a node, than the same thing is repeated for that node as the root. This is repeated until ns is reached for each path. The worst-case complexity of this algorithm is (κ × |E|), where κ is the number of near-critical paths and |E| is the number of edges in the graph.

Finally, probabilistic analysis is done for all the near-critical paths, one path at a time, to find the delay PDF of each. In the end, we rank the paths based on their PDFs by some confidence point. We determine the probabilistic critical path and visualize the change in rank of paths, going from deterministic to statistical analysis. It is important to note that the means of the delay PDFs of the paths are not the same as their deterministic delays, because the inter-delay is not a linear function of RVs and therefore the expected value of the delay, is not the delay of the expected values. The intra-delay calculations for all the paths yield a complexity of (κ × Ω × QUALITYintra 2 ), where Ω is the number of layer-RVs in the path. However, the bottleneck complexity of the methodology comes from the inter-delay PDFs, which have a worst-case complexity of (κ × QUALITYinter R ).

Results and Discussion

We implemented a program that reads the circuit-description as a Design Exchange Format (DEF) file and gets the variability information from the user, in order to generate and rank the delay-PDFs of critical paths. It applies all the steps of the methodology described in the previous section and depicted in Fig. 1. The considered RVs were: tox, L ef f , V dd , VTn and VTp. Their PDFs were assumed to be Gaussian, truncated at their 6σ points. Typical standard-deviations were taken from [START_REF] Nassif | Delay Variability: Sources, Impacts and Trends[END_REF]. The (x, y) coordinates of the gates were extracted from the DEF files in order to account for spatial correlations. We used a 4 layer model along with a fifth random layer, and we divided the total variances equally over all layers. We tested our program on the ISCAS85 benchmark circuits, using 130nm technology nominal values. The computations were performed using PDF discretizations of QUALITYintra = 100 and QUALITYinter = 50. We shall illustrate how these values were chosen as an optimal trade-off between run-time and solution accuracy. Table 2 summarizes the results.

The confidence constant C provides a mean for controlling the number of near-critical paths to consider. We started with a minimum of C = 0.05 (except for c6288) and we increased C until we converged to a critical path that is not changing. In Table 2, we used the minimum value of C that found the correct probabilistic critical path. However, for benchmark c6288, even for a C = 0.005, the number of near-critical paths was more than a hundred thousand, which is unacceptable both in terms of run-time and memory. So we used C = 0.001, which still yielded about 900 paths. Column 7 shows the number of critical paths for each circuit. This number depends heavily on the choice of the confidence C, because all paths with a delay within CσC of the nominal critical delay are considered near-critical, where σC is the deterministic critical path's standard deviation. In addition, the type of benchmark (structure of the graph), as well as the variability and correlations of parameters, affect σC , and therefore the number of critical paths. To see this, one can observe the large difference in critical path numbers, between circuits for which we used the same C, such as benchmarks c880 and c1355: For C = 0.05, the former had 3 near-critical paths, while the latter had 1596. Fig. 3 shows how close the delay PDFs of the 1 st , 798 th and 1596 th paths are in c1355. It should be noted that Column 8 shows the mean of the probabilistic critical path PDF which is very close but not equal to the nominal critical path delay, because the delay is non-linear. Columns 9 and 5 depict the 3σ points for the critical path delay PDFs and the % overestimation of the worst-case analysis over the 3σ points of the probabilistic critical path delay. On average, an overestimation of 55% is found, which demonstrates how unacceptable and persistent the conservatism resulting from worst-case analysis can get. Fig. 4 shows the critical path's intra-and inter-delay PDFs for c432, as well as the convolution of the two, which is the total delay PDF. In addition, the 3sigma point of the total probabilistic (statistical) PDF compared to the worst-case deterministic analysis.

In order to characterize the different effects of inter-and intravariations on the delay, we performed probabilistic timing analysis of c432, for 3 different scenarios of inter-and intra-variances, for the same total RV variabilities. Table 3 demonstrates the results. One can clearly see that the larger the inter-variability is, the larger will be the path-delay standard-deviation. In fact, large inter-variability increases the possibility that all the gates in a path are subject to worst-case variations. Even the number of near-critical paths increases because of the increase of σC .

Table 3. inter-and intra-variations

Column 11 in Table 2 demonstrates the new rank of deterministic critical path after applying probabilistic timing analysis. Some paths remained unchanged (such as c432, c880, c7552), while others were nominally faster paths (notably in c1355 -what used to be the 40th slowest deterministic path is now the critical path in the probabilistic analysis). This is because of the spatial correlations' impact on c1355's topology, increasing the variability in a path, causing their 3σ delays to become very big. The larger the variances and correlations, the more deterministic and probabilistic ranks will tend to differ from one another. Fig. 5 draws the probabilistic rank of the first 100 paths of benchmark c1355 versus their deterministic ranks, when around 1600 near-critical paths were analyzed. Fig. 6 illustrates the same plot for the first 100 paths of benchmark c7552 with the same number of analyzed near-critical paths. We observe that the first plot is considerably far from the first bisector, which means that distances between paths are very small for c1355 compared to the existing amount of variability. On the other hand, in the plot for c7552, we can see that the ranks' changes are minor. This can be explained by the fact that the graph of c1355 is more bushy than that of c7552, therefore paths are much closer in terms of their delays. Thus, path delays are very vulnerable to change ranks for the c1355 case (as depicted in Fig. 5). This is not the case for c7552 where the delay of the paths are more distinctive. Moreover, spatial correlations can be accounted for larger delay variances in c1355, due to its DEF circuit description. In conclusion, it is the topology and placement of the circuit that usually determine changes in critical path ranks.

Finally, the last column in Table 2 shows the run-times of the whole program, from the input-file parsing, to finding the deterministic critical path, to computing their delay PDFs and ranking them. For a path-based approach such as ours, the path delay complexity is directly proportional to the number of near-critical paths to consider, hence those run-times are very strong functions of C. Moreover, the run-time also varies with the graph structure: very close deterministic delays will imply a considerable number of nearcritical paths. Lastly, run-time is a strong function of the discretization points QUALITYintra and QUALITYinter.

In order to measure the trade-off between accuracy and runtime, as far as QUALITYintra and QUALITYinter are concerned, we calculated the delay PDF of c499, using a series of discretization combinations. The most accurate would be the one done with the most discretization points. As an optimal trade-off between accuracy and run-time, we chose the point (QUALITYintra = 100, QUALITYinter = 50) because it yielded an accuracy within 0.009% of the 3σ point with the highest discretization, with a runtime of 0.4 seconds. Thus, the work presented in this paper uses QUALITYintra = 100, QUALITYinter = 50.

One last observation involves the typical minimum value of C. Looking at Table 2, we can see that the biggest C we needed, to find the absolute critical path, was 0.1. This means that, for typical circuits, the probabilistic critical path is within 10% of a typical path-delay standard-deviation from the deterministic critical path. This is a vital advantage for path-based statistical analysis, because it means that typically the number of near-critical paths to consider is acceptable.

Conclusion

We presented a framework for performing statistical timing analysis. Five random variables have been accounted for in spacial correlations. The work presented has a polynomial run time. Our findings confirmed that the worst-case analysis overestimates path delays by more than 50%. Moreover, depending on the circuit's topology and placement information, a path's probabilistic rank with respect to delay could be very different from its deterministic rank.

Figure 1 .

 1 Figure 1. Methodology Flowchart

Figure 2 .

 2 Figure 2. Finding near-critical paths

Figure 3 .

 3 Figure 3. Delay PDFs of the 1 st , 798 th and 1596 th paths in c1355

  50% inter-die, 50% intradie 267.074 35.577 32.674 14.076 54 75% inter-die, 25% intradie 266.889 41.388 39.960 10.778 76

  (a) Intra-delay PDF (b) Inter-delay PDF (c) Total path delay PDF for c432

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Inter-, Intra-, and Total PDF for c432

1 .

 1 Initialization ROOT = n f LABEL ROOT = D -Cσ C FIRST ROOT = true 2. Find Near-Critical Paths (ROOT, LABEL ROOT ,FIRST ROOT ) For all fan-in nodes ni of ROOT If (LABELn i >= LABEL ROOT -Wi) If (FIRST ROOT = true) Add ni to the path Update delay and node count of path Find Near-Critical Paths(ni, LABEL ROOT -Wi, true)

	Else
	Create new path
	Find all previous nodes to ROOT in previous path
	Add previous nodes, ROOT and ni to new path
	Set delay and node count of new path
	Find Near-Critical Paths(ni, LABEL ROOT -Wi, true)
	End If
	End If
	End For

Table 2 .

 2 Results

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)


Marie Dehestru
email: marie.dehestru@erametgroup.com

Moncef Stambouli

Morgan Gohin

Dominique Pareau

Yves Le Quesne

Study of the aggregation of Cyanex 272 and Ionquest 290 by Vapor Pressure Osmometry

The extractant molecule bis-2,4,4-trimethylpentyl phosphoric acid (BTMPP), is the active component of Cyanex 272 and Ionquest 290. BTMPP is known to form aggregates, especially dimers, in nonpolar organic diluents. The aggregation results from strong intermolecular hydrogen bonding between the P=O and P-OH groups. The knowledge of this aggregation is important, as it influences solvent extraction equilibriums by reducing the amount of extractant molecules available for the extraction.

In the present work, aggregation of BTMPP was studied by Vapor Pressure Osmometry (VPO). The average numbers of aggregation were calculated for several solutions containing various concentrations of extractant for both Cyanex 272 and Ionquest 290. The average numbers of aggregation are around 2 -2.2 at 50°C, depending on the concentration of BTMPP. It was concluded that BTMPP forms monomers, dimers and trimers. Then the equilibrium constants of dimerization and trimerization were estimated by a numerical optimization procedure developed on Matlab. The logarithms of the constants of dimerization and trimerization are respectively estimated at 3.0±0.2 and 4.8±0.2 at 50°C (constants in L/mol and L²/mol² resp.). VPO measurements were also conducted at 43°C and 60°C. Only small differences were observed for the average numbers of aggregation and the equilibrium constants calculated at 43°C, 50°C and 60°C. It was concluded that the temperature dependence of the aggregation of BTMPP is low from 43°C to 60°C.

Introduction

The extractant molecule bis-2,4,4-trimethylpentyl phosphoric acid (BTMPP), is the active component of Cyanex 272 and Ionquest 290 (CX272 and IQ290 resp.). BTMPP is widely studied for its hydrometallurgical applications of metals separation, especially cobalt and nickel. It is generally accepted that monobasic organophosphorus acid extractants like BTMPP form dimers in diluents of low polarity. The dimer arises from strong intermolecular hydrogen bonding between the P=O and P-OH groups as shown Figure 1. [START_REF] Ritcey | Solvent extraction, principles and applications to process metallurgy[END_REF], [START_REF] Tait | [END_REF].

The knowledge of the aggregation of an extractant is important, as it influences the extraction equilibriums by reducing the amount of extractant molecules available for the extraction. Most of the studies on solvent extraction by BTMPP consider that BTMPP exists only as dimer in the organic phase. The predominance of the dimer form was demonstrated in toluene by Vapor Pressure Osmometry (VPO) [START_REF] Tait | [END_REF], [3]. The existence of monomers and of other more aggregated forms such as trimer is also suggested [4]. Indirect methods involving biphasic equilibriums were used to estimate the value of the equilibrium constant of dimerization (Table 1). The dimerization constant is optimized together with the constants of the other equilibriums involved (distribution, dissociation…), in order to fit measured data such as pH, distribution of BTMPP or conductivity [5]- [START_REF] Sella | Diphasic acido-basic properties of organophosphorus acids[END_REF].

Table 1: review of the dimerization constants of CX272 or pure BTMPP available in the litterature, with the conditions and the methods used for their estimation.

It was pointed out that dimerization constant cannot be determined precisely from these biphasic methods, because of the too low impact of the variation of the constant on the measured data [7], [START_REF] Sella | Diphasic acido-basic properties of organophosphorus acids[END_REF]. This could explain the dispersion observed between the constants reported in Table 2. Furthermore, some studies are conducted with BTMPP purified from CX272 whereas others use non-purified CX272. Sastre et al. estimated the dimerization constant of BTMPP by VPO [10].

In the present work, the aggregation of BTMPP is studied by VPO at 50°C, temperature commonly chosen for industrial applications. The aggregation constants are calculated by least square optimization method. In order to know the impact of temperature, tests were also conducted at 43°C and 60°C.

Material and methods

Reagents

The extractants Cyanex 272 (Cytec) and Ionquest 290 (Rhodia) are used without further purification. The diluents chosen for the VPO tests are n-hexane (Sigma-Aldrich) for the measurements at 43°C and 50°C, and n-heptane (Panreac) at 60°C, because their boiling points (69°C and 98°C resp.) are adapted to VPO at such temperatures and they have dielectric properties similar to kerosene. We assumed that the length of chain of the alkane diluent does not modify the aggregation of BTMPP. Tri-octylamine (TOA) (Sigma-Aldrich) is used as standard because it is known to be monomeric and to dissolve well in alkane diluents [11]. All the solutions are prepared by weighting.

VPO technique

VPO consists in measuring the difference of resistance ∆R between two thermistors: one containing a drop of pure diluent suspended, and the other a drop of solution (diluent + solute). The two thermistors are located in an atmosphere saturated with diluent vapor. ∆R is directly proportional to the difference of temperature between the thermistors. It is also a function of the difference in vapor pressure of the two drops, and so directly related to the activity of the diluent, or to the apparent molar fraction of solutes in the solution (X app ). X app is the molar fraction of all the solute species present in the solution:

The simplification of the equation of Clausius Clapeyron leads to a proportional relation between ∆R and X app [12].

∆R = K.X app (constant K determined by calibration)
This proportional relation, valid for low X app , can also be used with molalities instead of molar fractions for solutions of low molalities. Although most of the VPO studies work with molalities for convenience, we chose to work with molar fractions that led to better calibrations than molalities because of the high concentrations of our solutions.

VPO experiments are carried out at 43°C, 50°C and 60°C using a Knauer K-7000 Vapor pressure osmometer. The drop size is kept as constant as possible. Due to the time dependence of the signal, the measurements are taken when the maximum of ∆R is attained and stabilized for a short period of time. Three measurements are done for each sample showing a good repeatability. The average of the three measurements is calculated.

Calculation method

∆R are measured by VPO for several solutions containing extractant at various analytical molar fractions, called X tot .

X app are calculated from the ∆R and K. The average numbers of aggregation (N agg ) are obtained from the couples X tot and X app : Nagg = Xtot / Xapp Then an optimization procedure developed on Matlab based on the least square method is used to estimate the equilibrium constants of aggregation from the experimental average number of aggregation. For convenient use of the constants, the molar fractions are converted into molar concentrations.

Results and discussion

Presence of impurities

CX272 and IQ290 are respectively 88% and 95% pure. The presence of impurities influences the estimation of N agg obtained by VPO because they are included in X app and X tot . The main impurity is tri(2,4,4-trimethylpentyl) phosphine oxide. For simplification, we assumed that the phosphine oxide impurities are monomeric, although mixed aggregates phosphine-oxide / BTMPP could be formed, as it was suggested for TOPO and D2EHPA [13]. Other assumptions were tried: similar aggregation behavior of the impurities and BTMPP, formation of mixed aggregates… However these assumptions (results not presented here) led to different values of aggregation constants between CX272 and IQ290.

Comparison of Cyanex 272 and Ionquest 290 at 50°C

Depending on the concentration of BTMPP, its average number of aggregation varies from 1.8 to 2.2 (Figure 2), so it was concluded that BTMPP forms monomers (AH), dimers (A 2 H 2 ) and trimmers (A 3 H 3 ):

K dim K trim
The procedure of optimization developed on Matlab is used to estimate the aggregation constants (for calculations in concentrations): CX272:

Log(K dim ) = 3.1 ± 0.2 Log(K trim ) = 4.9 ± 0.2 (at 50°C) IQ290:

Log(K dim ) = 3.3 ± 0.2 Log(K trim ) = 4.9 ± 0.2 (at 50°C) The optimized aggregation constants allow calculating the experimental average numbers of aggregation with good agreement. However, we observe that the method is very sensitive to small variations of the experimental average numbers of aggregation. Moreover, the aggregation constants are only approximations due to the presence of impurities and the assumption on their aggregation. Therefore the incertitude on the logarithm of the optimized aggregation constant is around 0.2.

Influence of temperature

Experimental tests were also conducted at 43°C (in n-hexane) and 60°C (in n-heptane) with CX272, in order to see the influence of such temperature variations (Figure 3). The differences obtained between the three couples of constants are low, comparable to the incertitude of the method (Table 2). So it is concluded that temperature variations from 43 to 60°C have small influence on the aggregation. At first approximation, it can be considered that aggregation is independent on temperature from 43°C to 60°C, taking the average of the constants calculated at each temperature:

TOA, 50°C CX272, 50°C IQ290, 50°C 1 1,2 1,4 1,6 1,8 2 2,2 2,4 2,6
Log(K dim ) = 3.0 ± 0.2
Log(K trim ) = 4.8 ± 0.2 Knowing the equilibrium constants of dimerization and trimerization, the distribution of BTMPP between monomers, dimers and trimers can be calculated at various analitycal concentrations of BTMPP (Figure 4).

Contrarily to the common assumption that BTMPP is only dimeric, the proportion of BTMPP present in the form of monomer or trimer can be elevated. For concentrations higher than 0.15mol/L, the proportion of BTMPP present in the dimers is lower than 60%, and higher than 35% in the trimers.

TOA

Conclusion

The aggregation of BTMPP has been studied by VPO. The average numbers of aggregation of CX272 and IQ290 are measured for various concentrations of extractant at 43°C, 50°C, 60°C. As the average number of aggregation reaches 2-2.2, it was concluded that BTMPP exists as monomer, dimers and trimers. The equilibrium constants of dimerization and trimerization are estimated by least square optimization method. Nearly similar values of equilibrium constants were obtained for both CX272 and IQ290 when assuming impurities of the commercial extractants are monomeric. No significant influence of temperature has been detected on the average aggregation number in the range 43 -60°C. Log(K dim ) and Log(K trim ) are respectively estimated at 3.0 ± 0.2 and 4.8 ± 0.2. At low concentrations of BTMPP (0.005-0.1mol/L), the extractant is mainly present as dimer. However at higher concentrations the fraction of BTMPP that is present as trimer becomes elevated.

Figure 1 :

 1 Figure 1: Schema of dimer BTMPP

Figure 2 :

 2 Figure 2: Calibration curve with TOA, VPO signals and average numbers of aggregation obtained at 50°C with CX272 and IQ290 at various concentration of extractant

Figure 3 :

 3 Figure 3: Calibration curves with TOA, VPO signals and average numbers of aggregation obtained at 43°C, 50°C and 60°C with CX272 at various concentration of extractant

Figure 4 :

 4 Figure 4: fraction of BTMPP that is present in the form of monomer, dimer or trimer, as a function of the analytical concentration of BTMPP.

Table 2

 2

	TOA
	CX272
	and IQ 290
	: Logarithm of BTMPP equilibrium
	constants of aggregation estimated with
	CX272 at 43°C, 50°C and 60°C.