Design of Passive Analog Electronic Circuits using Hybrid Modified UMDA Algorithm

Hybrid evolutionary passive analog circuits synthesis method based on modified Univariate Marginal Distribution Algorithm (UMDA) and a local search algorithm is proposed in the paper. The modification of the UMDA algorithm which allows to specify the maximum number of the nodes and the maximum number of the components of the synthesized circuit is proposed. The proposed hybrid approach efficiently reduces the number of the objective function evaluations. The modified UMDA algorithm is used for synthesis of the topology and the local search algorithm is used for determination of the parameters of the components of the designed circuit. As an example the proposed method is applied to a problem of synthesis of the fractional capacitor circuit.


Introduction
In many areas of the engineering applications Estimation of Distribution Algorithms (EDA) have proved excellent capabilities of dealing with various kinds of problems of different nature.However, in the area of the evolutionary electronics utilization of the EDA algorithms has not been sufficient.Several papers focused on the evolutionary design of the whole analog circuit (the topology and the parameters of the components) have been published.In [2] Zinchenko published method of synthesis of the passive analog circuits using Univariate Marginal Distribution Algorithm (UMDA) algorithm however the method has two drawbacks.Parallel connection of the components in the branches of the circuits is not possible and with increasing number of the nodes of the designed circuit the number of the components is increasing more than necessary.Another paper focused on the evolutionary synthesis of the passive analog circuits using EDA algorithms was published by Torres [3] however the method employs UMDA algorithm only for determination of the number of the resistors, capacitors and inductors in the designed circuit.Minimal Switching Graph Problem solved using hybrid EDA algorithm method was presented in [4].The method employs the UMDA to sample start search points and a hill-climbing algorithm to find local optimum of the search space.
The proposed EDA algorithm solves the drawbacks of the method [2] mentioned above.The maximum number of the nodes and the maximum number of the components of the synthesized circuit can be set independently.Also encoding of the parallel components is possible.Furthermore the proposed method employs efficient hybrid approach which significantly reduces the number of the objective function evaluations.For the same optimization problem the number of the evaluations of the objective function of the proposed method is almost four times lower than for simulated annealing algorithm [12].The flowchart of the proposed method is presented in Fig. 1.Population P is formed of binary vectors of length 135 bits which are initialized randomly with seeding of 10 bits (n c = 10).Parameters storage e ps is formed of a vector of length 135 consisting of real numbers in the range 0,1 .Parameters storage e ps is dynamically optimized during the whole synthesis process and it is adapted to the selected topologies in the selection phase of the algorithm.The vector includes a component value for every single admittance of the used expanded fully connected ad-mittance network (see Sec. 6).During initialization phase e ps is set randomly with uniform distribution.

Hybrid Modified UMDA Algorithm
In the next step selected population P s is formed of the good individuals of the previous population P.
Probabilistic model M of selected population P s is built.For more information on building of the probabilistic model in UMDA algorithm please refer to [5].
Probabilistic model M built in the previous step is used for generation of new samples of solutions P g .The new samples are generated using Stochastic Universal Sampling method (SUS) and are repaired using the repairing method described in Sec. 3.
The generated samples are evaluated using the topological information stored in the individuals of P g and the parameters of the components stored in e ps .If the condition of execution of the local search algorithm is fulfilled, the local search algorithm tries to optimize the parameters of the current solution.If the accuracy of the current solution is improved, then storage of the parameters e ps is updated according to the results of the local search algorithm.Detailed description of the cost evaluation phase and the local search algorithm is presented in Sec. 4.
In the next step new population P is formed of the best individuals of P g and selected population P s .The described process is repeated until one of the termination criteria of the algorithm is met.

