## Derivation of System Equations

The electrical network of example 2 is shown above. The elaborated network consists of $$l = 34$$ branches and $$k = 9$$ nodes.

The proper tree corresponding to the electrical network is shown above. Since capacitors $$C_{TG1}$$, $$C_{OUTN4}$$, $$C_{INP4}$$, $$C_{OUTN1}$$, $$C_{INP3}$$ and $$C_{INP2}$$ are members of the proper tree, the state vector may be written as:

$\mathbf{x} = \begin{bmatrix} v_{C_{TG1}} \\ v_{C_{OUTN4}} \\ v_{C_{INP4}} \\ v_{C_{OUTN1}} \\ v_{C_{INP3}} \\ v_{C_{INP2}} \\ \end{bmatrix}$

Determining the proper tree becomes increasingly complex as the number of branches and nodes are increased. As the network is comprised of 9 nodes, there are 8 basic cuts. To ease the process of determining all 8 basic cuts, it is often more convenient to form the $$\mathbf{Q_f}$$ matrix (where rows 0 to 7 correspond to the basic cuts containing $$v_{in}$$, $$v_{S}$$, $$C_{TG1}$$, $$C_{OUTN4}$$, $$C_{INP4}$$, $$C_{OUTN1}$$, $$C_{INP3}$$ and $$C_{INP2}$$ respectively:)

$\mathbf{Q_f} = \begin{bmatrix} -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & \\ 0 & 1 & -1 & 0 & 0 & -1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & 1 & 0 & -1 & -1 & -1 & 1 & 1 & 0 & -1 & -1 & -1 & 0 & -1 & 0 & \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & 1 & 0 & -1 & -1 & -1 & 0 & 0 & -1 & 0 & 0 & -1 & 1 & -1 & 0 & \\ 0 & 0 & 0 & 1 & 1 & 1 & -1 & -1 & 0 & 1 & 1 & 0 & -1 & -1 & -1 & -1 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & -1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & -1 & -1 & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ \end{bmatrix} \begin{bmatrix} i_{IN} \\ i_{TG1} \\ i_{C_{TG1}} \\ i_{P1} \\ i_{C_{OUTP1}} \\ i_{C_{INP1}} \\ i_{N1} \\ i_{C_{OUTN1}} \\ i_{C_{INN1}} \\ i_{P2} \\ i_{C_{OUTP2}} \\ i_{C_{INP2}} \\ i_{N2} \\ i_{C_{OUTN2}} \\ i_{C_{INN2}} \\ i_{TG2} \\ i_{C_{TG2}} \\ i_{TG3} \\ i_{C_{TG3}} \\ i_{P3} \\ i_{C_{OUTP3}} \\ i_{C_{INP3}} \\ i_{N3} \\ i_{C_{OUTN3}} \\ i_{C_{INN3}} \\ i_{P4} \\ i_{C_{OUTP4}} \\ i_{C_{INP4}} \\ i_{N4} \\ i_{C_{OUTN4}} \\ i_{C_{INN4}} \\ i_{TG4} \\ i_{C_{TG4}} \\ i_S \end{bmatrix} = \mathbf{0} \label{eq:ex2_qf}$

The differential equation for $$v_{C_{TG1}}$$ is obtained from the basic cut involving capacitor $$C_{TG1}$$. From row 2 of equation \eqref{eq:ex2_qf}:

$$i_{TG1} - i_{C_{TG1}} - i_{C_{INP1}} - i_{C_{INN1}} + i_{TG2} - i_{C_{TG2}} = 0 \\ \label{eq:ex2_kcl1}$$

Equation \eqref{eq:ex2_kcl1} may be rewritten in terms of the system state variables by noting that $$v_{C_{INP1}} = v_{C_{TG1}} - v_{C_{OUTN1}}$$, $$v_{C_{INN1}} = v_{C_{TG1}}$$, $$v_{C_{TG2}} = v_{C_{TG1}}$$ and $$v_{TG2} = v_{C_{OUTN1}} - v_{C_{INP2}} - v_{C_{TG1}}$$:

