Simulation

1. Simulation • Objective • Given a circuit, to find out whether it behaves in the desired manner or not • How accurate do we want to be? • Device-level simulation • Circuit-level simulation (SPICE) • Timing simulation/logic simulation • etc.

2. Circuit simulation • Given a netlist of elements • transistors, resistors, capacitors, transmission lines, etc. • Find currents and voltages at all points of interest • Example R=1 V=1 R=1

3. More formally... • Write the network equations • Topological equations • Kirchoff’s Current Law (KCL): current entering each node is 0 • Kirchoff’s Voltage Law (KVL): voltage around each loop is 0 • Device equations: describe the behavior of the device, e.g., V= IR • Solve the system of equations to find V, I everywhere

4. 2 3 1 b2 b3 b5 b4 b1 4 The incidence matrix • Properties • nxb matrix (n = number of nodes, b = number of branches) • Each column has exactly one +1 and one -1 • Rank = n-1 for a connected graph

5. b2 b3 2 3 1 b5 b4 b1 4 KCL and KVL • KCL: A I = 0 (A: n x b,I: b x1) • KVL: vb= ATvn (vb: b x 1, vn: n x 1)

6. b2 b3 2 3 1 b5 b4 b1 4 Device equations • General form: • Ib = Y vb+ s • If b1 is a current source, others are resistors

7. Circuit equations • (Use reduced incidence matrix now – remove ground node) • Set of equations describing the circuit • KCL: AIb = 0 • KVL: vb = ATvn • DE: Ib= Yvb + s • Some simple algebra: • DE + KVL: Ib= YATvn + s • DE + KVL + KCL: A (YATvn + s) = 0 AYATvn = - As • This is the nodal formulation of circuit equations

8. Modified nodal formulation • Nodal formulation cannot handle all device types (e.g., voltage sources) • DE’s can be of two types • Type 1 branches: Ib1= Yvb1 + s1 • Type 2 branches: Z Ib2+ Gvb2 = s2 • Rewrite KCL and KVL slightly • KCL: • KVL:

9. Modified nodal formulation (Contd.) • Set of equations describing the circuit • KCL: A1Ib1+ A2Ib2 = 0 • KVL: vb1 = A1Tvn vb2 = A2Tvn • DE: Ib1= Yvb1 + s1 Z Ib2+ Gvb2 = s2 • Simple algebra again leads to

10. Writing the circuit equations • Multiplying matrices can be messy! • Can write circuit equations by inspection by creating element “stamps”; need not multiply any matrices • Basic circuit elements • Conductance: i = g v • Current source: i = k • Voltage source: v = k • Controlled sources • CCCS: i = k i’ • CCVS: v = k i’ • VCCS: i = k v’ • VCVS: v = k v’

11. n1 n2 g Stamp for a conductance • Nodal formulation • DE + KVL: Ib= YATvn + s • DE + KVL + KCL: A (YATvn + s) = 0, i.e., AYATvn = - As • KCL << DE << KVL • Conductance • KCL: • KCL+DE: • KCL+DE+KVL:

12. Element Stamps: Admittances (Type 1) • Conductance g • For AC analysis only • Single frequency excitation • Inductance: admittance of 1/j L • Capacitance: admittance of j  C • Similar stamp: replace g by 1/j L or j  C

13. Independent current source (Type 1) • Current source • KCL: • KCL+DE: • iJ has been eliminated entirely • Does not impact KVL since Ib1= Yvb1 + s1has zero entries in Y • Alternative interpretation: voltage across a current source can take on any value n1 n2 J

14. Independent voltage source (Type 2) • KCL: A1Ib1+ A2Ib2 = 0 • KVL: vb1 = A1Tvn vb2 = A2Tvn • DE: Ib1= Yvb1 + s1 Z Ib2+ Gvb2 = s2 • Current through voltage source is a variable that is not eliminated • Assume: value of the voltage source = K volts • KCL KVL+DE

15. VCCS Device equations: Ijk = 0 Ilm = g Vjk • Type I element: can write as I = YV + s • Stamp j l + v gv - k m Col j Col k

