**Analog CAD**

## **Memetic Method for Passive Filters Design**

Tomasz Golonek and Jantos Piotr

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/53716

## **1. Introduction**

The design of analog passive filters with specialized (not typical) frequency responses is not a trivial problem. The presence of finite load impedances for filter sections and limited qual‐ ity factors of coils are just two of many concerns which a design engineer has to take into account. Additionally, classical techniques of filters synthesis require assuming of the ap‐ proximation type (e.g. Butterworth, Chebyshev) before calculating the filter transfer func‐ tion's poles and zeroes. This choice is frequently a challenge itself.

One of the methods allowing for elimination of the mentioned problems is the use of evolu‐ tionary computations (EC). Evolutionary techniques are a well known and frequently used tool of global optimization [1-3]. This kind of the optimization imitates natural processes of individuals' competition as candidates for reproduction. Better fitted individuals have high‐ er survival probability and their genetic material is preferred. During the recombination process some parts of parents' genotypes are exchanged and offspring individuals are creat‐ ed. A new generation collected after the succession procedure conserves the features consist‐ ed in the previous genotypes. Besides, to assure a system resistance for a stagnation effect, mutation operations are applied to EC. The most popular sorts of EC approaches are: genet‐ ic algorithm (GA), genetic programming (GP), evolutionary strategies (ES), differential evo‐ lution (DE) and gene expression programming (GEP).

The main drawback of evolutionary approaches is an ineffective and insufficient local opti‐ mization. This property and significant computational efforts necessary for a huge genera‐ tion processing predispose the EC applicability especially to the NP hard global searching problems [4-10]. This chapter describes the passive filters synthesis method by means of EC. The process of circuits' automated designing is a very complex issue. A wide area of solu‐ tions should be probed during the early stage of computations and its local parameters should be finally optimized. In contrary to the alternative systems [4-6], the method present‐

© 2013 Golonek and Piotr; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Golonek and Piotr; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

ed in this chapter is based on an application of a hybrid system - a synergy of genetic pro‐ gramming (GP) (used for the purpose of determining an optimal network of a passive filter circuit) and a deterministic local search by the means of Hooke and Jeeves method (HJM), which enables the system to find accurate values of the filter's elements. The proposed de‐ sign system allows for obtaining the desired frequency response and, optionally, production yield optimization.

Section 2 explains the general algorithm of the proposed system, Sections from 3 to 5 present the descriptions of the important details of the algorithm. Next, in Section 6, the exemplary results of an automated circuit design are placed. Finally, in Section 7, some considerations of the method future development and final conclusions are presented.

## **2. Optimization process overview**

The process is initialized with the desired filter specifications. Additional algorithms' pa‐ rameters, e.g. population size, number of Monte Carlo analyses, and the like, are assumed. The use of GP and HJM is briefly presented in following paragraphs.

In the presented research both filters topology and circuit parameters values are being opti‐ mized. As far as GP has been proven an effective tool of circuits networks determination, it is not as efficient in adjusting resistors, coils, or capacitors values. The latter has been solved by the use of a deterministic, non-gradient local search algorithm - Hooke and Jeeves direct search method [11]. A synergy of evolutionary global optimization and local optimization algorithm is called a memetic algorithm [12-16]. In the presented research the proposition of the memetic genetic programming (MGP) introduction for the purpose of the analog filters design is described.

The block diagram of the optimization process has been presented in Fig. 1. After system running, the first, primary generation of the *Nmx* individuals is created randomly. It is very important to assure the possible wide range of dispersion for the starting solutions, so the diversity of this population is extremely desired. In the proposed system, the uniform prob‐ ability of the primary individuals' randomization with the maximal allowed size limitation is applied. To evaluate the actual, random solutions, the fitness function is determined for each individual from population. The distances between the parameters of the evaluated phenotypes and the target specifications are checked during this process. Next, the GP sys‐ tem is executed. During reproduction, the mating pool is collected with the reproduction method that prefers more fitted individuals. An adequate strength of selection pressure has to be kept during this process, and it has crucial impact for a system convergence. During recombination process, pairs selected from an intermediate pool are crossed with assumed probability. Besides, offspring genotypes can be mutated and it assures adequate veracity of the population and allows for achieving the new regions of the searching space. Detailed in‐ formation about the genetic operations and the fitness function of the GP part of the pro‐ posed memetic system are included in Section 4.

**Figure 1.** The process overview

ed in this chapter is based on an application of a hybrid system - a synergy of genetic pro‐ gramming (GP) (used for the purpose of determining an optimal network of a passive filter circuit) and a deterministic local search by the means of Hooke and Jeeves method (HJM), which enables the system to find accurate values of the filter's elements. The proposed de‐ sign system allows for obtaining the desired frequency response and, optionally, production

Section 2 explains the general algorithm of the proposed system, Sections from 3 to 5 present the descriptions of the important details of the algorithm. Next, in Section 6, the exemplary results of an automated circuit design are placed. Finally, in Section 7, some considerations

The process is initialized with the desired filter specifications. Additional algorithms' pa‐ rameters, e.g. population size, number of Monte Carlo analyses, and the like, are assumed.

In the presented research both filters topology and circuit parameters values are being opti‐ mized. As far as GP has been proven an effective tool of circuits networks determination, it is not as efficient in adjusting resistors, coils, or capacitors values. The latter has been solved by the use of a deterministic, non-gradient local search algorithm - Hooke and Jeeves direct search method [11]. A synergy of evolutionary global optimization and local optimization algorithm is called a memetic algorithm [12-16]. In the presented research the proposition of the memetic genetic programming (MGP) introduction for the purpose of the analog filters

The block diagram of the optimization process has been presented in Fig. 1. After system running, the first, primary generation of the *Nmx* individuals is created randomly. It is very important to assure the possible wide range of dispersion for the starting solutions, so the diversity of this population is extremely desired. In the proposed system, the uniform prob‐ ability of the primary individuals' randomization with the maximal allowed size limitation is applied. To evaluate the actual, random solutions, the fitness function is determined for each individual from population. The distances between the parameters of the evaluated phenotypes and the target specifications are checked during this process. Next, the GP sys‐ tem is executed. During reproduction, the mating pool is collected with the reproduction method that prefers more fitted individuals. An adequate strength of selection pressure has to be kept during this process, and it has crucial impact for a system convergence. During recombination process, pairs selected from an intermediate pool are crossed with assumed probability. Besides, offspring genotypes can be mutated and it assures adequate veracity of the population and allows for achieving the new regions of the searching space. Detailed in‐ formation about the genetic operations and the fitness function of the GP part of the pro‐

of the method future development and final conclusions are presented.

The use of GP and HJM is briefly presented in following paragraphs.

**2. Optimization process overview**

posed memetic system are included in Section 4.

yield optimization.

52 Analog Circuits

design is described.

This heuristic stage impacts on the coded circuits topology especially (i.e. the values of the filter resistances, inductances and capacitances are optimized not effectively), however the values of its elements are adjusted on the next stage during local deterministic optimization.

The newly created population of solutions represents a sort of circuits with non optimal values of its elements and now they are determined by means of HJM (pattern search) algorithm.

The differential of the fitness function is unknown. Hence, the local search algorithm could not have required its gradient computation. Among various optimization methods that meet this requirement, the pattern search method is characterized by its simplicity. The HJM method consists of two moves, i.e. exploratory and pattern, repeated sequentially until stop conditions are fulfilled.

Evolutionary algorithms are a trade-off between a global and local optimization. As far as they provide a good solution it never is but a sub-optimal one. The aim of the presented research was to achieve the best possible solution. Hence, it has been decided to use local search algo‐ rithms. The application of a local search algorithm after the evolutionary optimization is com‐ pleted would improve the performance of the last chosen individual. Though, it would not affect the optimization process itself. Memetic solutions aim to improve the overall optimiza‐ tion process. The local search algorithm is applied within the evolutionary algorithm's main loop. It allows for faster convergence of the process. Applying a full local search process within each of the evolutionary iterations would influence the required computation time significant‐ ly and negatively. Therefore, only a few iterations of the HJM were allowed. For the purpose of improving the final result a full HJM cycle is applied after the evolutionary cycle is over.

In the presented approach the Hooke and Jeeves method had to be modified so the search is carried out in the assumed search space (E24 sequence of resistances, capacitances and in‐ ductances). The detailed description of the algorithm is presented in Section 5.

Next, after the topology and the values of elements determination, if the production yield optimization is required, the Monte Carlo analysis can be applied. The algorithm is termi‐ nated if all specification (and yield) requirements are met or if the last allowed generation *Gmx* was reached.

## **3. Circuit structure coding**

The filter circuit structure is coded with the use of a binary tree structure. An exemplary tree is illustrated in Fig. 2. Its inner nodes contain functions, however terminals keep their argu‐ ments. The length from a root node to the most remote leaf is defined by the maximal depth *Dmx* and it reaches the value *Dmx*=4 for the illustrated case. Finally, the tree presented in Fig.2 can be decoded as a general expression given below:

$$y = f\_1\left(f\_2\left(t\_1, f\_3\left(t\_2, t\_3\right)\right), f\_4\left(t\_4, t\_5\right)\right). \tag{1}$$

**Figure 2.** The example of tree structure

This kind of structure allows for defining the circuit topology and values of its elements in the flexible way with the assumed type of connection functions defined in **F** and the set of basic terminal blocks defined in **T**. Besides, it enables an easy way to the genetic operations implementation in the GP system. The set F contains unique symbols for connection types and for the described system it was defined as follows:

$$\mathbf{F} = \{ \underline{\qquad}, \; \underline{\qquad} \qquad \| \, \| \, \underline{\qquad} \tag{2}$$

where contained symbols denote a serial, a cascade and a parallel configuration respec‐ tively.

ly and negatively. Therefore, only a few iterations of the HJM were allowed. For the purpose of improving the final result a full HJM cycle is applied after the evolutionary cycle is over.

In the presented approach the Hooke and Jeeves method had to be modified so the search is carried out in the assumed search space (E24 sequence of resistances, capacitances and in‐

Next, after the topology and the values of elements determination, if the production yield optimization is required, the Monte Carlo analysis can be applied. The algorithm is termi‐ nated if all specification (and yield) requirements are met or if the last allowed generation

The filter circuit structure is coded with the use of a binary tree structure. An exemplary tree is illustrated in Fig. 2. Its inner nodes contain functions, however terminals keep their argu‐ ments. The length from a root node to the most remote leaf is defined by the maximal depth *Dmx* and it reaches the value *Dmx*=4 for the illustrated case. Finally, the tree presented in Fig.2

This kind of structure allows for defining the circuit topology and values of its elements in the flexible way with the assumed type of connection functions defined in **F** and the set of basic terminal blocks defined in **T**. Besides, it enables an easy way to the genetic operations implementation in the GP system. The set F contains unique symbols for connection types

*y f ft ftt ftt* = 1 21 3 23 4 45 ( ( ) , ,, ,. ( ) ( )) (1)

**F** = -- >> { } , ,|| , (2)

ductances). The detailed description of the algorithm is presented in Section 5.

*Gmx* was reached.

54 Analog Circuits

**3. Circuit structure coding**

**Figure 2.** The example of tree structure

can be decoded as a general expression given below:

and for the described system it was defined as follows:

Generally, passive filters circuits are constructed by the use of resistors, capacitors and mag‐ netic coils (inductors), so the T-typed three-terminal (two-port) blocks which contain these elements (one in each branch) were chosen as basic ones:

$$\mathbf{T} = \left\{ \mathbf{R}\boldsymbol{\upsilon}\_{e}\boldsymbol{o}\_{u}\mathbf{L}\boldsymbol{\upsilon}\_{e}\boldsymbol{o}\_{u}\mathbf{C}\boldsymbol{\upsilon}\_{e}\boldsymbol{o}\_{u'}\mathbf{L}\boldsymbol{\upsilon}\_{e}\boldsymbol{o}\_{u}\mathbf{C}\boldsymbol{\upsilon}\_{e}\boldsymbol{o}\_{u}\mathbf{R}\boldsymbol{\upsilon}\_{e}\boldsymbol{o}\_{u'}\mathbf{C}\boldsymbol{\upsilon}\_{e}\boldsymbol{o}\_{u'}\mathbf{R}\boldsymbol{\upsilon}\_{e}\boldsymbol{o}\_{u}\mathbf{L}\boldsymbol{\upsilon}\_{e}\boldsymbol{o}\_{u}\right\},\tag{3}$$

where integer numbers *ve* for *e*=(1,..,*E*) and *ou* for *u*=(1,..,*U*) determine the value and its order for the element identified by the antecedent descriptor (i.e. *R* for a resistor, *L* for a inductor and *C* for a capacitor). In the proposed system, discrete values *vi* are selected from a set **V** that contains elements from a widely known practical series E24 (*E*=24) prepared for the components manufactured with tolerance equal to *δtol*=5%:

$$\mathbf{V} = \begin{pmatrix} 10, 11, 12, 13, 15, 16, 18, 20, 22, 24, 27, 30, \\ 33, 36, 39, 43, 47, 51, 56, 62, 68, 75, 82, 91 \end{pmatrix} \text{.} \tag{4}$$

Order values are included in the set **O** and they should enable the desired range of coding, and consist of *U* integer ciphers ( 0 ≤ *ou* ≤ 9 ) which determine exponents of decimal multipliers:

$$\mathbf{O} = \left\langle o\_{0\prime}, \dots, o\_{\{\mathbf{U}\cdot\mathbf{1}\}} \right\rangle \; . \tag{5}$$

Finally, values of resistance *R*, inductance *L* and capacitance *C* are calculated from the equa‐ tions (6) ÷ (8) respectively:

$$\mathcal{R} = \upsilon\_{\varepsilon} \cdot 10^{\alpha\_{\mu}} \{\Omega\}\_{\nu} \tag{6}$$

$$L = v\_e \cdot 10^9 \text{\AA}\text{\AA}\text{\AA}\tag{7}$$

$$\mathbf{C} = \boldsymbol{\upsilon}\_e \cdot \mathbf{1} \mathbf{0}^{\boldsymbol{\theta}\_u} \text{[pF]}.\tag{8}$$

The idea of the terminals' T (basic circuit blocks) coding is illustrated in Fig.3, where the three exemplary two-ports and adequate coding strings are presented. Impedances of ele‐ ments from two-ports branches of these blocks are calculated from:

$$Z\_{\mathbb{R}} = \mathbb{R},\tag{9}$$

$$Z\_L = \mathbb{Z}\,\pi f \mathbf{L},\tag{10}$$

$$Z\_{\mathbb{C}} = \frac{1}{j2\pi f \mathbb{C}}.\tag{11}$$

where *f* is a signal frequency and *j* denotes imaginary unit. Impedances *Z*1, *Z*2, *Z*3 for termi‐ nal illustrated in Fig.3a can be determined according to equations (9), (10), (11) and its ele‐ ment sequence (i.e. *Z*1=*ZR*, *Z*2=*ZL*, *Z*3=*ZC* for RLC-typed leaf, *Z*1=*ZL*, *Z*2=*ZC*, *Z*3=*ZR* for LCRtyped leaf and *Z*1=*ZC*, *Z*2=*ZR*, *Z*3=*ZL* for CRL-typed leaf).

**Figure 3.** Terminal block structures: a) general, b) exemplary two-ports coded by the respective above strings

Finally, basing on elements from (2) and (3), the trees which code passive filters circuits can be constructed and some example is presented in Fig.4.

**Figure 4.** Exemplary structures of: a) genotype, b) phenotype