$$i_{TG1} - C_{TG1}\frac{dv_{C_{TG1}}}{dt} - C_{INP1}\left(\frac{dv_{C_{TG1}}}{dt} - \frac{dv_{C_{OUTN1}}}{dt}\right) - C_{INN1}\frac{dv_{C_{TG1}}}{dt} + i_{TG2} - C_{TG2}\left(\frac{dv_{C_{TG1}}}{dt}\right) = 0 \\ \label{eq:ex2_kcl2}$$

Note that the current sources $$i_{TG1}$$ and $$i_{TG2}$$ are functions of their respective branch voltages, which may be expressed in terms of the state vector. If a similar approach is adopted for rows 3 - 7 of equation \eqref{eq:ex2_qf}, the following system of equations, expressed in matrix form, may be derived:

$\begin{bmatrix} C_{TG1}+C_{INP1}+C_{INN1}+C_{TG2} & 0 & 0 & -C_{INP1} & 0 & 0 \\ 0 & C_{OUTN4}+C_{TG3}+C_{OUTP3}+C_{OUTN3}+C_{INN3}+C_{OUTP4}+C_{INN4}+C_{TG4} & C_{TG3}+C_{OUTP3}+C_{OUTN3}+C_{INN3}+C_{INN4}+C_{TG4} & 0 & C_{TG3}+C_{INN3}+C_{TG4} & 0 \\ 0 & C_{TG3}+C_{OUTP3}+C_{OUTN3}+C_{INN3}+C_{INN4}+C_{TG4} & C_{INP4}+C_{TG3}+C_{OUTP3}+C_{OUTN3}+C_{INN3}+C_{INN4}+C_{TG4} & 0 & C_{TG3}+C_{INN3}+C_{TG4} & 0 \\ -C_{INP1} & 0 & 0 & C_{OUTN1}+C_{OUTP1}+C_{INP1}+C_{OUTP2}+C_{OUTN2}+C_{INN2} & 0 & -C_{OUTP2}-C_{OUTN2} \\ 0 & C_{TG3}+C_{INN3}+C_{TG4} & C_{TG3}+C_{INN3}+C_{TG4} & 0 & C_{INP3}+C_{TG3}+C_{INN3}+C_{TG4} & 0 \\ 0 & 0 & 0 & -C_{OUTP2}-C_{OUTN2} & 0 & C_{INP2}+C_{OUTP2}+C_{OUTN2} \\ \end{bmatrix} \begin{bmatrix} \frac{dv_{C_{TG1}}}{dt} \\ \frac{dv_{C_{OUTN4}}}{dt} \\ \frac{dv_{C_{INP4}}}{dt} \\ \frac{dv_{C_{OUTN1}}}{dt} \\ \frac{dv_{C_{INP3}}}{dt} \\ \frac{dv_{C_{INP2}}}{dt} \\ \end{bmatrix} = \begin{bmatrix} i_{TG1} + i_{TG2} \\ i_{TG3} + i_{P3} + C_{OUTP3}\frac{dv_S}{dt} - i_{N3} + i_{P4} + C_{OUTP4}\frac{dv_S}{dt} - i_{N4} \\ i_{TG3} + i_{P3} + C_{OUTP3}\frac{dv_S}{dt} - i_{N3} + i_{TG4} \\ i_{P1} + C_{OUTP1}\frac{dv_S}{dt} - i_{N1} + i_{P2} + C_{OUTP2}\frac{dv_S}{dt} - i_{N2} - i_{TG2} - i_{TG3} \\ i_{TG3} + i_{TG4} \\ i_{TG2} - i_{P2} - C_{OUTP2}\frac{dv_S}{dt} + i_{N2} \\ \end{bmatrix} \label{eq:ex2_eqns}$

Notice that the right hand side of equation \eqref{eq:ex2_eqns} contains references to the time derivative of the supply voltage $$\left(\frac{dv_S}{dt}\right)$$. These terms arise whenever a co-tree branch forms a loop with the $$v_S$$ branch. For example, $$v_{C_{OUTP1}} = v_S - v_{C_{OUTN1}}$$ and therefore the time derivative may be written as:

$$\frac{dv_{C_{OUTP1}}}{dt} = \frac{dv_S}{dt} - \frac{dv_{C_{OUTN1}}}{dt} \\$$

