Transcript A x

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.
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
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
The incidence matrix
1
b2
2
b3
b4
b1
4
3
b5
1 1 1 0 0 0 
2  0  1 1 1 0 
A
3  0 0 1 0 1 


4  1 0 0  1  1
• 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
KCL and KVL
1
• KCL: A I = 0 (A: n x b, I: b x 1)
 i (b1) 
0 0 0 
1 1
 0 
i
(
b
2
)


 0 
0 
 0 1 1 1
i (b3)    
 0 0  1 0 1 
 0 
i
(
b
4
)


  
 1 0 0  1  1 i (b5)  0


• KVL: vb = AT vn (vb: b x 1, vn: n x 1)
 v(b1)  1 0 0 1

 
  v1 
v
(
b
2
)
1

1
0
0

 
 v 2 
 v(b3)   0 1  1 0   

 
  v3
v
(
b
4
)
0
1
0

1

 
 
v(b5)  0 0 1  1 v 4

 

b1
b2
b3
2
b4
4
3
b5
Device equations
1
• General form:
– Ib = Y v b + s
b2
b1
b3
2
b4
4
• If b1 is a current source, others are resistors
0
0
0   v(b1)   I 
 i (b1)  0 0

 

  
0
0  v(b 2) 0
i (b 2) 0 g 2 0
i (b3)   0 0 g 3 0
0   v(b3)   0

 

  
i
(
b
4
)
0
0
0
g
4
0
v
(
b
4
)

 

 0 
i (b5)  0 0
0
0 g 5  v(b5)  0

 
3
b5
Circuit equations
• (Use reduced incidence matrix now – remove ground node)
• Set of equations describing the circuit
– KCL:
– KVL:
– DE:
A Ib = 0
vb = AT vn
Ib = Y vb + s
• Some simple algebra:
– DE + KVL:
Ib = Y AT vn + s
– DE + KVL + KCL: A (Y AT vn + s) = 0
AY AT vn = - A s
• This is the nodal formulation of circuit equations
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 = Y vb1 + s1
– Type 2 branches: Z Ib2 + G vb2 = s2
• Rewrite KCL and KVL slightly
– KCL:
A1
– KVL:
 v b1   A1T 
v    T  vn
 b2   A2 
I 
A2  b1   0
Ib2 
Modified nodal formulation (Contd.)
• Set of equations describing the circuit
– KCL:
– KVL:
– DE:
A1 Ib1 + A2 Ib2 = 0
vb1 = A1T vn
vb2 = A2T vn
Ib1 = Y vb1 + s1
Z Ib2 + G vb2 = s2
• Simple algebra again leads to
 A YAT
 1 T1
 GA2
A2   v n   A1s1 
   

Z   Ib 2   s 2 
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’
Stamp for a conductance
• Nodal formulation
– DE + KVL:
Ib = Y AT vn + s
– DE + KVL + KCL: A (Y AT vn + s) = 0, i.e., AY AT vn = - A s
– KCL << DE << KVL
• Conductance
– KCL:
– KCL+DE:
– KCL+DE+KVL:
n1 ..  1 ..
.. .. ..

n 2 ..  1 ..

.. .. ..
..  : 
.. ig 
0
..  : 
 
..  : 
n1 ..  g ..
.. .. ..

n 2 ..  g ..

.. .. ..
n1 ..  g ..
.. .. ..

n 2 ..  g ..

.. .. ..
..  : 
.. v g 
0
..  : 
 
..  : 
 : 
 g ..  
v
.. ..  n1 
 : 0
 g ..  
 v
.. ..  n 2 
 : 
n1
g
n2
Element Stamps: Admittances (Type 1)
• Conductance g
n1 ..  g ..  g ..
.. .. .. .. ..


n 2 ..  g ..  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
Independent current source (Type 1)
• Current source
– KCL:
– KCL+DE:
n1 ..  1 ..
.. .. ..

n2 ..  1 ..

.. .. ..
n1 ..
..

n 2 ..

..
..
..
..
..
..
..
..
..
..  : 
.. iJ 
0
..  : 
 
..  : 
.. :  J 
.. :  

.. :  J 
   