#### **3.1. Analysis of the blocks connections**

An important capability of the technique described is that it assures galvanic (direct) con‐ nection of the ground line between input and output ports and it is very desired for practi‐ cally implemented filter circuits. Besides, as can be seen in the analysis below, each type of the basic blocks configurations allows for defining an equivalent T-typed circuit and it uni‐ fies the interpretation of genotypes of any shapes and sizes.

**Figure 5.** Basic blocks connected in series and its equivalent circuit

2 , *Z fL <sup>L</sup>* = p

<sup>1</sup> , <sup>2</sup> *ZC j fC* p

where *f* is a signal frequency and *j* denotes imaginary unit. Impedances *Z*1, *Z*2, *Z*3 for termi‐ nal illustrated in Fig.3a can be determined according to equations (9), (10), (11) and its ele‐ ment sequence (i.e. *Z*1=*ZR*, *Z*2=*ZL*, *Z*3=*ZC* for RLC-typed leaf, *Z*1=*ZL*, *Z*2=*ZC*, *Z*3=*ZR* for LCR-

**Figure 3.** Terminal block structures: a) general, b) exemplary two-ports coded by the respective above strings

Finally, basing on elements from (2) and (3), the trees which code passive filters circuits can

An important capability of the technique described is that it assures galvanic (direct) con‐ nection of the ground line between input and output ports and it is very desired for practi‐ cally implemented filter circuits. Besides, as can be seen in the analysis below, each type of

typed leaf and *Z*1=*ZC*, *Z*2=*ZR*, *Z*3=*ZL* for CRL-typed leaf).

56 Analog Circuits

be constructed and some example is presented in Fig.4.

**Figure 4.** Exemplary structures of: a) genotype, b) phenotype

**3.1. Analysis of the blocks connections**

(10)

<sup>=</sup> (11)

Two blocks configured in a series (genotype node denoted as '--') are illustrated in Fig.5. A simple analysis of this circuit leads to equations defining the impedances for the equivalent circuits:

$$Z\_{1s} = Z\_{1a} + Z\_{1b},\tag{12}$$

$$Z\_{\mathfrak{Z}s} = Z\_{\mathfrak{Z}a} + Z\_{\mathfrak{Z}b},\tag{13}$$

$$\mathbf{Z}\_{\mathfrak{B}^{\mathsf{a}}} = \mathbf{Z}\_{\mathfrak{A}^{\mathsf{a}}} + \mathbf{Z}\_{\mathfrak{A}^{\mathsf{a}}}.\tag{14}$$

**Figure 6.** Basic blocks connected in cascade and its equivalent circuit

As can be seen in Fig.6, for a cascade connection (genotype node denoted as '>>') a conver‐ sion from a Π to a T-typed structure can be used for an equivalent circuit determination:

$$Z\_{1c} = Z\_{1a} + Z\_{p1} = Z\_{1a} + \frac{\left(Z\_{2a} + Z\_{1b}\right)Z\_{3a}}{Z\_{\text{Tlab}}},\tag{15}$$

$$Z\_{2\varepsilon} = Z\_{2b} + Z\_{r2} = Z\_{2b} + \frac{(Z\_{2a} + Z\_{1b})Z\_{3b}}{Z\_{\text{Tulb}}},\tag{16}$$

$$Z\_{3c} = \frac{Z\_{3d} Z\_{3b}}{Z\_{\text{T1ab}}},\tag{17}$$

where:

$$Z\_{\rm THo} = Z\_{\rm 2a} + Z\_{\rm 3a} + Z\_{\rm 1b} + Z\_{\rm 3b} \,. \tag{18}$$

The parallel configuration (genotype node denoted as '||') is the last kind of connection as‐ sumed in (2) and it can be easily modeled by an equivalent circuit after a few transformations ex‐ plained in Fig.7. Finally, impedances for an equivalent circuit can be calculated from relations:

$$Z\_{1p} = \frac{Z\_{1xy} Z\_{3xy}}{Z\_{1xy}},\tag{19}$$

$$Z\_{2p} = \frac{Z\_{2\text{ny}} Z\_{3\text{ny}}}{Z\_{\text{Tray}}},\tag{20}$$

$$Z\_{3p} = \frac{Z\_{1xy} Z\_{2xy}}{Z\_{1xy}},$$

where:

$$\frac{1}{Z\_{1xy}} = \frac{1}{Z\_{1x}} + \frac{1}{Z\_{1y}} = \frac{1}{Z\_{1x} + Z\_{3x} + \frac{Z\_{1y}Z\_{3x}}{Z\_{1x}}} + \frac{1}{Z\_{1b} + Z\_{3b} + \frac{Z\_{1b}Z\_{3b}}{Z\_{1b}}},\tag{22}$$

$$\frac{1}{Z\_{2xy}} = \frac{1}{Z\_{2x}} + \frac{1}{Z\_{2y}} = \frac{1}{Z\_{24} + Z\_{34} + \frac{Z\_{24}Z\_{34}}{Z\_{7x}}} + \frac{1}{Z\_{2b} + Z\_{3b} + \frac{Z\_{2b}Z\_{3b}}{Z\_{7b}}},\tag{23}$$

#### Memetic Method for Passive Filters Design http://dx.doi.org/10.5772/53716 59

$$\frac{1}{Z\_{3xy}} = \frac{1}{Z\_{3x}} + \frac{1}{Z\_{3y}} = \frac{1}{Z\_{14} + Z\_{24} + \frac{Z\_{14}Z\_{24}}{Z\_{14}}} + \frac{1}{Z\_{1b} + Z\_{2b} + \frac{Z\_{1b}Z\_{2b}}{Z\_{1b}}},\tag{24}$$

$$Z\_{\rm T1xy} = Z\_{\rm 1xy} + Z\_{\rm 1xy} + Z\_{\rm 1xy}.\tag{25}$$

**Figure 7.** Basic blocks connected in parallel and the equivalent circuit

As can be seen in Fig.6, for a cascade connection (genotype node denoted as '>>') a conver‐ sion from a Π to a T-typed structure can be used for an equivalent circuit determination:

1 1 11 , *a ba*

2 2 22 , *a bb*

3 3 <sup>3</sup> , *a b c*

*Z*P

*ab Z Z*

The parallel configuration (genotype node denoted as '||') is the last kind of connection as‐ sumed in (2) and it can be easily modeled by an equivalent circuit after a few transformations ex‐ plained in Fig.7. Finally, impedances for an equivalent circuit can be calculated from relations:

> 1 3 <sup>1</sup> , *xy xy p*

2 3 <sup>2</sup> , *xy xy p*

*Z*P

1 2 <sup>3</sup> , *xy xy p*

*Z*P

1 11 1 3 1 3 1 3 1 3 1 11 1 <sup>1</sup> , *xy x y a a b b a a b b*

2 22 2 3 2 3 2 3 2 3 1 11 1 <sup>1</sup> , *xy x y a a b b a a b b*

*Z ZZ Z Z Z Z ZZ ZZ*

*Z ZZ Z Z Z Z ZZ ZZ*

*Ta Tb*

*Ta Tb*

*Z Z*

*Z Z*

*xy Z Z*

*xy Z Z*

*Z*P

*Z*

*Z*

*Z*

=+= +

=+= +

*xy Z Z*

*c at a*

*c bt b*

*Z*

*ZZZZ*

where:

58 Analog Circuits

where:

*ZZZZ*

( ) 2 13

*Z*P

( ) 2 13

*Z*P

*ab Z ZZ*

*ab Z ZZ*

<sup>+</sup> =+=+ (15)

<sup>+</sup> =+=+ (16)

<sup>2313</sup> . *Z ZZZZ* <sup>P</sup>*ab a a b b* =+++ (18)

= (17)

= (19)

= (20)

= (21)

++ ++ (22)

++ ++ (23)

The transformations described above allow to obtain the resultant impedances *Z*1*<sup>r</sup>*, *Z*2*<sup>r</sup>* and *Z*3*<sup>r</sup>* of the filter circuit (Fig.8) coded by a genotype tree. Finally, a frequency response of the filter loaded by the impedance *Zo* can be calculated from:

$$K = \frac{\mathcal{U}\_{out}}{\mathcal{U}\_{\text{ht}}} = \frac{Z\_o}{\left(Z\_o + Z\_{2r}\right)\left(1 + \frac{Z\_{1r}}{Z\_{3r}}\right) + Z\_{1r}}.\tag{26}$$

**Figure 8.** The resultant equivalent filter circuit

#### **4. Genetic operations and fitness calculation**

#### **4.1. Reproduction and crossover**

The first genetic operation executed after primary generation evaluation is reproduction and it completes a mating pool. A rang method of reproduction is applied in the proposed solu‐ tion. This kind of operation assures the minimization of the evolutionary system's tendency to promote the average fitted phenotypes (selection pressure regulation). The probability *Psel* of an individual selection for an intermediate pool of candidates for recombination depends on rang *r* value:

$$P\_{sol} = P\_{mn} + \left(1 - P\_{mn}\right)\frac{r}{N\_{max}},\tag{27}$$

where *Pmn* denotes an assumed minimal value of the resultant probability (to avoid the cur‐ rently most fitted genotypes domination), and *Nmx* is a total number of individuals in the population. The rang *r* is equal to the position of an individual in population ordered by a fitness value (i.e. for the worst *r*=1 and for the best evaluated *r*=*Nmx*).

The crossover (recombination of the genetic material) process is carried out with probability *Pcr* for two parental genotypes selected randomly (with uniform probability) from a mating pool and its idea is illustrated in Fig.9. After two crossover points CP1 and CP2 random designation for genotypes (individually for each one), adequate sub-trees are exchanged between them. This process allows for preserving and propagating the genetic information of the well fitted individuals inside the generation and this operation promises to obtain better evaluated off‐ spring genotypes. Additionally, for the proposed system, the protection against too much tree growing was applied to crossover process. Recombination is accepted only if for each offspring individual the maximal depth does not exceed the assumed *Dmx* value.

**Figure 9.** The idea of the crossover process

#### **4.2. Genotypes mutations**

**Figure 8.** The resultant equivalent filter circuit

**4.1. Reproduction and crossover**

on rang *r* value:

60 Analog Circuits

**4. Genetic operations and fitness calculation**

The first genetic operation executed after primary generation evaluation is reproduction and it completes a mating pool. A rang method of reproduction is applied in the proposed solu‐ tion. This kind of operation assures the minimization of the evolutionary system's tendency to promote the average fitted phenotypes (selection pressure regulation). The probability *Psel* of an individual selection for an intermediate pool of candidates for recombination depends

( ) 1 , *sel mn mn*

where *Pmn* denotes an assumed minimal value of the resultant probability (to avoid the cur‐ rently most fitted genotypes domination), and *Nmx* is a total number of individuals in the population. The rang *r* is equal to the position of an individual in population ordered by a

The crossover (recombination of the genetic material) process is carried out with probability *Pcr* for two parental genotypes selected randomly (with uniform probability) from a mating pool and its idea is illustrated in Fig.9. After two crossover points CP1 and CP2 random designation for genotypes (individually for each one), adequate sub-trees are exchanged between them. This process allows for preserving and propagating the genetic information of the well fitted individuals inside the generation and this operation promises to obtain better evaluated off‐ spring genotypes. Additionally, for the proposed system, the protection against too much tree growing was applied to crossover process. Recombination is accepted only if for each offspring

*<sup>r</sup> PP P*

fitness value (i.e. for the worst *r*=1 and for the best evaluated *r*=*Nmx*).

individual the maximal depth does not exceed the assumed *Dmx* value.

*mx*

*<sup>N</sup>* = +- (27)

The main goal of the mutation operations is the protection against the optimization process stagnation in the local area of the searching space. Four types of mutation procedure are used in the system described:


**Figure 10.** Genotype mutations illustration

