Synthesizing Electrically Equivalent Circuits for Use in Electrochemical Impedance Spectroscopy through Grammatical Evolution

: Electrochemical impedance spectroscopy (EIS) is an important electrochemical technique that is used to detect changes and ongoing processes in a given material. The main challenge of EIS is interpreting the collected measurements, which can be performed in several ways. This article focuses on the electrical equivalent circuit (EEC) approach and uses grammatical evolution to automatically construct an EEC that produces an AC response that corresponds to one obtained by the measured electrochemical process(es). For ﬁtting purposes, synthetic measurements and data from measurements in a realistic environment were used. In order to be able to faithfully ﬁt realistic data from measurements, a new circuit element (ZARC) had to be implemented and integrated into the SPICE simulator, which was used for evaluating EECs. Not only is the presented approach able to automatically (i.e., with almost no user input) produce a more than satisfactory EEC for each of the datasets, but it also can also generate completely new EEC conﬁgurations. These new conﬁgurations may help researchers to ﬁnd some new, previously overlooked ongoing electrochemical processes.


Introduction
Electrochemical impedance spectroscopy (EIS) [1] is an important electrochemical technique that is used to detect changes and ongoing processes in a given material. It can be applied in fields such as solid oxide fuel cell (SOFC) development, monitoring protective coatings properties, measuring the effectiveness of electrodes and more. At its core, EIS consists of measuring an AC signal and presenting it as a function of frequency at a constant amplitude [2]. The results of an EIS study offer a large amount of additional data such as the distribution of relaxation times (DRT) [3][4][5] and the distribution of diffusion times (DDT) [6,7], which can be used to gain deeper insights into ongoing electrochemical processes. The main challenge of EIS is interpreting the collected measurements, which can be conducted in several ways. For example, one can use Electrical Equivalent Circuits (EECs) or apply the DRT and DDT approaches (Scheme 1). The EEC application is custom in EIS study [8,9], and there are also modern articles on how to improve and modify EEC analysis [10][11][12][13][14][15] EEC analysis.
This article focuses on the EEC approach and proposes an innovative method that automatically constructs an EEC that produces an AC response, which has the same impedance characteristic(s) as the measured electrochemical process(es). The resulting circuit also allows a glimpse into the type of material that was measured (conductor, conductor with an oxide layer, an SOFC electrode, etc.). Finding an EEC manually is a difficult task that requires in-depth knowledge of electrical circuit elements (and system under the investigation) and comprises several steps. One must first select the appropriate EEC topology and parameters and then tune the parameters by using an optimization algorithm such as, for example, the Levenberg-Marquardt algorithm (LMA) [16][17][18] or the Nelder-Mead method (NMA) [19,20]. Once this optimization is complete, the optimized EEC is selected as the representation of the measured system. This process can be somewhat simplified with the use of sets of pre-determined circuits, but it is unlikely that such limited sets will be able to fit all the possible signals that one usually encounters in EIS. Scheme 1. Different approaches commonly applied in EIS study.
In this paper, a novel approach is presented that aims to automate all of the steps in the EEC approach (Scheme 1) by autonomously generating and optimizing a matching EEC. Expert knowledge of circuits and circuit simulations is combined with evolutionary techniques in order to avoid the severe limitations of a conventional EEC approach, which looks for a solution only from a limited number of EECs with predetermined topology. The presented approach [21] is based on Grammatical Evolution (GE), which is a well known Evolutionary Computation technique and allows the user to pick any number of electrical components, sources, and even custom-created elements, such as a ZARC element (see, e.g., [1,4,5]). An important advantage of the proposed approach is that once the setup is complete-which consists only of the list of desirable circuit elements, the objective function and the EIS data to be fitted-the fitting process is autonomous and requires no further intervention on the part of the user. Note that the proposed approach lowers the required skill level of the practitioner; once a working setup has been found (i.e., the right set of elements for the current measurement set), all the subsequent fittings require no additional user input except for the new set of EIS values to be fitted.
The next section briefly overviews grammatical evolution and explains the grammar that was used to describe EECs and the slightly adapted genetic operators of mutation and crossover. The used genetic parameters, objective function, and the sets of data used for fitting the EECs are described in Section 3. Finally, the results of fitting the sets of data are presented, and the paper is concluded with a short discussion of the value of the presented contribution.