.. :  
n1
n2
J
• iJ has been eliminated entirely
• Does not impact KVL since Ib1 = Y vb1 + s1has zero entries in Y
• Alternative interpretation: voltage across a current source can
take on any value
Independent voltage source (Type 2)
– KCL:
– KVL:
–
–
–
–
A1 Ib1 + A2 Ib2 = 0
vb1 = A1T vn
vb2 = A2T vn
DE:
Ib1 = Y vb1 + s1
Z Ib2 + G vb2 = s2
Current through voltage source is a variable that is not
eliminated
Assume: value of the voltage source = K volts
n1 .. .. .. .. .. 1   vn1 
:

 
KCL
KVL+DE
 
.. .. .. .. .. ..   : 
n1 ..

..
n 2 ..

..
..
..
..
..
..
..
..
..
..
..
..
..
.. 1   : 

.. ..   : 
 0
..  1  : 

.. ..   : 
 
iV 
n 2 ..

..
..

1
.. ..
.. ..
.. ..
..  1
..
..
..
..
..
..
..
..
:
 
:
 1 vn 2   : 
    
..   :   : 
..   :   : 
   
..   iV   K 
VCCS
j
k
+
v
-
Device equations:
Ijk = 0
Ilm = g Vjk
l
gv
m
• Type I element: can write as I = Y V + s
j  ..
.. .. .. .. ..  v 

• Stamp
k  ..
.. .. .. .. ..  v 
j
k
:  ..

l  g
m  g

:  ..
Col j
..
g
g
..
..
..
..
..
..
..
..
..
Col k
..
..
..
..
..  : 
   0
..  vl 
.. vm 
 
..  : 
Using the stamps
• b1 is a current source J, others are conductances
g2..g5
b2
b3
2
3
1
b1
b4
b5
0
• All are Type 1 elements
• KCL+KVL+DE gives us
1  g2
 g2
0   v1   J 

   
2  g 2 g 2  g 3  g 4
 g 3  v2    0 
3  0
 g3
g 3  g 5  v3   0 
Another example
1
2
R=1
V=1
R=1
0
• Equations:
1  1  1 1  v1  0

   
2  1 2 0 v2   0
3  1
0 0  iv  1
• Solution (easy to check) v1 = 1, v2 = 0.5, iv = -0.5
Solution of nonlinear equations
• Solution by iteration
I
– Find linear approximation around
a guess solution
– Solve resulting system of linear eqs.
– Iterate until convergence
• Example:
V
– Diode equation: I = f(V) = Io (eV/Vt – 1)
Local linear
approximation
– Assume a kth iteration guess (Ik, Vk)
– Taylor series expansion
I = Ik + [df/dV]k (V-Vk) + higher order terms
(neglected)
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
Solution of a circuit containing diodes
• Example:
conductance ak
+
V
-
conductance g
+
V
-
conductance g
– Linearize the circuit as shown earlier
bk
– 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)
“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)
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
Jk x =Jk xk – F(xk)
where J = Jacobian matrix = analog of the derivative
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
F1 
 F1 F1 F1

 x
x2 x3
xn 
1


 F2 F2 F2  F2 
J   x1 x2 x3
xn 
 



 
 F
Fm Fm
Fm 
m



xn 
 x1 x2 x3
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
1. Linearize F(x) = 0 as Jk x =Jk xk – F(xk) to get n-m linear
eqs
2. Solve these and A x = b – together n linear eqs in n vars.
3. Call this solution xnewk+1; map to xk+1 and increment k
4. If converged, quit; else go to step 1
•
Convergence criterion
||xk+1 – xk|| <  and ||Fk+1 – Fk|| < 
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)
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
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
Solving these equations
• Given a system of linear equations
A x = 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
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
• A x = b becomes L U x = b
• Substitute U x = y where y is an intermediate
variable (vector)
• Solve L y = b (easy) to find y
• Then solve U x = y to find x
Gaussian Elimination (GE)
3 1 1  x1  5

   
1
3
1

  x2   5
1 1 3  x3  5
1   x1   5 
3 1

  

0
8
/
3
2
/
3
x

10
/
3
2

  

0 2 / 3 8 / 3  x3  10 / 3
1   x1   5 
3 1

  