The ideas of all these genotype modifications are illustrated in Fig.10. In case of node muta‐ tion, it is replaced by a random one (function or terminal adequately) selected from the ini‐ tially assumed set (2) or (3). During the sub-tree mutation process it is deleted or replaced by a randomly created one. Sub-tree cutting allows for simplifying the phenotype, so it is de‐ sired to the filter circuit size minimization. On the other hand, the new sub-tree adding can lead to a genotype growing, so this process is controlled and only results with the maximal depth up to *Dmx* are accepted.

#### **4.3. Fitness value calculation**

The quality of the phenotype is evaluated by means of fitness value calculation:

$$Q = \sqrt{\frac{\sum\_{m=1}^{M} \left( \left| \mathbb{K}\_{m} \right| - \left| \mathbb{K}\_{0m} \right| \right)^{2}}{w\_{a} \Big|\_{a=1}^{A}}} + \max \left( \frac{\left( \left| \mathbb{K}\_{m} \right| - \left| \mathbb{K}\_{0m} \right| \right) \Big|\_{m=1}^{M}}{w\_{a} \Big|\_{a=1}^{A}} \right). \tag{28}$$

Equation (28) is a sum of two components. The first one is a mean square distance of the |*Km*| amplitude response of the tested phenotype from the desired region of the filter de‐ sign specifications |*K*0*<sup>m</sup>*|and it is calculated in *M* points totally of frequency response for the assumed range of analysis. Besides, to assure the same impact on optimization for each frequency band *a* defined in filter specifications, its distance is averaged by division by the number *wa* of frequency points included in the actually analyzed band. The sec‐ ond one is the average value of the maximal deviation of the frequency point for bands and it additionally prevents the system from stagnation at average solutions. The memet‐ ic optimization system minimizes *Q* during evolutionary cycles and for the best pheno‐ type this distance reaches zero.

After offspring individuals creation its terminal nodes (i.e. filter elements values) are searched with the use of HJM in the way described in the next section and after collecting the new generation it replaces the previous one during the succession process. Finally, the best fitted genotype codes the structure of the filter circuit.

### **5. Local search**

#### **5.1. Hooke and Jeeves local search algorithm**

The pattern search optimization algorithm consists of two types of moves, i.e. exploratory move and pattern move.

The purpose of the first of them is to find and utilize information about the optimized function values around the current base point **b(***k=0***)** (individual). In increasing order of indexes, each of the following circuit parameters' values is given a small increment (first‐ ly in positive and then, if required, in the negative direction). The value of the fitness function is checked and if there is a progress noticed, the value of this variable will be kept as a new **b(***k+1***)** vector. This step is being repeated with reduced range of the param‐ eters values change. When no fitness function value improvement is possible, a pattern move, starting from the current point, is made.

The aim of the pattern move is to speed up the search by information gained about the optimized function and to find the best search direction. A move from the current **b(***k+1***)** in the direction given with **b(***k+1***)** *-***b***k* is made. The fitness function value is calculated in the point given by

$$\mathbf{p}\_k = \mathbf{b}\_{\binom{k+1}{k+1}} - \mathbf{2}\mathbf{b}\_k. \tag{29}$$

The procedure is continued with a new sequence of exploratory moves starting from the point **p***k*. If the achieved fitness function's value is lower than the one around point **b***k*, then a new base point **b(***k+2***)** has been found. In this case a new pattern move (29) is carried out. Otherwise, the pattern move from **b(***k+1***)** is dropped and the new exploratory sequence starts from **b(***k+1***)** . The procedure is continued until the length of the step for each of the variables has been reduced to the value assumed previously.

#### **5.2. The implementation of the Hooke-Jeeves method**

The ideas of all these genotype modifications are illustrated in Fig.10. In case of node muta‐ tion, it is replaced by a random one (function or terminal adequately) selected from the ini‐ tially assumed set (2) or (3). During the sub-tree mutation process it is deleted or replaced by a randomly created one. Sub-tree cutting allows for simplifying the phenotype, so it is de‐ sired to the filter circuit size minimization. On the other hand, the new sub-tree adding can lead to a genotype growing, so this process is controlled and only results with the maximal

The quality of the phenotype is evaluated by means of fitness value calculation:

*Q max*

å

best fitted genotype codes the structure of the filter circuit.

**5.1. Hooke and Jeeves local search algorithm**

( ) ( ) <sup>2</sup> <sup>0</sup> <sup>0</sup> 1 1 1 1

Equation (28) is a sum of two components. The first one is a mean square distance of the |*Km*| amplitude response of the tested phenotype from the desired region of the filter de‐ sign specifications |*K*0*<sup>m</sup>*|and it is calculated in *M* points totally of frequency response for the assumed range of analysis. Besides, to assure the same impact on optimization for each frequency band *a* defined in filter specifications, its distance is averaged by division by the number *wa* of frequency points included in the actually analyzed band. The sec‐ ond one is the average value of the maximal deviation of the frequency point for bands and it additionally prevents the system from stagnation at average solutions. The memet‐ ic optimization system minimizes *Q* during evolutionary cycles and for the best pheno‐

After offspring individuals creation its terminal nodes (i.e. filter elements values) are searched with the use of HJM in the way described in the next section and after collecting the new generation it replaces the previous one during the succession process. Finally, the

The pattern search optimization algorithm consists of two types of moves, i.e. exploratory

The purpose of the first of them is to find and utilize information about the optimized function values around the current base point **b(***k=0***)** (individual). In increasing order of indexes, each of the following circuit parameters' values is given a small increment (first‐ ly in positive and then, if required, in the negative direction). The value of the fitness

*K K K K*


*w w* = = = =

*<sup>M</sup> <sup>M</sup> m m m m m m A A a a a a*

.

(28)

è ø

depth up to *Dmx* are accepted.

62 Analog Circuits

**4.3. Fitness value calculation**

type this distance reaches zero.

**5. Local search**

move and pattern move.

The search space in the presented research is not only discrete but also limited by the E24 sequence. Considering it several correcting procedures needed to be implemented.

#### *5.2.1. Increasing and decreasing values of the circuit parameters*

The first problem to be addressed was to properly increase the values of the circuit parame‐ ters. It has been achieved by translating the **V** and **O** sets into two vectors of integers:

$$\mathbf{V}\_{\mathbf{V}} = \{1, \dots, 24\},\tag{50}$$

$$\mathbf{V\_{O}} = \{1, \dots, 7\},\tag{31}$$

where: **Vv**(1) denotes v1=10, **Vv**(2) denotes v2=11, etc. **Vo**(1) denotes *o*1, **Vo**(2) denotes o2, etc.

Let stepk be the current size of the variable *vari* change. The *i-th* variable is given with:

$$
tau\_l = \left(\upsilon\_v^\vee, \upsilon\_o^\vee\right),\tag{32}$$

where *vv i* is the translated value of ve and *vo i* is the translated value of ou (see (3)-(5)).

The increment of the vari is carried out according to the presented algorithm:


The decrement procedure is carried in a similar way, i.e.:


if *vo i* <1 then *vv i* =11∧*vo i* =0

*5.2.2. The "in-loop" local search*

The local search procedure implemented within the evolutionary loop has been implement‐ ed in a way that the computation time was not affected significantly.

First of all, there has been a strict limit of pattern searches number set. Secondly, the initial step size has been small. It allows for exploring only a small area around the base point giv‐ en with following individuals.

### **6. Exemplary Results**

For the efficiency of the technique presentation, an automated design of the specialized filter was realized. The assumed, desired amplitude response specifications for a filter are as fol‐ lows (*A*=7 bands):

$$\begin{array}{lcll} & \left| \left[ K\_{0} \right] \in \{ -\infty, -20 \} \text{dB} & \wedge & f = \{ 1e-7, 0.1 \} \text{MHz} \\ \left| K\_{0} \right| \in \{ -\infty, 0 \} \text{dB} & \wedge & f = \{ 0.1, 0.325 \} \text{MHz} \\ \left| K\_{0} \right| \in \{ -3, 0 \} \text{dB} & \wedge & f = \{ 0.325, 0.375 \} \text{MHz} \\ \left| K\_{0} \right| \in \{ -\infty, 0 \} \text{dB} & \wedge & f = \{ 0.375, 0.6 \} \text{MHz} \\ \left| K\_{0} \right| \in \{ -30, -20 \} \text{dB} & \wedge & f = \{ 0.6, 0.7 \} \text{MHz} \\ \left| K\_{0} \right| \in \{ -\infty, 0 \} \text{dB} & \wedge & f = \{ 0.7, 0.9 \} \text{MHz} \\ \left| K\_{0} \right| \in \{ -\infty, -40 \} \text{dB} & \wedge & f = \{ 0.9, 1 \} \text{MHz} \end{array} \tag{33}$$

for an assumed load resistance *Z*0=*R*0=150. The MGP system was executed with the initial parameters:


**•** *vv i* =*vv*

64 Analog Circuits

**•** if *vv i*

**•** *vv i* =*vv i* −24

**•** *vo i* =*vo <sup>i</sup>* + 1

**•** if *vo i*

**•** *vv i* =*vv i* −*ste pk*

**•** if *vv i*

**•** *vv i* =*vv <sup>i</sup>* + 24

**•** *vo i* =*vo i* −1

if *vo i* *<sup>i</sup>* <sup>+</sup> *ste pk*

>24 then:

>7 then *vv*

<1 then:

<1 then *vv*

*i* =11∧*vo i* =0

*5.2.2. The "in-loop" local search*

en with following individuals.

**6. Exemplary Results**

lows (*A*=7 bands):

*i* =24∧*vo i* =7

The decrement procedure is carried in a similar way, i.e.:

The local search procedure implemented within the evolutionary loop has been implement‐

First of all, there has been a strict limit of pattern searches number set. Secondly, the initial step size has been small. It allows for exploring only a small area around the base point giv‐

For the efficiency of the technique presentation, an automated design of the specialized filter was realized. The assumed, desired amplitude response specifications for a filter are as fol‐