The Unitation Constraints
Generally desired specifications of an analog circuit are easier to reach using an analog circuit of higher complexity.Due to the fact the evolutionary analog circuits synthesis methods tend to evolve analog circuits with complexity as large as possible.Without restriction of the number of the components of the evolved circuit its complexity becomes higher than necessary.
Therefore the number of the components of the evolved circuit should be restricted to a user define value.Since in the proposed encoding method (section 6) the number of the components of the encoded circuit is determined by the number of the "ones" of the binary characteristic vector c the restriction of the number of the components leads to a problem with unitation constraints [6].
Then the unitation value of x is defined as Value of unitation function u(x) depends only on the number of the "ones" in an input binary vector x.The unitation values of two vectors with the same numbers of "ones" are equal.
A problem with unitation constraints is defined as solution e in which unitation value u(e) (the number of the "ones" in solution e) is restricted to a defined number [6].
As was described above the analog circuit synthesis problem has to be viewed as a problem with unitation constraints [6].Modification of Factorized Distribution Algorithm (FDA) [7] which enables solving problems with unitation constraints was described in [6].Modification of the UMDA algorithm which is able to handle the problems with unitation constraints is proposed in the text bellow.Pseudocode of the original UMDA algorithm [5] is presented in Fig. 2.
step0: Set k = 1.Generate n i 0 points randomly.step1: Select n s ≤ n i points.Compute the marginal frequencies r k;i (x i ) of the selected set.step2: Generate n i new points according to the distribution q k+1 (x) = ∏ n i=1 r k;i (x i ).Set k = k + 1. step3: If not terminated, go to step1.Fig. 2. Pseudo-code of the original UMDA algorithm [1].
In [6] to handle unitation constraints problems only generation phase (sampling) of FDA algorithm was modified.The presented approach was adopted also in the proposed modification of the UMDA algorithm.The UMDA algorithm was implemented using toolbox MATEDA 2.0 [8].Pseudo-code of the modified UMDA algorithm is presented in Fig. 3. step0: Set k = 1.Generate n i 0 points randomly.step1: Select n s ≤ n i points.Compute the marginal frequencies r k;i (x i ) of the selected set.step2: Generate n i new points according to the distribution q k+1 (x) = ∏ n i=1 r k;i (x i ).Set k = k + 1. step3: With regard to unitation constraints repair generated points.step4: If not terminated, go to step1.Fig. 3. Pseudo-code of the modified UMDA algorithm [1].
One additional step (step3) was added to the original UMDA algorithm.In step3 the generated samples are repaired to satisfy the desired unitation constraints.
The repairing function is applied to every single individual of the population of the generated samples in step2.Only n c "ones" with the highest marginal frequencies r k;i of every generated sample are accepted.The rest of the "ones" of the samples are set to zero.If the number of the "ones" of the sample is equal or lower than n c then the sample is accepted without any modification and no repairing is performed.This way the number of the "ones" (which corresponds to the number of the components of the encoded analog circuit) of every generated sample never exceeds n c .Note that the repairing function is applied only for the part of the encoding vector e which encodes the topology of the solution (characteristic vector c in Fig. 13).

The Local Search Algorithm
Evaluation of the cost values and using of the local search algorithm will be discussed in the section.Principal flowchart of this phase is presented in Fig. 4. The following procedure which is described bellow is performed for every single individual P g (i) of the population of the generated samples P g .
Based on the topology information stored in the binary vector of individual P g (i) appropriate set of parameters P 1 is loaded from parameters storage e ps and cost value c 1 of individual P g (i) is evaluated.
If the condition of execution of the local search algorithm (LSA) is fulfilled (rand < P LSA ) the LSA tries to improve accuracy of individual P g (i).The probability of execution of the LSA is set to P LSA = 0.02.The initial point of the search of the LSA is set to parameters set P 1 .The number of the objective function evaluations of the LSA is set to MaxFunEvals = 100.The results of the LSA are the cost value of the optimized solution c 2 and set of the optimized parameters P 2 .
If the LSA is successful in improving of the accuracy of individual P g (i) and its cost value was improved (c 2 < c 1 ) then value cost of individual P g (i) is set to c 2 (cost := c 2 ) and the appropriate parameters of parameters storage e ps are updated according to parameters set P 2 .
If the condition of execution of the LSA was not fulfilled (rand > P LSA ) or the LSA was not able to improve the cost value of individual G(i) (c 2 ≥ c 1 ) the resulting value cost of individual P g (i) is set to c 1 (cost := c 1 ).After performing of the described process cost value cost of individual P g (i) and updated parameters storage e ps are obtained.During the run of the algorithm the parameters stored in e ps are adapted to the topological information of the good individuals selected in the selection phase of the algorithm (Fig. 1).This way the information about the topology stored in P g and the information about the parameters stored in e ps are mutually optimized and the whole synthesis process is directed towards the promising areas of the solution space.