16. b2 b3 2 3 1 b5 b4 b1 0 Using the stamps • b1 is a current source J, others are conductances g2..g5 • All are Type 1 elements • KCL+KVL+DE gives us

17. R=1 V=1 Another example 1 2 • Equations: • Solution (easy to check) v1 = 1, v2 = 0.5, iv = -0.5 R=1 0

18. Solution of nonlinear equations • Solution by iteration • Find linear approximation around a guess solution • Solve resulting system of linear eqs. • Iterate until convergence • Example: • Diode equation: I = f(V) = Io (eV/Vt – 1) • Assume a kthiteration guess (Ik, Vk) • Taylor series expansion I = Ik + [df/dV]k (V-Vk) + higher order terms I V Local linear approximation (neglected)

19. Diode example (contd.) • Linearizing the diode equation • I = f(V) = Io (eV/Vt – 1); df/dV = (I0/Vt) eV/Vt • First order Taylor series expansion I = Ik + [df/dV]k (V-Vk) = Ik + (I0/Vt) exp(Vk/Vt) (V-Vk) = ak V + bk where ak = (I0/Vt) exp(Vk/Vt); bk = Ik - (I0/Vt) exp(Vk/Vt) Vk • Circuit interpretation: • conductance of ak in parallel with a current source of bk conductance ak bk

20. conductance ak + V - + V - conductance g conductance g bk Solution of a circuit containing diodes • Example: • Linearize the circuit as shown earlier • Solve the linear circuit;get (Inewk+1,Vnewk+1) • This lies on the linear approx., not on the I-V characteristic! • Choose (Ik+1, Vk+1) on diode characteristic • Continue iteratively until convergence I (Ik+1,Vk+1) (Ik,Vk) V (Inewk+1,Vnewk+1)

21. “Current and voltage iterations” • While mapping (Inewk+1,Vnewk+1) to (Ik,Vk), need to get reasonable values for both • In previous example: fixed Vk+1= Vnewk+1; found corresponding Ik+1 from diode equation • Fixing Ik+1=Inewk+1 would have been a bad idea – no Vk+1! • Heuristic for a diode: • If Vnewk+1  Vcutin • set Ik+1= Inewk+1; find Vk+1 = Vt log(Ik+1/I0 + 1) • If Vnewk+1 < Vcutin • set Vk+1= Vnewk+1; find Ik+1 = I0 (exp(Vk+1/Vt) – 1)

22. General procedure for nonlinear elements • Intuition: Consider one function F of one variable x F(x) = 0 • First order Taylor series about a guess xk F(x) = F(xk) + [dF/dx]k (x – xk) = 0 • x = xk – {[dF/dx]k}-1 F(xk) (Newton-Raphson update formula) • Can also write this as {[dF/dx]k} x = {[dF/dx]k} xk – F(xk) • Analogy: for a vector of m functions F = [F1 F2 … Fm]T of variables x = [x1 x2 … xn]T Jkx =Jkxk – F(xk) where J = Jacobian matrix = analog of the derivative

23. Definition of Jacobian • For a system of m nonlinear equations F = [F1 F2 F3 … Fm]T in n variables x = [x1 x2 x3 … xn]T, the Jacobian J is defined as

24. Solution of a set of linear+nonlinear equations • Given a set of equations in n variables • A x = b (m linear equations) • F(x) = 0 (n-m nonlinear equations) • Choose initial guess xk, k = 0 • Linearize F(x) = 0 as Jkx =Jkxk – F(xk) to get n-m linear eqs • Solve these and A x = b – together n linear eqs in n vars. • Call this solution xnewk+1; map to xk+1 and increment k • If converged, quit; else go to step 1 • Convergence criterion ||xk+1 – xk|| <  and ||Fk+1 – Fk|| < 

25. Simple example • Solve x + y – 10 = 0, xy – 25 = 0 • (We know that the solution is at x=5, y=5) • Choose a guess (x0,y0) = (6,4) • Linearize nonlinear equation, f(x,y) = xy – 25 = 0 • Jk = [{df/dx}k {df/dy}k] = [yk xk] • Linearized equation: • [yk xk] [x y]T = [yk xk] [xk yk]T – f(xk,yk), i.e., yk x + xk y = xk yk + 25 • Iteration 1: Solve x+y=10;4x+6y=49  x1 = 5.5, y1 = 4.5 • Iteration 2: Solve x+y=10;4.5x+5.5y=49.75  x2= 5.25, y1 = 4.75 • and so on… • Converges towards (5,5)