> ( ( ( ( ( (

(33)

, 20 dB 1 7 ,0.1 MHz ,0 dB 0.1,0.325 MHz 3,0 dB 0.325,0.375 MHz ,0 dB 0.375,0.6 MHz , 30, 20 dB 0.6 ,0.7 MHz ,0 dB 0.7 ,0.9 MHz , 40 dB 0.9 ,1 MHz

*K f e K f K f K f K f K f K f*

<sup>ì</sup> Î -¥ - Ù = - <sup>ï</sup> ï Î -¥ Ù = <sup>ï</sup> Î - Ù = <sup>ï</sup>

í Î -¥ Ù = <sup>ï</sup> Î- - Ù = <sup>ï</sup> <sup>ï</sup> Î -¥ Ù = <sup>ï</sup> <sup>ï</sup> Î -¥ - Ù = <sup>î</sup>

ed in a way that the computation time was not affected significantly.

ï


and for the orders (5) of values for searched elements:

$$\mathbf{O} = \{0, 1, \dots, 6\}\ \ . \tag{34}$$

The best found genotype evaluated for the specifications (33) is illustrated in Fig.11 and it codes the temporary version of the filter circuit placed in Fig.12a. Next, due to radically small or high values for some elements and due to the obtained kinds of connections (e.g. shortened or opened branches, the same types of elements connected in series or in parallel), unnecessary elements from this circuit can be removed or simplified. This stage of the circuit synthesis can be supported by a simulation software tool. In Fig.12b the final filter circuit obtained for the considered example is presented. In comparison to the temporary one, its size was radically reduced without amplitude response degradation.

**Figure 11.** The resultant genotype

**Figure 12.** The structures of the designed filters: a) temporary version, b) final version

The quality of the specifications (33) keeping can be seen in Fig.13. The rectangles defining the allowed region for the filter amplitude response and simulation results obtained for the resul‐ tant circuit from Fig.12b are placed there. Only one restriction from (33) (the last frequency band) was a little violated, but the other ones are fully fulfilled. It should be emphasized, that due to discretizaton of the available values of components to a practical E24 series (4) not all theoretical shapes of frequency responses assumed on filter specifications defining stage will be practically reachable and it is necessary to conciliate with some inaccuracies. Besides, too high tolerance decreasing leads to the undesirable growing of production costs.

**Figure 13.** The amplitude response of the designed filter

## **7. Conclusions**

The automated system for a passive filter circuits design was presented in this chapter. Any initial information about the filter circuit structure is not necessary for filter synthesis, only desired specifications should be defined. The circuit's topology as well as its elements val‐ ues are optimized together in the MGP system. Thanks to the deterministic algorithm of the local searching engaging (HJM), the speed of convergence to the well evaluated solutions during the evolutionary computations grows significantly and the values of the filter's ele‐ ments are adjusted to the most fitted ones for an actual circuit topology. Finally, after redun‐ dant elements elimination, the filter circuit is obtained with components selected from a practical production series (i.e. for practically reachable nominal values), and it makes the circuit realization easier. Besides, the genotypes sub-trees deletion applied in the system en‐ ables the complexity of the circuit minimization.

For the future system improving, the additional criteria for fitness calculation (28) can be added. For example, the group delay value and the kind of circuit elements (the quantity of inductors minimization) may be optimized. Besides, the applying of the models of real ele‐ ments should be considered (to the parasitic parameters taken into account). Additionally, the proposed system can be adapted to the active circuits automated synthesis.

## **Author details**

**Figure 12.** The structures of the designed filters: a) temporary version, b) final version

66 Analog Circuits

high tolerance decreasing leads to the undesirable growing of production costs.

**Figure 13.** The amplitude response of the designed filter

The quality of the specifications (33) keeping can be seen in Fig.13. The rectangles defining the allowed region for the filter amplitude response and simulation results obtained for the resul‐ tant circuit from Fig.12b are placed there. Only one restriction from (33) (the last frequency band) was a little violated, but the other ones are fully fulfilled. It should be emphasized, that due to discretizaton of the available values of components to a practical E24 series (4) not all theoretical shapes of frequency responses assumed on filter specifications defining stage will be practically reachable and it is necessary to conciliate with some inaccuracies. Besides, too

Tomasz Golonek\* and Jantos Piotr

\*Address all correspondence to: tgolonek@polsl.pl

Silesian University of Technology, Poland

## **References**


## **Interval Methods for Analog Circuits**

## Zygmunt Garczarczyk

[5] Tathagato Rai Dastidar, P. P. Chakrabarti, Partha Ray, "A Synthesis System for Ana‐ log Circuits Based on Evolutionary Search and Topological Reuse", IEEE Transac‐

[6] Budzisz H., "Evolutional searching for circuit structures", Electronic Letters, Vol. 34,

[7] Golonek T., Jantos P., Rutkowski J., Stimulus with limited band optimization for ana‐ logue circuit testing, Metrology and Measurement Systems, Vol. XIX, No. 1, pp.

[8] Jantos P., Rutkowski J., Evolutionary methods to analogue electronic circuits yield optimisation, Bulletin of the Polish Academy of Sciences - Technical Sciences, Vol. 56,

[9] Golonek T., Rutkowski J., Genetic-Algorithm-Based Method for Optimal Analog Test Points Selection, IEEE Trans. on Cir. and Syst.-II., Vol.54, No.2, 2007, pp. 117-121.

[10] Golonek T., Grzechca D., Rutkowski J., Application of Genetic Programming to Edge Decoder Design, Proc. of the Inter. Symposium on Cir. and Sys., ISCAS 2006, Greece,

[11] Hooke R. , Jeeves T. A., Direct search solution of numerical and statistical problems,

[12] Moscato P., "On Evolution, Search, Optimization, Genetic Algorithms and Martial Arts: Towards Memetic Algorithms", Caltech Concurrent Computation Program (re‐

[13] Land M. W. S., Evolutionary Algorithms with Local Search for Combinatorial Opti‐

[14] Ozcan E., "Memes, Self-generation and Nurse Rostering". Lecture Notes in Computer Science. Lecture Notes in Computer Science, Springer-Verlag, 3867, 2007, 85–104.

[15] Chen X., Ong Y. S., Lim M. H., Tan K. C., "A Multi-Facet Survey on Memetic Compu‐ tation", Evolutionary Computation, IEEE Transactions on , vol.15, no.5, 2011, pp.

[16] Ong Y. S., Lim M. H., Zhu N., Wong K. W., "Classification of adaptive memetic algo‐ rithms: a comparative study," Systems, Man, and Cybernetics, Part B: Cybernetics,

Assoc. Computing Machinery J., Vol.8 , No.2, 1960, 212-229.

IEEE Transactions on , vol.36, no.1, 2006, pp. 141-152.

tions on Evolutionary Computation, Vol. 9, No. 2, 2005, 211-224.

No. 16, 1998, 1543-1545.

Issue 1, 2008, pp. 9-16.

73-84, 2012.

68 Analog Circuits

4683-4686.

port 826), 1989.

mization, 1998.

591-607.

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/54848

## **1. Introduction**

Concepts of interval analysis are suitable for solving some problems for linear and nonlinear analog circuits. The effects of circuit parameter uncertainties on circuit performance are always of great concern to the system designer. It is desirable to know a priori estimates of the circuit response, subject to parameter uncertainties. In the first section of this chapter a problem of calculating of the operating regions (solutions) for linear circuits with parameters done as interval numbers is considered. Such circuits are described by linear interval equations. The set of all possible operating points of the circuit may have a very complicated structure and it is usually impractical to calculate it. Applying ideas of interval analysis it is possible to calculate multidimensional rectangular region bounding the set of operating points. In the paper an algorithm of iterative evaluation of the bounds of operating regions is applied. The second section deals with the finding DC solutions of nonlinear, inertialess circuits. Using idea of a continuation method one can solve nonlinear equations describing circuit by tracing so-called solution path. A version of the predictor-corrector method for computing points of continua‐ tion path of a nonlinear equation is presented. The target of this approach is to control the corrector step in such a manner that the predictor step is sufficiently large but that the corrector step does not jump to another continuation path. Additionally interval analysis offers possi‐ bility to solve that problem by applying generalized bisection method based on some interval operator associated with nonlinear equation. In the section Krawczyk operator is used in ndimensional box-searching of all solutions. For both sections numerical studies are reported in order to illustrate and verify presented approach.

## **2. Analysis of linear circuits with interval data**

The changes in the performance of a linear circuit due to the variations in circuit parameters are of great practical importance in engineering analysis and design. Such perturbations of

© 2013 Garczarczyk; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Garczarczyk; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

parameters may model the effects of uncertainties in manufacturing tolerances, ageing of components, environmental causes and the like. because of these uncertainties, the value of the parameters of a given linear circuit may frequently be treated as belonging to suitable intervals.

From this point of view, the aim of this paper is to consider solutions of linear circuits using paradigm of interval computations [1]. This problem is closely related to the tolerance analysis of electric and electronic circuits especially the worst-case analysis [3].

To explain problem formulation let us consider following simple example. We are interested in computation of nodal voltages in the R-ladder circuit of Fig.1.

Parameters of the circuit are done as interval number, viz. E∈[5.67, 6.93], R1, R2, R4, R6∈[0.09, 0.11], R3, R5∈[1.8, 2.2]. Applying nodal analysis we obtain following system of equations:

$$
\begin{bmatrix}
\text{[1.98, 2.42]} & \text{[-2.2, -1.8]} & \text{[0, 0]} \\
\text{[-2.2, -1.8]} & \text{[3.69, 4.51]} & \text{[-2.2, -1.8]} \\
\text{[0, 0]} & \text{[-2.2, -1.8]} & \text{[1.89, 2.31]} \\
\end{bmatrix}
\begin{bmatrix}
V\_1 \\ V\_2 \\ V\_3 \\
\end{bmatrix} = \begin{bmatrix}
\text{[0.51, 0.76]} \\
\text{[0, 0]} \\
\text{[0, 0]} \\
\end{bmatrix}
$$

**Figure 1.** R-ladder circuit

It's seen that the coefficients and right-hand sides of the system are not determined exactly, but are only known to lie within some real intervals.

Denote a system of linear interval equations by

$$\mathbf{A} \cdot \mathbf{x} = \mathbf{b}.\tag{1}$$

The coefficients and right-hand sides of the system are not determined exactly, but are only known to lie within some real intervals. Such a system of linear interval equations represents a family of ordinary linear systems which can be obtained from it by fixing coefficients and right-hand-side values in the prescribed intervals. Each of these systems, under the assumption that each A∈**A** is nonsingular, has a unique solution, and all these solutions constitute a socalled solution set S. The solution set of Eq.(1) can be expressed as

$$\mathbf{S} = \{ \mathbf{x}; \mathbf{A} \mathbf{x} \mathbf{=} \mathbf{b}, \ A \in \mathbf{A}, \ \mathbf{b} \in \mathbf{b} \}. \tag{2}$$

It forms some multidimensional operating region of a circuit.

It the next points, we introduce some notations and ideas concerned with a structures of operating regions, drawbacks of some circuit descriptions, and method of rectangular (interval) evaluating of the operating regions.

As the examples, we give the calculation of variations of the frequency response and nodal voltages of the RC circuits with interval data. We also test our approach when the widths of matrix elements in Eq.(1) are large.

#### **2.1. Notations and preliminaries**

parameters may model the effects of uncertainties in manufacturing tolerances, ageing of components, environmental causes and the like. because of these uncertainties, the value of the parameters of a given linear circuit may frequently be treated as belonging to suitable

From this point of view, the aim of this paper is to consider solutions of linear circuits using paradigm of interval computations [1]. This problem is closely related to the tolerance analysis

To explain problem formulation let us consider following simple example. We are interested

Parameters of the circuit are done as interval number, viz. E∈[5.67, 6.93], R1, R2, R4, R6∈[0.09, 0.11], R3, R5∈[1.8, 2.2]. Applying nodal analysis we obtain following system of equations:

It's seen that the coefficients and right-hand sides of the system are not determined exactly,

The coefficients and right-hand sides of the system are not determined exactly, but are only known to lie within some real intervals. Such a system of linear interval equations represents a family of ordinary linear systems which can be obtained from it by fixing coefficients and right-hand-side values in the prescribed intervals. Each of these systems, under the assumption that each A∈**A** is nonsingular, has a unique solution, and all these solutions constitute a so-

<sup>R</sup> <sup>3</sup> <sup>V</sup> <sup>R</sup> <sup>5</sup> <sup>2</sup>

R <sup>2</sup> R <sup>4</sup> R <sup>6</sup>

*V*1 *V*2 *V*3 =

V <sup>3</sup>

**Ax b** . = (1)

S = x:Ax=b, A A, { Î Î }b b . (2)

0.51, 0.76 0, 0 0, 0

of electric and electronic circuits especially the worst-case analysis [3].

1.98, 2.42 −2.2, −1.8 0, 0 −2.2, −1.8 3.69, 4.51 −2.2, −1.8 0, 0 −2.2, −1.8 1.89, 2.31

<sup>R</sup> <sup>1</sup> <sup>V</sup> <sup>1</sup>

in computation of nodal voltages in the R-ladder circuit of Fig.1.

E

but are only known to lie within some real intervals.

called solution set S. The solution set of Eq.(1) can be expressed as

Denote a system of linear interval equations by

**Figure 1.** R-ladder circuit

intervals.

70 Analog Circuits

All vectors in this paper will have n components and all matrices will be of size n×n. The sets of real vectors, real matrices, interval vectors, and interval matrices are represented by the lower case, upper case, lower case bold, and upper case bold letters respectively. A real, scalar interval **x** is given by *x* ¯ , *x*¯ , where the endpoints of an interval are *x* ¯ <sup>≤</sup> *<sup>x</sup>*¯. We denote the center of an interval **x** by *m*(**x**)=(*x* ¯ <sup>+</sup> *<sup>x</sup>*¯)/2, and width of an interval by *w*(**x**)= *<sup>x</sup>*¯ <sup>−</sup> *<sup>x</sup>* ¯ . For interval vectors and matrices these concepts are defined via the elements. If **A**=(**a**ij) is an interval matrix, then m(**A**) = (m(**a**ij)), see for example [2]. We say that a real vector x, is contained in an interval vector **x**, and we write x∈**x**, if *x* ¯*<sup>i</sup>* <sup>≤</sup> *xi* <sup>≤</sup> *<sup>x</sup>*¯*<sup>i</sup>* for all i = 1,...,n. A real matrix A is contained in an interval matrix **A**, and we write A∈**A**, if *a* ¯*ij* <sup>≤</sup>*aij* <sup>≤</sup>*a*¯*ij* for all i, j = 1,...,n. An interval matrix **A** is called regular if detA ≠ 0 for each A∈**A**. Let | *A*| =( |*aij* | ) denote matrix with absolute values of elements and A ≥ 0 (and similar relations) are meant componentwise (i.e. aij ≥ 0). An interval matrix **A** is called inverse stable if | *A*−<sup>1</sup> | >0 for each A∈**A**, i.e. if each inverse matrix element is nonzero [12]. Subsequently we introduce some concepts of matrix theory [11]. First, we use the (componentwise) natural partial ordering on sets of real rectors and matrices:

$$A \succeq B \iff a\_{ij} \ge b\_{ij} \text{ for } \text{ i. } \ j = 1, \dots, m$$

For inverse stable matrix A the signature matrix Z is defined by

$$z\_{ij} = \begin{cases} 1 & \text{if } \quad a\_{ji} \ge 0 \\ -1 & \text{if } \quad a\_{ji} \le 0 \end{cases} \quad i, \quad j = 1, \dots, m$$

(notice the transposition of indices). For two n×n matrices A = (aij) and B = (bij) their compo‐ nentwise product is defined as A\*B=(aijbij). Further, diag A denote the diagonal vector of A, i.e. diag A = (a11,...,ann) T and *ρ*(*A*)denotes spectral radius of A.

#### **2.2. Structure of an operating region**

In the first section we introduce a linear interval **Ax** = **b** with an interval matrix of coefficients **A** = [A-Δ, A+Δ] and a right-hand side interval vector **b** = [b-δ, b+δ] describing the linear circuit. Here A = m(**A**), Δ = w(**A**)/2 and b = m(**b**), δ = w(**b**)/2. The set S of all possible operating points (solutions of Eq.(1)) of linear circuit may have a very complicated structure. This set is not generally an interval vector.

A number of authors have described ways to compute an x satisfying (2). See [12] for reference. In their pioneer work Oettli and Prager [8] showed that the solution set of a linear interval equation with regular matrix A is described following

$$\mathbf{x} \in \mathbb{S} \iff \left| A\mathbf{x} - b \right| \le \Delta \left| \mathbf{x} \right| + \delta \tag{3}$$

It is known [12] that S can be represented as a union of at most 2n convex polyhedra. The intersection of S with each orthant of Rn is convex. In general, S itself is not convex.

EXAMPLE 1. To illustrate that the solution set S given by (3) is not simple, let consider a linear system with input-output relationship written as y(jω) = K(jω,p) x(jω).

Frequency response K(jω,p) vary with vector p of some parameters ranging in known intervals. It can be shown [9] that for a fixed frequency changes of K(jω,p) are described by the system of two interval equations

$$
\begin{bmatrix} \begin{bmatrix} a \ \ b \end{bmatrix} & -\begin{bmatrix} c \ \ d \end{bmatrix} \end{bmatrix} \begin{bmatrix} y\_1 \\ y\_2 \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} 1 \ \end{bmatrix} & \begin{bmatrix} 1 \end{bmatrix} \\\\ \begin{bmatrix} a \ \end{bmatrix} & \begin{bmatrix} a \end{bmatrix} \end{bmatrix} \begin{bmatrix} y\_1 \\ y\_2 \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} 1 \ \end{bmatrix} & \begin{bmatrix} 1 \\ 0 \ \end{bmatrix} \end{bmatrix}
$$