Grammatical Evolution
The presented approach of automatic circuit generation is based on grammatical evolution [22], which is an evolutionary computation technique. Basically, evolutionary computation is a family of algorithms for global optimization inspired by biological evolution. A typical evolutionary computation algorithm generates an initial set of candidate solutions and updates them iteratively. Each new generation is produced by applying selection and mutation operators. As a result, the overall fitness (specified by a special fitness function) of the population will gradually increase until an acceptable solution is obtained.
Unlike many conventional EA algorithms, grammatical evolution uses populations of variable-length binary string genomes (also known as chromosomes, where each codon is a consecutive group of eight bits, representing an integer value) to select production rules in a Backus-Naur form (BNF) grammar definition. Thus, GE does not perform the evolution process on the actual circuits but rather on binary stings, which increases the flexibility of the algorithm. By using appropriate BNF grammar, one can easily integrate the expert/domain knowledge in order to speed up the process (such as limiting the ZARC values as shown in [4,5,23]).
The GE approach was implemented in the Python programming environment with the addition of the PyOpus package to access the SPICE circuit simulator [24] and is described in-depth in articles [21,25]. Basically, one step of the GE algorithm produces a generation of binary strings, each of which is then used to select production rules in a BNF grammar to build an EEC. Thus, each of the obtained EECs is then evaluated in the SPICE simulator for its fitness.

The Grammar
The Backus-Naur form is a notation for specifying a language grammar in the form of so-called production rules [26]. Production rules include terminals, which are elements that appear in the language (e.g., 0, 1, +, input, etc.), and nonterminals, which must be expanded in one or more terminals and/or nonterminals. BNF grammar is often expressed by the tuple {N, T, P, S}, where N and T stand for the sets of terminals and nonterminals, respectively, P is a set of production rules mapping nonterminals to terminals, and S is a start symbol, which must be a nonterminal. When there is more than a single production that can be applied to a certain nonterminal, the different production rules are delimited by a vertical bar.
In order to be able to evaluate an EEC produced in the GE process, one needs to build a SPICE netlist from a binary string, which is then evaluated in the SPICE simulator. Standard components (i.e., resistors and capacitors) as well as the EIS specific ZARC element are employed to build EECs. Each component has at least two parameters-the numbers of the connecting ports and the element value(s) (i.e., resistances, capaticances, etc.). This is the BNF used to produce a SPICE netlist: and P is summarized in Table 1. The first row in the table contains a production rule for the <netlist> nonterminal. The rule ensures that a netlist is composed of up to twelve circuit elements (or parts), separated by spaces. This was performed to prevent bloating (i.e., the uncontrolled growth of the size of a circuit due to unnecessary parallel/serial elements). The second row shows that each of the twelve elements (parts) can become either one of the three actual elements (i.e., resistor, capacitor, or ZARC) or a nonexistent element (None). It can be observed from the next three rows how the three elements are specified. They contain two connecting ports (determined by the <gpair> rule) and the value(s) appropriate for a specific element. A ZARC element contains specialized numeric values (zexp and znum). These are obtained according to expert knowledge (as described in [4,5,23]) in order to narrow down the solution search space. For example, as we are especially interested in the capacitive characteristic of electrochemical processes, we only allow ZARC elements with n between 0.50 and 0.99. Table 2 shows the value ranges of the elements that can be produced using the rules from Table 1.
The last row of Table 1 contains another hard-wired optimization. Instead of specifying each connection port separately, the <gpair> production rule was created which holds all the feasible port combinations (thus, producing some so-called composite terminals made up from terminals from (2)). In that way, two benefits were achieved-a higher chance of creating a working circuit (i.e., where one actually receives a signal from input to output) and a significantly lower chance of creating an illegal circuit with closed loops. Note that there is still a certain (small) probability that the algorithm produces an invalid circuit when using the above rules-in some conditions, the process of mapping the codons to the actual circuit can result in an element that is connected to itself. Such a circuit can create serious problems during the simulation (i.e., the simulator could either crash or become stuck in an endless loop). In order to avoid such situations, the system was augmented with a special subroutine that detects such circuits and simply eliminates them from the evolutionary process even before they become evaluated for their fitness. Some may argue that this is not the best possible strategy since it may result in a loss of important genetic material at the beginning of the run. That being said, this was still the most efficient solution. Devising a set of production rules that would prevent this problem is simply not feasible.