26. Another problem in this example! • At the solution (5,5), the linearized equation is • 5x+5y = 50, or x+y=10 – identical to the other equation! • Solution will be anywhere along this line • In this case, this happened because we got (un)lucky • However, this illustrates a problem of ill-conditioned equations • Example: x+y = 10 and (1+) x + y = 10 • A solution exists, but suffers from numerical problems for small epsilon • As we will see later: this can cause pivots to become very small in LU factorization, leading to inaccuracies under finite machine precision

27. Typical application to circuits • Linearize elements one by one based on guess • Generate stamps for these linearized elements • For the diode, this is the stamp for conductance ak and current source bk • Combine with stamps for other (type 1 and type 2) elements • Solve equations • Find values for next iterate and repeat until convergence

28. Solving these equations • Given a system of linear equations Ax = b (A: n x n; x,b: n x 1) Solve for x • Simple approaches • Cramer’s rule - requires determinant computation: expensive • Gaussian elimination • Perform row transformations on the A matrix to make it upper triangular • Perform the same row transformations on the RHS vector b • Solve upper triangular system

29. LU factorization: outline • Perform LU decomposition • Write A = L . U • L: lower triangular matrix • U: upper triangular matrix • A is written as a product of L and U • Ax = b becomes L U x = b • Substitute U x = y where y is an intermediate variable (vector) • Solve Ly = b (easy) to find y • Then solve Ux = y to find x

30. Gaussian Elimination (GE) • Step 1 • Row2 – 1/3 Row1 • Row3 – 1/3 Row1 • In other words, m21 = 1/3, m31 = 1/3 • Step 2 • Row3 –1/4 Row2 • i.e., m32= ¼ • Step 3: Backward sub. • x3 = x2 = x1 = 1

31. GE in general • Given A x = b, or • Convert to U x = y, i.e., where i.e., (L y = b)

32. Relation to LU factorization • Perform LU decomposition • Write A = L . U • L: lower triangular matrix • U: upper triangular matrix • A is written as a product of L and U • Ax = b becomes L U x = b • Substitute U x = y where y is an intermediate variable (vector) • Solve Ly = b (easy) to find y • Then solve Ux = y to find x • Previous slide tells us how to find L, U !!

33. Versions of LU factorization • Gauss • Store multipliers in L, find U • 1’s on the diagonals of L • While processing multipliers mk*, update as (k,k) Updated submatrix

34. Computational cost • An update at (k,k) involves a constant number of arithmetic operations on (n – k)(n – k + 1) elements • Cost  k=1 to n-1 (n – k)(n – k + 1) = O(n3)

35. Other versions of LU • Same operations in a different order • Doolittle: 1’s on diagonal of L • Crout: 1’s on diagonal of U • Order of operations: while processing (k,k) (k,k) Updated entries (L, U stored in the same matrix here)

36. (k,k) Doolittle (Need to calculate ukj values before lik values – ukkrequired for lik) Complexity: same as Gauss (same operations, different order) Crout: similar update formulae

37. Pivoting • Cannot find m32– divide by 0 • Does not mean solution does not exist • Can overcome by reordering rows 2 and 3 to get a nonzero on the diagonal – pivoting • (Coincidentally, in this case, that also means that GE is over!) m21= 1 m31 = 1/3

38. Sparsity • Consider a matrix structure • x implies non-zero element • After LU factorization, get • All sparsity lost • Many fill-in’s (0  x) • Can use pivoting to reduce fill-in’s

39. (k,k) When is a fill-in created? • Computations during LU factorization • Fill-in created if is a zero and any one pair of are both nonzeros

40. Pivoting for sparsity - example • Reorder the matrix • Row 1  Row 4 • Column 1  Column 4 becomes • LU factors have the structure