where the intervals [a,b] and [c,d] represent the ranges of values of Re{*K*( *jω*, *p*) <sup>−</sup>1} and Im{*K*( *jω*, *p*) <sup>−</sup>1}, respectively, and K(jω,p) = y1 + jy2. In Fig.2 is presented the region of frequency response changes caused by variations of system parameters done by intervals [a,b] = [1,2] and [c,d] = [-0.5,0.5]. Coordinates of the points fixing the shape of the region are following: A(4/9, 2/9), B(4/3, 2/3), C(1, 0), D(4/3, -2/3), E(4/9, -2/9), and F(1/2, 0).

**Figure 2.** The region of changes of frequency response for a fixed frequency

Because S is generally so complicated in shape, it is usually impractical to try to use it. Instead, it is common practice to seek the interval vector x containing S which has the narrowest possible interval components. We say we "solve" the problem when we find x.

#### **2.3. Rectangular evaluation of an operating region**

A number of authors have described ways to compute an x satisfying (2). See [12] for reference. In their pioneer work Oettli and Prager [8] showed that the solution set of a linear interval

It is known [12] that S can be represented as a union of at most 2n convex polyhedra. The

EXAMPLE 1. To illustrate that the solution set S given by (3) is not simple, let consider a linear

Frequency response K(jω,p) vary with vector p of some parameters ranging in known intervals. It can be shown [9] that for a fixed frequency changes of K(jω,p) are described by the system

> *y*1 *y*2

<sup>−</sup>1}, respectively, and K(jω,p) = y1 + jy2. In Fig.2 is presented the region of frequency

B

C

<sup>F</sup> <sup>y</sup> <sup>1</sup>

Because S is generally so complicated in shape, it is usually impractical to try to use it. Instead, it is common practice to seek the interval vector x containing S which has the narrowest possible

D

<sup>=</sup> 1, 1 0

d

(3)

<sup>−</sup>1} and

*x S Ax b x* Î Û - £D +

intersection of S with each orthant of Rn is convex. In general, S itself is not convex.

*a*, *b* − *c*, *d c*, *d a*, *b*

where the intervals [a,b] and [c,d] represent the ranges of values of Re{*K*( *jω*, *p*)

response changes caused by variations of system parameters done by intervals [a,b] = [1,2] and [c,d] = [-0.5,0.5]. Coordinates of the points fixing the shape of the region are following: A(4/9,

system with input-output relationship written as y(jω) = K(jω,p) x(jω).

2/9), B(4/3, 2/3), C(1, 0), D(4/3, -2/3), E(4/9, -2/9), and F(1/2, 0).

y 2

A

E

interval components. We say we "solve" the problem when we find x.

**Figure 2.** The region of changes of frequency response for a fixed frequency

equation with regular matrix A is described following

of two interval equations

Im{*K*( *jω*, *p*)

72 Analog Circuits

From the standpoint of view of interval computations the formulation of circuit equations based of mesh analysis or modal analysis are not necessarily suitable for the networks with interval parameters. They result in a system of linear interval equations whose coefficients are not independent. We have the wider range of the coefficients since the elements of meshimpedance matrix or mode-admittance matrix are given by the linear combinations of the interval parameters of resistors, capacitors, etc. Hence the width of the elements of the matrix coefficients and right-hand side vector of Eq.(1) becomes resultantly larger. This makes the evaluation of interval solution worse. Secondly owing to the interval dependency the linear combination of interval numbers gives us possibility to have the meaningless combination of the parameters. In order to avoid the drawbacks mentioned above the formulation by hybrid equation is well suited for solving linear interval systems [4], [10].

The problem treated in this section is how to compute the interval hull of S, i.e. the interval vector x with components

$$\begin{array}{l} \underline{\mathbf{x}}\_{i} = \min \{ \mathbf{x}\_{i}, \mathbf{x} \in \mathcal{S} \} \\ \underline{\mathbf{x}}\_{i} = \max \{ \mathbf{x}\_{i}, \mathbf{x} \in \mathcal{S} \} \end{array} \qquad i = 1, \ldots, n \tag{4}$$

describing the exact range of components of the solution if coefficients and right-hand sides of Eq.(1) are allowed to vary in the given intervals.

Method of computing the vector **x** = *x* ¯ , *x*¯ described by (4) is based on the results of Rohn [5], [6] slightly modified. Namely, the vectors of left endpoints and right endpoints of the interval (rectangular) evaluation of the interval hull are computed iteratively by

$$\underline{X}^{k+1} = A^{-1} \left( B - Z \ast \left( \Lambda \left| \underline{X}^k \right| + D \right) \right), \quad \underline{\mathbf{x}} = \operatorname{diag} \underline{X}$$

$$\overline{X}^{k+1} = A^{-1} \left( B + Z \ast \left( \Lambda \left| \underline{X}^k \right| + D \right) \right), \quad \overline{\mathbf{x}} = \operatorname{diag} \overline{X} \tag{5}$$

$$k = 0, 1, \ldots$$

Where B=beT, D=δeT, e=(1,1,...,1)T∈Rn. Z is a signature matrix of A-1. Recommended initial guess for the sequence (5) is

$$
\underline{\mathbf{X}}^o = \overline{\mathbf{X}}^o = \text{diag } \mathbf{A}^{-1} \mathbf{B}.\tag{6}
$$

It is assumed that coefficient matrix **A** is regular and inverse stable. It is assured if following conditions are satisfied

$$\rho\left(\left|A^{-1}\right|\Lambda\right)<1\tag{7}$$

and

$$
\Delta \left( I - \left| A^{-1} \right| \Delta \right)^{-1} \left| A^{-1} \right| < I \tag{8}
$$

(I is the unit matrix).

Conditions (7) and (8) under which the method works are satisfied if all coefficients of A-1 are nonzero and Δ is sufficiently small. Under this conditions the sequences {*x* ¯ *<sup>k</sup>* } , {*x*¯ *<sup>k</sup>* } converge from starting point (6). Note that each column of *X* ¯ *<sup>o</sup>* , *<sup>X</sup>*¯*<sup>o</sup>* is equal to an approximate solution to Ax=b.

#### **2.4. Numerical examples**

EXAMPLE 2. Consider the simple RC voltage divider circuit of Fig.3, where the resistor and the capacitor values are allowed to vary from 4.4Ω to 4.6Ω and from 520 μF to 580 μF respec‐ tively.

We want to calculate amplitude variations of the frequency response due to the changes of

resistance and capacitance. For a fixed frequency, amplitude of the frequency response *K*(*ω*)= |*U*2( *jω*)/ *U*1( *jω*)| is obtained as a function of variations of R and C by solving fol‐ lowing two linear interval equations

$$
\begin{bmatrix}
1 & -\omega \mathsf{L} \underline{R} \underline{C} \underline{\mathsf{C}} \underline{\mathsf{C}} \begin{bmatrix} \mathsf{L}^{(\mathsf{R})} \\ 1 \end{bmatrix} \begin{bmatrix} \boldsymbol{U}\_{2}^{(\mathsf{R})} \\ \boldsymbol{I} \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}
\end{bmatrix}
$$

where *R* = *R* ¯ , *<sup>R</sup>*¯ , *<sup>C</sup>* <sup>=</sup> *<sup>C</sup>* ¯ , *<sup>C</sup>*¯ and U1=1, *<sup>U</sup>*<sup>2</sup> <sup>=</sup>*U*<sup>2</sup> (*R*) <sup>+</sup> *jU*<sup>2</sup> (*I* ) .

In Fig.4, we have plotted K(ω) as a function of ω. Hatched region show the result by Oettli-Prager condition. Outer bounds were computed with use of the sequences (5). Both results are in a fairly good agreement.

**Figure 4.** Variations of the frequency response of a voltage divider circuit

EXAMPLE 3. Consider the linear RC circuit illustrated in Fig.5 [9]. The interval parameter with the center value m and the width 2r is denoted as (m, r). The numerical data are following: Yk = Gk + jBk, Gk=(1, ε), Bk=(0.2, 0.2ε), J=10

**Figure 5.** The circuit with complex admittances done as interval numbers.

In spite of remarks in previous section on drawbacks of nodal analysis we applied here nodal equations. It allowed us to fulfill assumptions of Rohn's algorithm although the elements of node-admittance matrix are dependent. In order to carry the sequences (5) working with Eq.(1) represented by real interval numbers we replaced the interval complex nodal equations

*Y V* = *J*

by the real interval equation

$$\mathbf{A} \cdot \mathbf{x} = b$$

where

( ) <sup>1</sup>

Conditions (7) and (8) under which the method works are satisfied if all coefficients of A-1 are

EXAMPLE 2. Consider the simple RC voltage divider circuit of Fig.3, where the resistor and the capacitor values are allowed to vary from 4.4Ω to 4.6Ω and from 520 μF to 580 μF respec‐

R

We want to calculate amplitude variations of the frequency response due to the changes of resistance and capacitance. For a fixed frequency, amplitude of the frequency response *K*(*ω*)= |*U*2( *jω*)/ *U*1( *jω*)| is obtained as a function of variations of R and C by solving fol‐ lowing two linear interval equations

> ¯ *C* ¯ , *<sup>R</sup>*¯*C*¯

In Fig.4, we have plotted K(ω) as a function of ω. Hatched region show the result by Oettli-Prager condition. Outer bounds were computed with use of the sequences (5). Both results are

(*R*) <sup>+</sup> *jU*<sup>2</sup>

1 −*ω R*

, *<sup>C</sup>*¯ and U1=1, *<sup>U</sup>*<sup>2</sup> <sup>=</sup>*U*<sup>2</sup>

, *<sup>R</sup>*¯*C*¯ <sup>1</sup>

*ω R* ¯ *C* ¯

U <sup>1</sup>

¯

C

U <sup>2</sup>

*U*2 (*R*)

(*<sup>I</sup>* ) <sup>=</sup> <sup>1</sup> 0

*U*2

(*I* ) .

*A* 1 - D < (7)

( ) <sup>1</sup> 1 1 *IA A I* - - - D- D < (8)

¯

*<sup>o</sup>* , *<sup>X</sup>*¯*<sup>o</sup>* is equal to an approximate solution

*<sup>k</sup>* } , {*x*¯ *<sup>k</sup>* } converge

r

nonzero and Δ is sufficiently small. Under this conditions the sequences {*x*

from starting point (6). Note that each column of *X*

and

74 Analog Circuits

to Ax=b.

tively.

(I is the unit matrix).

**2.4. Numerical examples**

**Figure 3.** RC voltage divider

where *R* = *R*

¯

in a fairly good agreement.

, *<sup>R</sup>*¯ , *<sup>C</sup>* <sup>=</sup> *<sup>C</sup>*

¯

$$A = \begin{bmatrix} G & \cdot B \\ B & G \end{bmatrix} \quad \text{x} = \begin{bmatrix} V\_R \\ V\_I \end{bmatrix} \quad b = \begin{bmatrix} J\_R \\ J\_I \end{bmatrix}$$

$$J\_R = \text{Re}(f), \quad J\_I = \text{Im}(f)$$

$$G = \text{Re}(Y), \quad B = \text{Im}(Y), \quad V\_R = \text{Re}(V), \quad V\_I = \text{Im}(V)$$


The results are tabulated below

**Table 1.** Ranges for nodal voltages

Results are essentially the same as in [10] although the widths of nodal voltages are greater. It confirms the influence of coefficient dependence on evaluation of interval solutions.

EXAMPLE 4. To test properties of presented approach let consider the active two-port shown in Fig.6.

Current-controlled description of the two-port is following

$$
\begin{bmatrix} \mu\_1\\ \mu\_2 \end{bmatrix} = \begin{bmatrix} R\_1 + R\_3 - \operatorname{g} R\_1 R\_3 & R\_3 + \alpha - \operatorname{g} R\_1 R\_3\\ R\_3 - \operatorname{g} + \operatorname{g} R\_2 R\_3 & R\_1 + R\_3 + \operatorname{g} R\_1 R\_3 \end{bmatrix} \begin{bmatrix} i\_1\\ i\_2 \end{bmatrix}
$$

For numerical data

$$\mathbf{R}\_1 = \mathbf{R}\_2 = \mathbf{R}\_3 = 1 \quad , \quad \alpha = \begin{bmatrix} -1, 1 \end{bmatrix} \text{.} \quad \beta = 2 \text{ .} \ g = \begin{bmatrix} 0, 1 \end{bmatrix}$$

**Figure 6.** The active two-port with interval parameters.

interval representation is following

$$
\begin{bmatrix} \mu\_1\\ \mu\_2 \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} 1 \ \ddots \ 2 \end{bmatrix} & \begin{bmatrix} 0, \ 1 \end{bmatrix} \begin{bmatrix} i\_1\\ \vdots \end{bmatrix} \\\ \begin{bmatrix} -1, \ 0 \end{bmatrix} & \begin{bmatrix} 2, \ 3 \end{bmatrix} \begin{bmatrix} i\_2\\ i\_2 \end{bmatrix} \end{bmatrix}
$$

If excitations are determined as U1 ∈ [2,3], U2 ∈ [0,1] the solutions set S has form shown in Fig.7.

**Figure 7.** The solutions set of i1 and i2

The results are tabulated below

76 Analog Circuits

**Table 1.** Ranges for nodal voltages

in Fig.6.

For numerical data

**nodal voltage ε = 0.05**

**nodal voltage ε = 0.1**

Current-controlled description of the two-port is following

*u*1 *u*2 =

u1

**Figure 6.** The active two-port with interval parameters.

interval representation is following

**Re Jm**

**Re Jm**