0
8
/
3
2
/
3
x

10
/
3
2

  

0 0 5 / 2  x3   5 / 2 
• 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
GE in general
• Given A x = b, or
• Convert to
U x = y, i.e.,
 a11 a12

a21 a22
 


 an1 an 2
 a1n   x1   b1 
   
 a2n   x2  b2 

       
   
 ann   xn  bn 
u11 u12

 0 u22
 


0
 0
 u1n   x1   b1' 
   
 u2n   x2  b2' 

       
   ' 
 unn   xn  bn 
b1'  b1
where
b2'  b2  m21b1'
b3'  b3  m31b1'  m32b2'
i.e.,

bn'  bn  mn1b1'  mn 2b2'    mn,n 1bn' 1
0
 1

m21 1
 m31 m32


 
m
 n1 mn 2
0
0
1

mn3





0  x1   b1 
   
0  x2  b2 
0  x3   b3 
   
      
1  xn  bn 
(L y = b)
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
• A x = b becomes L U x = b
• Substitute U x = y where y is an intermediate
variable (vector)
• Solve L y = b (easy) to find y
• Then solve U x = y to find x
• Previous slide tells us how to find L, U !!
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
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)
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)
Doolittle
k 1
ukj  akj   lkpu pj ; j  k , k  1,, n
p 1
(k,k)
k 1
aik   lipu pk
lik 
p 1
ukk
; i  k , k  1,, n
(Need to calculate ukj values before lik values – ukk required for lik)
Complexity: same as Gauss (same operations, different order)
Crout: similar update formulae
Pivoting
3 1 1   x1  5

   
3
1
2

  x2   6
1 1 3  x3  5
m21 = 1
1   x1  5
3 1

   
0
0
1

  x2   6
0 2 / 3 8 / 3  x3  5
m31 = 1/3
• 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!)
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
x x x

x x 0
x 0 x

 x 0 0
x

0
0

x 
x

x
x

 x
x

x
x

x 
x
x
x
x
x
x
x
x
When is a fill-in created?
• Computations during LU factorization
(k,k)
• Fill-in created if is a zero and any one
pair of
are both nonzeros
Pivoting for sparsity - example
• Reorder the matrix
x x x

x x 0
x 0 x

 x 0 0
– Row 1  Row 4
– Column 1  Column 4
x x x

x x 0
x 0 x

 x 0 0
x

0
0

x 
becomes
x

0
0

 x
x

0
0

x 
0 0 x

x 0 x
0 x x

x x x 
– LU factors have the structure
x

0
0

 x
0 0 x

x 0 x
0 x x

x x x 
Iterative solution of linear equations
• For a system of equations
• For an
a11x1  a12x2    a1n xn  b1

an1x1  an2 x2    annxn  bn
initial guess x1(0), …, xn(0),
can write
x1(1)  b1   a12x2(0)    a1n xn(0)  / a11



x2(1)  b2   a21x1(0)  a23x3(0)   a2n xn(0)  / a22




x1(1)  b1   a12x2(0)    a1n xn(0)  / a11



x2(1)  b2   a21x1(1)  a23x3(0)   a2n xn(0)  / a22




xn(1)  bn   an1x1(0)    an,n 1xn(0)1  / ann



xn(1)  bn   an1x1(1)    an,n 1xn(1)1  / ann



