Helsinki University of Technology Department of Electrical and Telecommunications Engineering HOMOTOPY METHODS IN DC CIRCUIT ANALYSIS

Helsinki University of Technology Department of Electrical and Telecommunications Engineering Vesa Linja-aho HOMOTOPY METHODS IN DC CIRCUIT ANALYSIS...
Author: Gavin Cook
0 downloads 1 Views 266KB Size
Helsinki University of Technology Department of Electrical and Telecommunications Engineering

Vesa Linja-aho

HOMOTOPY METHODS IN DC CIRCUIT ANALYSIS

Thesis submitted for examination for the degree of Master of Science in Technology Espoo 18.9.2006

Thesis supervisor: Prof. Martti Valtonen Thesis instructor: Lic.Sc. (Tech.) Mikko Honkala

Helsinki University of Technology

Abstract of master’s thesis

Author: Vesa Linja-aho Title: Homotopy Methods in DC Circuit Analysis Date: 18.9.2006

Language: English

Number of pages: 37

Department: Electrical and telecommunications engineering Professorship: Circuit theory

Code: S-55

Supervisor: Prof. Martti Valtonen Instructor: Lic.Sc. (Tech.) Mikko Honkala In this thesis, a comprehensive literature study and simulation-based comparison for using homotopy methods in DC circuit analysis was made. Homotopy methods are a subclass of continuation methods, which have been successfully used to solve nonlinear sets of equations. Fixed-point, variable stimulus, and modified variable stimulus homotopy functions with plain Newton-Raphson, ordinary differential equation-based, and pseudo-arclength-based solver algorithms for finding DC operating points for electric circuits were tested. Plain Newton-Raphson and pseudo-arclength methods were implemented in MATLAB and tested with the MATLAB-APLAC interface. The simulation results are discussed.

Keywords: Circuit simulation, homotopy methods, DC analysis

¨n Diplomityo ¨ tiivistelma

Teknillinen korkeakoulu

Tekij¨a: Vesa Linja-aho Ty¨on nimi: Homotopiamenetelm¨at virtapiirien DC-analyysiss¨a P¨aiv¨am¨a¨ar¨a: 18.9.2006

Kieli: Englanti

Sivum¨a¨ar¨a: 37

Osasto: S¨ahk¨o- ja tietoliikennetekniikka Professuuri: Circuit theory

Koodi: S-55

Valvoja: Prof. Martti Valtonen Ohjaaja: TkL Mikko Honkala T¨ass¨a diplomity¨oss¨a suoritettiin kattava kirjallisuustutkimus ja simulaatiovertailu homotopiamenetelmist¨a tasavirtapiirien analyysiss¨a. Homotopiamenetelm¨at ovat kontinuaatiomenetelmien alaluokka. Homotopiamenetelmi¨a on menestyksellisesti k¨aytetty ep¨alineaaristen yht¨al¨oryhmien ratkaisemiseen. Kiintopistehomotopiafunktio, muuttuvan her¨atteen homotopiafunktio ja modifioitu muuttuvan her¨atteen homotopiafunktio testattiin virtapiiriyht¨al¨oiden ratkaisemisessa. Ratkaisualgoritmeina k¨aytettiin Newton-Raphson-kontinuaatiota, pseudokaarenpituusalgoritmia ja differentiaaliyht¨al¨oihin pohjautuvaa menetelm¨a¨a. Newton-Raphson-kontinuaatio ja pseudokaarenpituusalgoritmi implementoitiin MATLABiin ja testattiin MATLAB-APLAC-rajapinnalla. Simulointituloksia tarkasteltiin.

Avainsanat: Piirisimulointi, homotopiamenetelm¨at, DC-analyysi

Preface ”I devoted myself to study and to explore by wisdom all that is done under heaven. What a heavy burden God has laid on men!” - The Holy Bible, NIV, Eccl. 1:13

I wish to thank professor Martti Valtonen and my instructor Mikko Honkala for patiently supporting me with this endless-seeming project. I deeply appreciate support from all my colleagues at the Circuit theory laboratory, especially Luis Costa for helping me with the LATEX typesetting program. The great work atmosphere, especially the fun discussions at the coffee table, helped me withstand the moments when this project felt frustrating. Financial support received through MOSAICS project is acknowledged. This work could not have been finished without support from my lovely wife Minttu — thank you for your encouragement and help. I dedicate this thesis to my dear parents Hannele, who died this year, and Martti, who died last year. You always supported me in my studies — so sad you did not see me graduate. And finally, to anyone who is starting his or her master’s thesis and reads this: choose to do your thesis as full-time employee, not part-time as I did. The other work assignments will wait you, no matter how well-paid and interesting they are.

Otaniemi, September 18, 2006 Vesa Linja-aho iv

Contents Abstract

ii

Abstract (in Finnish)

iii

Preface

iv

Contents

v

Symbols and notational conventions

vii

1 Introduction

1

2 Background on DC Circuit Analysis

1

2.1

Circuit Simulation . . . . . . . . . . . . . . . . . . . . . . . . .

1

2.2

APLAC Circuit Simulation and Design Tool . . . . . . . . . . .

2

2.3

DC Analysis and Newton-Raphson Method . . . . . . . . . . . .

2

2.4

Convergence-Aiding Techniques . . . . . . . . . . . . . . . . . .

2

2.4.1

Model Damping . . . . . . . . . . . . . . . . . . . . . . .

3

2.4.2

Source Stepping . . . . . . . . . . . . . . . . . . . . . . .

3

2.4.3

Voltage Limiting . . . . . . . . . . . . . . . . . . . . . .

4

2.4.4

Diode Damping . . . . . . . . . . . . . . . . . . . . . . .

4

2.4.5

Norm Reduction . . . . . . . . . . . . . . . . . . . . . .

4

2.4.6

Temperature Sweeping . . . . . . . . . . . . . . . . . . .

5

3 Homotopy Methods

6

3.1

Bifurcation and Folds . . . . . . . . . . . . . . . . . . . . . . . .

7

3.2

Homotopy Function . . . . . . . . . . . . . . . . . . . . . . . . .

9

3.2.1

Fixed-Point Homotopy Function . . . . . . . . . . . . . .

9

3.2.2

Variable Stimulus Homotopy Function . . . . . . . . . . 10

3.2.3

Modified Variable Stimulus Homotopy Function . . . . . 11

3.2.4

Combinations and Other Homotopy Functions . . . . . . 11

3.3

Solving the Homotopy Equations . . . . . . . . . . . . . . . . . 12 v

3.4

3.3.1

Standard Newton–Raphson Solver . . . . . . . . . . . . . 12

3.3.2

ODE-based Solver

3.3.3

Pseudo-arclength Solver . . . . . . . . . . . . . . . . . . 13

3.3.4

HOMPACK – A Suite of Solver Algorithms . . . . . . . 15

. . . . . . . . . . . . . . . . . . . . . 13

Practical Viewpoints . . . . . . . . . . . . . . . . . . . . . . . . 17 3.4.1

Selecting Parameters . . . . . . . . . . . . . . . . . . . . 17

3.4.2

Ensuring the Convergence . . . . . . . . . . . . . . . . . 17

4 MATLAB Tests 4.1

19

Experiment setup . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.1.1

MATLAB–APLAC Interface . . . . . . . . . . . . . . . . 19

4.1.2

Implementation and Benchmark Circuits . . . . . . . . . 19

4.2

A Simple Simulation Example . . . . . . . . . . . . . . . . . . . 20

4.3

Standard Newton-Raphson Solver . . . . . . . . . . . . . . . . . 22

4.4

4.3.1

The Fixed-Point Homotopy Function . . . . . . . . . . . 23

4.3.2

The Variable Stimulus Homotopy Function . . . . . . . . 24

4.3.3

The Modified Variable Stimulus Homotopy Function . . 25

4.3.4

Test for NR-Continuation with Equal Parameters . . . . 25