V1 [5.7635, 6.2724] [-1.2652, -11447] V2 [2.3055, 2.5090] [-0.5061, -0,4579] V3 [1.1527, 1.2545] [-0.2530, -0.2289]

V1 [5.5336, 6.5529] [-1.3355, -10927] V2 [2.2135, 2.6212] [-0.5342, -0.4371] V3 [1.1067, 1.3106] [-0.2671, -0.2185]

Results are essentially the same as in [10] although the widths of nodal voltages are greater. It

EXAMPLE 4. To test properties of presented approach let consider the active two-port shown

*R*<sup>1</sup> + *R*<sup>3</sup> − *gR*1*R*<sup>3</sup> *R*<sup>3</sup> + *α* − *gR*1*R*<sup>3</sup> *R*<sup>3</sup> −*β* + *gR*2*R*<sup>3</sup> *R*<sup>1</sup> + *R*<sup>3</sup> + *gR*1*R*<sup>3</sup>

R1= R2= R3= 1 , α= -1,1 , β= 2 , g = 0,1

R <sup>1</sup> R <sup>2</sup>

i <sup>1</sup> i <sup>2</sup> <sup>i</sup> <sup>2</sup> <sup>i</sup> <sup>1</sup>

g u0

<sup>u</sup> <sup>R</sup> <sup>3</sup> u2 <sup>0</sup>

*i* 1 *i* 2

confirms the influence of coefficient dependence on evaluation of interval solutions.

Computing bounds on operating region for i1 and i2 with recommendation given by (6) and slightly different starting intervals following results were obtained:

$$i\_1 = \begin{bmatrix} i\_{1\prime} & \dot{i}\_1 \mathbf{J} & i\_1 \in \mathbb{I} - 1.5 \ \text{ 0.3} \mathbf{J} & \dot{i}\_1 \in \mathbb{I} \text{3.9} \ \text{5.2} \end{bmatrix}$$

$$i\_2 = \begin{bmatrix} i\_{2\prime} & \dot{i}\_2 \mathbf{J} & i\_2 \in \mathbb{I} - 2.0 \ \text{ -1.2} \mathbf{J} & \dot{i}\_2 \in \mathbb{I} \text{2.5} \ \text{4.0} \end{bmatrix}$$

It's seen that this approach give poor results if widths of matrix coefficients are large.

#### **3. On finding operating points of nonlinear, inertialess circuits**

A continuation method is a method of dc circuit analysis that under proper conditions, is guaranteed to find a circuit's operating point [26], [27]. It can be also used in calculating a circuit's periodic steady-state response [24], [25]. Finding the dc operating point is essential for circuit simulation: both steady-state and transient analyses require a priori knowledge of circuit's dc operating point. The idea behind a continuation method is to embed in the circuit's equations

$$\mathbf{F}(\mathbf{x}) = \mathbf{0} \tag{9}$$

( nodal or hybrid ) an additional parameter λ such that the circuit that corresponds to λ = 0 is trivial or easy to find: the circuit that corresponds to λ = 1 is the original circuit whose solution is desired. Hence, a mapping (homotopy) H: Rn+1 → Rn is constructed, where H(x, 0) = 0 gives the solution to the initial ( start ) circuit and H(x, 1) = 0 gives the solution to the circuit being simulated. To find this solution, one generate a sequence of points {(*x <sup>k</sup>* , *λ <sup>k</sup>* )}*<sup>k</sup>* =0 *<sup>N</sup>* , where xk is on or near the path and x0 is a known solution of H(x,0) = 0. The path in Rn+1 begins at λ = 0 and ends at λ = 1. An important problem in following homotopy path is a control of the step from a point xk to xk+1. The design philosophies of a step size control are numerous and varied. See, e.g. [28] - [31].

Most of the methods on the following continuation paths are based on approximate step control. Approximate step control methods are successful and fast when the path is smooth and isolated, but problems arise when there are many paths near some points. In that case, algorithms based on approximate step control methods may jump from one path to another. Also, if rapid changes in curvature occur along the path, the method based on approximate step control sometimes even erroneously reverses orientation. However, appropriate use of interval analysis gives us guarantee that the predictor algorithm will not jump from one path to another, or, indeed, jump over different legs of the same path. We present an interval step control for tracing continuation paths which assures that the predictor-corrector iterations will not jump across paths, and each predictor step is as large as possible, subject to verification that the path is unique with the given interval extension.

#### **3.1. Step size control**

To trace a path from a known solution (*x* <sup>0</sup> , *λ* <sup>0</sup> ), we first predict the solution for *λ* =*λ* <sup>0</sup> + *Δλ* and then correct the prediction using Newton method with λ fixed. Specifically, for small Δx and Δλ, the Taylor expansion for H gives

$$H(\mathbf{x}^{0} + \Delta \mathbf{x}, \boldsymbol{\lambda}^{0} + \Delta \boldsymbol{\lambda}) \approx H(\mathbf{x}^{0}, \boldsymbol{\lambda}^{0}) + \boldsymbol{I}\_{\boldsymbol{x}}(\mathbf{x}^{0}, \boldsymbol{\lambda}^{0}) \Delta \mathbf{x} + \boldsymbol{I}\_{\boldsymbol{\lambda}}(\mathbf{x}^{0}, \boldsymbol{\lambda}^{0}) \Delta \mathbf{x} \tag{10}$$

where *Jx*(*<sup>x</sup>* <sup>0</sup> , *λ* <sup>0</sup> ) and *Jλ*(*<sup>x</sup>* <sup>0</sup> , *λ* <sup>0</sup> ) are the Jacobian matrices of H with respect to x and λ. For the prediction step, we have *H* (*x* <sup>0</sup> , *λ* <sup>0</sup> )=0, so setting *H* (*x* <sup>0</sup> + *Δx*, *λ* <sup>0</sup> + *Δλ*)=0 gives

$$
\Delta \mathbf{x} = -\sf{J}\_{\mathbf{x}} (\mathbf{x}^{0}, \mathscr{X}^{0})^{-1} \sf{J}\_{\mathbf{x}} (\mathbf{x}^{0}, \mathscr{X}^{0}) \Delta \mathbf{y}. \tag{11}
$$

Since this is only approximate, we correct the solution at the new value of *λ* =*λ* <sup>0</sup> + *Δλ* by Newton-Raphson algorithm starting with initial guess *x*<sup>0</sup> <sup>=</sup> *<sup>x</sup>* <sup>0</sup> <sup>+</sup> *<sup>Δ</sup>x*. A common numerical practice is to stop the Newton iteration whenever the distance between two iterates is less than a given tolerance, i.e., when

$$\left\|\mathbf{x}\_{k+1} - \mathbf{x}\_k\right\| < \varepsilon \tag{12}$$

However, just the fact that (12) is satisfied does not guarantee the existence of a solution. We try to overcome this difficulty by first performing the three Newton steps and using them to compute an interval vector that is very likely to contain a solution of the correction step. We adopt here results presented in [22]. Consider Newton's method

$$\mathbf{x}\_{k+1} = \mathbf{x}\_k - \mathbf{J}\_x(\mathbf{x}\_{k'}\mathcal{A})^{-1} H(\mathbf{x}\_{k'}\mathcal{A}), \quad \mathbf{x}\_0 = \mathbf{x}^0 + \Delta \mathbf{x}, \quad k = 0, 1, \tag{13}$$

Let

is desired. Hence, a mapping (homotopy) H: Rn+1 → Rn is constructed, where H(x, 0) = 0 gives the solution to the initial ( start ) circuit and H(x, 1) = 0 gives the solution to the circuit being

or near the path and x0 is a known solution of H(x,0) = 0. The path in Rn+1 begins at λ = 0 and ends at λ = 1. An important problem in following homotopy path is a control of the step from a point xk to xk+1. The design philosophies of a step size control are numerous and varied. See,

Most of the methods on the following continuation paths are based on approximate step control. Approximate step control methods are successful and fast when the path is smooth and isolated, but problems arise when there are many paths near some points. In that case, algorithms based on approximate step control methods may jump from one path to another. Also, if rapid changes in curvature occur along the path, the method based on approximate step control sometimes even erroneously reverses orientation. However, appropriate use of interval analysis gives us guarantee that the predictor algorithm will not jump from one path to another, or, indeed, jump over different legs of the same path. We present an interval step control for tracing continuation paths which assures that the predictor-corrector iterations will not jump across paths, and each predictor step is as large as possible, subject to verification

, *λ* <sup>0</sup>

+D +D » + D + D

 l

 ll

0 0 00 00 00 ( , ) (,) (,) (,) *Hx x Hx J x x J x <sup>x</sup>*

, *λ* <sup>0</sup>

l

*k k* <sup>1</sup> *x x*

0 01 0 0 (,) (,) . *<sup>x</sup> x Jx Jx* l

Since this is only approximate, we correct the solution at the new value of *λ* =*λ* <sup>0</sup> + *Δλ* by Newton-Raphson algorithm starting with initial guess *x*<sup>0</sup> <sup>=</sup> *<sup>x</sup>* <sup>0</sup> <sup>+</sup> *<sup>Δ</sup>x*. A common numerical practice is to stop the Newton iteration whenever the distance between two iterates is less than a given tolerance, i.e., when

e

 l

and then correct the prediction using Newton method with λ fixed. Specifically, for small Δx and Δλ, the Taylor expansion for H gives

, *λ <sup>k</sup>* )}*<sup>k</sup>* =0

), we first predict the solution for *λ* =*λ* <sup>0</sup> + *Δλ*

(10)

l

) are the Jacobian matrices of H with respect to x and λ. For


 ll

)=0, so setting *H* (*x* <sup>0</sup> + *Δx*, *λ* <sup>0</sup> + *Δλ*)=0 gives

<sup>+</sup> - < (12)

*<sup>N</sup>* , where xk is on