(Gauss-Jacobi)
(Gauss-Seidel)
In matrix form
• Write A = L + D + U
– L: lower  matrix; U: upper  matrix; D: diagonal matrix
• G-J:
implies
i.e.,
• G-S
i.e.,
(L+D+U) x = b
D x(k+1) = b – (L+U) x(k)
x(k+1) = D-1 b – D-1 (L+U) x(k)
(L+D) x(k+1) = b – U x(k)
x(k+1) = (L+D)-1 b – (L+D)-1 U x(k)
• For both, the update formula has the form
x(k+1) = M x(k) + c
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
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)
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.
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
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
– 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
0
h
2h
3h
t
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
• Trapezoidal rule
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.
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)
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
FE stability (contd.)
• FE formula for test equation: xk = (1+ h )k x0
• Possible values of 
– If  > 0, xk   as k   provided |1+ h | > 1 (always true)
– If  = 0, xk = x0 for all t
– If  < 0, xk  0 as k   provided |1+ h | < 1
• This is true provided –1 < 1+ h  < 1, i.e., if h < 2/||
• In other words, there is an upper bound on the maximum step
size in the physically meaningful region!
– Region of stability
h plane
Radius = 1
(-1,0)
Backward Euler stability
• BE: xn+1 = xn + h dxn+1/dt, i.e., xn+1 = xn + h  xn+1, i.e., xn+1 = xn/(1-h)
– x1 = x0 /(1-h); x2 = x1 /(1-h) = x0/(1-h)2; …xk = x0/(1-h)k
• Stability:
– If  > 0, xk   as k   provided |1- h| < 1
• This is true provided –1 < 1 - h  < 1, i.e., if h < 2/
• Here, the upper bound on the maximum step size in the physically nonmeaningful region!
– If  = 0, xk = x0 for all t
– If  < 0, xk  0 as k   provided |1-h | > 1 (always true)
• Region of stability
h plane
Radius = 1
(1,0)
• Conclusion: BE is a better formula for physical systems
Trapezoidal rule stability
• TR: xn+1 = xn + h/2 (dxn/dt + dxn+1/dt), i.e., xn+1 = xn + h/2 (xn + xn+1)
– i.e., xn+1 = xn (1+h/2) /(1-h/2)
– Can infer that xk = x0 [(1+h/2)/(1-h/2)]k
• Stability:
– If  > 0, xk   as k   provided |(1+h/2)/(1- h/2)| > 1
• Always true!
– If  = 0, xk = x0 for all t
– If  < 0, xk  0 as k   provided |(1+h/2)/(1- h/2)| < 1
• Always true!
h plane
• Region of stability
– Entire h plane!
• WARNING!! Stability guaranteed, not accuracy!
Stability and accuracy
• Stability is a necessary but not sufficient condition
• Need accuracy too!
• If h is chosen so that h = -4 then
– xk = x0 [(1+h/2)/(1-h/2)]k = x0 (-1/3)k
– Oscillates between position and negative values – not true
for the actual solution!
• Similarly h  = 2 means that the solution shoots to
infinity at the first time step
Linear multistep formulae
• Try to fit a kth order polynomial
– FE/BE fit a linear function and are order 1 LMS formulae
– For dx/dt = f(x), require k values of function f or derivative
df/dx at previous time points
• Explicit formulae (e.g., FE)
– Evaluation of xn+1 does not depend on dxn+1/dt
• Implicit formulae (e.g., BE, TR)
– Evaluation of xn+1 depends on dxn+1/dt
• No details in this class – for further reference, see
Chua and Lin’s book
Gear’s formulae
• LMS formulae that belong to the family
xn+1 = (j=1 to k j xn-j ) + h  dxn/dt
(the above is a kth order Gear’s formula)
• Specific formulae
– First order (=BE):
xn+1 = xn + h dxn+1/dt
– Second order:
xn+1 = 4/3 xn – 1/3 xn-1 + 2/3 h dxn+1/dt
– Third order:
xn+1 = 18/11 xn – 9/11 xn-1 + 2/11 xn-2 + 6/11 h dxn+1/dt
• Good stability properties (“stiffly stable”) and used for circuit
simulation applications
Application to circuit elements
• Capacitor
– I = C dV/dt
– Work instead with charge: q = C V
(better formulation for nonlinear caps)
– Objective: to find a “companion model” converting a differential
equation to a linear equation (perhaps going through a nonlinear
equation on the way)
• Forward Euler
–
–
–
–
qn+1 = qn + h dqn/dt  (CV)n+1 = (CV)n + h In
Linear capacitor: C is constant, so that Vn+1 = Vn + h/C In
Since Vn+1 is the unknown, can rewrite as V = Vn + h/C In
Capacitor is replaced by a constant voltage source at this time step!
Problem with forward Euler
• Loop of capacitors becomes a loop of V sources
+
V1
-
+ V2 -
+
V3
-
• No guarantee that KVL will hold around the loop
– (There are some ways of getting around this)
• This, plus stability problems, imply that FE is not used in circuit
simulation
• Also: possible numerical problems if C is very small
• In general: avoid the use of any explicit formula
Backward Euler – Companion models
• Linear capacitor
– qn+1 = qn + h dqn+1/dt  (CV)n+1 = (CV)n + h In+1
– Linear capacitor: C is constant, so that
(C/h)Vn+1 = (C/h) Vn + In+1
– Since Vn+1, In+1 are unknowns, can rewrite as
I = (C/h) V – (C/h) Vn
– Capacitor is replaced by a type 1 element: constant current
source in parallel with a conductance - at this time step!
• Linear inductor: work with flux  = LI; V = d/dt
– n+1 = n + h dn+1/dt  (LI)n+1= (LI)n + h Vn+1
 V = L/h I – L/h In