Time derivatives of system inputs can often present problems to the numerical solution of the system. The inputs can often change rapidly with respect to time, resulting in derivatives large in magnitude and causing the equations to become stiff. In such situations, small time steps and stringent error tolerances are often required to prevent the solution from becoming unstable. In this example, since the derivative is of the supply voltage, assumed to be constant for all time, the derivative evaluates to zero.

## Solving the System of Equations

Octave's function lsode can be used to solve Ordinary Differential Equations (ODEs) of the form:

$$\frac {d\mathbf{x}}{dt} = f(\mathbf{x}, t) \\ \label{eq:oct1}$$

The ODE solver requires a function to compute the vector of derivatives of equation \eqref{eq:oct1}. For this example, the derivatives are evaluated by the function dypc_example2.

	  function dy = dypc_example2(y,t)

V_CTG1 = y(1);
V_COUTN4 = y(2);
V_CINP4 = y(3);
V_COUTN1 = y(4);
V_CINP3 = y(5);
V_CINP2 = y(6);

dy = zeros(6,1);

% supply voltage and time derivative
V_VS = 1;
dVS = 0;

% setup din input voltage and clock sources.
% din is a 50/50 duty cycle pulse with a period of 25ns.
% clock is a 50/50 duty cycle pulse with a period of 10ns.
V_VDIN = pulse(25e-9, 12.5e-9, 1e-12, 1e-12, t);
V_VCIN = pulse(10e-9, 5e-9, 1e-12, 1e-12, t+5e-9);
V_VCINB = pulse(10e-9, 5e-9, 1e-12, 1e-12, t);

% instantiate the CMOS transistor models and calculate the input/output capacitance
% and drain current.

% CMOS model parameters
Ln = 10;
Wn = 10;
Lp = 10;
Wp = 20;
scale = 50e-9;
Cox_dash = 25e-15;

%N1
Vgs = y(1); % in parallel with CTG1
Vds = y(4);
F = nmos_digital(Vgs, Vds, Ln, Wn, scale, Cox_dash);
CINN1 = F(1);
COUTN1 = F(2);
i_IN1 = F(3);

%P1
V_COUTP1 = -V_COUTN1 + V_VS;
V_CINP1 = V_CTG1 - V_COUTN1;
Vsg = V_CINP1 - V_COUTP1;
Vsd = -V_COUTN1 + V_VS;
F = pmos_digital(Vsg, Vsd, Lp, Wp, scale, Cox_dash);
CINP1 = F(1);
COUTP1 = F(2);
i_IP1 = F(3);

%N2
V_CINN2 = V_COUTN1;
Vgs = V_CINN2;
Vds = V_COUTN1 - V_CINP2;
F = nmos_digital(Vgs, Vds, Ln, Wn, scale, Cox_dash);
CINN2 = F(1);
COUTN2 = F(2);
i_IN2 = F(3);

%P2
V_COUTP2 = -V_COUTN1 + V_CINP2 + V_VS;
Vsg = V_CINP2 - V_COUTP2;
Vsd = -V_COUTN1 + V_CINP2 + V_VS;
F = pmos_digital(Vsg, Vsd, Lp, Wp, scale, Cox_dash);
CINP2 = F(1);
COUTP2 = F(2);
i_IP2 = F(3);

%N3
V_CINN3 = V_CINP3 + V_CINP4 + V_COUTN4;
Vgs = V_CINN3;
Vds = V_CINP4 + V_COUTN4;
F = nmos_digital(Vgs, Vds, Ln, Wn, scale, Cox_dash);
CINN3 = F(1);
COUTN3 = F(2);
i_IN3 = F(3);

%P3
V_COUTP3 = -V_CINP4 - V_COUTN4 + V_VS;
Vsg = V_CINP3 - V_COUTP3;
Vsd = -V_CINP4 - V_COUTN4 + V_VS;
F = pmos_digital(Vsg, Vsd, Lp, Wp, scale, Cox_dash);
CINP3 = F(1);
COUTP3 = F(2);
i_IP3 = F(3);

%N4
V_CINN4 = V_CINP4 + V_COUTN4;
Vgs = V_CINN4;
Vds = y(2);
F = nmos_digital(Vgs, Vds, Ln, Wn, scale, Cox_dash);
CINN4 = F(1);
COUTN4 = F(2);
i_IN4 = F(3);