ODE-Based Curve Tracing . . . . . . . . . . . . . . . . . . . . . 25 4.4.1

The Fixed-Point Homotopy . . . . . . . . . . . . . . . . 25

4.5

Pseudo-Arclength Solver with the Planar Corrector . . . . . . . 27

4.6

Pseudo-Arclength Solver with the Spherical Algorithm . . . . . 27

5 Discussion

29

6 Conclusion

31

References

32

Appendix A

35

Appendix B

36

vi

Symbols and notational conventions Symbolic Notations x x d d∂ t

scalar node voltage vector derivative with respect to t partial derivative with respect to t ∂t f (x) set of nodal equations for a circuit J Jacobian matrix Gmin the minimum conductance in model damping technique VT thermal voltage of a semiconductor diode IS saturation current of a semiconductor diode ǫ current imbalance in norm reduction technique λ homotopy parameter H(x, λ) homotopy function a initial guess vector for homotopy function G conductance scaling diagonal matrix for homotopy function u unit vector in λ − x hyperplane ∆s step length in pseudo-arclength methods

Abbreviations AC APLAC BJT DC JFET MATLAB MOSFET NR ODE SPICE

alternating current circuit simulator and design tool (originally Analysis Program for Linear Active Circuits) bipolar junction transistor direct current junction field effect transistor Matrix Laboratory (a computer program) metal-oxide semiconductor field effect transistor Newton–Rhapson ordinary differential equation Simulation Program with Integrated Circuit Emphasis

vii

1

Introduction

In computer-based circuit simulation, finding the DC operating point of the circuit is usually the first step of all other analyses. Therefore, if the simulator cannot find the operating point, also all other simulations are impossible to run. Determining the DC operating point requires solving a set of nonlinear equations. The most usual methods are based on Newton–Raphson iteration with a bunch of convergence-aiding techniques. Strongly nonlinear circuits and circuits which have multiple DC operating points are very difficult for conventional solver algorithms. One way to improve convergence is to use so called homotopy methods. Homotopy methods can help to find the DC operating point for a troublesome circuit. For this thesis, homotopy methods have been studied and tested with troublesome circuits. The practical testing of homotopy algorithms was carried out using APLAC circuit simulation and design tool [1, 2] and MATLAB software. Before doing this, a comprehensive literature study on homotopy methods was carried out. The main goal of the research was testing the homotopy methods with practical circuits. The practical implementation of a homotopy method consists of two parts: deciding which homotopy equation to use, and solving the equation.

2 2.1

Background on DC Circuit Analysis Circuit Simulation

Circuit simulation is an important part of electronic circuit design at present. By simulating the operation of the circuit beforehand, the need for expensive prototypes is avoided. Using circuit simulation, the design process is accelerated and designers can test circuits which can be practically impossible to test in real world [3]. For example, it can be simulated how a certain improvement in characteristics of semiconductor component will affect the performance of a device. The simulation will help in deciding how much effort is reasonable to take in improving the certain characteristic. Circuit simulation can be used to optimize the yield of the manufacturing process, thus reducing the manufacturing costs. Dealing with component value tolerances is practically impossible without circuit simulation, especially when designing large circuits.

1

2.2

APLAC Circuit Simulation and Design Tool

APLAC is a flexible circuit simulator suite originally created by professor Martti Valtonen. Nowadays, the development is made under the AWR-APLAC Corporation.

2.3

DC Analysis and Newton-Raphson Method

The problem of finding the DC operating point can be described as a set of nonlinear equations f (x) = 0 (1) where x is a node voltage vector. The basic method for numerically solving the circuit equations is the Newton– Raphson method [3, ch. 4.2]. The equation 1 is solved iteratively with the formula x(j+1) = x(j) − [J (xj )]−1 f (x(j) ),

(2)

where J represents the Jacobian matrix of the function f (x) and j is the running index for the iteration count [6, ch. 5.4.3]. The Newton-Raphson method can be guaranteed to be globally convergent for single-BJT circuits [16]. In certain conditions described in [16], the convergence can be guaranteed for most of the practical BJT-circuits. The DC analysis in APLAC circuit simulator is carried out by solving modified nodal equations of the circuit. The set of modified nodal equations is formed via the gyrator transform. [2]

2.4

Convergence-Aiding Techniques

In DC analysis, Newton–Raphson algorithm converges well if the initial guess for the operating point is good. For a computer program, it is very difficult to make a good guess for the operating point. Additionally, almost all practical circuits are strongly nonlinear. An ordinary semiconductor diode, for example, has a very steep exponential characteristic function. Because the convergence of the plain Newton–Raphson algorithm cannot be guaranteed for all practical circuits, a number of convergence-aiding methods have been developed. Model damping, source stepping, voltage limiting, diode damping, norm reduction, piecewise-linear analysis [4] and dogleg method [5] are in use in the APLAC circuit simulator. 2

No method guarantees convergence alone. APLAC uses different combinations of convergence-aiding techniques to find the operating point. Usually the source stepping method is the last resort when all the other methods fail [2]. Some of the convergence aiding techniques are introduced in the following sections. They are described in [2]. One of the methods, the source stepping method, is a so called continuation method or homotopy method. Homotopy methods have been efficiently used to solve sets of nonlinear equations. 2.4.1

Model Damping

In model damping, every nonlinear static source in the circuit is shunted with a conductance having the value G = 10m Gmin,

(3)

where Gmin is the minimum conductance in the circuit and m is an integer. Conductance G is also added between every node and ground.1 The default values for Gmin and m are 10−12 and 6, respectively. User can alter both of the values. The iteration of the DC operating point is started with the m given, and after the iteration converges, m is reduced by one. This is repeated, until m = 0. 2.4.2

Source Stepping

In the source stepping method, all the voltage sources in the circuit are multiplied by a constant k ∈ [0, 1]. When making a DC analysis, the iteration is started with k = 1 (the original circuit). If convergence is not achieved, the constant k is decreased and a new iteration is attempted. When the iteration converges, k is increased gradually, until k = 1. The operating point found with the previous convergent iteration is used as an initial guess for the next iteration. Source stepping method has an easy physical interpretation — the circuit is “powered up” gradually by increasing the source voltages. Source stepping is a very efficient method when dealing with strongly nonlinear circuits. 1

A popular convergence aiding technique which consists only of adding small conductances between all nodes and ground is called the Gmin stepping method [10].

3

2.4.3

Voltage Limiting

Voltage limiting is a simple method for preventing semiconductor components from switching on and off repeatedly during the iteration. The user can define a constant Umax , which describes the maximum allowed voltage change between iteration cycles. If the voltage change between two iteration cycles exceeds Umax , all the voltages are multiplied with a constant β so that the maximum voltage change is Umax : where

ˆ n+1 = un + β(un+1 − un ) u

(4)

Umax . max |un+1 − un | If none of the voltage changes exceeds Umax , β = 1. β=

2.4.4

(5)

Diode Damping