– Resistance in series with a voltage source
Resistance L/h
- L/h In +
Choice of timestep and LTE
• LTE = Local truncation error
• FE/BE can be thought of as truncated Taylor series
– xn+1 = xn + h dxn/dt + h2/2 d2xn/dt2 + higher order terms
– Truncation error = h2/2! d2xn/dt2 + higher order terms, which is
approximated as h2/2! d2xn/dt2
• Choose a time step so that
LTE = h2/2 d2xn/dt2 < 
i.e., h < sqrt(2/xn”)
xn” is estimated by divided differences
• For BE,
LTE = h2/2 d2xn+1/dt2 <   h < sqrt(2/xn+1”)
(Can use FE to estimate xn+1” and then apply divided differences)
Sensitivity calculation
The direct method
• Finds the sensitivity of all responses to a single parameter
change
• System equation (e.g., MNA): A x = b
• Perturbation: A  A + A ; b  b + b
– x  x + x
• After perturbation
– (A + A)(x + x) = (b + b)
– Substituting A x = b and neglecting products of  terms, we get
A x = b - A x
– Solve this system of equations to get the sensitivity change of all
parameters x wrt a single parameter perturbation
• Same LHS as the original system; LU factors available
Adjoint sensitivities: Tellegen’s theorem
• Finds the sensitivity of a single response to all parameter changes
• For any pair of circuits with the same topology (i.e., incidence matrix)
(branch_voltage_vector)^T (branch_current_vector) = 0
• Consider three circuits with the same topology
Circuit 1
The original circuit
Voltage = V
Current = I
Circuit 2
The perturbed circuit
Voltage = V + V
Current = I + I
Circuit 3
The adjoint circuit
(to be defined)
Voltage = 
Current = 
– Circuits 1 and 3: VT = 0; TI = 0
•
[Also,  T  = VT I = 0 – conservation of power]
– Circuits 2 and 3: T (V + V) = 0; T(I + I) = 0
– Combining the two, TV = 0 and TI = 0
– Subtracting the second equation from the first,
(TV - TI) = 0, i.e., all branches (j Vj - j Ij ) = 0
Perturbations in circuits
• Consider circuits with R, V, J only
• Rewrite all branches (j Vj - j Ij ) = 0 as
R(j Vj - j Ij) + V(j Vj - j Ij) + J (j Vj - j Ij) = 0
• Perturbations:
– Increase the value of a voltage source to vj  vj + vj
Vj = vj ; Ij is determined by the rest of the circuit
– Increase the value of a current source Jj  Jj + Jj
Ij = Jj ; Vj is determined by the rest of the circuit
– Increase the value of a resistor Rj  Rj + Rj
Vj = IjR
(Vj+Vj) = (Ij+Ij)(R Rj)
Therefore (neglecting products of  terms)
Vj = Ij Rj + Rj Ij
Perturbations (contd.)
• Therefore, for a resistor,
R(j Vj - j Ij ) = R(j Ij Rj + j Rj Ij - j Ij )
• At this point, we have no specs on the adjoint circuit
but its topology: we are free to choose the  - 
relations as we want
– Choose j = j Rj
– Why?
• This cancels out the last two terms with circuit-dependent Ij
terms
• R(j Vj - j Ij ) = R(j Ij Rj + j Rj Ij - j Ij ) = Rj Ij Rj
• This term now depends only on the perturbation, not on the
circuit
Adjoint sensitivity computation
R(j Vj - j Ij) + V(j Vj - j Ij) + J (j Vj - j Ij) = 0
Rewrite as
V j Ij - J j Vj = R j Ij Rj + V j Vj + J j Ij
• LHS contains  terms corresponding to circuit-determined
variations
• RHS contains  terms corresponding to parameter variations
– To find sensitivities of voltages with respect to parameters
• Set j = 0; x = -1 for the voltage of interest, 0 for all others
• Vx = R j Ij Rj + V j Vj + J j Ij
• Vx/Rj = j Ij ; Vx/Vj = j ; Vx/Ij = j
Adjoint sensitivities
• To find sensitivities of currents with respect to multiple
parameters
– Set x = 1 for the current of interest, 0 for all others; j = 0
– Ix = R j Ij Rj + V j Vj + J j Ij
– Ix/Rj = j Ij ; Ix/Vj = j ; Ix/Ij = j
• In general: to find sensitivities of f(Vx,Ix) with respect to multiple
parameters
– f = f/Ix Ix + f/Vx Vx
– Compare with LHS of
V j Ij - J j Vj = R j Ij Rj + V j Vj + J j Ij
– Set
• x = f/Ix ; j = 0 for all others,
• x = - f/Vx j = 0 for all others
– f/Rj = j Ij ; f/Vj = j ; f/Ij = j
Example: Adjoint sensitivity computation
• Response of interest: voltage Vx
Original circuit
Modified circuit
V1
Vx
V1
R3
V=1
Vx
R3
R1
R2
gnd
V=1
R1
I=0
R2
gnd
Add a fake current source (open ckt) since the response
must be V across a current source (or I through a V source)
Adjoint circuit for Vx sensitivity
 = 0 for all voltage sources