Application of the Method
A problem of synthesis of a fractional capacitor circuit was adopted from [9] and will be used for demonstration of the synthesis capabilities of the proposed method.The goal is to synthesize a circuit with input impedance ( 2) In Fig. 5 there is a circuit realization of function (2) as presented in [9].For the rest of the section, the circuit will be called original approximation circuit.Comparison of the magnitude and the phase of Z in of the original approximation circuit and (2) is presented in Fig. 6 and Fig. 7 respectively.Deviation of the magnitude and the phase of Z in of the original approximation circuit and (2) is presented in Fig. 8 and Fig. 9 respectively.The highest deviations of the magnitude and the phase of Z in of the original approximation circuit are presented in Tab. 1.The following sections will be focused on synthesis of the fractional capacitor circuit problem using the proposed hybrid modified UMDA method.

The Encoding Method
The used encoding method is based on the idea of fully connected admittance network.For chosen number of the nodes the fully connected admittance network is formed by connecting the admittances between all combinations of the nodes of the network.The number of the admittances of the fully connected admittance network with n n nodes can be calculated according to ( 3) Every single admittance of the fully connected admittance network can be replaced by resistor, capacitor, inductor or their parallel combination.Therefore the largest circuit which can be for chosen number of the nodes n n obtained is the circuit where every single admittance of the fully connected admittance network is replaced by parallel combination of resistor, capacitor and inductor.In this paper such circuit is denoted as expanded fully connected admittance network and includes 3 n adm components.Example of the expanded fully connected admittance network is presented in Fig. 10.The expanded fully connected admittance network can be represented using complete multigraph with three multiple edges at the most [10].Complete multigraph G c corresponding to expanded fully connected admittance network N c is presented in Fig. 11.Nodes n 0 to n 3 of network N c correspond to vertices v 0 to v 3 of complete multigraph G c .Branches of network N c correspond to edges of complete multigraph G c .For example edges e 5 (1), e 5 (2), e 5 (3) (on complete multigraph G c ) correspond to components L 5 , R 5 ,C 5 (in network N c ) respectively.Then the problem of searching of the topology of the analog RLC circuits can be defined as searching of subgraph G s on complete multigraph G c [10].
(1 ) Subgraph G s can be encoded using binary characteristic vector c of length 3 n adm .Every single bit of characteristic vector c represents including or not including of the corresponding edge of the complete multigraph G c in subgraph G s .For example complete multigraph G c is encoded using characteristic vector c of length 18 bits where c(i) = 1 for i ∈ {1, 2, . . ., 18}.Example of subgraph G s , corresponding analog circuit and its characteristic vector c are presented in Fig. 12.
e 6 (2) Encoding vector e of every single solution is formed of two parts.The first one represents the topology and is encoded using characteristic graph c approach as described above.The number of the nodes of the expanded fully connected admittance network was experimentally chosen n n = 10.To enable direct comparison of the accuracy of the circuits synthesized using the proposed method and the original approximation circuit presented in [9] the maximum number of the components of the synthesized circuit was set to the number of the components of the original approximation circuit (n c = 10).According to (3) the topology is encoded using binary characteristic vector c of length 135 bits (3 n adm ).
The second part of encoding vector e represents information about the parameters of the components and is represented by vector of real numbers p of length n c .Schematic of the used encoding vector e is presented in Fig. 13.The topological information is encoded using binary characteristic vector c = {b1, b2, b3, ..., b135} and the parameters (the values of the components) are encoded using vector of real numbers p = {dbl1, dbl2, db3, ..., dbl10}.Based on the components selected in the topological part of the information (characteristic vector c) corresponding values of the parameters are loaded from storage of the parameters e ps and copied to vector p.For example let's assume that the topology is encoded using characteristic vector c( j) = 1 for j ∈ {1, 3, 15, 18, 28, 52, 78, 92, 107, 115}.Based on the information in characteristic vector c corresponding parameters of storage of the parameters e ps will be loaded to vector of the parameters p as follows: p(k) = e ps ( j) for k ∈ {1, 2, . . ., N c } and j ∈ {1, 3, 15, 18, 28, 52, 78, 92, 107, 115} .
Based on the parameters in vector p the values of the components are calculated using formula (4) where r are values of the parameters loaded from vector of the parameters p. Formula (4) was formed to map the values of the parameters stored in vector p to suitable range of the values of the components.Since the parameters in vector p are set in the range <0,1> corresponding values of the components for the lowest r = 0 and for the highest r = 1 are v min = 0.0061 and v max = 7.3685 × 10 3 respectively.Note that formula (4) is used for all three types of components RLC.Since the used angular frequency range is from 0.01 rad/s to 100 rad/s, the values of the components are set to non-realistic values.