Mapping Process Examples
For convenience, we now present two examples of mapping an individual produced during the evolutionary process to an actual electrical circuit.

A Single Element Example
Let us start with a simple example of mapping a chromosome to a single element. Assume one has a chromosome {12, 127, 209, 21, 76} that they want to map to a single element using the grammar from Section 2.1. Note that, as the chromosome is to be mapped to a single element, rather than a complete netlist, the start symbol will be S = {part}. The goal is to map the start symbol onto terminals by reading codons (i.e., integer values) from a chromosome from which a corresponding production rule is selected by using the following modulo operation.
(rule number) = (codon integer value) mod (number of rules for the current nonterminal) As the start symbol has four possible rules as shown in the second row of Table 1, the codon value of 12 selects the first production rule (i.e., 12 mod 4 = 0), hence choosing a resistor. The remaining four codons are then used to map the four nonterminals contained in the production rule for <res> to the corresponding four terminals. The entire mapping process is summarized in Table 3. The first column lists the complete chromosome, while the right operands and results of the modulo operations are given in the third column. The last column provides the actual selected terminals by using the corresponding production rule on nonterminals in the second row. The resulting element is a resistor connected to ports 1 and 4 with a resistance value of 2 kΩ (i.e., rXX (1 4) 02e3).

Genetic Operations and Circuits
Evolutionary computation features several operators used to manipulate and combine individual solutions. For the purposes of working with EEC, these operators have to be adjusted accordingly, most noticeably mutation and crossover.

Circuit Mutation
Mutation results in a change in the chromosome of the individual, which is reflected in the final interpretation of the chromosome. The change can be minor or major depending on which codon is selected to mutate. A minor change occurs when, for example, a chromosome from Section 2.2.1 changes from {12, 127, 209, 21, 76} to {12, 127, 209, 70, 76}. The effect of this change is that the resistor changes its value from two kiloohms to one kiloohm, as shown in Figure 3.  Table 4.

Before Mutation
Rxx (1 4) 2e3 Ω An even more drastic mutation occurs if the same element changes to {119, 127, 209, 21, 76}. This string now represents two empty elements (119 and 127 both map to rule number 3, None), and the next codon (209) is used to determine the next element in the circuit, which will be a capacitor connected to ports 1 and 3 (mapped from <gpair> using the codon with the value of 21).

Circuit Crossover
While mutation targets a single circuit, crossover takes two different circuits (parents) and swaps their elements in order to create new circuits (children) that might be better suited to solving the given problem.
The system first selects two appropriate circuits (usually the circuits that are amongst the better performing ones). In the example (shown in the top row of Figure 5), there are two circuits with identical topology but different component values. In the first step, a random element from Parent 1 is selected-the 13 Ω resistor drawn in solid black. The algorithm then searches for a similar element in the second circuit (Parent 2) and finds the 400 Ω resistor. Since two elements of the same type could be found, the crossover can proceed and the two elements are swapped. The chromosomes are adjusted accordingly so as to reflect the change. The resulting circuits (show in the bottom row of Figure 5) are then used as potential candidates in the next generation and (hopefully) perform better than their parents.  Note that we used a slightly modified version of crossover, which focuses on the phenotype instead of considering all available codons. Usually, crossover selects two completely random codons (or substrings of codons) and swaps them. This proved to be quite destructive for the circuits as it usually resulted in the circuit becoming nonfunctional either due to disconnected elements, impossible values, short circuits, or other similar drastic changes. Therefore, we limited the crossover to elements and values in order to preserve circuit integrity and increase the chance of converging towards a working solution. Table 5 shows the GE parameters that were used in the experiments. The parameters were chosen based on the previous experience with automatic creation of electronic circuits using this method [21]. For this research, we selected Sheppard's objective [27] function as it is commonly used by researchers for EEC fitting purposes (see, for example, [11]), which makes it easier to compare the results with those in the literature.