x = -1 for the current source related to Vx ; 0 else
Solve adjoint circuit for all branch , 
 =0
Vx = R j Ij Rj + V j Vj + J Ij
Vx/Rj = j Ij ; Vx/Vj = j ; Vx/Ij = j
1
x
R3
R1
R2
 = -1
gnd
[Pillage, Rohrer and Visweswariah}
Transient adjoint sensitivities
• So far, all discussions on linear/linearized nonlinear
elements
• Handling circuits with capacitors (for example)
• Starting point: integrate all branches (j Vj - j Ij ) = 0
Numerous minor notational errors in this slide – struggles with MS Equation Editor.
Please see the handout on the web page for correct info.
Transient sensitivity (contd.)
T
  ( ( ) (t)  ψ ( ) (t)) dt
j
j
j
j
all capacitors 0



I j (t)  C Vc  δI j (t)  Vj δC  Cδ Vj
Each summation term becomes
T


 ( ( ) (t)  ψ ( ) V (t)δt
j
j
j
 ψ j ( )Cδ Vj (t)) dt
j
0
b
b
b
du
Integrate by parts   uvdt  u  vdt - 
dt
a
a
a
T

 ψ ( )Cδ V dt ψ ( )CδC
j
j
j
T
j 0
0
 vdt dt 

T 
  ψ j ( )CδCjdt
0
The entire integral becomes
T


- ψ j ( )CδCj   ( j ( ) j (t)  ψ j ( ) Vj δC  ψ j ( )CδCj ) dt
T
0
0
Numerous minor notational errors in this slide – struggles with MS Equation Editor.
Please see the handout on the web page for correct info.
Transient sensitivity (contd.)

(Temporari ly) select  j ( )  - Cψ j ( ), ψ j (  0)  0
[This is a " negative" capacitor in the adjoint circuit, with initial voltage 0]
The integral becomes
T

  ψ ( ) V (t)δt
j
j
dt
0
Negative capacitor? ? Now set   T  t (i.e., d  -dt)

Replace the temporary selection by j ( )  Cψ j ( ), ψ j (  0)  0 where  is as defined above
Choose current, voltage sources in the adjoint circuit to select the I, V, or f of interest

f
0 C dt  0  ψ j ( ) Vj (t) dt
T
T
To get LHS 
Vx (t c )
, select f  Vx (t) (t - t c )
C