The Objective Function
Cost value cost is according to (7) computed as weighted summation of the magnitude and the phase differences.Difference of magnitude ∆ m is according to (5) calculated as weighted absolute value of differences of desired magnitude function f md and magnitude of current solution f mc over m = 101 frequency points in the range 0.01 rad/s to 100 rad/s.Similarly difference of phase ∆ p is according to (6) calculated as weighted absolute value of differences of desired phase function f pd and phase of current solution f pc .
Weights w cm and w cp were set to 1 and 2 respectively.Setting of weights w dm , w d p is presented in Tab.Frequency responses of the current solution f mc and f pc are obtained using nodal analysis method implemented in Matlab.

Settings of the Proposed Algorithm
The goal of the synthesis is to design a circuit which approximates function (2).The only information supplied to the system is desired magnitude and phase characteristics (2), maximum number of the used components (n c = 10), maximal number of the nodes (n n = 10) and types of the used components (resistors, capacitors, inductors).

Experiments and The Solutions
Comparison of the magnitude and the phase characteristics of Z in of the best found approximation circuit and desired function (2) are presented in Fig. 19 and Fig. 20 respectively.As can be seen in Fig. 19 to Fig. 22 the maximum deviations of the magnitude and phase responses are located at the boundaries of the used frequency range.At these frequencies the unoptimized areas of the frequency response (0 to 10 −2 rad/s and 10 2 to ∞ rad/s) affect behavior of the circuit in the area where the optimization was performed.The zeros and poles diagram of the synthesized circuit is presented in Fig. 23.All coefficients of the denominator of the approximation function of Z in are positive therefore stability of the synthesized circuits is guaranteed.Although the probability of using of all three component types (resistors, capacitors, inductors) was equal during the synthesis process, none of the circuits presented in Fig. 15 to Fig. 18 include any inductors.As the proposed algorithm was constrained to use only n c = 10 components, it seems that using only capacitors and resistors allows the method to reach lower cost values than in solutions where inductors are included.

Comparison of The Results
In the section the best synthesized approximation circuit obtained using the proposed algorithm will be compared to the original approximation circuit designed in [9] by a classical method of the analog circuits design.
Since the proposed evolutionary synthesis method was configured to use the same circuit complexity (10 components at the most) as the original approximation circuit, accuracy of both circuits can be directly compared.
Except the deviations at the boundaries of the used frequency range (as was commented above) for the original approximation circuit the highest deviation of the magnitude is ∆ m = 0.61 dB at angular frequency 0.33 rad/s.For the best solution of the proposed method the highest deviation of the magnitude is ∆ m = 0.27 dB at angular frequency 0.30 rad/s.Thus in terms of deviation of the magnitude the accuracy of the synthesized circuit is more than twice better than the original approximation circuit.Comparison of the deviations of the magnitude of Z in for both circuits is presented in Tab.The highest phase deviation inside the used frequency range is for original circuit ∆ p = 5.2 • at angular frequency 0.14 rad/s and for the synthesized circuit it is ∆ p = 1.5 • at angular frequency 0.58 rad/s.Thus phase accuracy of the best synthesized circuit is more than three times better than the original approximation circuit.Comparison of the maximum deviations of the phase of Z in is presented in Tab.In [12] the same problem of synthesis of the fractional capacitor circuit was solved using simulated annealing method.The method was able to reach solutions of the same accuracy however the number of required evaluations of the objective function was almost four times higher.Comparison of accuracy of the best solutions and numbers of evaluations of the objective function of the proposed EDA method and simulated annealing method [12] is presented in Tab. 8. simulated annealing best cost evals method [12] 2.45 440e3 the proposed best cost evals EDA method 2.43 120e3 Tab. 8. Comparison of the proposed EDA method and simulated annealing method presented in [12].