Objective Function and Evaluation
An objective function is used to calculate the cost or fitness of a solution. The cost of an individual reflects how well the individual solves a given problem-the lower the cost the closer the solution to the desired behavior. How the objective function is selected and/or constructed depends strongly on the problem at hand. In some cases, the objective function represents a simple numerical comparison of points between two curves while it can be a complex multi-objective function in others that compares different aspects of the solution such as cut-off frequency, damping, gain, and combines them into a single measure. It was already mentioned that, in the experiments, we employed the objective function proposed by Sheppard et al. [27], which is commonly used for solving Complex Nonlinear Least Squares (CNLS) problems (e.g., [8,11]).
Here, m, y exp j , y com j , and w j are the number of data points, the jth value of experimental impedance data to be fitted, the jth value of impedance data computed by the proposed algorithm, and the weighting modulus factor [28] associated with the jth data point. The goal of the evolutionary algorithm is to minimise the value of S.

Data Sets Used for Evaluation
The experimental impedance data to be fitted were obtained either via mathematical models (so-called synthetic data) or from measurements in a realistic environment [5].

Randles Circuits
The first two sets of synthetic data were generated using two circuits that contain only resistors and capacitors, without any constant phase elements. The first of the two circuits contains three elements (an equivalent of one Randles circuit [23]), and the second one contains five elements (an equivalent of two serially linked Randles circuits). The impedance characteristics and the corresponding Randles circuits can be observed in Figures 6 and 7.

Cole-Cole Model
A more complex Cole-Cole (ZARC) [29] model was employed for the next set of synthetic data. This model is widely used in Distribution of Relaxation Time (DRT) studies [3][4][5] and can also be applied when developing new algorithms for solving CNLS problems [13]. The following equation was used to compute single (K = 1) and double (K = 2) ZARC element data: where ω j is an angular frequency associated with the jth impedance data point, i is the imaginary unit, R k is the kth resistance, n k is a factor related to a constant phase element [11], and τ 0,k is the kth time constant. ZARC data in this study were computed by using (6) and values in Table 6. An example of such synthetic data can be observed in Figure 8 where a double ZARC (i.e., k = 2) was used to generate the impedance characteristic. This characteristic is defined by two suppressed semi-circuits that can be frequently found in EIS study [5,30]. Table 6. Parameters used to compute single and double ZARC data.

Solid Oxide Fuel Cell Data
Realistic data were obtained from the measurements of solid oxide fuel cells, which were already used in the previous work [5]. The measured cells were industrial-sized and characterized by an active surface of 80 cm 2 and an operating temperature of 800 • C. Humidified hydrogen was used as a fuel for a fuel electrode, and air was used as fuel for an air electrode. More details regarding the experimental setup can be found in [31]. The impedance characteristic obtained from the measurements is shown in Figure 9. Note that, compared to synthetic data, this set features a lot less comparison points since it was obtained from real measurements. We were, therefore, faced with two optionseither to evaluate the model using only the original data points or try to increase the amount of points using interpolation. We opted for interpolation since it permits better accuracy with respect to the evaluation of the results (more data points help to match the curve with better accuracy). We used the linear interpolation included in the Numpy Python package. Incidentally, without interpolation, the algorithm was unable to evolve an appropriate EEC since almost any random circuit fitted in the original twenty data points quite well. Using interpolation, every single run produced an acceptable solution.

Results
This section presents the results of running the grammatical evolution algorithm in trying to approximate the data provided in the previous section.

Randles Impedance Characteristic Matching
After initial fiddling with system settings and parameters (maximum number of elements and a few others), the system was able to consistently produce matching circuits with a good fit with the Sheppard's value of 0.29 or lower. What is more, at least half of the population (i.e., 150 circuits) produced viable (and different) solutions, which means that the system is able to produce several candidates in a single run. Figures 10 and 11 show that the matched impedance characteristics (dashed orange) overlaid over the original (solid blue) characteristic and layouts of the best matching circuits for each case, respectively.
Note that the presented approach generates some additional elements that can be later removed with a bit of expert knowledge. An example of such a simplification is shown in Figure 12, where the elements that have no effect since they are too large to receive any current were removed. Thus, the capacitor with 29 F was removed since it is effectively an open circuit, and no current passes through it at measured frequencies. The entire bottom subcircuit (the two resistors of 290 MΩ and 2.8 MΩ and the 560 µF capacitor) can also be removed since there will be no current over these elements as there is a path with much lower resistance (i.e., via the 13.3 Ω and 40 Ω resistors).