Diode damping is a convergence improvement method for dealing with the steep exponential characteristic curve of diodes. It is closely related to the voltage limiting method. Every diode monitors its voltage, and the damping factor is calculated ! VT β = 1, if un+1 ≤ VT ln √ (6) 2IS and ! −un + 1) VT ln( un+1 VT VT , if un+1 > VT ln √ , (7) β= un+1 − un 2IS where VT and IS are the thermal voltage and the saturation current of the diode, respectively. The new voltage is then calculated from the equation (4). 2.4.5

Norm Reduction

Norm reduction is practically equivalent to line search optimization. On each iteration cycle, the current imbalance ε between the feeding and loading circuits is monitored. If the current imbalance ε=

max I



∂i − i(u) − u ∂u u=un

(8)

increases, a damping factor β = 0.5 is applied. If the norm still increases, parabolic fitting is used in order to find a damping factor which reduces the error. The procedure is repeated until an error-reducing step length is found or the maximum error criteria is satisfied. 4

2.4.6

Temperature Sweeping

In the temperature-sweeping [8] method the temperature of the circuit is swept from a certain value (usually zero) to the final value, which is the operating temperature of the circuit. The method requires that the component models work in the whole temperature range. For example, if a model is verified only to work at temperatures near room temperature, the temperature sweep has no purpose. The temperature sweeping method is not implemented in APLAC simulator.

5

3

Homotopy Methods

“A continuous transformation from one function to another” is called a homotopy [7]. For this thesis, a number of scientific publications concerning homotopy methods were examined. Homotopy methods are a subset of continuation methods. In a continuation method, a solution of an algebraic nonlinear problem is treated as an instant of a dynamic problem [3]. For instance, the source stepping method is a continuation method. Sweeping the source voltages from small values to the final value is like powering up the circuit. In a homotopy method, an auxiliary equation is constructed for solving the original problem. If the original problem is f (x) = 0,

(9)

one example of a homotopy equation for the problem is H(x, λ) = λf (x) + (1 − λ)g(x) = 0.

(10)

The idea of the homotopy method is the following: an easy function is chosen as g(x). After that, the parameter λ is swept from 0 to 1, so that the easy problem is continuously deformed into the original problem. The homotopy function H(x, λ) may be constructed and solved in numerous ways. The concept of using a homotopy method for solving a set of nonlinear equations consists of 1. the formation of the homotopy function 2. choosing an efficient solver 3. making a working implementation of them. If one of the three elements has problems, the whole homotopy method may not work (Figure 1). For normal circuit equations, the solutions for the homotopy function when λ is swept form a continuous path in λ–x –plane. For example in [9], the homotopy method is defined to be “the class of continuation methods where the homotopy parameter does not necessarily vary monotonically as the path from the original solution to the final solution is followed”. According to the definition, the source stepping, for example, is not a homotopy method. Homotopy methods can be divided into two groups. In natural parameter homotopy methods, the homotopy parameter is a parameter with an obvious physical interpretation. That is, for example, if we multiply the supply voltage 6

Figure 1: The ingredients of a well-working homotopy method. of the circuit by λ [12]. The opposite of the natural parameter homotopy methods are artificial parameter homotopy methods. The homotopy in Equation 10 is a good example of an artificial parameter homotopy method. Homotopy methods have been used successfully for improving convergence in DC circuit analysis [11, 14]. The main idea of homotopy methods is to convert the set of circuit equations into a simple equation, and convert it gradually into the original equation set. The homotopy methods are less sensitive to the initial guess than the bare Newton–Raphson method. Some homotopy methods can, with certain conditions, even be globally convergent [8]. The main drawback of homotopy methods is their computational intensity [11]. For example, curve–tracing methods which use predictor and corrector need a lot of computing power. Other difficulties with practical use of homotopy methods pertain to the implementation of homotopy methods [13]. The algorithm for solving the homotopy equations has to be carefully implemented. Solving the DC operating point with a homotopy method is usually done by tracing the zero curve on the λ–x -plane. It is also possible to just step the embedded parameter forward. Numerical problems can be met if the zero curve has sharp turns, or the homotopy curve bifurcates.

3.1

Bifurcation and Folds

The main purpose of homotopy methods is to find the DC operating point for a circuit, when all conventional methods fail. Homotopy methods can guarantee convergence, if implemented properly. The two major problems homotopy solver algorithm may face are bifurcation 7

and folds in the homotopy curve. Bifurcation means that the homotopy curve in λ–x -plane splits into two or more curves, as in Figure 2. 1

0.9

0.8

0.7 Bifurcation point

x

0.6

0.5

0.4

0.3

0.2

0.1

0

0

0.2

0.4

0.6 λ

0.8

1

1.2

Figure 2: An example of bifurcation. The bifurcation points can be troublesome for homotopy solvers. The solver can diverge especially if the new branches of the curve diverge steeply from each other. Folds in the homotopy curve may cause convergence problems too. If the curve folds back steeply, as in Figure 3 while λ ≈ 10.5, the curve-tracing algorithm may get stuck in the steep turn. If a monotonic continuation, where λ always increases, is used, folds cause the algorithm to fail. The reason for this is that if the λ always increases, the numerical algorithm cannot follow the curve to the final solution. In Figure 3, when λ is stepped third time, the algorithm seeks the solution near the rightmost blue dot. Because λ cannot be decreased and the correct solution for that value of lambda is too far away, the algorithm will fail. The source stepping, for example, can suffer from this problem. The problems caused by folds can be dealt with the so called curve tracing methods. Instead of just stepping the parameter λ forward, the zero curve of the homotopy function is followed by a numerical algorithm. The simplest way to do this is to take a step on the curve to the direction of the tangent vector of the curve, and then use a corrector algorithm to get back on the original curve. 8

Very steep folds may cause problems even for curve tracing methods. One way for avoiding steep folds is to use complex parameter homotopies [21].

Figure 3: Monotonic continuation cannot handle folds.

3.2

Homotopy Function

In this thesis, the fixed-point, variable stimulus and modified variable stimulus homotopy functions described in e.g., [8] are reviewed and tested with different solvers. 3.2.1

Fixed-Point Homotopy Function

In a fixed-point homotopy, the original set of nodal functions f (x)

(11)

is transformed into homotopy function H(x, λ) = (1 − λ)G(x − a) + λf (x), 9

(12)

where x is the node voltage vector, λ ∈ [0, 1] is the homotopy parameter and a is a randomly chosen vector. G ∈ Rn ×Rn is an embedded conductance scaling parameter. When the homotopy parameter λ = 0, the function represents a simple circuit having node voltages described by a. When λ = 1, the original function is obtained. By varying the parameter λ from 0 to 1, the simple circuit is gradually transformed into the original (difficult) circuit. The fixed-point homotopy has a simple physical representation. All the nodes in the circuit are connected to ground via a branch, which consists of a conλ Gk in series with a voltage source having voltage ak , where k is ductance 1−λ the index of the node. The fixed-point homotopy can be considered as an improved version of the Gmin stepping method. The voltage sources in series with the conductances guarantee that the homotopy does not bifurcate. In a circuit with two or more operating points, the vector a defines which of the operating points is found. Bifurcation is described more accurately in Section 3.1. Homotopy curve may also fold back, which results in finding multiple operating points, if solver can handle a folding curve. Because of the simple physical interpretation, the fixed-point homotopy is relatively easy to implement into a circuit simulator. Tampering of the existing models is not required. 3.2.2

Variable Stimulus Homotopy Function

Variable stimulus homotopy function is H(x, λ) = (1 − λ)G(x − a) + f (x, λ),

(13)

where the function f (x) is modified so that the node voltages of nonlinear elements are multiplied by the factor λ. Thus, the solution of the function at the starting point λ = 0 is the solution of a linear circuit. Then the circuit is gradually converted to a more and more nonlinear circuit. One disadvantage of the variable stimulus homotopy is that the simulator core and/or component models have to be tampered in order to multiply the voltages of nonlinear circuit elements. According to [8], one advantage of variable stimulus homotopy is a remarkable improvement of numerical stability when dealing with p-n junctions of transistors and diodes.

10

3.2.3

Modified Variable Stimulus Homotopy Function

The modified variable stimulus homotopy is closely related to the variable stimulus homotopy. The equation is slightly different H(x, λ) = (1 − λ)G(x − a) + f (λx).

(14)

The main advantage of the modified variable stimulus homotopy is its simplicity. The circuit element models need no modifications, because all of the node voltages are multiplied with the same factor λ. The modified variable-stimulus homotopy has been tested earlier [9] with ADVICE–HOMPACK -interface. ADVICE is a circuit simulator and HOMPACK is a suite of homotopy methods on Fortran language. The results were promising. Two circuits which did not converge with the ordinary simulator, converged successfully with HOMPACK homotopy solver. For more detailed description of HOMPACK, see Chapter 3.3.4. 3.2.4

Combinations and Other Homotopy Functions

The homotopy function can be constructed in practically infinite number of ways. However, certain simplicity should be preserved. In addition to the three homotopies presented, there are several others to mention. In the variable gain homotopy the transistor forward current gains are multiplied with the homotopy parameter H(x, λ) = (1 − λ)G(x − a) + f (x, λα),

(15)

where the α is a column vector consisting of the forward current gains of the transistors. According to [8], the variable gain homotopy is the fastest converging homotopy for bipolar circuits. The implementation of the variable gain homotopy requires tampering the bipolar transistor models, which makes it harder to test. In the hybrid homotopy introduced in [8], the starting point for variable gain homotopy is calculated with the variable stimulus homotopy. In the same paper, an update homotopy, which is useful in determining how the DC operating point varies over a range of circuit parameters (such as temperature or some component values) H(x, λ) = S(λ)G(x − a) + λf b (x) + (1 − λ)f a (x)

(16)

is introduced. The function S(λ) is an any smooth function with S(0) = S(1) = 0, f a (x) is the nodal function for circuit just being calculated (for 11

example at temperature of 10 ◦ C), and f b (x) describes the circuit that should be solved with a new parameter value (for example at temperature of 15 ◦ C). In [14], a two-parameter homotopy algorithm called BLHOM for large-scale MOSFET circuits, is presented. The method utilizes two parameters, λ1 and λ2 , which control the drain–source characteristics and transfer characteristics, respectively. According to the publication, the homotopy converges on virtually all practical circuits and has led to significant savings in the AT&T company. Two interesting categories of homotopy methods to mention are complex parameter homotopy methods and piecewise-linear homotopy methods. According to [21], complex parameter homotopy methods are very effective when dealing with folds in the homotopy curve. Piecewise-linear homotopy methods [23] are considered less efficient than conventional predictor-corrector methods, but require no smoothness of the equations to be solved. Homotopy methods have also been successfully used in fault diagnosis of piecewise-linear circuits [20].

3.3

Solving the Homotopy Equations

A robust solver is a necessity in order to utilize homotopy methods in convergence aiding. For this thesis, three different solvers were tested: a simple Newton-Raphson continuation, an ordinary differential equation (ODE) -based solver, and pseudo-arclength continuation (with both spherical and planar corrector). 3.3.1

Standard Newton–Raphson Solver

The simplest way to utilize the homotopy equations is just stepping the homotopy parameter, λ, forward and using the standard Newton-Raphson iteration. Starting from λ = 0, the homotopy parameter is increased a bit and the equation is solved with Newton-Raphson iteration. This procedure is repeated until λ = 1 and therefore the original problem solved. The result of previous iteration is used as an initial guess for next iteration cycle. The algorithm goes as follows: 1. Start with λ = 0 and solve H(x, λ) = 0 2. Step λ forward. 3. Solve H(x, λ) = 0 with Newton–Raphson algorithm 4. If λ < 1, go to step 2. 5. Return the final solution 12

3.3.2

ODE-based Solver

The ODE-based method is an algorithm for tracing the zero curve of the homotopy function. In the method, the homotopy function is parameterized with respect to the arc length of the zero curve. Then, a differential equation generating that trajectory is derived and solved. In the ODE-based approach, the zero curve of the homotopy function is treated as a solution of an ordinary differential equation. The ODE is formed by parameterizing both x and λ by the arc length s H(λ, x) = 0

(17)

H(λ(s), x(s)) = 0.

(18)

The problem can be treated as the solution of differential equation d H(λ(s), x(s)) = 0. ds

(19)

From the chain rule, dλ/ds DH(λ(s), x(s)) dx/ds

!

=0

(20)

where DH(λ(s), x(s))

(21)

is the Jacobian of the homotopy function. With initial conditions λ(0) = 0,

x(0) = a,



dλ dx , ds ds

!

=1

(22)

2

the differential equation (19) can be solved. The more detailed description is found in [30] and [11]. 3.3.3

Pseudo-arclength Solver

The pseudo-arclength method described in [19] is based on tracing the zero curve of the homotopy function H(x, λ). Because of working with the actual length of the zero curve is difficult, a so called pseudo-arclength method is used. The basic idea of the pseudo-arclength method is to calculate the tangent vector for the zero curve of H(x, λ), take a step to the direction of the tangent and then use a corrector iteration to get back on the zero curve (Figure 4). The pseudo-arclength method belongs to the class of so called predictor-corrector algorithms — the curve is followed by taking a predictor step towards the tangent, and then a corrector algorithm will be executed to get back on the zero curve. 13

Figure 4: A predictor-corrector method. Consider the zero curve of H(x, λ) in the λ−x hyperplane. Let the unit vector be u = (ux, uλ ).

(23)

The tangent vector can be computed by solving a system of equations     

dH(x,λ) ux dx

+

dH(x,λ) uλ dλ

kuxk2 +

u2λ

= 0 .

(24)

= 1

The solution of equation 24 is    

where

  

uλ = ± √

1 1+kvk22

,

(25)

ux = u λ v

v=

−dH(x,λ) dλ . dH(x,λ) dx

(26)

When the unit vectors ux and uλ have been calculated, the predictor step is taken ( xi+1 = xi + ux∆s , (27) λi+1 = λi + uλ ∆s

14

where ∆s is the chosen step length. Because equation (25) has two possible solutions (the curve has two tangent vectors pointing in opposite directions), a test must be made in order to prevent the algorithm reverse the direction on the curve. A simple way to do this is to evaluate, which of the selections takes us further from the previous point on the curve (if the wrong sign is selected, the new predicted point is very close to the previous predicted point). Then, a corrector algorithm (Newton-Raphson) is executed in order to get back on the zero curve. In order to get back on the solution curve, the corrector algorithm needs one extra equation added to the system of original equations. The original system consists of n modified nodal equations and n unknown voltages, and we have added the homotopy parameter λ. To invoke the corrector algorithm, one extra equation g(x, λ) is added. The most common choice is to use equation g(x, λ) = u · ((xi+1 , λi+1 ) − (xi , λi )) − ∆s = 0

(28)

which means that the next point lies on a hyperplane perpendicular to the unit tangent vector u. Using a hyperplane as the additional equation may cause problems if the zero curve of H(x, λ) has very steep turns on it. For example, if the curve makes a small-radius 180 degrees turn, the hyperplane doesn’t intersect the zero curve at all and the corrector will fail. Other choice is to use a hyperball instead of a hyperplane [22]. A ball will always intersect the zero curve at two points. The equation for the hyperball condition is g(x, λ) = (x1 − c1 )2 + (x2 − c2 )2 + . . . + (xn − cn )2 + (λ − λn )2 − r 2 = 0, (29) where cn corresponds to the previous point on the zero curve. The hyperball will always intersect the zero curve in two points, regardless of the shape of the curve (Figure 5). After selecting the additional equation, the system of n+1 unknowns and n+1 equations is to be solved, to get back to the zero curve. The most common algorithm is the Newton-Raphson, introduced in the Chapter 2.3: "

3.3.4

dH(x,λ) dx dg(x,λ) dx

dH(x,λ) dλ dg(x,λ) dλ

#"

xj+1 − xj λj+1 − λj

#

=−

"

H(x, λ) g(x, λ)

#

.

(30)

HOMPACK – A Suite of Solver Algorithms

HOMPACK [29] is a suite of program codes in FORTRAN language. In HOMPACK, three different methods for solving systems of nonlinear equations 15

Figure 5: The hyperball condition for corrector. are implemented: ODE-based, normal flow and augmented Jacobian matrix method. [30]. The HOMPACK algorithms have been successfully used to find DC operating points in a circuit simulator [34]. HOMPACK is designed for artificial parameter homotopies and it cannot handle homotopy curves which have bifurcations. The ODE-based solver utilizes the same technique introduced in Chapter 3.3.2 – the zero curve of the homotopy function is treated as a solution of an ordinary differential equation. The normal flow algorithm is based on a concept called the Davidenko flow: when parameter a of the homotopy function is varied, the zero curve of the homotopy function H(x, λ) also varies. This family of zero curves is known as Davidenko flow. In the normal flow algorithm, the corrector iterates back on the zero curve along flow, which is normal to the Davidenko flow. The principle of normal flow method is illustrated in [31]. The augmented Jacobian algorithm uses a Hermite cubic predictor and a quasiNewton corrector, and the tangent vector of the zero curve is augmented into the Jacobian matrix of the homotopy map. The augmented Jacobian algorithm is a modified version of an algorithm known as PITCON [32].

16

3.4 3.4.1

Practical Viewpoints Selecting Parameters

For the fixed-point homotopy, the variable stimulus homotopy and the modified variable stimulus homotopy, the formation of homotopy function has two circuit-specific parameters: the initial voltage vector a and the conductance scaling matrix G. Determining values for these parameters has a significant effect in ensuring fast convergence. The randomness of a also guarantees a bifurcation-free path. It will also prevent finding the non-stable operating point of the circuit. One difficulty with handling the homotopy function is the choice of the impedance scaling matrix G. It may be hard to automatically calculate reasonable values. In [9], the matrix G was chosen to be a diagonal matrix with constant entries. The value for elements should correlate with the average impedance level of the circuit. For example, in a circuit dealing with power electronics and low impedances like less than ten ohms, adding conductances of 1 mS may not linearize the circuit and aid convergence adequately. Because it is hard to define good separate values for diagonal elements of G, choosing all of them to be the same is a rational choice. In the simulation examples, the value of the diagonal elements is denoted with G. The choice of vector a has also an impact on convergence – it defines the starting point of the homotopy path. For fixed–point, variable stimulus, and modified variable stimulus homotopies, the solution for H(x, λ) = 0 for λ = 0 is x = a. Selection of the parameter defines also, which operating point is found, if the circuit has multiple operating points. If a curve tracing method which can find multiple operating points is used, the vector a will determine the starting point for the homotopy curve and affect the shape of the homotopy curve. Therefore, a determines, which of the operating points are found. For example, if the circuit has 11 operating points, three of them can be found with a certain value of a and two of them with a different value of a. How the value of a affects, depends on the circuit topology. 3.4.2

Ensuring the Convergence

The main advantage of homotopy methods is their robustness. The fixedpoint, variable stimulus, and modified variable stimulus homotopies can be guaranteed to be globally convergent with certain conditions described in [15]. The homotopy method has to be implemented carefully in order to ensure convergence. In [15] it has been proven that the three homotopy functions studied in this thesis are globally convergent for a circuit satisfying passivity and nogain properties. Almost all practical circuits fulfill these requirements. All 17

circuits consisting of non-negative resistors, diodes, tunnel diodes, thyristors, BJTs, MOSFETs and JFETs are passive and possess the no-gain property [15].

18

4

MATLAB Tests

Some of the homotopy methods reviewed were tested on actual circuits with the MATLAB–APLAC interface. The fixed-point homotopy function, variable stimulus homotopy function, and modified variable stimulus function were tested with a Newton–Raphson solver. Additionally, the fixed-point homotopy function was tested with a pseudo-arclength (both with planar and spherical corrector) and the ODE-based algorithm.

4.1 4.1.1

Experiment setup MATLAB–APLAC Interface

The tests for the homotopy methods were run on a specific MATLAB–APLAC –interface developed in the Circuit theory laboratory of the Helsinki University of Technology . The interface allows for using the MATLAB language to test numerical methods with circuits described in the APLAC input language. The interface takes node voltages as an input, and returns the values of the circuit equations as well as the Jacobian matrix. Programming and testing numerical methods with MATLAB is more convenient in comparison with Clanguage programming, although the interface sets certain limitations for the maximum size of the benchmark circuit. The interface also slows down the computing. The availability of source code of the simulator is an essential thing in testing of the homotopy methods. For example it is also possible to import circuit description files of the popular SPICE circuit simulator, which makes it possible to run tests with different kind of circuits. Benchmarking circuit description files for SPICE are widely available. 4.1.2

Implementation and Benchmark Circuits

The algorithms for testing a simple Newton-Raphson and pseudo-arclength method were implemented. The source codes are presented in appendices A and B, respectively. Tests for ODE-based method were made with a ready MATLAB source code by professor Heath Hoffmann [33]. The source code was modified to use sparse matrices. The solver utilizes a variable step 4th order predictor-corrector method to trace the solution curve. Combinations of homotopy functions and solvers were tested with 12 circuits. The circuits used for testing and benchmarking the methods are summarized in Table 1. The circuit schmitt.i is a simple circuit described in [15] and the circuit ref6.i is a circuit with 11 DC operating points, described in [17]. 19

Table 1: Summary of benchmark circuits Name palikka.i simdc.i toronto.i voter.i flipflop.i globtest.i dcprob.i ssjudc.i schmitt.i ref6.i dc demo.i dcprob2.i

Nodes 2 3 26 1710 5 113 475 653 6 17 129 487

Diode 1 55 3 3 -

BJT 1 2 17 104 151 2 6 17 104

MOS 58 1710 42 41 41

OpAmp 4 4

The following simulation sets were made: • A simple simulation example to demonstrate the effect of initial guess. • A test set for simple NR continuation to qualitatively test their capabilities • A test set for simple NR continuation with fixed set of parameters • A test set with Hoffmann’s ODE solver • A test set with pseudo arclength solver, with planar and spherical correctors

4.2

A Simple Simulation Example

In order to better understand the homotopy methods, a test for a simple nonlinear set of equations was carried out with MATLAB. Consider the set of equations ( x2 + y 2 − 1 = 0 , (31) −x2 + y = 0 which describes a unity circle and parabola crossing each other in two points y=



5−1 ≈ 0.61803, 2

s√

x=± 20

5−1 ≈ ±0.78615 2

(32)

2 x y 1.5

1

x,y

0.5

0

−0.5

−1

−1.5

−2

0

1

2

3

4

5

λ

6

7

8

9

10

11

Figure 6: Fixed-point homotopy curves of Equation 31.

x y

1 x=0.7862

0.5

y=0.6180

x,y

0

−0.5 x=−0.7862

−1

0.92

0.94

0.96

0.98

1 λ

1.02

1.04

1.06

Figure 7: The solution is found at λ = 1.

21

1.08

30 x y

20

x,y

10

0

−10

−20

−30

0

0.1

0.2

λ

0.3

0.4

Figure 8: Diverging homotopy. The ODE-based solver by Heath Hoffmann [33] was used in tracking the zero curve. The fixed-point homotopy was used, resulting in function (

(1 − λ)(x − a1 ) + λ(x2 + y 2 − 1 ) = 0 . (1 − λ)(x − a2 ) + λ(−x2 + y )= 0

(33)

The initial guess is important for finding the solutions. For example,if the initial guess is [a1 = 0.1, a2 = −0.1], the algorithm will find both the solutions. If the initial guess is [a1 = 3, a2 = 4], the homotopy curve will fold back before λ = 1 and none of the solutions is found (see Figure 8). Results for different values of a are collected in Table 2.

4.3

Standard Newton-Raphson Solver

The standard Newton-Raphson combined with homotopy functions is not an efficient way to solve the circuit equations. The main drawback of the algorithm is that it is not able to deal with foldbacks of the curve. Three homotopy equations were tested: the fixed-point homotopy, the variable stimulus homotopy, and the modified variable stimulus homotopy. Newton-Raphson continuation was tested using following values: the step for λ = 0.05, the conductance G = 10−2 and the criterion for terminating the corrector iteration is that the norm of the function is smaller than 10−5 . The number of iterations depends on initial guess a. The solver was run five times and the average is written on Table 3. The range for initial guess was a ∈ [0, 1]. 22

Table 2: The effect of initial guess a1 0.1 0.0 0.1 0.0 -0.1 -0.1 0.1 -0.1 1.2 1 3

a2 0.0 0.1 0.1 -0.1 0.0 -0.1 -0.1 0.1 1 1.2 4

Solutions found Positive Positive Positive Both Positive Both Both Both Both Both None

Table 3: The standard Newton-Raphson continuation - the number of Jacobian evaluations. Fixed-point Variable stimulus palikka.i 50 154 simdc.i 82 58 (1 toronto.i 63 57 voter.i fail fail flipflop.i 70 72 globtest.i fail fail dcprob.i fail fail ssjudc.i fail fail schmitt.i 72 90 ref6.i 143(1 452(1 dc demo.i fail fail dcprob2.i fail fail 1 Converged only after modifications to

4.3.1

Mod. variable stimulus 369 289 367 fail 331 fail fail fail 1071 994 fail fail certain parameters.

The Fixed-Point Homotopy Function

Toronto.i required setting G = 10−4 in order to converge. Voter.i caused a floating point error near λ = 1. Both of the stable operating points of flipflop.i are found. The operating point which is found depends on the initial guess of a. Globtest.i has a nearly singular Jacobian at λ = 0.95. Dcprob.i will fail at 23

λ = 0.65 for the same reason. Ssjudc.i is even more difficult, it will fail in a nearly singular Jacobian matrix when λ = 0.05. MATLAB detects the (nearly) singularity of Jacobian matrix from the condition number of the matrix. For a well-conditioned matrix, the condition number is close to 1. For a nearly singular matrix, the condition number is very large, 10100 for example. The most common reason for the singular Jacobian is that the iteration diverges. The schmitt.i circuit has two operating points, of which just one is found with the parameters described. However, if the range of the initial guess is extended to a ∈ [0, 2], the other operating point is found, too.

The circuit ref6.i has numerical problems at λ = 1. The Jacobian remains non-singular but the iteration does not converge. The problem can be worked around by raising the norm condition (the criteria for stopping the iteration based on the norm of the homotopy function) from 10−5 → 10−4 and changing the range of initial guess to a ∈ [0, 10]. In this manner, at least four different operating points (of 11 possible) are found. The dc demo.i first converges well but at λ = 1 a floating point error occurs due to nearly singular Jacobian matrix. The circuit dcprob2.i has similar problems at λ = 0.9. 4.3.2

The Variable Stimulus Homotopy Function

Palikka.i and simdc.i converged without problems. Toronto.i converged without the modifications needed in fixed-point homotopy. Some of the initial guesses resulted in numerical failure due to nearly singular Jacobian matrix. Voter.i has a nearly singular Jacobian near λ = 1, and fails to converge. For flipflop.i, both of the stable operating points are found. Simulation of circuits globtest.i, dcprob.i and ssjudc.i fail in similar manner as in fixed-point homotopy. Both of the operating points of the Schmitt trigger are found when the initial voltage vector range is extended to a ∈ [0, 2].

The circuit ref6.i needs the initial voltage vector values to be a ∈ [0, 10] in order to converge. Unlike in fixed-point homotopy, the iteration stop criteria needs no changes. Unfortunately, only one operating point is found. Some iterations do not converge. If the stopping criteria for function norm is altered to 10−4 , three more operating points are found. The circuits dc demo.i and dcprob2.i fail. The iteration does not converge at λ = 1.

The variable stimulus homotopy doesn’t seem to have any remarkable advantages when compared to the fixed-point homotopy. 24

4.3.3

The Modified Variable Stimulus Homotopy Function

The modified variable stimulus homotopy seems to converge very poorly, if G is few decades smaller than 1. For example, palikka.i did not converge after 10000 iterations. Therefore, value 1 was used. The simple circuits palikka.i and simdc.i did converge well. The voter.i did converge extremely slowly and resulted in a numerical failure due to closely singular Jacobian matrix. Both of the operating points of flipflop.i are found. Globtest.i will fail at λ = 0.95 and the circuits dcprob.i and ssjudc.i fail too, all due to closely singular Jacobian matrix. For schmitt.i the initial voltage range must be a ∈ [0, 4] in order to find both of the operating points. The circuit ref6.i converges, but only one operating point is found. Dc demo.i and dcprob2.i both fail due to nearly singular Jacobian.

Combination of the Newton-Raphson solver and the modified variable stimulus homotopy seem to be inefficient if compared to the other two methods, because more iterations are needed, see Table 3. 4.3.4

Test for NR-Continuation with Equal Parameters

Some methods, which do not converge, can be made convergent by adjusting the error criteria, step length, or the conductance scaling factor. However, to compare objectively the efficiency of the methods, a test with constant parameters is made. The results are shown in Table 4. The values are an average of five different solver runs. The source code for the method is included in Appendix A.

4.4

ODE-Based Curve Tracing

The ODE-based homotopy solver was used with default tolerance options with the fixed-point homotopy function and two different values for G. The summary of Jacobian evaluations needed for finding the DC operating point (or the first DC operating point, if multiple operating points exist) is presented in Table 5. The values are calculated as an average of five different solver runs. 4.4.1

The Fixed-Point Homotopy

The operating points of palikka.i and simdc.i are found easily. The operating point of toronto.i is found easily, but the solver experiences numerical oscillation at λ = 1. That is, the zero curve is traced many times over the line λ = 1, so that the solver finds the same operating point about ten times. The 25

Table 4: The standard Newton-Raphson continuation with different G and step length. The number of Jacobian evaluations. Circuit Step G palikka.i simdc.i toronto.i voter.i flipflop.i globtest.i dcprob.i ssjudc.i schmitt.i ref6.i dc demo.i dcprob2.i

Fixed-point 0.05 0.001 10−1 10−3 10−1 10−3 46 47 2001 1092 53 55 2001 1126 46 fail fail 1103 fail fail fail fail 62 59 2013 1059 fail fail fail 1321 fail fail fail fail fail fail fail fail 70 56 2001 1106 fail fail fail 1198 fail fail fail 2649 fail fail fail fail

Variable stim. 0.05 0.001 10−1 10−3 10−1 10−3 57 339 2001 4072 73 55 2227 1007 44 64 2217 1003 fail fail fail fail 69 71 2021 2001 fail fail fail 2568 fail fail fail fail fail fail fail fail 151 159 2734 2524 fail fail fail 2001 fail fail fail fail fail fail fail fail

Mod. var. stim 0.05 0.001 10−1 10−3 10−1 10−3 49 147 2001 2001 55 55 2012 2007 51 263 2007 2012 fail fail fail fail 60 107 2022 1319 fail fail fail 3287 fail fail fail fail fail fail fail fail 55 82 2270 1430 fail 178 3954 2192 fail fail fail fail fail fail fail fail

problem can be worked around by setting the lambda tolerance (parameter that defines, how close to λ = 1 the solver must be, in order to accept the solution) parameter to 10−15 . The voter.i circuit can not be simulated with this method due to lack of memory. The matrices used in the algorithm are too large to fit in the buffer of MATLAB-APLAC –interface. All the operating points (two stable and one unstable) of circuit flipflop.i are found, if the lambda tolerance is dropped to 10−15 . If the default value is used, only one of the operating points is found at a time. The initial guess defines which operating point is found. The fixed-point homotopy and ODE-based curve tracing algorithm is the first homotopy–algorithm combination that finds the DC solution of globtest.i. The solver has same kind of oscillation near λ = 1 like with toronto.i, but the problem is fixed if the lambda tolerance is changed to 10−15 . Solving the circuit dcprob.i is too large to be tested with ODE algorithm. After 30 minutes of CPU time the homotopy parameter λ has progressed to 0.005. At that speed, the test run would take at least 100 hours. Additionally, the node voltages begin to rise very steeply, indicating that the algorithm will diverge. Circuit ssjudc.i has similar problems, the MATLAB-APLAC -interface crashes at the start of the ODE solver. All of the operating points of schmitt.i are found, if the lambda tolerance parameter is changed to 10−15 . The number of Jacobian evaluations needed for finding all the operating points is 510. The ODE solver finds only one operating point of ref6.i. The circuit dc demo.i diverges, the voltages approach infinity when λ → 1. The same problem occurs with dcprob2.i. 26

Table 5: The ODE-based homotopy method with different values of G. palikka.i simdc.i toronto.i voter.i flipflop.i globtest.i dcprob.i ssjudc.i schmitt.i ref6.i dc demo.i dcprob2.i

4.5

G = 0.1 83 207 477 fail (memory) 197 1372 fail (too slow to be tested) fail (program crash) 158 344 fail fail

G = 0.001 115 260 470 184 1334

175 310 fail fail

Pseudo-Arclength Solver with the Planar Corrector

The pseudo-arclength solver introduced in Chapter 3.3.3 was tested with different step size and the conductance scaling factor G. The results are in Table 6. The method seems to be computationally efficient, and choice of G has no significant influence in the count of Jacobian evaluations. The source code for the solver is included in Appendix B.

4.6

Pseudo-Arclength Solver with the Spherical Algorithm

A spherical algorithm is a predictor-corrector method for curve tracing. The outlines of the algorithm are presented at the end of Chapter 3.3.3. The tests were run for two different values for G and sphere radius. The results are presented in Table 7.

27

Table 6: The arclength algorithm for fixed-point homotopy function with different G and step length. Circuit Step palikka.i simdc.i toronto.i voter.i flipflop.i globtest.i dcprob.i ssjudc.i schmitt.i ref6.i dc demo.i dcprob2.i

G = 0.1 0.2 0.05 78 306 93 247 96 552 fail fail 65 317 fail fail fail fail fail fail 106 583 215 1011 fail fail fail fail

G = 0.001 0.2 0.05 84 314 71 77 151 626 fail fail 87 301 307 fail fail fail fail fail 153 581 239 1235 fail fail fail fail

Table 7: The Spherical algorithm and fixed-point homotopy function with different radius of sphere and different values of G. Circuit Radius palikka.i simdc.i toronto.i voter.i flipflop.i globtest.i dcprob.i ssjudc.i schmitt.i ref6.i dc demo.i dcprob2.i

G = 0.1 0.01 0.005 1383 2942 1543 2397 2605 6082 fail fail 1162 2294 fail fail fail fail fail fail 3297 6459 4633 8223 fail fail fail fail

28

G = 0.01 0.01 0.005 fail 4553 fail fail fail fail fail fail fail fail fail fail fail fail fail fail fail 7758 fail fail fail fail fail fail

5

Discussion

In the literature, the homotopy methods are mentioned to be globally convergent. In theory, that is the fact. In practice, while tracing the zero curve of the homotopy function, the usual numerical problems are encountered. The Jacobian may be nearly singular near λ = 1. The zero curve has very steep folds which cause the method to jump back on a point on the curve already traced and therefore the algorithm will end up in an endless loop. As with the plain Newton–Raphson method without homotopy implementation or any other convergence-aiding technique, the homotopy methods seem to be very sensitive to the initial guess a. A good starting point is essential for fast convergence of the algorithm [25, 27]. For example, the circuits schmitt.i and ref6.i failed to converge on some initial guesses. The most common problem with the curve tracing methods tested, was singularity or nearly singularity (condition number > 1090 ) of the Jacobian matrix near the final solution λ = 1. The actual reasons of singularity are hard to find out for a set of hundreds of very complex nonlinear equations. One possible reason may be high-impedance nodes in the circuit [26]. One goal for this thesis was to try to find out which homotopy method (or methods) would be good enough to implement into the APLAC circuit simulator. Based on this thesis, the pseudo-arclength looks to be the most robust method for such implementation. An efficient algorithm for step length control is a necessity when implementing homotopy methods into a production version of a circuit simulator. In this thesis, the step length was kept constant, but in order to save computation time and improve accuracy at problem points on the zero curve, an efficient implementation of step control algorithm should be reviewed. A good starting point for such a study would be [24]. The main problems confronted with testing the homotopy methods were numerical problems while tracing the zero curve of the homotopy function. In the MATLAB–APLAC -interface, all the other convergence-aiding techniques implemented in APLAC were not in use. By implementing a homotopy method directly into the simulator core (written in C programming language), the convergence-aiding techniques already implemented would be used to improve the numerical stability of the homotopy algorithm. In [21] it is stated that complex parameter homotopy methods can find all solutions of circuit without passing through folds or singular points. The complex parameter homotopy methods would be worth deeper study and implementation tests. Another interesting kind of homotopy method for automatically finding all the DC operating points for circuits is described in [28]. 29

The disadvantage of the method is that implementation requires tampering the component models of the nonlinear elements in the circuits. It may be also worth to test calling HOMPACK or PITCON source codes directly from the simulator.

30

6

Conclusion

In this thesis, improving the convergence of DC analysis of a circuit simulator by using homotopy methods was studied. A literature review about homotopy methods in finding DC operating points of electric circuits was made. Three methods of constructing homotopy functions and three methods for solving them were reviewed more closely and tested with APLAC circuit simulator utilizing the MATLAB–APLAC -interface developed in the Circuit Theory laboratory of the Helsinki University of Technology. Two solver algorithms were implemented in MATLAB language and one ready algorithm was tested. The results from the MATLAB tests were reviewed. According to MATLAB tests, choosing the homotopy function is not a crucial part of a successful homotopy method. The solver algorithm, which is capable of dealing with problematic points on the zero curve, is a necessity for an efficient homotopy method. The succession of the homotopy methods tested depends on the initial guess a of the homotopy path, especially with large and complex circuits. The fixed-point homotopy method was found to be robust, easy to implement, and computationally efficient. All the homotopy methods were clearly helpful in improving convergence — the circuit simdc.i, for example, does not converge with ordinary Newton–Raphson iteration. The arclength continuation method with planar corrector turned out to be efficient and stable when tracing the zero curve of the fixed-point homotopy function. The homotopy methods were tested without any other convergence-aiding methods. Homotopy methods combined with other convergence-aiding techniques have potential to improve the DC analysis convergence drastically in APLAC circuit simulator.

31

References [1] Aplac Solutions Corporation, Finland, APLAC 8.00 Reference Manual and User’s manual, 2004, http://www.aplac.com. [2] M. Valtonen, P. Heikkil¨a, H. Jokinen, and T. Veijola, “APLAC – objectoriented circuit simulator and design tool,” in Low-power HF Microelectronics: a Unifed Approach (G. A. S. Machado, ed.), pp. 333-372, IEE, 1996. [3] J. Ogrodzki, Circuit Simulation Methods and Algorithms, CRC Press Inc., 1994. [4] J. Roos, “Improving the Speed and Convergence of DC Analysis by Means of Self-Generating Lookup Tables and Piecewise-Linear Analysis,” Doctoral thesis, Helsinki University of Technology, Espoo, 1999. [5] M. Honkala, J. Roos, and V. Karanko, “On Nonlinear Iteration Methods for DC Analysis of Industrial Circuits”, in Mathematics in Industry 8: Progress in Industrial Mathematics at ECMI 2004 (A. Di Bucchianico, R. M. M. Mattheij, and M. A. Peletier, ed.), pp. 144-148, Springer, Berlin/Heidelberg, 2006. [6] L.O. Chua and P.-M. Lin, Computer-Aided Analysis of Electronic Circuits: Algorithms & Computational Techniques, Prentice & Hall Inc., 1975. [7] E. W. Weisstein, CRC Concise Encyclopedia of Mathematics, CRC Press, 1999. [8] L. Trajkovic, R. C. Melville and S. C. Fang, “Finding dc operating points of transistor circuits using homotopy methods,” Proc. IEEE Int. Symp. Circuits and Systems, Singapore, June 1991, pp. 758-761. [9] L. Trajkovic, R. C. Melville and S. C. Fang, “Improving dc convergence in a circuit simulator using a homotopy method,” Proc. IEEE Custom Integrated Circuits Conference, San Diego, CA, May 1991, pp. 8.1.1-8.1.4. [10] Cadence Design Systems, USA, PSpice Convergence Guide, 2003, http://www.orcad.com [11] A. Dyess, E. Chan, H. Hofmann, W. Horia, and L. Trajkovic, “Simple implementations of homotopy algorithms for finding dc solutions of nonlinear circuits,” Proc. IEEE Int. Symp. Circuits and Systems, Orlando, FL, June 1999, vol. VI, pp. 290–293.

32

[12] R. C. Melville, L. Trajkovic, S. C. Fang, and L. T. Watson, “Artificial parameter homotopy methods for the dc operating point problem,” IEEE Trans. on CAD, vol. 12, pp. 861–877, June 1993. [13] L. Watson, “Globally Convergent Homotopy Methods,” SIAG/OPT Views-and-News, vol. 11, pp. 1–6, April 2000. [14] J. S. Roychowdhury, R. C. Melville, “Homotopy techniques for obtaining a DC solution of large-scale MOS circuits,” Proceedings of the 33rd Annual Conference on Design Automation, United States, pp. 286-291, June 1996. [15] L. Trajkovic, R. C. Melville and S. C. Fang, “Passivity and no-gain properties establish global convergence of a homotopy method for dc operating points,” Proc. IEEE Int. Symp. Circuits and Systems, New Orleans, LA, May 1990, pp. 914-917. [16] C. P. Reames and A.N. Willson, Jr., “Global Convergence of Newton’s Method in the DC analysis of Single-Transistor Networks,” Electronics Letters, vol. 18, no. 12, pp. 519-520, June 1982. [17] K. Yamamura, T. Sekiguchi and Y. Inoue, “A Fixed-Point Homotopy Method for Solving Modified Nodal Equations,” IEEE Transactions on Circuits and Systems, vol. 46, no. 6, June 1999. [18] M. M. Green, “An Efficient Continuation Method for Use in Globally Convergent DC Circuit Simulation,” ISSSE ’05 Proceedings, pp. 497-500, October 1995. [19] P. Rodrigues, Computer-aided analysis of nonlinear microwave circuits, Artech House, 1998. [20] A. Robotycki, R. Zielonko, “Fault Diagnosis of Analog Piecewise Linear Circuits Based on Homotopy,” IEEE Transactions on Instrumentation and Measurement, vol. 51, pp. 876-881, no. 4, August 2002. [21] D. M. Wolf, S. R. Sanders, “Multiparameter Homotopy Methods for Finding DC Operating Points of Nonlinear Circuits,” IEEE Transactions on Circuits and Systems–I: Fundamental Theory and Applications, vol. 43, no. 10, October 1996. [22] K. Yamamura, “Simple algorithms for tracing solution curves,” IEEE Transactions on Circuits and Systems–I: Fundamental Theory and Applications, vol.40, no.8, pp. 537-541, August 1993. [23] E. L. Allgower, K. Georg, Numerical Continuation Methods – An Introduction, Springer-Verlag, 1990. 33

[24] R. Kearfott, Z. Xing, “An Interval Step Control for Continuation Methods”, SIAM J. Numer. Anal., vol. 31, pp. 892–914, 1994. [25] L. Trajkovic, Wiley Encyclopedia of Electrical and Electronics Engineering, vol. 9, pp. 171–176, John Wiley & Sons, Inc., 1999. [26] L. Goldgeisser, M. Green, “Using Continuation Methods to Improve Convergence of Circuits with High Impedance Nodes”, IEEE International Symposium on Circuits and Systems, May 28–31, 2000, Geneva. [27] L. Trajkovic, W. Mathis, “Parameter embedding methods for finding dc operating points: formulation and implementation”, Proc. NOLTA ’95, Las Vegas, NV, Dec. 1995, pp. 1159–1164. [28] L. Goldgeisser, M. Green, “A Method for Automatically Finding Multiple Operating Points in Nonlinear Circuits”, IEEE Transactions on Circuits and Systems, vol. 52, no. 4, April 2005. [29] L.T. Watson, “A Globally Convergent Algorithm for Computing Fixed points of C 2 Maps”, Appl. Math. Comput. 5, 1979, pp. 297-311. [30] L.T. Watson, S.C. Billups, A.P. Morgan, “ALGORITHM 652: HOMPACK: a suite of codes for globally convergent homotopy algorithms”, ACM Transactions on Mathematical Software (TOMS), vol. 13, no. 3, pp. 281-310, September 1987. [31] S. A. Ragon, Z. Grdal, L.T. Watson, “A comparison of three algorithms for tracing nonlinear equilibrium paths of structural systems”, Internat. J. Solids Structures, vol. 39, pp. 689–698, 2002. [32] W. Rheinboldt and J. Burkardt, “Algorithm 596: A Program for a Locally Parameterized Continuation Process”, ACM Trans. Math. Softw., vol. 9, no. 2, pp. 236–241, 1983. [33] Heath Hoffmann, ODE Based Homotopy power.eecs.berkeley.edu/heath/matlab.html

Solver,

http://www-

[34] L. Trajkovic, E. Fung, S. Sanders, “HomSPICE: simulator with homotopy algorithms for finding DC and steady-state solutions of nonlinear circuits”, Proc. ISCAS ’98, Monterey, CA, May 1998, pp. 227–231.

34

Appendix A % A source code for testing simple NR continuation % with MATLAB-APLAC interface. clear all; file_name=’/u1/vesa/erikoistyo/ajot/toronto/toronto.i’; % Select the homotopy function hom=3 % 1 = fixed-point, 2 = modified var.stimulus and 3 = var.stimulus

[h,n,n2] = startaplac(’../aplac.run’,file_name); % Start APLAC a=rand(n,1); % Initialize the random vector. x=a; l=0.0; gg=0.001; step=0.001; [ff,J,f,dl] = acalcf(x,h,n,l,a,hom,n2,gg); % Calculate function nstep=1; % values and Jacobian. niter=1; while l1e-4 % Solve the function with NR method dx = J\-f; x = x+dx; [ff,J,f,dl] = acalcf(x,h,n,l,a,hom,n2,gg); normf=norm(f) niter=niter+1; % Count the Jacobian evaluations end xx(nstep)=x(1); ll(nstep)=l; nstep=nstep+1; end endaplac(h); % Close APLAC. niter % Show number of Jacobian evaluations.

35

Appendix B % File for testing arc-length method for homotopies. close all; % Ikkunat sulki clear all; % Siivousoperaatio hom=1; % 1 is fixed-point, 2 is mod.var.stim, 3 is var.stim %file_name=’/u1/vesa/erikoistyo/ajot/toronto/toronto.i’; [h,n,n2] = startaplac(’../aplac.run’,file_name); a = rand(n,1); l = 0.0; % Homotopy parameter lambda ds=0.05; % The step length nn=1; rst=0; itertotal=0; x=a; gg=0.1; xp=0; lp=0; change=1; lastnorm=1; %predictor while l1e-6 iter=iter+1; 36

[f,J,hh]=acalcf(xp,h,n,lp,a,hom,n2,gg); g=[ux;ul]’*[xp-x;lp-l]-ds; dxl=[J,f+gg*(a-xp);ux’,ul]\-[hh;g]; newxl=[xp;lp]+dxl; lp=newxl(n+1); % l -> lp xp=newxl(1:n); itertotal=itertotal+1; change=norm(hh)-lastnorm; lastnorm=norm(hh); end x=xp; l=lp; % Store the points for drawing curves. lvec(nn)=l; x1vec(nn)=x(1); x2vec(nn)=x(2); nn=nn+1; itervec(nn)=iter; itercnt=0; itertotal=itertotal+1; end

endaplac(h); % closes APLAC plot(lvec,x1vec,lvec,x2vec) % Plot some node voltage curves. itertotal % Show the total number of iterations.

37

Suggest Documents