41. Iterative solution of linear equations • For a system of equations • For an initial guess x1(0), …, xn(0), can write (Gauss-Jacobi) (Gauss-Seidel)

42. In matrix form • Write A = L + D + U • L: lower  matrix; U: upper  matrix; D: diagonal matrix • G-J: (L+D+U) x = b implies D x(k+1) = b – (L+U) x(k) i.e., x(k+1) = D-1b – D-1 (L+U) x(k) • G-S (L+D) x(k+1) = b – U x(k) i.e., x(k+1) = (L+D)-1b – (L+D)-1 U x(k) • For both, the update formula has the form x(k+1)= M x(k) + c

43. Convergence of Iterative Methods x(k+1)= M x(k) + c converges when x(k+1) = x(k) = x*, i.e., when x* = (I – M)-1 c • Progress of iterations: • x(1) = M x(0) + c • x(2) = M x(1) + c = M2 x(0) + (M+I) c • x(3) = M x(2) + c = M3 x(0) + (M2+M+I) c • : • x(k) = M x(k-1) + c = Mk x(0) + (Mk-1+…+M+I) c • x(k) converges to x*for large k if • Mk 0 • (Mk-1+…+M+I)  (I-M)-1

44. Convergence (Contd.) • This happens if (M) < 1 where • (M) is the spectral radius of M • (M) is defined as the magnitude of the largest eigenvalue of M • In reality can have convergence if max eigenvalue magnitude is 1 and the eigenvalue is simple (i.e., not a multiple root of the characteristic equation)

45. Back to circuit equations.. • Diagonal dominance of a matrix A |akk|  j=1 to n, j  k |akj| • Diagonal dominance with at least one inequality being strict for a connected circuit implies positive definiteness of A • Resistive circuits with current sources and at least one resistance to ground satisfy this.

46. Transient analysis • Requires the solution of differential equations • Elements such as capacitors • I = C dV/dt : C is constant for a linear capacitor, can be a function of voltage for a nonlinear capacitor • V = L dI/dt : L is constant for a linear inductor, can be a function of current for a nonlinear inductor • First consider the basic problem of solving a differential equation

47. x h t 2h 3h 0 Intuitive way of numerically solving a differential equation • Consider dx/dt = f(x); x(t=0) = x0 • Start an x vs. t plot from given initial value x(t=0) = x0 • Derivative at x0 is f(x0) • In the plot, find x1 at time t=h based on a linear extrapolation using the derivative – thjs is reasonable if h is small • Having obtained x1, find derivative at t=h as f(x1); extrapolate linearly to get x2 and so on. • Update formula • x(t+(n+1)h) = x(t+nh) + h [dx/dt](t+nh) • Or using simpler notation, xn+1 = xn + h dxn/dt • Basically, all this says is dxn/dt = (xn+1-xn)/h! • This is the Forward Euler method

48. A few simple numerical integration methods • Forward Euler method: xn+1 = xn + h dxn/dt • Backward Euler method: xn+1 = xn + h dxn+1/dt • Trapezoidalrule xn+1 = xn + h/2 (dxn/dt + dxn+1/dt) • These are all based on local linear approximations • Can also use more derivatives and get formulas based on approximations that are locally quadratic, cubic, etc.

49. Numerical Stability • Test equation: dx/dt =  x, x(t=0) = x0 • Behavior of this equation: • If  > 0, x   as t   (real systems don’t do this!) • If  = 0, x = x0 for all t • If  < 0, x  0 as t   • (The last two are the meaningful regions) • Necessary condition for any numerical integration formula: must satisfy the above limiting conditions. This is a criterion for stability of the formula. • (Other definitions of stability also exist)

50. Numerical Stability: Forward Euler • Test equation: dx/dt =  x, x(t=0) = x0 • FE: • xn+1 = xn + h dxn/dt, i.e., xn+1 = xn + h  xn = (1+ h ) xn • For a constant time step h, • x1 = (1+ h ) x0 • x2 = (1+ h ) x1 = (1+ h )2 x0 • x3 = (1+ h ) x2 = (1+ h )3 x0 • : • xk = (1+ h )k x0 • As t  , k   and xk must satisfy stability conditions