Harmonic Analysis and Design
Of Nonlinear Electronic Circuits
By
DAVID ANTHONY WAYNE
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLHl T
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1971 .
To
Sharon, my wife
and
Trudy and Ted, my children:
I lock forward to getting to know them again.
ACKNOWLEDGEMENTS
I wish to express sincere appreciation to
Dr. S.W. Director, Chairman of my Supervisory Committee,
for his guidance and enthusiasm during this research.
Thanks are due also to the members of the Supervisory
Committee, Dr. F.A. Lindholm, Dr. E.R. Chenette,
Dr. A.J. Brodersen, and Dr. A.W. Westerberg, for reviewing
and commenting on this work.
Finally, I wish to express loving thanks to my wife,
Sharon, for her patience and understanding over the past
three years and for her excellent typing of this manuscript.
Financial support from an NDEA Title IV Fellowship
and an NSF Traineeship is gratefully acknowledged.
iii
TABLE OF CONTENTS
Page
ACKNOWLEDGEMENTS . . . . . . . . . iii
LIST OF TABLES . . . . . . . . vi
LIST OF FIGURES . . . . . . .... . . vii
ABSTRACT. . . . . . . . . .. .. x
CHAPTER 1
INTRODUCTION. .... . . . . . . . . 1
1.1 Background . . . . . . . . 1
1.2 TimeDomain Design . . . . . . 4
1.3 References . . . . . . ... . 7
CHAPTER 2
CALCULATION OF STEADYSTATE RESPONSE . . ... 12
2.1 Background . . . . . . . . 12
2.2 The State VariableNewton Method
of SteadyState Calculation. . . . .. 34
2.3 The Adjoint NetworkNewton
Method of SteadyState Calculation .. .. 37
2.4 References . . . . . . . .. 48
CHAPTER 3
FREQUENCYDOMAIN ASPECTS OF NONLINEAR
CIRCUIT DESIGN. .. . . . . . . . .53
3.1 The Design Algorithm . . . . . .. 55
3.2 Design Examples. . . . . . . .. 64
3.3 Discussion . . . . . . . . 78
3.4 References . . . . ... . . . 80
CHAPTER 4
CONSIDERATIONS FOR OSCTILTLTOR DESIGN .. 82
4.1 Determination of Frequency of
Oscillation. . . . .... .. . . .. 82
4.2 Calculation of Gradient Components
for Automated Oscillator Design ..... 88
4.3 References . . . . . . . . 105
CHAPTER 5
SUMMARY AND SUGGESTIONS FOR FUTURE RESEARCH . 106
Page
APPENDIX A
DERIVATION OF GRADIENT COMPONENTS . . .. 109
References . . . . . . . . . 122
APPENDIX B
FLOW CHART OF DESIGN ALGORITHM . . . . 123
BIOGRAPHICAL SKETCH. . . . . . . . 131
LIST OF T;.ELE3
Table P age
21 Results of Application of Galerkin's
Method to Circuit of Fig. 21. . . .. 18
22 Results of Application of Adjoint
NetworkOptimization Method to
Circuit of Fig. 21. . . . ... . 32
23 Results of Application of Adjoint
NetworkOptimization Method to
Circuit of Fig. 26 . . . . .. 35
24 Results of Application of Steady
State Determination Methods to
Circuit of Fig. 26. . . . . . ... 50
25 Results of Application of Steady
State Determination Methods "o
Circuit of Fig. 27 .. . . ..... 52
31 Transistor Parameters ........... 65
32 Results of Application of Design
Algorithm to Circuit of Fig. 33 . . .. 73
33 Results of Application of Design
Algorithm to Circuit of Fig. 37 . . .. 79
LIST OF FIGURES
Figure
2 1 Circuit used in examples of steady
state determination methods. . . .
2 2 Flow of operations for implementation
of Galerkin's method ... . . . .
2 3 Capacitor and inductor initial
condition models
(a) Linear timeinvariant. . . .
(b) Nonlinear. . . . . . .
2 4 The original and adjoint networks used
ir gradient calculation. . . . .
Page
15
. . 17
. . 20
. . 20
. . 23
2 5 Flow of operations for adjoint network
optimization method of steadystate
calculation. . . . . . . . ... .31
2 6 Voltage doubler circuit of Lhe
second example . . . . . . . . 33
2 7 Flow of operations for adjoint network
Newton method of steadystate calculation; . 47
2 8 Circuit for adjoint networkNewton method
example. . . . . . . . . 51
3 1 Flow of operations for automated
design algorithm . . . . . . .
3 2 Arrangement of adjoint network to
obtain gradient components . . . .
3 3 The shuntseries feedback pair
amplifier of the first example . . .
3 4 Initial configuration for example
of Fig. 33
(a) Output waveform. . . . . .
(b) Frequencydomain components
of output signal . . . . .
(c) Adjoint network excitation. . ..
. 60
S. 62
. 66
. 67
. 67
. 67
Figure
3 5 Circuit configuration for example
of Fig. 33 after one design
iteration
(a) Output waveform. . . . . ... 69
(b) Frequencydomain components
of output signal . . . ... 69
(c) Adjoint network excitation .. ... 69
3 6 Final circuit configuration for
example of Fig. 33
(a) Output waveform. . . . . . 71
(b) Frequencydomain components
of output signal . . . . .. 71
(c) Adjoint network excitation . .. 71
3 7 The largesignal audio power amplifier
of the second example . . . . ... 74
3 8 Initial circuit configuration for
example of Fig. 37
(a) Output waveform. . . . . ... 75
(b) Frequencydomain components
of output signal . . . ... 75
(c) Adjoint network excitation ..... 75
3 9 Circuit configuration for example
of Fig. 37 after one design iteration
(a) Output waveform. . . . . ... 76
(b) Frequencydomain components
of output signal . . . ... 76
(c) Adjoint network excitation . .. 76
310 Final circuit configuration for
example of Fig. 37
(a) Output waveform. . . . . ... 77
(b) Frequencydomain components
of output signal . . . ... 77
(c) Adjoint network excitation . . 77
4 1 Oscillator circuit
(a) Original network . . . . .. 83
(b) Adjoint network. . . . . ... 83
4 2 Output waveforms of network of Fig. 41
(a) Original network output voltage v 84
(b) Adjoint network output voltage 0 84
4 3 Phase plane diagram for an oscillator . 86
viii
Page
Figure
4 4 (a) Negativeresistance oscillator.
(b) Characteri'stics of negative
resistance device . . . .
4 5 Original and adjoint networks used
to calculate required partial
derivatives. . . . . . . .
4 6 RC feedback oscillator network . .
4 7 An oscillator circuit which may not be
partitioned. . . . . . . .
A 1 Networks used to obtain the summation
terms in the gradient component
expressions. . .' . . . . .
B 1 Flow chart of main program . . .
B 2 Flow chart of SEARCH routine .. ....
B 3 Flow chart of FUNCT routine . . .
B 4 Flow chart of DCANAL routine . . .
B 5 Flow chart of FOURR routine. . . .
Page
. . 90
. . 90
. . 96
. . 103
. . 104
. . 119
. . 124
. . 125
. . 126
* . 129
* . 130
Abstract of Dissertation Presented to the Graduate Council
of the University of Florida in Partial Fulfillment
of the Requirements for the Degree of Doctor of Philosophy
HARMONIC ANALYSIS AND DESIGN
OF NONLINEAR ELECTRONIC CIRCUITS
By
David Anthony Wayne
December, 1971
Chairman: Dr. Stephen W. Director
Major Department: Electrical Engineering
Methods are developed for the automated analysis and
design of nonlinear networks whose performance is charac
terized in frequencydomain terms. Such circuits are
generally found in communications systems and include large
signal amplifiers, mixers, frequency multipliers, oscillators,
and others. These design procedures utilize the efficient
"adjoint network" method, sparse matrix techniques and the
fast Fourier transform.
The problem of determination of steadystate response of
a nonlinear network to a periodic input is also treated. Two
methods are developed which calculate this response efficiently.
Special considerations for oscillator design are
discussed. The adjoint network corresponding to an
oscillator network may exhibit an unstable response: one
which grows without bound. Methods for dealing with this
difficulty are presented. Techniques are also developed
for calculating frequency of oscillation.
Several computer programs have been written to utilize
the algorithms presented. A number of problems using these
programs are presented to demonstrate the effectiveness of
the procedures.
CHAPTER 1
INTRODUCTION
The purpose of this work is the development of methods
for the automated design of nonlinear electronic networks
whose response is characterized in frequencydomain terms.
These circuits are generally found in communications systems
and include largesignal amplifiers, mixers, frequency
multipliers and dividers, oscillators and others.
1.1 Background
The analysis and design of nonlinear networks by
conventional analytic or graphical means are difficult and
tedious. In addition, these methods are generally applicable
only in situations in which one or two nonlinearities dominate
the behavior of the circuit. The resulting solution may only
be a firstorder approximation and detailed information about
the circuit found only through the construction of prototypes,
an expensive luxury in_';:t j sted circuit production.
Highspeed digital computers have had a profound impact
on the methods of circuit analysis and design. Computer
methods make practical the numerical simulation of large
nonlinear circuits without unrealistically limiting the
complexity of the nonlinearities. Development of automated
2
circuit analysis. algorithms began in the early 1960's. Programs
such as ECAP [1] and SCEPTRE [2] allowed the circuit design
engineer to accurately analyze a network. The designer then
modified the circuit, using his experience and intuition.
Again the computer was used to perform the analysis task. This
procedure was repeated iteratively until the design goals were
met (or until the designer gave up).
Increasing circuit complexity has been the incentive for
the development of fully automated design methods. In these
procedures, the computer is given an active role in the design
cycle. Moreover, in addition to analyzing the network, the
computer determines how the design parameters must be adjusted
in order to achieve the design goals. Techniques that may lead
to fully automated design have evolved through the efforts of
many researchers [3,4,5]. Hachtel and Rohrer [6,7] first
presented in 1967 a highly efficient scheme for automated
design which used parameter optimization. In their approach,
variational calculus methods were applied to the statevariable
formulation of the network equations to obtain the required
parameterspace error gradient. Director and Rohrer [8]
developed even more practical algorithms based on the "adjoint
network" concept. The adjoint network technique has been
successfully applied to dc [91 and smallsignal ac [10] net
work design. Simultaneous design for dc and smallsignal dc
performance has also been investigated [11,12]. Timedomain
design of switching circuits has been described by several
References will appear at the end of each chapter.
authors [13,14,15]. Timedomain automated design provides
the basis for the present work.
For a timedomain design scheme to be practical, the
required numerical integration must be done efficiently.
Computational effort required for numerical integration is
roughly proportional to the number of increments, or "steps,"
into which the time interval is divided. The "order" of a
numerical integration method may be defined as the highest
order of a polynomial in the time variable to which the method
would give an exact solution. Higher order methods may be
seen to provide a better approximation where the responses
are rapidly varying. Gear [16,17] has proposed a scheme in
which order of integration and time step are both adjusted to
minimize computational effort. Calahan [18] has developed
techniques for the application of Gear's algorithm to time
domain network analysis.
The system of equations which describes a network is
frequently "sparse"; that is, the coefficient matrix contains
relatively few nonzero entries. Exploitation of this sparse
ness has two advantages: First, computer memory requirements
are reduced by saving only the nonzero coefficients. Second,
only the nonzero entries need to be processed during the
solutionresulting in a time savings. Techniques which take
advantage of sparseness: innetwork equations have been
developed by Gustavson et al. [19], Tinney and Walker [20],
and others. A novel design technique based on sparse matrix
manipulation methodscalled the "tableau" approachhas been
presented by Hachtel et al. [21,22].
1.2 TimeDomain Design
Timedomain automated design begins with specification of
desired time responses v(t) and i(t), as well as an initial
network configuration. A performance function which reflects
the design specifications is then formulated. A suitable
performance function is the leastsquares function:
tf 1 2
(x) = W/ {v(t)[v(t,x) v(t)
o0 (1.1)
+ wi(t)[i(t,x) i(t)]2}dt,
where x represents the vector of designable parameters and
w and wi are nonnegative weighting functions which allow for
design flexibility. At its minimum, E(x) yields the optimum
set of network parameters x. This minimum is found by a
technique such as the method of steepest descent, the Fletcher
Powell method, or another function minimization scheme [23,24,25].
These procedures minimize the performance function iteratively
by adjusting the parameter vector in a direction dependent
on the gradient
SxT
~= i'xl' "x2 X (1.2)
(where superscript T indicates transposition).
For these optimization schemes to be useful, the
performance function e and the gradient Vc must be found in
an efficient manner. A successful technique for calculation
of the gradient components is the adjoint network approach
mentioned previously. This method evolves from application
of Tellegen's theorem [26] to the network of interest and a
topologically equivalent adjoint network.
A brief description of the method to calculate the
gradient components follows:
1. Analyze the original network. Retain branch
voltages and currents of designable elements
(vB(t) and iB(t)).
2. Analyze the adjoint network in reverse time
(T = t + tf t). Retain branch voltages
and currents of designable elements (B (T)
and (BT))
3. Compute the gradient components:
t
S iR(t) R(t + tf t)]dt (1.3)
to
for resistor branches,
s tf dvC(t)
c / dt )c(t0 + tf t)]dt (1.4)
t
for capacitor branches, and
t di (t)
= / dt L (t0 + t t)Idt (1.5)
to
for inductor branches.
S(These expressions are derived in Appendix A.)
Notice that only two complete network analyses are
required to obtain all gradient components. There are,
however, drawbacks in the method. Due to the time reversal
in the adjoint network analysis, values of branch voltage
and current must be stored for all time points during the
analysis of the original circuit. The large number of
values to be saved may cause storage problems. Also, since
the adjoint network excitation may not be a smooth time
function, the adjoint network must be analyzed with care.
In some nonlinear transient applications, design specifi
cations may be more meaningful if expressed in terms of
frequency rather than time. In the following chapters, we
present techniques for dealing with several aspects of non
linear network design for frequencydomain specifications.
The fast Fourier transform is used to obtain harmonic infor
mation about the network response. A performance function
is then formulated which reflects the frequencydomain aspects
of the design requirements. The parameterspace gradient of
this performance function is found by use of the adjoint
network method. The performance function is minimized to
realize the design requirements. In Chapter 2, we discuss
techniques for calculation of the steadystate response of a
nonlinear network driven by a periodic source. Use of the
fast Fourier transform to obtain harmonic information requires
that the network be in the steadystate condition. Network
examples illustrate the effectiveness of the steadystate
calculation. The algorithm for harmonic automated design is
developed in Chapter 3. Two practical design examples are
included. In Chapter 4, we consider special techniques
required for automated oscillator design. The determination
of oscillation frequency is one of the problems discussed.
Chapter 5 contains a summary and recommendations for future
research.
1. 3 References
1. "IBM 1620 Electronic circuit analysis program," IBM
Corporation, White Plains, New York, Application Program
1620EE02X, 1965.
2. H.W. Mathers, S.R. Sedore, and J.R. Sents, "Automated
digital computer program for determining responses of
electronic circuits to transient nuclear radiation
(SCEPTRE)," IBM Corporation, Oswego, New York, Technical
Report No. AXWLTR66126, Sceptre Users Manual, February,
1967.
3. D.A. Calahan, "Computer solution of the network reali
zation problem," Proceedings of the Second Annual Allerton
Conference on Circuit and System Theory, pp. 175193, 1964.
4. P.O. Scheibe and E.A. Huber, "The application of Carroll's
optimization technique:.to:,.network synthesis," Proceedings
of the Third Annual Allerton Conference on Circuit and
System Theory, pp. 182191, 1965.
5. D.A. Calahan, "Computer design of linear frequency
selective networks," Proceedings of the IEEE, Vol. 53,
No. 11, pp. 17011706, November, 1965.
6. R.A. Rohrer, "Fully automated network design by digital
computer: preliminary considerations," Proceedings of
the IEEE, Vol. 55, No. 11, pp. 19291939, November, 1967.
7. G.D. Hachtel and R.A. Rohrer, "Methods for the optimal
design and synthesis of switching circuits and other
nonlinear networks," Proceedings of the IEEE, Vol. 55,
No. 1 pp. 18641877, November, 1967.
8. S.W. Director and R.A. Rohrer, "The generalized adjoint
network and network sensitivities," IEEE Transactions
on Circuit Theory, Vol. CT16, No. 3, pp. 318322,
August, 1969.
9. S.W. Director and R.A. Rohrer, "On the design of resistance
nport networks by digital computer," IEEE Transactions
on Circuit Theory, Vol. CT16, No. 3, pp. 337346,
August, 1969.
10. S.W. Director and R.A. Rohrer, "Automated network design
the frequencydomain case," IEEE Transactions on Circuit
Theory, Vol. CT16, No. 3, pp. 330337, August, 1969.
11. A.J. Brodersen, S.W. Director, and W.A. Bristol, "Simul
taneous automated ac and dc design of linear integrated
circuit amplifiers," IEEE Transactions on Circuit Theory,
Vol. CT18, No. 1, pp. 5058, January, 1971.
12. B.A. Wooley, "Automated design of dccoupled monolithic
broadband amplifiers," IEEE Journal of SolidState
Circuits, Vol. SC6, No. 1, pp. 2434, February, 1971.
13. D.A. Calahan, "Optimization of switching circuits,"
Proceedings of the Second Biennial Cornell Electrical
Engineering Conference, pp. 282288, October, 1969.
14. P.M. Russo and R.A. Rohrer, "Computer optimization of
the transient response of an ECL gate," IEEE Transactions
on Circuit Theory, Vol. CT18, No. 1, pp. 197199,
January, 1971.
15. C.W. Ho, "Timedomain sensitivity computation for networks
containing transmission lines," IEEE Transactions on
Circuit Theory, Vol. CT18, No. 1, pp. 114122, January,
1971.
16. C.W. Gear, "Simultaneous numerical solution of differential
algebraic equations," IEEE Transactions on Circuit Theory,
Vol. CT18, No. 1, pp. 8994, January, 1971.
17. C.W. Gear, "The automatic integration of ordinary
differential equations," Communications of the ACM, Vol.14,
No. 3, pp. 176179, March, 1971.
18. D.A. Calahan, "Numerical considerations for implementation
of a nonlinear transient circuit analysis program,"
IEEE Transactions on Circuit Theory, Vol. CT18, No. 1,
pp. 6673, January, 1971.
19. F.G. Gustavson, W. Liniger, and R. Willoughby, "Symbolic
generation of an optimal Crout algorithm for sparse
systems of linear equations," Journal of the ACM, Vol. 17,
No. 1, pp. 87109, January, 1970.
20. W.F. Tinney and J.W. Walker, "Direct solutions of sparse
network equations by optimally ordered triangular factori
zation," Proceedings of the IEEE, Vol. 55, No. 11, pp. 1801
1809, November, 1967.
21. G.D. Hachtel, F.G. Gustavson, R. Brayton, and T. Grapes,
"A sparse matrix approach to network analysis," Proceedings
of the Second Biennial Cornell Electrical Engineering
Conference, pp. 6882, October, 1969.
22. G.D. Hachtel, R.K. Brayton, and F.G. Gustavson, "The sparse
tableau approach to network analysis and design," IEEE
Transactions on Circuit Theory, Vol. CT18, No. 1, pp. 101
113, January, 1971.
23. R. Fletcher and C.M. Reeves, "Function minimization by
conjugate gradients," Computer Journal, Vol. 7, No. 2,
pp. 149154, July, 1964.
24. R. Fletcher and M.J.D. Powell, "A rapidly convergent
descent method for minimization," Computer Journal,
Vol. 6, No. 2, pp. 163168, July, 1963.
11
25. D.G. Wilde and C.S. Beightler, Foundations of Optimi
zation, Englewood Cliffs, New Jersey: PrenticeIIall,
1967.
26. C.A. Desoer and E.S. Kuh, Basic Circuit Theory, New York:
McGrawHill, 1969.
CHAPTER 2
CALCULATION OF PERIODIC
STEEADYSTATE RESPONSE
A rapid method for periodic steadystate analysis is
of utmost importance in the harmonic design of nonlinear
networks. A conventional technique is to merely
numerically integrate the system equations until the response
becomes periodic. In lightly damped systems the initial
transient may be significant for a length of time much
greater than the time of one period. In these cases, the
method of continuously integrating the differential
equations until the transient part of the solution becomes
negligible is unacceptable, since a considerable amount of
computer time is wasted on calculation of the initial
transient. In this chapter, a method is presented which
uses a Newton iteration to obtain the periodic steadystate
response of a nonlinear network.
2.1 Background
Several methods have been described for the determin
ation of the periodic steadystate response. One such
method proposes the use of a "residual" which indicates
how well an assumed solution satisfies the system
differential equations. This technique is referred to as
Galerkin's method [1], and was used by Baily [2,3] for
steadystate analysis of nonlinear circuits. The solution
of the nonlinear differential equation
x(t) f[x(t), t] = 0 (2.1)
is assumed to be a linear combination of suitable functions.
Since a periodic solution is sought, the solution is
expressed as a Fourier series:
L
x(t) t) Xt) = a + a a sin(nt) + b cos(nt), (2.2)
n=l
where L, the number of harmonics, is chosen to yield suitable
accuracy. Since the solution X(t) is not exact, if it is
substituted into the original differential equation a non
zero residual results:
e(t) = X(t) f[X(t), t] (2.3)
This residual is a measure of how well the assumed solution
satisfies the original equation. The square of the residual
integrated over the time interval of interest is then
minimized by suitable adjustment of the coefficients a
n
and b Any of the multidimensional search routines such
n
as FletcherPowell [4] or conjugate gradients [5] may be
used for this purpose. These algorithms require calculation
of the gradient of the residual with respect to the
coefficients a and b
n n
In Baily's formulation of Galerkin's method, the state
equations of the network are written and the above approxi
mation is made for each state variable. An alternate
method would be to use nodal analysis. The above approxi
mation would then be used for each node voltage.
Galerkin's method is easily illustrated in terms of
an example. Consider the network of Fig. 21. The equation
describing this system ip.
Gv(t) + Cv(t) + f(v) A sin t = 0, (2.4)
where
f(v) = I [exp(Xv) 1]. (2.5)
The assumed solution is
L
V(t) = a0 + [an sin(nt) + b cos(nt)], (2.6)
n=l
and the residual is
e(t) = GV(t) + CV(t) + f(V) A sin t. (2.7)
Since we wish to minimize
t +T
E = f/ [e(t)]2 dt, (2.8)
to
the gradient components , c ... are required. These
a 1
gradient components are easily calculated:
t 0+T
t0+T
= = e[G + XI exp(XV)] dt, (2.9)
0 t, s
id=Is [exp (Xv) 1]
Fig. 21
Circuit used in examples
of steadystate determination methods
A sin t
t0+T
C / e[{G + XI exp(XV)}sin(nt) + n cos(nt)]dt,
n t0 (2.10)
t T
t0+T
= / e[{G + XI exp(XV)}cos(nt) n sin(nt)]dt,
n tO (2.11)
for n = 1, ***, L, where L is the number of harmonics to be
included. A FORTRAN computer program was written to test
Galerkin's method. A "package" conjugate gradients search
routine [6] was chosen. The algorithm is illustrated in the
flow chart of Fig. 22. Results of the analysis for several
values of L are given in Table 21.
A drawback of Galerkin's method is that if L harmonics
are required for suitable accuracy and there are m variables,
then m(2L+l) parameters must be determined. In practical
examples this number can be quite large. Moreover, optimi
zation methods tend to use excessive computer time,
especially if the number of optimization variables is large.
As may be seen from Table 21, even a very simple network
may require an excessive amount of computation.
A more efficient scheme emerges upon consideration
of what initial values of capacitor voltage and inductor
current yield periodic steadystate operation. First
we recognize that if a network operates in the periodic
steadystate, all voltages and currents in the network are
periodic. In particular, if the period is T then for all
SStart.
Spe ci fy
initial
estimate of
solution
Calculate error:
Approximate
t +T
S= /0 1 e2dt
t
by
= (h.1 2 1 2
2 2 n 1 2 n
= ( 1 2 en
(Trapezoidal Rule)
Small yesStop
Calculate gradients:
use Trapezoidal Rule
to numerically find
t +T
f =f / e[G+XI exp(Xv)]dt
o t
0
t +T
/ e[{G+XI exp(Av)}sin(nt)+n
t
0
t 0+T
f e[{G+Xs exp(Av) }cos(nt)n
t
0
cos (nt) ]dt
sin (nt) ] dt
for harmonics i = 1, N
Adjust parameters
Fig. 22
Flow of operations for implementation
of Galerkin's method
Da.
1
b=
I
Table 21
Results of Application of Galerkin's Method
to Circuit of Fig. 21
L 1 2 3 4
do 1.432 1.096 1.102 1.049
sin 1 1.688 1.817 1.836 1.833
2  0.5987 0.6016 0.6224
3   7.052E2 5.163E2
4    8.217E2
cos 1 1.219 1.342 1.298 1.308
2  0.1687 0.1018 0.1170
3   2.229E2 0.1439
4    6.112E2
Final error 6.765 3.911 3.506 2.029
Iterations 43. 121. 158. 240
CPU time 4.2 Sec. 34.3 sec. 62.6 sec. 177.2 sec.
capacitor voltages
vC(t) = vC(t + T), (2.12)
and for all inductor currents
iL (t) = i (t + T). (2.13)
At the initial time tO, we wish to determine a capacitance
voltage vC(t0) and inductance current iL(t ) such that
vC(t0 = V 0(t + T) (2.14)
and
iL(tO) = iL(t0 + T). (2.15)
It is convenient to model the capacitor with initial
voltage vC(t0) by a series combination of a capacitor with
zero initial conditions and a dc voltage source of value
ve = V(t 0). Inductors with initial current iL(t0) are
modeled by an inductor with zero initial conditions and a
parallel dc current source i = iL (t0) as shown in Fig. 23.
Then
vC(t) = vC' (t) + ve (2.16)
and
iL(t) = i '(t) + ij, (2.17)
where
V (t0) = 0
(2.18)
+ Lict)
vc' (t)
()ve = VC(tO)
+ C (t)
v (t) c
t+ i,(t)
vL (t)
J'Lit L
+ C(t)
vC (t
0f
qc(t) = f (vc(t),PC't)
d
iC(t) dt qC(t)
+ i (t)
vL (t
T
XL(t) = fL(iL(t),PLt)
vL(t) dt L(t)
+ T (t)
v,' (t) y
C(t) f (v (t)pt)
Ct) dt
+
vL (t)
= iL(t)
L(t) = fL' (i L'(t),PL't)
d
vL(t) dt (t)
Fig. 23
Capacitor and inductor initialcondition models
(a) Linear timeinvariant
(b) Nonlinear
+
vL (t)
ij = i(t0)
and
iL (t0) = 0. (2.19)
Thus if initial condition source values are chosen to cause
steadystate operation,
vC' (t0 + T) = 0 (2.20)
and
iL (t0 + T) = 0. (2.21)
The determination of these 'initial condition source values
may be accomplished in two ways: an optimization procedure
or a Newton iteration. In the optimization approach, an
error criterion is formulated, namely
E[= [ 1 [v(to + T)]2 + .1 [i' (t + T)]2. (2.22)
C L
This error function is zero if all initial capacitor voltages
and inductor currents are correct and greater than zero
otherwise. The object is then to minimize e by adjusting
the initial condition source values v and i.. All that is
e 3
required to use a .functio .miinimization algorithm to
perform this adjustment is the gradient of E with respect
to all the initial condition source values. The number of
optimization variables has thus been reduced to one per
energystorage element. The undesirable convergence
properties of optimization procedures (e.g., converging to
an incorrect local minimum) still remain. For this reason
a NewtonRaphson technique is preferred. However, in the
interest of completeness, the optimization method is
described in the following.
The adjoint network approach [7] yields an efficient
means for evaluating the necessary gradient components.
Consider the situation depicted in Fig. 24. We wish to
calculate the steadystate response of networki2, in which
capacitors and inductors have been modeled as discussed
earlier. The network is the topologically equivalent
adjoint network. Application of Tellegen's theorem [8] to
this situation (over one period) allows us to write
t +T
/ tT [vB(t)B(T) ig(t)PB(T)] dt = 0, (2.23)
t B
where vB and iB are the Bth branch voltage and current in
the original network, ~B and B are the Bth branch voltage
in the adjoint network /, and the indicated summation is
over all network branches. Upon an incremental change in
the parameters of the original network, (2.23) becomes
t +T
t0+T
/ [{vB(t) +Av (t)} B((T)
to B
{iB(t) + AiB(t) }B(T)] dt = 0, (2.24)
so that
t0+T
t B
ii v Li
VL V
vv (c
v v
used in gradient calculation
IvII c
R.
11.
i L_
v
n e
i\__L
Fig. 24
The original and adjoint networks
used in gradient calculation
This equation may now be written
t +T
/t
0
S [Av (t) (T) Ai (t)i (T)] dt
v
t0
t, i
t0
+ /
to
t0+T
+ /
t0+T
0
[Av.(t).i(T) Ai.(t) .i(T)] dt
[Av (t) j(T) Ai (t)4 (T)] dt
[Av. (t)(.(T) Ai .(t)$.(T)] dt
[ [AvN(t) (N T) AiN(t) N(T)] dt, (2.26)
N
where the summations over v, i, e, j, and N represent,
respectively, independent voltage source branches, independent
current source branches, capacitor initial condition voltage
sources, inductor initial current sources, and nonsource
branches. Since the values of independent sources do not
vary,
Av (t) = Aii(t) = 0.
(2.27)
The adjoint network branches which correspond to sources
in the original network are zerovalued sources. Thus
9 (T) = (.(T) = 0,
V 1
and
e (T) = 0.j(T) = 0.
Therefore, (2.26) becomes
tO"+
Ave (t)e) (T) dt +T
e to
t +T
= 
to
SAi (t) j(T) dt
j
S [AVN (t)N (T) AiN(t)N(T)] dt. (2.30)
N
The nonsource branches are now considered. Only nonlinear
resistor, capacitor, and inductor branches are assumed to
make up the original network. Other network element types
are considered in Ref. [7]. Equation (2.30) thus becomes
t0+T t+T
/ X Av e (T)dt / I Aijqj(T) dt
tO e to 3
t0+T
= /
t
0+T
tO
t0+T
0
tO
0
I/
to
S[AvR(t) R(T) AiR(t) R(T)] dt
R
I [Av' (t) C(T) Aic(t)lc(T)] dt
C
S[AvL (t) L(T) Ai' (t)L(T)] dt.
L
(2.28)
t +T
t0
(2.29)
(2.31)
By a proper interpretation of network variablesas in
Appendix Awe may write
t +T
t0+
to
SAv ee (T) dt 
e
t 0+T
j
to
}Ai i (T) dt
j JJ
t +T
= /t0+T
t 0
to
t +T
C to
f R(t)
SP R(T) R dt
S8f '(t)
C( dt pC aC] dt
Dfc (tf)
+ v (to) AvC '(tf)
8f (0 )
c t 0 (t)) Avy (tO)
C
+ t0+T
L t
d pf (t)
[L() dt PL AL] dt
L
Df (t )
 aiL (t0) AiL'(tf)
L
DfL (to)
+ i L( Ai L ) ,
L4
(2.32)
where tf = t0 + T and T = T t. Since energystorage
elements are modeled as described before,
qC(t) = f C(vc' (t),P',t) = fc(vc(tO) + VC'(t),PC't)
(2.33)
and
L(t) L I(iL(t)PL't) = fL(iL (t0) + iL (t),Pct)
(2.34)
Because element parameters pR', PC and pL do not vary,
e(T) dt Av /0
j t0
j (T) dt Aij
af'' (t ) a c '(t
C (t ) v f c c(tf)Avc (to)
C *C c
fL' (tf)
L iL L 0 L f
DfL (t )
iL L (tf) (t 0
(2.35)
Now, by letting
Av' (t0) = AiL' (t) = 0
v '(t )
C(t0) = fc'(t )
VC
iL (tf)
L(t0) = F (tf)
L
and dT = dt ,
we obtain
(2.36)
(2.37)
(2.38)
(2.39)
t +T
e (T) dTAve / +
j to
_.(T) dT i.
3 3
t0+T
e tO
t0+T
e tO
Sv (tf)v' (tf) + i (t'f) L' (tf)
C L
The variation of the error (2.22) is
AE = vc' (tf)AvC' (tf) + : iL' (tf)AiL' (tf)
C L
= : A v + i Ai. .i
e e j j 3
Upon consideration of (2.40) and (2.41) we may write
(2.40)
(2.41)
DC t +T
8v /
e tO
Th
j
 t0+
t0
8f (t ) Df '(tf)
ge(T) dT = avc e(tf) vc e(t0)
(2.42)
T 3fL (t0)
j (T) d= 3im j (tf)
(fL f)
L (t(t )
Di L 0
(2.43)
For linear timeinvariant capacitors and inductors, (2.42)
and (2.43) become
e
9v
e
C [ (t ) 9(t )]
e f e 0~)
and
D C
i .
j
Notice that two network analyses suffice to obtain the
required gradients.
(2.44)
(2.45)
and
 L [ j(tf) (to) .
The following algorithm emerges:
1. Estimate initial condition source values v
and i.. e
3
2. Analyze the original network over one period T.
3. Calculate the error c.
4. Analyze the adjoint network over one period in
reverse time (, = T t), letting initial
capacitor voltages be
vC' (tf)
'C (t0 (tf)
vL
or for linear timeinvariant capacitors
vC f
I c(to) = C
C 0 C '
and initial inductor current be
i L (tf)
S= fL (tf)
aVv
V C
or for linear timeinvariant inductors
iL' (t f)
5. Determine the gradient components by calculating
e FfC' (tf) fC' (t )
8v v e f ay v e 0
e C C
and
j 8fi (t ) 8f '(t )j(t)
= (t ) (t t ) ,
SLL
or for linear timeinvariant elements,
= C[e (tf) e(t0)]
e
j L[j.(tf) j(t]
6. Alter the initial source values according to
some gradient search technique (e.g., Fletcher
Powell or conjugate gradients).
A FORTRAN program was written to implement this
algorithm. The flow chart of the resulting program is shown
in Fig. 25. The results from application of this program
to the network of Fig. 21 are given in Table 22. Since
this technique yields the steadystate time response of
the network, the fast Fourier transform must be used if
harmonic information is desired. In the program example,
the IBM Scientific Subroutine Package fast Fourier transform
routine HARM [6] was used. Comparison of Tables 21 and
22 indicate the speed advantage of the adjoint network
optimization technique. However, when the adjoint network
optimization approach was applied to the circuit of
Fig. 26, a local minimum caused the method to yield an
incorrect solution. Results of this calculation are given
Start
Initial
estimate of
initial condition
Analyze network
over one
period
Calculate error e
c smal Yes
Analyze adjoint
network in reverse
time: T = t + tf t
Calculate gradient
components
Adjust initial condition
source values
Fig. 25
Flow of operations for adjoint network
optimization method of steadystate calculation
Table 22
Results of Application of Adjoint
NetworkOptimization Method to Circuit of Fig. 21
dc term
sin
coefficients
cos
coefficients
1.086
1.852
0.5405
2.204E2
1. 450E2
5.226E2
1.273
8.922E2
1. 591E1
5.296E2
1. 807E2
final error 4.5E13
iterations 3
CPU time 3.52 sec
i
if
10 sin(27Tt)
1 f
Fig. 26
Voltage doubler circuit of the second example
in Table 23.. The fact that the optimization approach was
unable to yield the correct solution for such a simple
network spurred the development of a more reliable method.
2.2 The State Variable Newton Method of SteadyState
Calculation
A Newton method for calculating the steadystate
response of a nonlinear network was described by Aprille
and Trick [9,101. In their approach, state variable
analysis was used to obtain the response over one cycle.
A state transition matri: is calculated for an "auxiliary"
system of linear differential equations with timevarying
coefficients. This state transition matrix is used to
compute the new set of initial condition values. The
periodic steadystate solution is sought for the system of
equations
x= f(x, t). (2.46)
If the backward Euler numerical integration formula is applied
to (2.46), we obtain
x(jh) = x[(j l)h + bh f[x(jh), jh] (2.47)
where h is the step size. Newton's method may be used to
solve (2.47):
i+l i i
x (jh) x (jh) {I hF[x (jh)]}l
{x (jh) x[(j l)h] hf[xi(jh)]} (2.48)
Table 23
Results of Application of Adjoint
NetworkOptimization Method to Circuit of Fig. 26
v (to) v (t )
Iteration CPU Time* 1 2 Error* *
24.9
35.4
43.0
52.4
66.7
70.5
77.3
0.394
3.108
4.729
4.415
4.349
4.348
4.34's
*Seconds IBM 360/65 CPU Time
A IC1 (to) + 4.3381
**Error = 33 +
4.338
VC2 (t) 9.8021
9.802
5.184
10.320
10.133
9.994
9.840
9.838
9.838
1.561
0.336
0.124
3. 734E2
6.412E3
5.978E3
5.978E3
where
F[x (jh)] = f[x(jh),jh]
~x(jh)
(2.49)
is the Jacobian of f. The "auxiliary" system of equations
is
Z = F[x(t)] Z .
(2.50)
Application of the backward Euler formula to (2.50) yields
Z(h) = Z(0) + hF[x(h) ] Z(h)
= {I hF[x(h)]}1 Z(0) ,
(2.51)
so that
k
k 1
Z(T) = Z(kh) = {I hF[x(ih)] } Z(0)
i=l
(2.52)
Thus, the state transition matrix 0 of (2.50) is simply
S= H {I hF[x(ih)]} .
=1
(2.53)
Define x0 to be the vector of initial condition values.
The Newton iteration to determine these values is
x = x [I T'(x [ T(x)]
~0 0 0 H 0 0
= [I T'(x)]1 [T(x ) T'(x )x ]
S[{T' (x )}. Ii [{T' (x )} T(x) x] (2.54)
0 e 
where
T(x3) x3(T, x3) (2.55)
0 ( 2 5
and
ax(T,x )
T'(x0) = (2.56)
~ ~ 0 3
x
I0
Aprille and Trick showed in their presentation that
k
T' (x ) = = {I. hF [x (ih)]} (2.57)
~ ~ ~ i=1 ~ ~ ~
and
S1 
{T'(x3)} = = H {I hF [x (ih)]}. (2.58)
01
Hence, from (2.54),
j+1 1 1 xj
1 [1 I] [1 x(T) x0] (2.59)
gives the new initial state. This "state variableNewton"
method must be implemented with a state variable analysis
scheme. Moreover, it is dependent upon the numerical
technique used to integrate the systems of equations (2.46)
and (2.50). In the following section, the adjoint network
concept is used to develop a Newton method for computation
of steadystate response which avoids these restrictions.
2.3 The Adjoint Network N Newton Method of SteadyState
Calculation
Define x0 to be the vector of initial condition
0
source values:
A T
x = (ve, ij) (2.60)
~0* u ev ~ i
and a function f :
f = f(x0) [V (tf), L' (tf)] (2.61)
which is zero for the proper choice of initial conditions
x0. Newton's method for determining the values of x
which cause f to be zero may be expressed as
F Ax = f (2.62)
where Dv '(t ) 3v '(tf)
Dv Di.
~e
F   (2.63)
iL (tf) L (tf)
Dv Di.
~e 3
is the Jacobian matrix. The "change vector" Ax is obtained
by solving the simultaneous linear equations (2.62). This
solution may be accomplished by a Gaussian elimination or
LU factorization procedure [11]. The change vector is
k
then added to the previous iterate x to obtain the new
iterate:
k+l k
x0 = x + Ax
When Ax is sufficiently small, the iteration is terminated.
In order to apply Newton's method to the steadystate
determination, the partial derivatives which make up the
Jacobian matrix must be determined. The adjoint network
concept may be employed to calculate these partial
derivatives efficiently. Once again, capacitors and
inductors are modeled as shown in Fig. 23. Upon appli
cation of Tellegen's theorem to the original networkIZ
and its mutually reciprocal adjoint network Zas shown
in Fig. 24we may write
t +T
/tT [vB(t)B)(T) iB(t)B(T)] dt = 0 (2.64)
t0 B
where vB and i are the Bth branch voltage and current in
the original networkZ, ,B' and )B are the Bth branch
voltage and current in the adjoin' network't, and the
indicated summation is over all network branches. Upon an
incremental change in the parameters of the original network,
(2.64) becomes
t +T
t0+T I [{vB(t) + AvB(t)} (B(T)
t B
{iB(t) + Ai(t)} B(T)] dt = 0 (2.65)
Subtraction of (2.64) from (2.65) yields
t +T
0 [AvB(t)cB(T) AiB(t)BB(T)] dt = 0. (2.66)
t B
Equation (2.66) may be rewritten as
t0+T
/ [Av (t)v (T) Ai (t)i (T)] dt
t v
0t+T
+ [Av.(t)4.(T) Ai.(t).i(T)] dt
0to
t +T
+ / [Av (t)A (T) Ai (t)4) (T)] dt
t e e e e
to
+ 0 [Avj(t) j(T) Ai (t)Mj(T) ] dt
t +
= / X [AvI(t)N(T) AiN(t)N(T)] dt ,
to N
0
(2.67)
where the summations over v, i, e, j, and N represent,
respectively, independent voltage source branches, independent
current source branches, capacitor initial condition voltage
sources, inductor initial condition current sources, and
nonsource branches. Since the values of the independent
sources are constant,
AV (t) = Ai (t) = 0 .
(2.68)
The adjoint network branches which correspond to sources in
the original network are zerovalued sources, so
v (T) = (T) = 0
V 1
(2.69)
(2.70)
( (T) = i (T) = 0 .
e J
Thus, (2.67) becomes
to+T t +T
SAV e e(T) dt /
t e t
SAi" j(T) dt
j
= /+T
to
I [AvN(t)N(T) AiN(t)pN(T)] at.
N
and
(2.71)
The nonsource branches are now considered. For brevity
in this presentation, only nonlinear resistor, capacitor,
and inductor branches are assumed to make up the original
network. Other elements are considered in Ref. [7].
Equation (2.71) thus becomes
t +T
t0
t +T
Av e) (T) at t
ee to
t+T
t0+T
to
t +T
t
Ai (T) dt
J j
I [A'v (tip(T) Ai (t)R(T)]dt
R
[AvC' (t)(C () Aic(t)lC(T)] at
C
[ [AvL(t) L (T) AiL'(t) L(T)] dt (2.72)
L
By a proper interpretation of the adjoint network
variablesas in Appendix Awe may write
t0+T
tO
t +T
SAv (T) dt /+T
e to
SAi j (T) dt
J
_ o0 8f (fR(t)
ft R i PR p R(T) Ap dt
t R ap R R
+ +T d fC' (t)
C+ [t o( Z) PC Apc] dt
fC' (t )
+ 3vcC c(t0) Av '(tf)
Sfc( (to)
v, c~(tf) Avc' (tO)
t0+T fL(t)
L+ t L Ltp L] dt
I f (tf)
iL L(t0)AiL f)
8f, (t )
+ iaiL L(tf) AiL'(t) (2.73)
where tf = t0 + T and T = T t. Since energystorage
elements are modeled as an element with zero initial
conditions and a source as shown in Fig. 23,
qC(t) = f' (vc'(t),Pc,t) = fC(vc(t0) + vc'(t),PCt)
and
L(t) = fL '(iL' t),PLt) = fL(iL(tO) + iL'(t),PLrt)
(2.74)
Because element parameters pR, PC' and PL do vary,
t0+T t+T
S/ 4e (T) dt Av x / 0 j(T) dt Ai.
e t e e jt0
v fc' (tf)
c '(t0)
C v C C 0 C f
Cf C (t ) c (t)AV '(t
S v C C f (to)C 0
43
f '(tf)
Lf
L L L L 0 L f)
f L' (t )
i C c(tf)Ai (Lo) (2.75)
By letting
Avc' (t) = AiL' (t0) = 0 (2.76)
Df (t f)
c (t0) = 1/ vc (2.77)
C. (to) = 0 \j i (2.78)
L (t0) = 0 (2.79)
and
dT = dt (2.80)
we obtain
t t
S/ c (T)dT Av J / t (T)dT Ai
e tO e j t 3
= Av (tf) (2.81)
1
Since
8v' (t ) 8v' (t )
C f CD f
Av' (tf) = 1 Av + Ai, (2.82)
i e e 3 3
the partial derivatives are seen to be
Dv (tf) t
S = e(T) dT
ei 0 k
44
CfC (t0) Cf (tf
vck ek(t) (t)
k k (2.83)
Dv (tf t
C. f tf
i. f (T) dT
jm to m
Df (t0) fm (t) t
mL m m
DiL (tf) DiL
m m (2.84)
DiL (tf)
Dvek= / (e(T) dT
av e
ek t k
k 0
af' (to) Dfk (tf)
t0) t Pe (t )
vCk eek(tf C ek(t0)
Dvk k k k (2.85)
and
iL. (tf) t
1 tF
Di. = / .jm (T) dT
3m t 0
DfL (t0) DfL (tf)
m m
= (t ) i (t) .
m mm (2.86)
For linear timeinvariant elements, these partial
derivatives are
v'. (tf)
v = C [ tf) ek(t0)] (2.87)
Sk k k
Dv. (t f)
L [ (t) ) j(t 0)] (2.88)
mm m
i' (t )
L. f
= [ek(tf ek (t0), (2.89)
ek k k
and
iL. f
i
S = L [L j (t, ) (t0)] (2.90)
Di. m (tf 0
3m m m
Thus, all that is required is to set one capacitor or
afC(tf)
inductor initial condition to the value 1/ or
f' (t ) C
1/ i respectively (or to 1/C or 1/L for
linear timeinvariant elements), in the adjoint network.
The adjoint network analysis is performed and the partial
derivatives of the final condition of that element with
respect to all initial conditions are calculated by
performing the subtractions indicated in (2.83) through
(2.86) or (2.87) through (2.90). These partial derivatives
occupy one row in the Jacobian (2.63). This procedure is
repeated, setting a different initial condition in the
adjoint network until all partial derivatives have been
calculated. A total of M adjoint network analyses are
thus required if there are M energystorage elements in the
network. The (M + 1) total analyses are the same number
as required to calculate the partial derivaties by
perturbation of the initial conditions. However, the
adjoint network method yields the exact partial derivatives
rather than the approximate value arrived at by the pertur
bation technique.
Once all the partial derivatives have been determined,
Newton's method (2.62) is applied to obtain the new initial
condition source values. The resulting algorithm is
summarized below:
1. Analyze the original network over one period.
2. Set initial condition in adjoint network and
perform analysis of adjoint network.
3. Repeat step 2 until all Jacobian entries have
been calculated.
4. Solve for change vector.
5. Calculate new initial conditions.
6. If change vector is greater than desired
go to step 1.
7. Stop.
A variation of the NewtonRaphson method may be used
to decrease the computation required. The Jacobian is
evaluated once and used until the error either stops
decreasing or increases. This procedure is known as the
"fixed tangent" method [11]. The algorithm which results
from incorporation of this technique is illustrated in
Fig. 27.
To establish the effectiveness of the algorithm, two
practical examples are presented. The first example is
the voltage doubler circuit depicted in Fig. 26. The
Fig. 27
Flow of operations for adjoint network
Newton method of steadystate calculation
steadystate solution was obtained by two methods:
continuous integration and the adjoint networkNewton
algorithm described earlier. Gear's numerical integration
procedure was used. Both methods were started at
V = VC = 0. The period was one second. Results of
I 2
this example are summarized in Table 24. The speed
advantage of the adjoint networkNewton procedure is
obvious.
The second example 'is the circuit shown in Fig. 28.
Again, continuous integration and the adjoint network
Newton method are compared. Initial conditions were
C = 2 = vC3 = i = 0. The period was 1/60 second.
c 2 3
A considerable saving in computer time is realized by use
of the adjoint networkNewton method as may be seen in
Table 25.
2.4 References
1. W.J. Cunningham, Nonlinear Analysis, New York:
McGrawHill, 1958.
2. E.M. Baily, "Steadys.tate harmonic analysis of nonlinear
networks,"Ph.D. The'sis, Stanford University, Stanford,
California, 1968.
3. E.M. Baily, "Steadystate harmonic analysis of nonlinear
networksa computeroriented approach," Proceedings
of the First Annual Houston Conference on Circuits,
Systems, and Computers, pp. 390399, 1969.
4. R. Fletcher and M.J.D. Powell, "A rapidly convergent
descent method for minimization," Computer Journal,
Vol. 6, No. 2, pp. 163168, July, 1963.
5. R. Fletcher and C.M. Reeves, "Function minimization
by conjugate gradients," Computer Journal, Vol. 7,
No. 2, pp. 149154, July, 1964.
6. "IBM System/360 scientific subroutine package," IBM
Corporation, White Plains, New York, 1966.
7. S.W. Director and R.A. Rohrer, "The generalized adjoint
network and network sensitivities," IEEE Transactions
on Circuit Theory, Vol. CT16, No. 3, pp. 318322,
August, 1969.
8. C.A. Desoer and E.S. Kuh, Basic Circuit Theory,
New York: McGrawHill, 1969.
9. T.J. Aprille and T.N. Trick, "Computeraided steadystate
analysis of nonlinear circuits," Proceedings of the 14th
Midwest Symposium on Circuit Theory, 1971.
10. T.J. Aprille and T.N. Trick, "Steadystate analysis of
nonlinear networks that have a periodic response,"
University of Illinois Coordinated Science Laboratory
Report R515, June, 1971.
11. C.E. Fr6berg, Introduction to Numerical Analysis,
Reading, Massachusetts: AddisonWesley, 1969.
Table 24
Results of Application of SteadyState Determination
Methods to Circuit of Fig. 26
Continuous Integration Adjoint NetworkNewton Method
v (t ) (t ) CP vC (t0) C (t0) CPU
Iteration 1 C2 Error* Time** 1 2 Error* Time**
1 0.3091 2.570 1.809 3.0 0.3091 2.570 1.809 11.8
2 .2518 4.160 1.518 4.2 4.094 9.067 0.131 12.7
3 .9192 5.301 1.247 5.4 4.310 9.705 1.6E2 13.5
4 1.462 6.155 1.035 6.2 4.330 9.801 2.1E3 14.4
5 1.950 6.821 0.895 7.1 4.338 9.802 0.0 15.3
6 2.361 7.360 0.705 7.9
7 2.703 7.797 0.581 8.0
8 3.015 8.128 0.476 9.6
9 3.239 8.425 0.394 10.4
10 3.424 8.671 0.326 11.2
20 4.194 9.646 4.91E2 18.8
40 4.334 9.792 1.94E3 33.3
100 4.335 9.803 7.94E4 75.2
Ivcl(to) + 4.3381
*Error = 4.338
1vC2(to) 9.802j
9.802
**Seconds IBM 360/65 CPU time
5Q 0.1 h
6 L0
10 6 L
+
+ + v _
L C v v 1 k
C C
V2 3
3
10 f LO 3f
10 sin(21r60t)
Fig. 28
Circuit for adjoint networkNewton method example
0
Table 25
Results of Application of SteadyState Determination Methods
to Circuit of Fic. 27
Continuous Integration Adjoint NetworkNewton Method
VC (to) i ) CPU vC (t) i (t) CPU
Iteration 1 L 0 Error* Time** 1 L 0 Error* Time**
1 1.970 0.3665 40.40 5.1 1.970 0.3665 40.40 36.5
2 3.365 1.160E3 1.76 8.4 7.311 2.054E2 1.470 38.7
3 7.635 1.364E2 0.67 11.6 8.343 2.863E2 2.252 41.0
4 6.222 0.1459 15.48 14.2 8.695 1.838E2 1.078 43.2
5 7.462 6.471E2 8.35 16.9 8.864 1.392E2 0.5648 48.4
10 8.287 9.572E2 11.69 29.1 9.053 9.320E3 3.424E2 56.2
14 9.119 2.587E2 1.87 38.6 9.066 9.024E3 0.0 64.8
20 8.930 2.546E2 1.84 52.4
30 9.108 1.231E2 0.37 75.0
40 9.078 6.620E3 0.27 98.0
50 9.059 8.523E3 5.63E2 120.9
60 9.064 9.392E3 4.11E2 143.6
70 9.068 9.098E3 8.42E3 167.2
80 9.067 8.963E3 6.87E3 190.2
90 9.066 9.009E3 1.66E3 213.2
100 9.066 9.028E3 4.43E4 236.5
Error =
I Cl (t0)
+ 9.066I liL(tO) 
+
9.024 1031
9.024 103
**Seconds IBM 360/65 CPU time
9.066
CHAPTER 3
FREQUENCYDOMAIN ASPECTS
OF NONLINEAR CIRCUIT DESIGN
In some instances, nonlinear networks are designed
to meet specifications in terms of frequencydomain, rather
than timedomain, behavior. The rationale for development
of automated design methods for these network types may be
established by several practical examples.
Largesignal amplifiers are often designed by methods
involving linearization or graphical constructions [1]. The
added complexity of integrated circuits decreases the
desirability of these approaches. Possible design require
ments for largesignal amplifiers might include specifications
concerning:
1. Power efficiency; i.e., the ratio of average
power delivered to the load to the power
consumed by the amplifier;
2. Harmonic distortion;
3. Maintenance of .. and 2. within some tolerance
over a specified frequency range.
Consider the specification of harmonic distortion.
Design for minimum distortion makes a great deal more sense
if specifications are formulated in frequencydomain terms.
There is no salient feature of a time waveform which reflects
distortion. Measuring distortion by taking the difference
between the response and a pure sinusoid is inadequate
because the two signals must have the same frequency and
phase for the subtraction to indicate waveform distortion.
Moreover, even if the two waves did have the same frequency
and phase, a loss of numerical accuracy results from
subtraction of two similar numbers. In the frequency domain,
distortion is found easily: it is simply the ratio of the
sum of the harmonics to the fundamental component. Thus, a
method which exploits the frequencydomain nature of this
problem is required.
Frequency multiplier circuits make use of the nonlinear
characteristics of some electronic circuit elements to
generate harmonics. Usually only one of these harmonics is
desired. Graphical methods are often used in the design of
frequency multiplier circuits [1]. Here again, an automated
design method based on frequencydomain informationthe
desired harmonic may be advantageously applied.
Finally, sinusoidal oscillators are generally designed
to meet specifications regarding
1. frequency of oscillation;
2. power output;
3. spectral quality (harmonic content of output);
4. efficiency;
5. frequency and amplitude stability.
In this chapter, general concepts are developed which
lead to an effective algorithm for the automated design of
nonlinear circuits for specified harmonic content. The
nature of oscillator operation dictates a need for special
methods. Consideration of these methods is undertaken in
Chapter 4.
3.1 The Design Algorithm
The algorithm for harmonic design is based on the time
domain automated design concepts discussed in Chapter 1.
The proposed method can be illustrated in terms of a large
signal amplifier. Let v(t) denote the actual output wave
form of the amplifier. Under reasonable assumptions (e.g.,
finite energy content [2])the Fourier transform of v(t)
exists:
Vfj J2nft A
V(f) = v(t) ejf dt I=[v(t)]. (3.1)
00
The inverse transform is
v fe j2'itft A 1
v(t) = V(f) e it df =2 [V(f)]. (3.2)
00
Only a discrete set of sampled values of v(t) is actually
stored in the computer... Therefore, in computer work the
discrete Fourier transform is used. Let v(k) represent the
sampled values of v(t) at the equally spaced times t = k At,
k = 0, 1, **, N 1, where N is the number of samples.
V(n) denotes the discrete values of V(f) at the equally
spaced frequency points fn = n Af, n = 0, 1, ".*, N 1.
The discrete Fourier transform is [3,4]:
N1
V(n) = [ v(k) ej2 nAf kAt
k=0 (3.3)
n = 0, 1, ***, N 1; and
Nl
v(k) = V(n) ej2" nAf kAt
n=O (3.4)
k = 0, 1, ***, N 1 .
Since At = C/N and Af = 1/J where 3 is the time interval
of interest,
Nl
N (k) j2nk/N A
V(n) = 1 v(k) e =aD [V(n)],
(3.5)
n = 0, 1, ***, N 1; and
N1
v(k) = 1 V(n) ej2nk/N [V(n)],
n=o (3.6)
k = 0, 1, ', N 1.
The Fourier series coefficients may be shown to be
1
D = V(n)
n N
(3.7)
n = 0, 1, "', N 1.
In practice, the calculation of the discrete Fourier
transform is accomplished by the fast Fourier transform
(FET), which is also called the CooleyTukey algorithm [5].
The output of a distortionless amplifier, v(t), is a
sinusoid of appropriate amplitude with period T = 1/fl,
where fl is the input signal frequency. The discrete Fourier
transform of v(t), denoted by V(n), has nonzero components
only at the fundamentalfrequency (n = 1) and possibly a dc
component (n = 0).
Consider the performance function
e(x) /T [e(t)]2 dt, (3.8)
2 0
where
A
e(t) = v(t) (t). (3.9)
The performance function is zero when the actual output is
distortionless and nonzero otherwise. At its minimum it
yields the optimum parameter values x. Since v(t) and v(t)
could be about equal in magnitude, the indicated subtraction
could give rise to numerical errors in e(t) and thus in E.
From the relationship [2]'
T
F2 = 1f2(t)dt (3.10)
S0
for the mean squared value of the signal
00
f(t) = D ej2Tnt/T
SD + 2 D ej2 nt/T (3.11)
n=l
we may write
0O
2 2 D2
F2 = D + 2 D (3.12)
n=l
where D = ID Equation (3.12) is a result of the
n n
orthogonality of the terms of the Fourier series. This
equation and (3.7) may now be used to evaluate the performance
function (3.8):
T
T ^ 2
e = / [v(t) v(t)] dt
0
L
T ^ 2 2
2 J(Do Do + 2 1 (Dn D) }
n=l
L
{[V(0) V(0)] + 2 1 [V(n) V(n)]2}. (3.13)
n=l
The summation is taken over the fist L harmonics. The value
of E is seen to be calculated directly from the results of
the fast Fourier transform: no integration over time is
necessary. For design flexibility, nonnegative weights may
be inserted:
L
T {W[V(0) V(0)]2 + 2Y Wn[V(n) V(n)]2}. (3.14)
The gradient of e with respect to all designable parameters
remains to be found. As indicated earlier, the adjoint
network concept is an extremely efficient method for finding
gradients. An additional concern is calculation of the
adjoint network excitation signals. The adjoint network
will be excited by the error signal
Q(t) = w(t) [v(t) v(t)]
(3.15)
or
c(k) = w(k) [v(k) v(k)]. (3.16)
Again, to avoid numerical inaccuracies which may arise from
taking differences of approximately equal numbers, the fast
Fourier transform is used:
1.
#(k) =D1 [Wn {V(n) V(n)}] (3.17)
The following design algorithm emerges for meeting
harmonic requirements:
1. Analyze the original net.rork. Retain branch
voltages and currents of designable elements
(vB(t) and iB(t)).
2. Use the fast Fourier transform to determine
the performance function and adjoint network
excitations in the frequency domain.
3. Use the inverse fast Fourier transform to
form the timedomain adjoint excitations.
4. Analyze the adjoint network in reverse time
(T = t0 + tf t), using the error signal
as excitation. Retain branch voltages and
currents of designable elements ( B(T) and
cB(T)).
5. Calculate the gradient components.
The performance function c and gradient Ve are used by a
function minimization algorithm to adjust the parameter
values. The gradient components are calculated as outlined
in Appendix A. A flow chart of the design algorithm is
shown in Fig. 31.
Initial
network
configuration
1I
Analyze
original
network
Perform
FFT
Form
error terms
and performance
function
Perform
inverse
FFT
Analyze
adjoint
network
Calculate
gradient
components
,v
Conjugate
gradients routine:
obtain new
parameter vector
Fig. 31
Flow of operations for automated design algorithm
dI
I
L
Done S
Stopj
For the design algorithm to be practical, efficient
techniques are mandatory for analysis of the original and
adjoint networks. Numerical integration methods of fixed
order are frequently used for transient analysis. Increased
efficiency is obtained by adjusting order of integration as
required to achieve the largest time step consistent with
allowable error at each time point. This is the basic
premise of Gear's integration algorithm [6,7]. Moreover,
an increase in speed and a decrease in the amount of storage
required is realized by using sparse matrix techniques. In
these methods, only the nonzero elements of the appropriate
matrices are stored and operated upon.
A FORTRAN program has been developed which uses the
'Gear integration method and sparse matrix techniques for
automated harmonic design. A detailed flowchart of this
program appears in Appendix B. In the analysis of the
adjoint network, the gradient components and network
equations are integrated concurrently. This integration is
done in a manner similar to that described by Hachtel et al.
[8], and may be visualized by appending to the network a
parallel combination of a unitvalued capacitor and a
current source. This concept is illustrated in Fig. 32.
The current source values are equal to the appropriate
products of branch voltages and currents so that, upon
integration, the node voltages correspond to the time
integral in the gradient component expressions. For
SR DC DL
in 1 1 1
nt
:k) I1 2 13
Fig. 32
Arrangement of, adjoint network
to obtain gradient components
example:
I = iR) R(T) (3.18)
dv (t)
2 = dt C (T) (3.19)
and
diL(t)
13 d= t L(T) (3.20)
are used to obtain gradient components for a resistor,
capacitor, and inductor, respectively.
The integration is carried out by the "companion model"
method described by Calahan [9]. In this approach, energy
storage elements are replacedat each time pointby a
model which represents the numerical integration formula
applied to the element. This technique allows a timedomain
analysis problem to be reduced to a repetitive dc analysis
problem.
The sparse matrix techniques used are based upon the
methods presented by Gustavson et al. [10]. Executable
machine code is generated for solution of the specific
circuit under consideration. The equations are ordered so
that the solution code which is generated reflects to best
advantage the original sparsity. The procedure utilized
corresponds to the second method of Tinney and Walker [11].
3.2 Design Examples
Two practical examples are presented to demonstrate
the effectiveness of the design algorithm. In order to
speed the simulation, chargestorage effects were neglected
in the transistor models used in the analysis. That is,
frequencies of interest were assumed to be so low that the
static (dc) EbersMoll transistor model was applicable.
This assumption was verified by a trial analysis in which
capacitive elements of the model were included. Comparison
with an analysis of the same network using the static model
revealed a negligible difference in response. Hence, the
original assumption was valid.
In the examples, function minimization was accomplished
by the conjugate gradients method [12]. Transistor parameters
used in the examples are shown in Table 31.
A shuntseries feedback pair amplifier, shown in Fig. 33,
is to be designed. This circuit is a Class A amplifier in
which the feedback and biasing elements are to be adjusted
so that operation in a linear region occurs. A desired gain
of 30 was specified and harmonic distortion was minimized.
Figures 34, 35, and 36 illustrate the progress of the
design. The initial output waveform, spectrum, and adjoint
excitation signal are shown in Fig. 34. Notice the
reflection in the adjoint excitation of the distortion of
the output signal. Figure 35 reveals the result of one
Table 31
Transistor Parameters
BF = 100.
R = 0.1
14
I = I = 1.0 x 10 A.
CS ES
+12v.
R >1K
r!10
600I
+ RR2
0.1 sin(21l000t) 1200
Fig. 33
The shuntseries feedback pair amplifier
of the first example
v0
volts 4
2
0
10
1
0.1
Vjn 10i2
5
10
4
i
amps
.2 .4 .6 .8
0 1 2 3 4 5 6 7 8 9 10
.2 .4
Fig. 34
Initial configuration for example of Fig. 33
(a) Output waveform
(b) Frequencydomain components of output signal
(c) Adjoint network excitation
, ms.
n
t, ms.
0
Fig. 35
Circuit configuration for example of Fig. 33
after one design iteration
(a) Output waveform
(b) Frequencydomain components
of output signal
(c) Adjoint network excitation
12
10
V0o 8
volts 6
4.
2
0
10 
.2 .4
.6 .8 1.
0.1.
Vn 10 2
3
10
4
10
5
10 n
0 1 2 3 4 5 6 7 8 9 10
(b)
0.3
0.2
amps
0.1
0o
t, ms.
.6 .8 1.
ms.
9
Fig. 36
Final circuit configuration for example of Fig. 33
(a) Output waveform
(b) Frequencydomain components
of output signal
(c) Adjoint network excitation
12
10
8
6
4
2
0 t, ms.
.2 .4 .6 .8 1.
(a)
10
1
0.1
102
103
10 3
10
105 I 1 I 1 >8n
0 1 2 3 4 5 6 7 8 9.10
(b)
.05
.04
.03
.02
.01
0 I\ ' 'A t, ms.
.01
.02
.b .0 iL.
design iteration. The output waveform is more nearly
sinusoidal and the reduction in the error between desired
and actual frequency components causes a much smaller
adjoint excitation signal. After four design iterations,
the situation is that shown in Fig. 36. Note the change
in the scale of Fig. 36(c). In circuits such as this,
gain and distortion are conflicting specifications. That
is, for a given input signal, increased gain results in
increased distortion. The design procedure accounts for
this conflict, yielding a design which compromises between
the two conflicting specification.. The availability of
weighting factors allows one or the other of the specifi
cations to be emphasized, if desired. The compromise is
seen by consideration of Table 32. Notice that although
distortion was quite small after the first iteration, the
fundamental component was too small. After the fourth
iteration, both gain and distortion were acceptable.
The largesignal Class B audio power amplifier of
Fig. 37 is to be designed for minimum distortion at a
power output of 20 watts (RMS). Four resistances are
adjusted by the program to achieve the design goals. The
design progress is depicted in Figs. 38, 39, and 310.
The initial output waveform, spectrum, and adjoint network
excitation signal are shown in Fig. 38. In this example,
the weighting factors of the harmonics were increased to
emphasize the requirement of minimum distortion. This
Results of Application of
Table 32
Design Algorithm
to Circuit of Fig. 33
Initial It. 1 It. 2 It. 3 It. 4 Desired
R1 5000 4747 4725 4708 3111
R2 5000 1982 2095 2110 2061 
Dist.(%) 28.2 1.38E2 1.46E2 1.47E2 1.26E2 0
v0 (dc) 5.99 5.71 5.70 5.70 5.51 5.50
v1 (fund) 5.81 2.86 3.01 3.03 2.96 3.00
v2 0.719 1.60E4 1.88E4 1.92E4 1.30E4 0
v3 0.492 6.55E5 7.41E5 7.53E5 7.01E5 0
v4 : 0.145 3.12E5 3.17E5 3.19E5 3.18E5 0
v5 2.08E2 1.82E5 1.89E5 1.87E5 1.87E5 0
V6 0.110 2.96E5 2.13E5 3.14E5 3.07E5 0
v7 4.74E2 2.26E5 2.38E5 2.40E5 2.33E5 0
v8 3.91E2 2.42E5 2.55E5 2.57E5 2.50E5 0
v9 3.94E2 2.35E5 2.47E5 2.49E5 2.43E5 0
vl0 2.26E2 2.07E5 2.18E5 2.20E5 2.13E5 0
Performance 4.45E3 3.30E5 2.04E5 2.02E5 6.92E7 0
function
0
+50v.
R3
R4
1001F=
S_390
20001
R O
)IF 2.7K
!v 10 390 390
'Hz K
100
pF
R1
Fig. 37
The largesignal audio power amplifier
of the second example
.6 .8
0 1 2
3 4 5
(b)
1L .1. I I I _t I 
6 7 8 9 10
Fig. 38
Initial circuit configuration for example of Fig. 37
(a) Output waveform
(b) Frequencydomain components of output signal
(c) Adjoint network excitation
V0
v0
volts
10
20
SVni
t, ms.
n
0.1
2
10
103
104
8
6
4
i 2
amps
 2
 4
6
 8
10
12
t, ms.
V0
volts
volts
0
10
20
10
1
0.1.
102
10
3
10
.6 .8 1.
10 4
0
1
_________________ L
1 2 3 4 5
(b)
7 8 9 10
.6 .8 1.
Fig. 39
Circuit configuration for example of Fig. 37
after one design iteration
(a) Output waveform
(b) Frequencydomain components of output signal
(c) Adjoint network excitation
I n
t, ms.
n
t, ms.
i
amps
77
20
10
volts
0  t, ms.
.2 .4 / .6 .8 1.
10
20 (a)
10
1
S 0.1
102
10 3
4 >I I I I n
0 1 2 3 4 5 6 7 8 9 10
(b)
0.4 4
0.3
0.2
amps
0.1
0 I t, ims.
0.1
0.2
0.3 (c)
Fig. 310
Final circuit configuratich for example of Fig. 37
(a) Output waveform
(b) Frequencydomain components of output signal
(c) Adjoint network excitation
emphasis is reflected in the adjoint excitation signal.
Figure 39 illustrates the results of the first design
iteration. The final situation, after three design
iterations, is depicted in Fig. 310. The final distortion
level is 0.3 per cent: an acceptable level for audio appli
cations. Results appear in Table 33.
3.3 Discussion
A method has been developed which is effective in
designing nonlinear networks for frequencydomain specifi
cations. This method makes use of three powerful concepts:
the adjoint network approach, Gear's numerical integration
algorithm, and the fast Fourier transform. Practical
examples using the FORTRAN program were presented to
demonstrate the effectiveness of the design procedure.
Several points need to be mentioned, however. The
examples were chosen so that no initial transient was
present in the simulation. In general, the simulation will
contain both transient and steadystate components. The
programto be applicable in these casesneeds to be
altered to allow for this situation. The steadystate
determination procedure of Chapter 2 should be incorporated.
A method of determining the transversal terms (see
Appendix A) must also be included.
Table 33
Results of Application of Design Algorithm
to Circuit of Fig. 37
Initial
It. 1
It. 2
Desired
1 50.K 47.7K 44.9K 
R2 130.K 129.K 129.K
3 1.50K 6.23K 5.81K 
R
4 6.80K 1.14K 1.10K 
Dist. (%) 8.08 0.299 0.303 0
V0 (dc) 0.145 5.48E2 5.09E2 0
V1 (fund) 21.3 20.9 20.0 20.0
V
2 0.379 4.18E3 3.67E3 0
V3 0.308 2.29E2 2.24E2 0
V
4 0.296 9.47E4 9.05E4 0
V5 0.267 1.41E2 1.36E2 0
V6 0.194 7.46E4 7.18E4 0
V7 0.124 1.02E2 9.90E3 0
V8 8.61E2 7.09E4 6.85E4 0
V9 5.47E2 8.14E3 7.89E3 0
V10 1.08E2 7.31E4 7.11E4 0
Performance
Function 1.85E2 4.63E4 2.96E5 0
3.4 References
1. F.K. Manasse, J.A. Ekiss, and C.R. Gray, Modern Transistor
Electronics Analysis and Design, Englewood Cliffs, New
Jersey: PrenticeHall, 1967.
2. M. Javid and E. Brenner, Analysis Transmission and
Filtering of Signals, New York: McGrawHill, 1963.
3. E.O. Brigham and R.E. Morrow, "The fast Fourier
transform" IEEE Spectrum, Vol. 4, No. 11, pp. 6370,
December, 1967.
4. G.O. Bergeland, "A guided tour of the fast Fourier
transform," IEEE Spectrum, Vol. 6, No. 7, pp. 4152,
July, 1969.
5. J.W. Cooley and J..W. Tukey, "An algorithm for the
machine calculation of complex Fourier series,"
Mathematics of Computation, Vol. 19, pp. 297301,
April, 1965.
6. C.W. Gear,"Simultaneous numerical solution of
differentialalgebraic equations," IEEE Transactions
on Circuit Theory, Vol. CT18, No. 1, pp. 8994,
January, 1971.
7. C.W. Gear, "The automatic integration of ordinary
differential equations," Communications of the ACM,
Vol. 14, No. 3, pp. 176179, March, 1971.
8. G.D. Hachtel, R.K. Brayton, and F.G. Gustavson, "The
sparse tableau approach to network analysis and design,"
IEEE Transactions on Circuit Theory, Vol. CT18, No. 1,
pp. 101113, January, 1971.
9. D.A. Calahan, "Numerical considerations for implemen
tation of a nonlinear transient circuit analysis
program," IEEE Transactions on Circuit Theory, Vol. CT18,
No. 1, pp. 66.73, January, 1971.
10. F.G. Gustavson, W. Liniger, and R.A. Willoughby, "Symbolic
generation of an optimal Crouc algorithm for sparse
systems of linear equations," Journal of the ACM, Vol. 17,
No. 1, pp. 87109, January, 1970.
11. W.F. Tinney and J.W. Walker, "Direct solutions of sparse
network equations by optimallyordered triangular
factorization," Proceedings of the IEEE, Vol. 55, No. 11,
pp. 18011809, November, 1967.
12. R. Fletcher and C.M. Reeves, "Function minimization by
conjugate gradients," Computer Journal, Vol. 7, No. 2,
pp. 149154, July, 1964.
CHAPTER 4
CONSIDERATIONS FOR OSCILLATOR DESIGN
Automated design of oscillator circuits requires
several special methods in addition to the techniques
developed in the preceding chapters. Automated network
design is based on iterative analysis of the network under
consideration and of the "adjoint" network. In order to
design oscillator circuits, therefore, it is necessary
to be able to determine the frequency of oscillation as
well as the steadystate waveform. Moreover, the steady
state determination algorithm described in Chapter 2
requires knowledge of the period of the waveform. Investi
gation of oscillator design reveals that the unstable
nature of oscillator circuits causes additional difficulties
in the analysis of the adjoint network. For example,
analysis of the oscillator network of Fig. 41 yielded the
results shown in Fig. 42. The adjoint network response
grows without bound. This chapter is concerned with
possible solutions to these problems.
4.1 Determination of Frequency of Oscillation
The frequency of oscillation must be accurately
determined in order to apply the harmonic design methods
i=f(v0)
i=f(v0)
gm =3.5
slope = 1
(1l)
slope = 100
(o .on)
0 r
r(T)
gm= 3.5
(b)
Fig. 41
Oscillator circuit
(a) Original network
(b) Adjoint network
v(t)
t
(a)
o(T)
(T=Tt)
(b)
Fig. 42
Output waveforms of networks of Fig. 41
(a) Original network output voltage v0
(b) Adjoint network output voltage 0
described in Chapter 3. A straightforward approach to this
determination may be based on the phase plane concepts
frequently used in control systems analysis [1,2]. The
phase plane is the plane of the variable under consideration
(for example, output voltage) and its time derivative.
Stable oscillations in this variable give rise to closed
loops in the phase plane referred to as "limit cycles."
One revolution around such a limit cycle represents one
period of the steadystate response. In the following
discussion, we refer to Fig. 43, which is assumed to be a
phase 'plane diagram representative of an oscillator circuit.
If the initial conditions of the circuit correspond to
point A, the "trajectory" is seen to eventually converge on
the 'limit cycle which represents steadystate operation.
One .revolution around the phase plane from point A to
point B may be taken as an estimate of the period of
oscillation, while the distance vA vB may be thought of
as an error which is zero when the circuit is operating in
the .steadystate condition. This error may be used in a
method which adjusts the initial condition to achieve steady
state operation.
One such method is the adjoint networkNewton method
of Chapter 2. This technique makes use of derivative
information to calculate a new initial condition which
reduces the error. Another revolution around the phase
VV
Limit
Cycle
Fig. 43
Phase plane diagram for an oscillator
plane is undertaken using the new initial condition value.
This process is repeated until the error decreases to an
acceptable amount. The period of oscillation is then
known: it is the time required to traverse one revolution
around the limit cycle.,
The application of the phase plane method just
described to a network".anialysis program is quite simple.
A phase plane plot is "constructed" for each capacitor
voltage and its derivative and each inductor current and
its derivative in the circuit. The initial conditions
determine an angle OC for each capacitor and 0L for each
inductor, where
S= ArctanC (t0) (4.1)
and
S= ArctanL i(0) (4.2)
If the initial capacitor voltage vC(t0) is given, the
current iC(t0) is found by replacing the capacitor by a
voltage source of value'.t; (t). A dc analysis is performed
and the current through this voltage source equals i(t0) .
Similarly, an inductor is replaced by a current source of
value iL(t0). A dc analysis is undertaken and the voltage
across this current source is vL(t0). A timedomain analysis
is then performed and the angles 0C and 8L are calculated
for each capacitor and inductor, respectively, at each step
of the analysis. When at least a full revolution has
occurred for each energystorage element, the analysis is
halted. Interpolation may be necessary to obtain this
angle precisely. The time of the estimated period of
oscillation is calculated and one iteration of the steady
state algorithm may then be applied. When periodic
steadystate operation has been achieved, the period is
known.
It should be noted that in !.ome instances, vC and VC
or iL and iL must be used because of the presence of a dc
component. The second derivative may be calculated by
differencing; i.e.,
v (t0+h) Vc(t )
v ~ h (4.3)
or
S L(t0+h) 1L(t0)
iL (4.4)
In some numerical integration methods which use higher
derivatives (such as Gear's algorithm [3]) this differencing
may not be necessary.
4.2 Calculation of Gradient Components for Automated
Oscillator Design
Proper operation of practical oscillator circuits
depends upon two factors: instability and nonlinearity.
The unstable nature of oscillator circuits insures the
growth of oscillations, while the nonlinear property limits
the growth and establishes a steady oscillation. Since the
adjoint network is tq be used in oscillator design, several
facts must be established about the adjoint network of an
oscillator. We demonstrated in Appendix A that nonlinear
elements in the original network become linear time
varying elements in the adjoint network. Thus, the
limiting of the growth of signals may not occur in the
adjoint network whichcorresponds to an oscillator. As
an example, consider the circuit shown in Fig. 44. In
this network, an increasing sinusoid is generated which
grows until excursions of the output reach the positive
segments of the nonlinear resistance characteristic. The
change in resistance limits these excursions and a steady
state condition is achieved. In the adjoint network
corresponding to this circuit, the resistor has a negative
value during the time intervals corresponding to time
segments in the original network analysis during which the
slope of the nonlinear resistance was negative. The signal
grows during this time interval. When the slope of the
resistance in the, ori inai>lnetwork was positive, the
adjoint network resistance also has a positive value.
During this time, the signal will decay. Situations thus
exist in which the effect of the growth is greater than
that of the decay.