Conclusion
Automated analog circuit synthesis approach based on hybrid evolutionary method employing modified UMDA algorithm and a local search algorithm was in the paper.Used hybrid approach enables to employ specialized methods for both sub problems of different nature (synthesis of the topology and determination of the parameters).
Synthesis of the topology which is combinational optimization problem was solved using modified UMDA algorithm.Determination of the parameters which is continuous optimization problem was solved using a local search algorithm.The principle of the method is based on mutual interaction of synthesis of the topology phase (modified UMDA algorithm) and determination of the parameters phase (the local search algorithm) of the desired solution.Modification of the UMDA algorithm which allows solving problems with unitation constrains was proposed in the paper.
The proposed method was verified on the problem of synthesis of the fractional capacitor circuit introduced in [9].Presented experiments have shown that the proposed algorithm is capable to synthesize solutions with accuracy overperforming solutions obtained using a classical method of the analog circuit design given in [9].Accuracy of the magnitude of Z in of the best obtained solution was more than twice better than the original approximation circuit.Accuracy of the phase of Z in of the best obtained solution was more than three times better than the original approximation circuit.
In [12] the problem of fractional capacitor circuit realization was solved using simulated annealing method (SA).The accuracy of the circuits synthesized using SA was the same as the solutions produced using the proposed EDA method.
However SA required much higher number of the objective function evaluations.While the proposed EDA method required 120e3 objective function evaluations for synthesis of the same problem SA required 440e3 objective function evaluations.
For demonstration purposes the presented algorithm was verified using nodal analysis circuit simulator implemented in Matlab.Another improvement of the efficiency of the algorithm can achieved using Model Order Reduction Techniques [13].

Load parameters P 1 2 c 2 < c 1 cost := c 2 BasedFig. 4 .
Fig. 4. Evaluation of the cost value and employing of the local search algorithm [1].

Fig. 6 .
Fig.6.Comparison of the magnitude characteristics of Z in regarding (2) and the original approximation circuit[1].

Fig. 7 .
Fig. 7. Comparison of the phase characteristics of Z in closer (2)and the original approximation circuit[1].

ωFig. 8 .
Fig. 8. Deviation of the magnitude of Z in of the original approximation circuit [1].

Fig. 9 .
Fig.9.Deviation of the phase of Z in of the original approximation circuit[1].

6 LFig. 12 .
Fig. 12. a) Example of subgraph G s b) analog circuit corresponding to G s c) characteristic vector c of G s [1].

Fig. 19 .
Fig. 19.Comparison of the magnitude of Z in of the best found approximation circuit and (2) [1].

Fig. 20 .
Fig. 20.Comparison of the phase of Z in of the best found approximation circuit and (2) [1].Absolute values of the deviations of the magnitude and the phase of Z in of the best synthesized circuit are presented in Fig.21.and Fig.22respectively.

Fig. 21 .Fig. 22 .
Fig. 21.Deviation of the magnitude of Z in of the best found approximation circuit [1].
Tab. 1.The highest deviations of the magnitude and the phase of Z in of the original approximation circuit.
2. All weight coefficients were set experimentally.
Tab. 2. Setting of weights w dm and w d p .
All parameters of the synthesized problem and settings of the proposed algorithm are summarized in Tab.3.Tab.3. Settings of the proposed algorithm.
There were executed 20 instances of the proposed algorithm.Average running time of a single execution was 11 min.The cost values of the solutions are presented in Tab. 4.
Three highest deviations of the magnitude and the phase characteristics of Z in of the best synthesized circuit are summarized in Tab. 5.The highest deviations of the magnitude and the phase of Z in of the best synthesized approximation circuit.
Comparison of the deviations of the magnitude of Z in of the original approximation circuit and the best synthesized circuit.
Tab. 7. Comparison of the deviations of the phase of Z in of the original approximation circuit and the best synthesized circuit.