simulated. To find this solution, one generate a sequence of points {(*x <sup>k</sup>*

that the path is unique with the given interval extension.

To trace a path from a known solution (*x* <sup>0</sup>

l

) and *Jλ*(*<sup>x</sup>* <sup>0</sup>

the prediction step, we have *H* (*x* <sup>0</sup>

 l

, *λ* <sup>0</sup>

e.g. [28] - [31].

78 Analog Circuits

**3.1. Step size control**

where *Jx*(*<sup>x</sup>* <sup>0</sup>

, *λ* <sup>0</sup>

$$\rho\_k = \left\| \boldsymbol{x}\_{k+1} - \boldsymbol{x}\_k \right\|\_{\boldsymbol{x}\_{\mathcal{I}}} \tag{14}$$

for some fixed k and

$$\mathbf{x} = \{ \mathbf{x} \in \mathbb{R}^n \: \left\| \mathbf{x}\_1 - \mathbf{x} \right\|\_{\infty} < \rho\_0 \}\tag{15}$$

Combining the properties of the Krawczyk operator and a corollary of the Newton-Kantoro‐ vich theorem it can be shown, that for three successive iterates the following relation holds

$$\begin{cases} 8\rho\_1^3 \\ \rho\_2 \Big|\_{\omega} \cdot \rho\_0^2 \Big(\|\mathbf{x}\_2\|\_{\omega} \cdot \rho\_0^2\Big)^{<\epsilon pps.} \end{cases} \tag{16}$$

In Fig. 8 there is a sketch of the idea of constructing the predictor - corrector step for twodimensional problems, i.e. x = (y, z).

We can compute interval (15) and using the Krawczyk operator K(x) we test whenever this interval contains a solution. It is true if the Krawczyk operator satisfy an inclusion

$$K(\mathbf{x}) = \mathbf{x}\_2 - \mathbf{C} \cdot H(\mathbf{x}\_2, \lambda) + (I - \mathbf{C} \cdot I\_x(\mathbf{x}, \lambda))(\mathbf{x} - \mathbf{x}\_2) \subseteq \mathbf{x} \tag{17}$$

then x contain solution of the correction step *<sup>x</sup>* <sup>∗</sup>. Here *Jx*(*x*, *<sup>λ</sup>*) denotes interval arithmetic evaluation of *Jx*(*x*, *λ*), *C* = *Jx*(*x*1, *λ*) <sup>−</sup>1, and I is identity matrix [20]. If additional condition

$$\left\| I - C \cdot f\_x(\mathbf{x}, \boldsymbol{\lambda}) \right\| < 1 \tag{18}$$

is satisfied [21], then the solution *x* <sup>∗</sup> is unique in x and the operator (17) generates a sequence of intervals *xr*, *r* =1, 2, …, *m*, that satisfy the relations

**Figure 8.** The predictor-corrector step for two dimensional problem

$$\mathbf{x}^\* \in \mathbf{x}\_{\mathbf{m}} \subseteq \mathbf{x}\_{\mathbf{m-l}} \dots \tag{19}$$

This nested sequence is stopped when the width of the interval xm becomes smaller than

μ > 0 and then *<sup>x</sup>* <sup>∗</sup> <sup>=</sup>*m*(*xm*). The correction step may be repeated several times before taking the next prediction step.

Step size control can be formalized in the following algorithm.

#### **Algorithm**

Unless otherwise indicated STEP {k+1} follows STEP {k}.

Given: Solution (*x* <sup>0</sup> , λ <sup>0</sup> ) of a homotopy*H*(*x*, λ) = 0.

STEP {1} Assume 0 < Δλ < 1;

STEP {2} Compute prediction step Δx according to (11);

STEP {3} Compute using Newton method (13) three successive iterates;

STEP {4} Verify the inequality (16) for some eps "/> 0;

STEP {5} Check the relation (17). If (17) is not satisfied take Δλ = Δλ/2 and go to {2};

STEP {6} Compute interval **x** according to (14) and (15);

STEP {7} Check the condition (18);

STEP {8} Generate a sequence (19) according to

```
x1=K(x), x2=K(x1)...xm=K(xm-1
```
until *w*(**x**m) <μ, μ > 0,

where *w*(**x**m) = *x*¯ *<sup>m</sup>* − *x* ¯ *<sup>m</sup>* ( component wise );

STEP {9} Compute

*x*\* =*m*(**x**m)=(*x* ¯ *<sup>m</sup>* <sup>+</sup> *<sup>x</sup>*¯ *<sup>m</sup>*) / 2 ( component wise );

> STEP {10} Assign *x* <sup>0</sup>← *x* ' = *x* <sup>∗</sup> and λ <sup>0</sup>←λ '

and go to {1} ( new prediction );

#### **3.2. Box — Bisection searching**

Introducing some interval operator ( e.g. Krawczyk operator, Hansen-Sengupta operator or some of their modifications [20]) for system of nonlinear equations we can formulate gener‐ alized bisection algorithm applicable to find solutions of equation (9).

Let denote Krawczyk operator associated with Eq. (9) as

$$K(\mathbf{x}, \mathbf{y}, \mathbf{F}) = y - \mathcal{Y}F(y) + R(\mathbf{x})(\mathbf{x} - \mathbf{y}) \tag{20}$$

Here *Y=*[*m(F'*(**x**))]-1 where *F'*(**x**) is an interval arithmetic evaluation of *F* '(*x*) (Jacobian matrix of Eq. (1)), y ∈ **x** ( e.g. y = m(**x**) ),

$$R(\mathbf{x}) = I - \mathcal{Y}F'(\mathbf{x})\tag{21}$$

and I is identity matrix.

Let T denotes the list of subregions yet to be tested. P is the list of subregions of B which may contain a solution to (9) but which are too small for further analysis, i.e. w(**x**) ≤ ε (ε - accuracy of searching).

Unless otherwise indicated STEP {k+1} follows STEP {k}.

Given: n-dimensionsl box (interval vector) B ⊆ D ⊆ Rn.

#### **Algorithm**

m m-1 *<sup>x</sup>*\* Î Í **x x** <sup>K</sup> (19)

+ Dx

r0

y

r0

r0

x0 = x0

This nested sequence is stopped when the width of the interval xm becomes smaller than

r0

r<sup>1</sup> r<sup>1</sup>

x2

x1

r1

(x',l')

r1

(x0 ,l0 )

**Figure 8.** The predictor-corrector step for two dimensional problem

z

80 Analog Circuits

, λ <sup>0</sup>

) of a homotopy*H*(*x*, λ) = 0.

Step size control can be formalized in the following algorithm.

Given: Solution (*x* <sup>0</sup>

Unless otherwise indicated STEP {k+1} follows STEP {k}.

STEP {2} Compute prediction step Δx according to (11);

STEP {4} Verify the inequality (16) for some eps "/> 0;

STEP {6} Compute interval **x** according to (14) and (15);

STEP {3} Compute using Newton method (13) three successive iterates;

STEP {5} Check the relation (17). If (17) is not satisfied take Δλ = Δλ/2 and go to {2};

**Algorithm**

STEP {1} Assume 0 < Δλ < 1;

μ > 0 and then *<sup>x</sup>* <sup>∗</sup> <sup>=</sup>*m*(*xm*). The correction step may be repeated several times before taking the next prediction step.

> STEP {1} Assign **x**← *B*; T ← ∅; P ← ∅; STEP {2} Compute F(**x**); STEP {3} If 0 ∉ F(**x**) go to STEP {11}; STEP {4} Compute Y, if Y is singular go to STEP {9};

STEP {5} Assuming y = m(**x**) compute K(**x**, y, F**)** and R(**x**);

STEP {6} If K(**x**, y, F**)** ∩ **x** = ∅ go to STEP {11};

STEP {7} If ||K(**x**, *y, F*)-**x**||≤*w*(**x**)/2|| , then **x** contains a solution, continue, else go to STEP {9};

STEP {8} ||R(**x**)<1|| If , then there is a unique solution to (1) in **x**, if w(**x**) ≤ ε terminate search, go to STEP {11};

STEP {9} If w(**x**) ≤ ε add **x** to the list P and go to STEP {11} else bisect **x** on **x**' and' **x**'' according to some rules of bisection;

STEP {10} Assign **x** ← **x**', add **x**'' to the head of the list T, go to STEP {2};

STEP {11} If list T is empty go to STEP {12}, otherwise take an interval vector **x** p at the head of list T, assign **x** ← **x** p, delete this vector from list T and go to STEP {2};

STEP {12} If list P is empty, terminate search with no solution in B, otherwise print list P and terminate search.

One can simply check that in the STEPS {5}, {7}, {8} for Algorithm with step size control and in the STEPS {2}, {4},{5} for box – bisection searching one needs to evaluate the ranges of nonlinear functions over some interval **x**. This problem can be solved by computing coefficients of Bernstein polynomials. For details see [34].

#### **3.3. Numerical examples**

As an illustration of the implementation of the approach outlined above, consider the following two examples.

EXAMPLE 5. Our first example concerns a simple circuit of Fig. 9 which consists of two tunnel diodes and two constant current sources and a linear resistor. The tunnel diode characteristics are defined following

$$h\_1(\mu\_1) = h\_2(\mu\_2) = h\_1(\mu) = 5.0\mu - 1.7\mu^2 + 0.15\mu^3$$

The network equations are obtained from the hybrid analysis as

$$F(\mathbf{x}) = F\left(\mu\_1, \mu\_2\right) = \begin{bmatrix} h\left(\mu\_1\right) \\ h\left(\mu\_2\right) \end{bmatrix} + \frac{1}{R} \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} \mu\_1 \\ \mu\_2 \end{bmatrix} - \begin{bmatrix} j\_1 \\ j\_2 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$

We introduce homotopy of the form

$$H(\mathbf{x}, \lambda) = \gamma (1 - \lambda) Q(\mathbf{x}) + \lambda F(\mathbf{x})$$

where λ ∈ [0,1] and γ is a random complex number, different from zero.

Start equations are following

$$Q(\mathbf{x}) = Q(\mu\_1, \mu\_2) = \begin{cases} \mu\_1^3 - \mathbf{1}, \\ \mu\_2^3 - \mathbf{1}, \end{cases}$$

**Figure 9.** The circuit of two tunnel diodes

STEP {5} Assuming y = m(**x**) compute K(**x**, y, F**)** and R(**x**);

STEP {10} Assign **x** ← **x**', add **x**'' to the head of the list T, go to STEP {2};

STEP {7} If ||K(**x**, *y, F*)-**x**||≤*w*(**x**)/2|| , then **x** contains a solution, continue, else go to STEP {9};

STEP {8} ||R(**x**)<1|| If , then there is a unique solution to (1) in **x**, if w(**x**) ≤ ε terminate search, go to STEP {11}; STEP {9} If w(**x**) ≤ ε add **x** to the list P and go to STEP {11} else bisect **x** on **x**' and' **x**'' according to some rules of

STEP {11} If list T is empty go to STEP {12}, otherwise take an interval vector **x** p at the head of list T, assign **x** ← **x** p,

One can simply check that in the STEPS {5}, {7}, {8} for Algorithm with step size control and in the STEPS {2}, {4},{5} for box – bisection searching one needs to evaluate the ranges of nonlinear functions over some interval **x**. This problem can be solved by computing coefficients of

As an illustration of the implementation of the approach outlined above, consider the following

EXAMPLE 5. Our first example concerns a simple circuit of Fig. 9 which consists of two tunnel diodes and two constant current sources and a linear resistor. The tunnel diode characteristics

*<sup>h</sup>*1(*u*1)=*h*2(*u*2)=*<sup>h</sup>* (*u*)=5.0*<sup>u</sup>* <sup>−</sup>1.7*<sup>u</sup>* <sup>2</sup> + 0.15*<sup>u</sup>* <sup>3</sup>

*H* (*x*, *λ*)=*γ*(1−*λ*)*Q*(*x*) + *λF* (*x*)

*Q*(*x*)=*Q*(*u*1, *u*2)={

> *u*1 <sup>3</sup> −1

> *u*2 <sup>3</sup> −1

*u*1 *u*2 − *j* 1 *j* 2 = 0 0

*h* (*u*1) *h* (*u*2) + 1 *R*

The network equations are obtained from the hybrid analysis as

*F* (*x*)= *F* (*u*1, *u*2)=

where λ ∈ [0,1] and γ is a random complex number, different from zero.

STEP {12} If list P is empty, terminate search with no solution in B, otherwise print list P and terminate search.

STEP {6} If K(**x**, y, F**)** ∩ **x** = ∅ go to STEP {11};

delete this vector from list T and go to STEP {2};

Bernstein polynomials. For details see [34].

**3.3. Numerical examples**

two examples.

are defined following

We introduce homotopy of the form

Start equations are following

bisection;

82 Analog Circuits

It's seen that the solutions of equation

$$H(\mathbf{x}, \ 0) = Q(\mathbf{x}) = 0$$

are distributed on unit circles and they form nine starting points (x0 , 0). For solving so prepared homotopy equations we splited them into real and imaginary parts.

Presented approach leads to the five real solutions

*Q*<sup>1</sup> =(1.805435, 7.325956), *Q*<sup>2</sup> =(3.609712, 6.937570), *Q*<sup>3</sup> =(6.103275, 6.103273), *Q*<sup>4</sup> =(6.937570, 3.609712), *Q*<sup>5</sup> =(7.325956, 1.805434).

Because above equation is equivalent polynomial of 9th degree there are also four complex solutions:

$$Q\_6 = \begin{pmatrix} 1.493997 + j2.483038, & 1.493997 + j9.932153 \end{pmatrix},$$

$$Q\_7 = \begin{pmatrix} 1.493997 - j2.483038, & 1.493997 - j9.932153 \end{pmatrix},$$

$$Q\_8 = \begin{pmatrix} 2.615029 + j2.812082, & 2.615029 + j1.687249 \end{pmatrix},$$

$$Q\_9 = \begin{pmatrix} 2.615029 - j2.812082, & 2.615029 - j1.687249 \end{pmatrix}.$$

EXAMPLE 6. Let consider circuit hybrid equation [17].

$$F(\mathbf{x}) = \begin{bmatrix} f\_1(\mathbf{x}\_1) \\ f\_2(\mathbf{x}\_1 \cdot \mathbf{x}\_2) \\ f\_3(\mathbf{x}\_2 \cdot \mathbf{x}\_3) \end{bmatrix} - \begin{bmatrix} 1 & 1 & 1 \\ 0 & -1 & 0 \\ 0 & 1 & 1 \end{bmatrix} \begin{bmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \mathbf{x}\_3 \end{bmatrix} - \begin{bmatrix} -5 \\ 5 \\ -5 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}$$

Nonlinear functions ( characteristics of nonlinear elements ) are following

$$\begin{aligned} \mathbf{f}\_1(\mathbf{x}\_1) &= \mathbf{x}\_1^3 \\ \mathbf{f}\_2(\mathbf{x}\_1, \mathbf{x}\_2) &= 0.25 \mathbf{x}\_1 \mathbf{x}\_2 \\ \mathbf{f}\_3(\mathbf{x}\_2, \mathbf{x}) &= 0.048(\mathbf{x} + \mathbf{x}) . \end{aligned}$$

For solving this polynomial system we have applied the same homotopy as in Example 5.

Using concept of m-homogenity [18],[19] we can introduce solutions that corresponds to λ=0 and the start system H(x, 0) = Q(x) = 0 look like this

$$\mathbf{Q}(\mathbf{x}) = \begin{cases} -0.509\mathbf{x}\_1^3 + 1.235\mathbf{x}\_1^2 \cdot 2.961\mathbf{x}\_1\mathbf{x}\_2 + 0.902\mathbf{x}\_1^2 \mathbf{x}\_2 \\ + 0.078\mathbf{x}\_1\mathbf{x}\_3 + 0.196\mathbf{x}\_1 + 2.294\mathbf{x}\_2 \cdot 0.235\mathbf{x}\_3 + 1 = 0 \\ -0.228\mathbf{x}\_1\mathbf{x}\_2 - 0.024\mathbf{x}\_1\mathbf{x}\_3 + 0.492\mathbf{x}\_1 + 0.443\mathbf{x}\_2 \\ + 0.069\mathbf{x}\_3 + 1 = 0 \\ -0.641\mathbf{x}\_2^2 \cdot 0.143\mathbf{x}\_3^2 + 0.642\mathbf{x}\_2\mathbf{x}\_3 + 1.603\mathbf{x}\_2 \\ -0.781\mathbf{x}\_3 + 1 = 0 \end{cases}$$

The start system respects the m-homogeneous structure [18] of the hybrid equation we want to solve and has eight start solutions. Namely

$$\begin{aligned} \mathbf{x}\_1^0 &= \text{(3.0, 2.0, 1.0)}, & \mathbf{x}\_2^0 &= \text{(1.5, 2.0, 1.0)},\\ \mathbf{x}\_3^0 &= \text{(2.0, 3.02.5, )}, & \mathbf{x}\_4^0 &= \text{(2.5, 1.5, 1.0)},\\ \mathbf{x}\_5^0 &= \text{(1.0, 2.0, 2.5)}, & \mathbf{x}\_6^0 &= \text{(2.0, 1.0, 3.0)},\\ \mathbf{x}\_7^0 &= \text{(1.5, 1.0, 2.0)}, & \mathbf{x}\_8^0 &= \text{(3.0, 2.0, 3.0)}.\end{aligned}$$

We applied the step size control previously described with tolerance eps = 0.01. We needed to trace only eight solution paths. Two pairs of them reached the same solutions, and so we obtained following six solutions of the hybrid equations:

$$\begin{aligned} \mathbf{x}\_A &= \text{(1.7155, 3.4992, 4.8340)}, \\ \mathbf{x}\_B &= \text{(2.1273, 3.2641, 9.2359)}, \\ \mathbf{x}\_C &= \text{(-0.8578+j1.0988, 5.6714-j1.9832, 2.6619+j1.9832)}, \\ \mathbf{x}\_D &= \text{(-0.8578-j1.0988, 5.6714+j1.9832, 2.6619-j1.9832)}, \\ \mathbf{x}\_E &= \text{(-1.9832+j1.5473, 5.3309-j2.8091, 7.1691+j2.2091)}, \\ \mathbf{x}\_F &= \text{(-1.9832-j1.5473, 5.3309+j2.8091, 7.1691-j2.2091)}. \end{aligned}$$

For circuit analysis only real solutions xA and xB are of interest.

EXAMPLE 7. In the last example an electric circuit depicted in Fig. 10 assembling some linear resistors, diodes and operational amplifier is considered (cf. [29]).

**Figure 10.** The circuit with operational amplifier and diodes

The voltage induced in the operational amplifier is modeled by:

$$\mu\_A(x) = 7.65 \arctan(1962x),$$

the current through the diodes satisfies:

f1(x1)=x1 3

and the start system H(x, 0) = Q(x) = 0 look like this

84 Analog Circuits


+0.069x3+1=0


<sup>0</sup> =(3.0, 2.0, 1.0), *<sup>x</sup>*<sup>2</sup>

<sup>0</sup> =(2.0, 3.02.5, ), *<sup>x</sup>*<sup>4</sup>

<sup>0</sup> =(1.0, 2.0, 2.5), *<sup>x</sup>*<sup>6</sup>

<sup>0</sup> =(1.5, 1.0, 2.0), *<sup>x</sup>*<sup>8</sup>


Q(x)={

to solve and has eight start solutions. Namely

*x*1

*x*3

*x*5

*x*7

obtained following six solutions of the hybrid equations:

*xA* =(1.7155, 3.4992, 4.8340), *xB* =(2.1273, 3.2641, 9.2359),

For circuit analysis only real solutions xA and xB are of interest.

resistors, diodes and operational amplifier is considered (cf. [29]).

f2(x1,x2)=0.25x1x2 f3(x2,x)=0.048(x+x).

For solving this polynomial system we have applied the same homotopy as in Example 5.

Using concept of m-homogenity [18],[19] we can introduce solutions that corresponds to λ=0


+0.642x2x3+1.603x2

<sup>0</sup> =(1.5, 2.0, 1.0),

<sup>0</sup> =(2.5, 1.5, 1.0),

<sup>0</sup> =(2.0, 1.0, 3.0),

<sup>0</sup> =(3.0, 2.0, 3.0).

+0.078x1x3+0.196x1+2.294x2-0.235x3+1=0 -0.228x1x2-0.024x1x3+0.492x1+0.443x2

The start system respects the m-homogeneous structure [18] of the hybrid equation we want

We applied the step size control previously described with tolerance eps = 0.01. We needed to trace only eight solution paths. Two pairs of them reached the same solutions, and so we

> *xC* =(−0.8578 + *j*1.0988, 5.6714− *j*1.9832, 2.6619 + *j*1.9832), *xD* =(−0.8578− *j*1.0988, 5.6714 + *j*1.9832, 2.6619− *j*1.9832), *xE* =(−1.9832 + *j*1.5473, 5.3309− *j*2.8091, 7.1691 + *j*2.2091), *xF* =(−1.9832− *j*1.5473, 5.3309 + *j*2.8091, 7.1691− *j*2.2091).

EXAMPLE 7. In the last example an electric circuit depicted in Fig. 10 assembling some linear

2 x2

$$i\_D(\mu) = 5.6 \cdot 10^{-8} (e^{25\mu} - 1);$$

resistances in ohms are following: R1 = 51, R2 = 39, R3 = 10, R4 = 104 , R5 = 0.201, R6 = 25.5, R7 = 0.62, R8 = 1, R9 = 13. Applying Kirchhoff's laws to the circuit leads to the system of nonlinear equations:

$$\begin{aligned} f\_1 &= \left(\upsilon\_1 - \upsilon\_3\right) / R\_4 + \left(\upsilon\_1 - \upsilon\_2\right) / R\_2 + \left(\upsilon\_1 + E\right) / R\_1 = 0 \\ f\_2 &= \left(\upsilon\_2 - \upsilon\_6\right) / R\_3 + i\_D(\upsilon\_2) + \left(\upsilon\_2 - \upsilon\_1\right) / R\_2 = 0 \\ f\_3 &= \left(\upsilon\_3 - \upsilon\_1\right) / R\_4 + \left(\upsilon\_3 - \upsilon\_4\right) / R\_6 = 0 \\ f\_4 &= \left(\upsilon\_4 - \upsilon\_3\right) / R\_6 + \upsilon\_4 / R\_7 + \left(\upsilon\_4 - \upsilon\_5\right) / R\_8 = 0 \\ f\_5 &= \left(\upsilon\_5 - \upsilon\_4\right) / R\_8 + \left(\upsilon\_5 - \upsilon\_6\right) / R\_9 + i\_D(\upsilon\_5) = 0 \\ f\_6 &= \left(\upsilon\_6 - \upsilon\_2\right) / R\_3 - \left[\mu\_A(\upsilon\_3 - \upsilon\_1) - \upsilon\_6\right] / R\_5 \\ &\qquad + \left(\upsilon\_6 - \upsilon\_5\right) / R\_9 = 0 \end{aligned}$$

Functions describing nonlinearities of diodes and operational amplifier are smooth functions so we have needed their expansion to compute range values. We have assumed Taylor expansion with degree = 4 and for Bernstein coefficients degree = 10 [34]. Applying box– bisection algorithm we have obtained solutions for two different ranges of the input voltage E. Results are tabulated as


**Table 2.** Solutions presented here are in good agreement with the results of [29].

## **4. Concluding remarks**

In first section, we applied an algorithm for bounding the solution set of a system of linear interval equation which works when the interval elements of **A** and the interval components of **b** are relatively narrow. In practice, however, we are often in situation like in Example 4 when the widths of these intervals are rather wide. This case creates difficulties. Presented here version of Rohn's algorithm doesn't work well. Gaussian elimination using interval arithmetic tends to give poor results. Interval Gaussian elimination can fail because of division by an interval containing zero. This can occur even when the interval coefficient matrix is regular i.e. does not contain a singular real matrix [2]. Maybe methods based on precondi‐ tioning of the linear interval equations [7] offer new possibilities. However, recently it was reported that the problem of computing the exact bounds of solution set for a system of linear interval equations is NP-hard (computationally intractable) [13]. Shortly speaking, NPhardness of a problem P means that if we are able to solve this problem in reasonable time, then we would be able to solve all problems from a very large class of complicated problems (called class NP) in reasonable time, and this possibility is widely believed to be impossible. Reasonable time means a time that does not exceed some polynomial of the length of the input. For exact definitions see, e.g., [15]. On the other hand it was also proved [14] that if the interval components of **A** and **b** are "thin" enough, then there exists a polynomial-time algorithm that computes the exact bounds for S in "almost all" cases ("almost all" in some reasonable sense). It means that solving linear interval equations isn't a hopeless task.

Second sections presents an approach to compute approximations to the solution path of a systems of nonlinear equations. The suggested method is based on predictor-corrector principle. We have combined common Newton-Raphson correction step with the properties of the interval Krawczyk operator. An interval step control for tracing continuation path assures that the predictor-corrector iterations were not jumping across paths. A robust and efficient path tracker was obtained by halving the homotopy parameter step whenever uniqueness of correction solution is not assured within a given tolerance of the continuation path and doubling it when the prediction step has been sufficiently accurate several consecu‐ tive times. Presented here path tracker was more effective than the tracing homotopy path via integration so-called basic differential equations [23]. Krawczyk operator was also useful in generalized box-bisection for searching all solutions of the circuit equations.

## **Author details**

E [ 0.320, 0.322 ] [ 0.599, 0.601 ] v1 [ 0.2325, 0.2354 ] [ 0.0439, 0.0493 ] v2 [ 0.6578, 0.6629 ] [ 0.5423, 0.5471 ] v3 [ 0.2303, 0.2375 ] [ 0.0464, 0.0492 ] v4 [ 0.2299, 0.2375 ] [ 0.0446, 0.0496 ] v5 [ 0.6155, 0.6208 ] [ 0.1232, 0.1293 ] v6 [ 9.5941, 9.6088 ] [ 1.1597, 1.1662 ]

In first section, we applied an algorithm for bounding the solution set of a system of linear interval equation which works when the interval elements of **A** and the interval components of **b** are relatively narrow. In practice, however, we are often in situation like in Example 4 when the widths of these intervals are rather wide. This case creates difficulties. Presented here version of Rohn's algorithm doesn't work well. Gaussian elimination using interval arithmetic tends to give poor results. Interval Gaussian elimination can fail because of division by an interval containing zero. This can occur even when the interval coefficient matrix is regular i.e. does not contain a singular real matrix [2]. Maybe methods based on precondi‐ tioning of the linear interval equations [7] offer new possibilities. However, recently it was reported that the problem of computing the exact bounds of solution set for a system of linear interval equations is NP-hard (computationally intractable) [13]. Shortly speaking, NPhardness of a problem P means that if we are able to solve this problem in reasonable time, then we would be able to solve all problems from a very large class of complicated problems (called class NP) in reasonable time, and this possibility is widely believed to be impossible. Reasonable time means a time that does not exceed some polynomial of the length of the input. For exact definitions see, e.g., [15]. On the other hand it was also proved [14] that if the interval components of **A** and **b** are "thin" enough, then there exists a polynomial-time algorithm that computes the exact bounds for S in "almost all" cases ("almost all" in some reasonable sense).

Second sections presents an approach to compute approximations to the solution path of a systems of nonlinear equations. The suggested method is based on predictor-corrector principle. We have combined common Newton-Raphson correction step with the properties of the interval Krawczyk operator. An interval step control for tracing continuation path assures that the predictor-corrector iterations were not jumping across paths. A robust and efficient path tracker was obtained by halving the homotopy parameter step whenever uniqueness of correction solution is not assured within a given tolerance of the continuation path and doubling it when the prediction step has been sufficiently accurate several consecu‐ tive times. Presented here path tracker was more effective than the tracing homotopy path via

**Table 2.** Solutions presented here are in good agreement with the results of [29].

It means that solving linear interval equations isn't a hopeless task.

**4. Concluding remarks**

86 Analog Circuits

Zygmunt Garczarczyk\*

Silesian University of Technology, Gliwice, Poland

## **References**


[30] Forster, W. Some computational methods for systems of nonlinear equations and systems of polynomial equations. J. Global Optimization (1992). , 2(4), 317-356.

[14] Rohn, J, & Kreinovich, V. Computing exact componentwise bounds on solution of linear systems with interval data is NP-hard. SIAM J. Matrix Anal. Appl. (1992). , 16(2),

[15] Lakeyev, A. V, & Kreinovich, V. If input intervals are small then interval computations are almost always easy. Reliable Computing 1995, Supplement (1995). , 134-139. [16] Garey, M. E, & Johnson, D. S. Computers and Intractability: A Guide to the Theory of

[17] Chua, L. O, & Lin, P. M. Computer-Aided Analysis of Electronic Circuits. Englewood

[18] Morgan, A, & Sommese, A. A homotopy for solving general polynomial systems that respects m-homogeneous structures. Appl.Math.Comput. (1987). , 24(2), 101-113. [19] Garczarczyk, Z. Circuit design problems via polynomial equations solving, In Procee‐ digs of the 12th European Conference on Circuit Theory and Design, ECCTD'August

[20] Moore, R. E, Kearfott, R. B, & Cloud, M. J. Introduction to Interval Analysis. Philadel‐

[21] Moore, R. E. A test for existence of solutions to nonlinear systems. SIAM J. Numer.

[22] Alefeld G, Gienger A, & Potra F. Efficient numerical validation of solutions of nonlinear

[23] Garcia, C. B, & Zangwill, W. I. Pathways to Solutions, Fixed Points and Equilibria.

[24] Lu, X, Li, Y, & Su, Y. Finding periodic solutions of ordinary differential equations via

[25] Trajkovic, L, Fung, E, & Sanders, S. HomSPICE: simulator with homotopy algorithms for finding dc and steady-state solutions of nonlinear circuits. In Proceedings of the 1998 International Symposium on Circuits and Systems, ISCAS'98, 31 May- 3 June

[26] Trajkovic, L, & Willson, A. N. Theory of dc operating points of transistor networks,

[27] Mathis, W, et al. Parameter embedding methods for finding dc operating points of transistor circuits. In Proceedings of the International Specialist Workshop on Nonlin‐

[28] Allgower, E. L, & Georg, K. Numerical Continuation Methods: An Introduction. New

[29] Seydel, R, & Hlavacek, V. Role of continuation in engineering analysis. Chemical

ear Dynamics of Electronics Systems, NDES'95, July (1995). Dublin, Ireland.

NP-Completeness. San Francisco: Freeman, (1979).

Cliffs, NJ: Prentice-Hall, Inc., (1975).

(1995). Istanbul, Turkey., 95, 27-31.

systems. SIAM J. Numer. Anal. (1994).

Englewood Clifs NJ: Prentice-Hall, Inc., (1981).

homotopy method. Appl. Math. Comput. (1996). , 78(1), 1-17.

Intern. J. Electronics and Comm. (1992). , 46(4), 228-240.

Anal. (1977). , 14(4), 611-615.

phia: SIAM, (2009).

(1998). Monterey, CA.

York: Springer Verlag; (1990).

Engineering Science (1987). , 42(6), 1281-1295.

415-420.

88 Analog Circuits


## **Chapter 5**