%P4
V_COUTP4 = -V_COUTN4 + V_VS;
Vsg = V_CINP4 - V_COUTP4;
Vsd = -V_COUTN4 + V_VS;
F = pmos_digital(Vsg, Vsd, Lp, Wp, scale, Cox_dash);
CINP4 = F(1);
COUTP4 = F(2);
i_IP4 = F(3);

%TG1
V_ITG1 = V_VDIN -V_CTG1;
F = tg(V_ITG1, V_VCIN, V_VCINB, Ln, Wn, Lp, Wp, scale, Cox_dash);
CTG1 = F(1);
i_ITG1 = F(2);

%TG2
V_ITG2 = -V_CTG1 + V_COUTN1 - V_CINP2;
F = tg(V_ITG2, V_VCINB, V_VCIN, Ln, Wn, Lp, Wp, scale, Cox_dash);
CTG2 = F(1);
i_ITG2 = F(2);

%TG3
V_ITG3 = V_COUTN1 - V_CINP3 - V_CINP4 - V_COUTN4;
F = tg(V_ITG3, V_VCINB, V_VCIN, Ln, Wn, Lp, Wp, scale, Cox_dash);
CTG3 = F(1);
i_ITG3 = F(2);

%TG4
V_ITG4 = -V_CINP3 - V_CINP4;
F = tg(V_ITG4, V_VCIN, V_VCINB, Ln, Wn, Lp, Wp, scale, Cox_dash);
CTG4 = F(1);
i_ITG4 = F(2);

% create the Qm and Tm matrices.
Qm = [(CTG1+CINP1+CINN1+CTG2) 0 0 -CINP1 0 0;
0 (COUTN4+CTG3+COUTP3+COUTN3+CINN3+COUTP4+CINN4+CTG4) (CTG3+COUTP3+COUTN3+CINN3+CINN4+CTG4) 0 (CTG3+CINN3+CTG4) 0;
0 (CTG3+COUTP3+COUTN3+CINN3+CINN4+CTG4) (CINP4+CTG3+COUTP3+COUTN3+CINN3+CINN4+CTG4) 0 (CTG3+CINN3+CTG4) 0;
-CINP1 0 0 (COUTN1+COUTP1+CINP1+COUTP2+COUTN2+CINN2) 0 (-COUTP2-COUTN2);
0 (CTG3+CINN3+CTG4) (CTG3+CINN3+CTG4) 0 (CINP3+CTG3+CINN3+CTG4) 0;
0 0 0 (-COUTP2-COUTN2) 0 (CINP2+COUTP2+COUTN2)];

Tm = [i_ITG1 + i_ITG2;
i_ITG3 + i_IP3 + COUTP3*dVS - i_IN3 + i_IP4 + COUTP4*dVS - i_IN4;
i_ITG3 + i_IP3 + COUTP3*dVS - i_IN3 + i_ITG4;
i_IP1 + COUTP1*dVS - i_IN1 + i_IP2 + COUTP2*dVS - i_IN2 - i_ITG2 - i_ITG3;
i_ITG3 + i_ITG4;
-i_IP2 - COUTP2*dVS + i_IN2 + i_ITG2];

% calculate the system derivatives by performing Gaussian elimination.
dy = Qm\Tm;


The function lsode uses a numerical integration algorithm to solve the system of equations. It must be linked to the dypc_example2 function and invoked with suitable arguments. For this example, the following code sequence is used to solve the system of equations for time $$t=0$$ to $$t=100 ns$$.

	  clear all;
abTol = 1e-3;

t = linspace (0, 100e-9, 30e3);
x0 = [0 0 0 0 0 0];

lsode_options("absolute tolerance",[abTol abTol abTol abTol abTol abTol]);
lsode_options("relative tolerance", 1e-4);
lsode_options("integration method", "stiff");
lsode_options("maximum step size", 0.5e-9);

tic
Y = lsode("dypc_example2", x0, t);
toc


Notice that in the above code fragment the "maximum step size" parameter is specified. This prevents the solver from using too large a time step such that important changes on the system inputs are "missed".