Cole-Cole Fitting
Cole-Cole data revealed a serious problem with the system-any simulation with a n k factor lower than 0.9 failed to produce a working circuit. However, this is not a surprise as such a factor requires a Constant Phase Element (CPE) in order to be included in a circuit and in order for the circuit to be able to produce the desired characteristic. Note that a CPE features frequency dependant values, which cannot be carried out using basic electronic components (i.e., resistors and capacitors) nor does such an element exist in the SPICE simulator. Fortunately, SPICE lends itself to implementing and integrating custom circuit elements; therefore, we created a model for the ZARC element based on (6) and added it to the Spice Opus circuit simulator as a built-in device. The model was implemented via the XSPICE extensions, which allow the user to load user-defined devices from dll files. However, because EIS only requires a small-signal AC analysis, a transient analysis model was omitted because of its undue complexity and the fact that it is not needed in the simulations. It is quite crucial that we built this model since it enabled us to continue to use the SPICE simulator. We believe this is the first implementation of the ZARC element in SPICE since we have not found it mentioned anywhere in the literature.
The results obtained with the updated engine are once again spot on, as observed in Figure 13.
The resulting circuit is shown in Figure 14 with ZARC values provided in Table 7. Again, the circuit was simplified (by removing unnecessary parallel/serial resistors and capacitors), but the resulting circuit still differs from the expected topology of two serial ZARC elements. The issue here lies in the fact that the presented system creates a circuit with a fitting impedance characteristic and matching time constants. It succeeds in this task since the results (extracted via Python Cole-Cole decomposition) indicate that the resulting circuit ( Figure 14) yields time constants values (0.007 and 0.00008) that are close to the ones (0.01 and 0.0001) used to prepare the original data (see Table 6). The final metric value is, therefore, slightly larger in the range of 10 due to a more complex problem (multiple ZARC elements), but the resulting circuit meets the expectations. We believe that this value could be further lowered with further fine tuning of the parameters, such as increasing the accuracy of element values (from two digits two four, using all possible exponent values, and increasing the number of generations in our experiment).

Solid Oxide Fuel Cell Fitting
As already noted in the last paragraph of Section 3, additional comparison points using interpolation had to be generated because the measured points were too few for the algorithm to be able to find a matching circuit. Taking more measured data points requires extra time, which prolongs DRT analysis that is vital for monitoring SOFC health. After the interpolation, the best obtained matching impedance characteristic is demonstrated in Figure 15, with the dashed line representing the fitted characteristic. Please note that the line appears smoother since it was generated via an AC sweep, while the original characteristic (solid line) comes from connecting the original (fewer) data points. Nevertheless, the characteristic was fitted with the Sheppards value of 300 and we believe that-if we had more measured data points available-the Sheppards value would be lower due to the disuse of interpolation. The resulting circuit (shown in Figure 16) is fairly simple given the complexity of input data. Regardless, a fitting circuit and DRFT characteristic was found without much trouble beyond interpolating the data points. This is significant since it opens many new possibilities of quickly fitting measured data and finding the relevant time constants.

Discussion
Grammatical evolution was used as an automatic method for creating a matching EEC circuit without much additional user input. Basically, the research started from previous work conducted on evolving electrical circuits [21], which turned out to be quite appropriate for the problem at hand. Amongst the most important additions that had to be implemented for this research was modified grammar, a new metric, and the implementation of a new element for the SPICE simulator.
In a series of experiments, we successfully created several solutions that fitted the supplied data-synthetic as well as real-with a sufficient level of accuracy. Some of the resulting circuits were similar to the known EEC solutions while others present a possible new standard model that should be further investigated by experts in the EIS field.
We believe that the presented approach offers at least two important benefits. Firstly, it automates the EEC creation process and renders the entire method accessible even to less experienced researchers who are only starting in this field. Secondly, it might produce new standard EEC circuits and help identify new, previously overlooked ongoing electrochemical processes.