Circuit Generation

In this notebook, we demonstrate how to generate a pool of equivalent circuit models from electrochemical impedance spectroscopy (EIS) measurements that best fit the data.

Set up the environment

AutoEIS relies on EquivalentCircuits.jl package to perform the EIS analysis. The package is not written in Python, so we need to install it first. AutoEIS ships with julia_helpers module that helps to install and manage Julia dependencies with minimal user interaction. For convenience, installing Julia and the required packages is done automatically when you import autoeis for the first time. If you have Julia installed already (discoverable in system PATH), it’ll get detected and used, otherwise, it’ll be installed automatically.

Note

If this is the first time you’re importing AutoEIS, executing the next cell will take a while, outputting a lot of logs. Re-run the cell to get rid of the logs.

import matplotlib.pyplot as plt

import autoeis as ae

ae.visualization.set_plot_style()

Load EIS data

Once the environment is set up, we can load the EIS data. You can use `pyimpspec <https://vyrjana.github.io/pyimpspec/guide_data.html>`__ to load EIS data from a variety of popular formats. Eventually, AutoEIS requires two arrays: Z and freq. Z is a complex impedance array, and freq is a frequency array. Both arrays must be 1D and have the same length. The impedance array must be in Ohms, and the frequency array must be in Hz.

For convenience, we provide a function load_test_dataset() in autoeis.io to load a test dataset. The function returns a tuple of freq and Z.

freq, Z = ae.io.load_test_dataset(preprocess=True)

[21:01:52] WARNING  10% of data filtered out.

Note

If your EIS data is stored as text, you can easily load them using numpy.loadtxt. See NumPy’s documentation for more details.

Let’s take a look at the test dataset before we proceed:

ae.visualization.plot_impedance_combo(freq, Z);

Generate equivalent circuits

Now that we have loaded the EIS data, we can generate a pool of candidate equivalent circuits using the generate_equivalent_circuits function. The function takes the impedance data and frequency as input and returns a list of equivalent circuits. It also takes many optional arguments to control the circuit generation process. The most important ones are:

• iters: Number of circuits to generate.

• complexity: Maximum number of elements in a circuit.

• terminals: Type of circuit components to use (e.g., R, C, L, or P).

• parallel: Whether to run the circuit generation in parallel.

• tol: Tolerance for accepting a circuit as a good fit.

The function uses a gene expression programming (GEP) algorithm to generate the circuits. The GEP algorithm is a genetic algorithm that evolves circuits by combining and mutating genes. The algorithm starts with a population of random circuits and evolves them over many generations to find the best circuits that fit the data. The following parameters control the GEP algorithm:

• generations: Number of generations to run the genetic algorithm.

• population_size: Number of circuits in the population for

The default values for these parameters are usually good enough for most cases. However, you can adjust them to get better results (e.g., increase both arguments if you’re not satisfied with the generated circuits, or decrease them if you want to speed up the process).

Note

Circuit generation is a lengthy process. It may take one minute per iteration on a modern CPU. We recommend generating at least 50 circuits to get a good pool of candidate circuits, which may take about an hour.

Note

By default, terminals is set to "RLP", i.e., the algorithm searches for circuits with resistors, inductors, and constant-phase elements. You can add capacitors by setting terminals="RCLP", but since a capacitor is a special case of a constant-phase element, it’s not necessary, and it makes the search less efficient.

kwargs = {
    "iters": 12,
    "complexity": 12,
    "population_size": 100,
    "generations": 30,
    "terminals": "RLP",
    "tol": 1e-3,
    "parallel": True,
}
circuits_unfiltered = ae.core.generate_equivalent_circuits(freq, Z, **kwargs)
circuits_unfiltered

	circuitstring	Parameters
0	[P1,P2]-R3	{'P1w': 1.6594864716689052e-06, 'P1n': 0.97132...
1	P1-R2-P3-[P4,[P5,R6]]	{'P1w': 0.0001709832771605446, 'P1n': 0.999990...
2	R1-L2-[P3,P4]	{'R1': 143.130151058375, 'L2': 8.9946264876557...
3	[P1,P2-L3]-P4-P5-R6	{'P1w': 1.7731744433018003e-06, 'P1n': 0.96184...
4	[[P1,R2],P3]-P4	{'P1w': 1.634327115205153e-06, 'P1n': 0.979494...
5	[R1,P2-[R3,P4-R5]]	{'R1': 6964972.938919069, 'P2w': 5.27530336220...
6	R1-[[[R2,P3],R4-L5],P6]	{'R1': 142.07487772062197, 'R2': 999016633.188...
7	[P1,[L2,R3]-P4]	{'P1w': 4.058893113496399e-07, 'P1n': 0.166711...
8	[[R1,P2],P3]-R4	{'R1': 966022225.0726572, 'P2w': 1.66878355310...
9	[P1,[P2,L3]-P4-R5-L6]	{'P1w': 5.019807735771084e-07, 'P1n': 0.212865...
10	R1-[L2-P3,P4]	{'R1': 142.06941780245384, 'L2': 5.0, 'P3w': 4...

Convergence

The Circuit generation algorithm is sensitive to the tol parameter, meaning that the order of magnitude of the tol needs to be proportional to the order of magnitude of the impedance data. There’s no one-size-fits-all value for tol, and we’re trying to make the algorithm tol-agnostic in future releases. For now, we’ve hacked a heuristic that internally scales the tol based on the impedance data. Nevertheless, you may still need to adjust the tol if you end up with no circuits (increase tol) or too many circuits (decrease tol). The default value is 1e-3. When increasing or decreasing tol, try doubling or halving the value to see if it helps.

Filter candidate equivalent circuits

Note that all these circuits generated by the GEP process probably fit the data well, but they may not be physically meaningful. Therefore, we need to filter them to find the ones that are most plausible. AutoEIS uses “statistical plausibility” as a proxy for gauging “physical plausibility”. To this end, AutoEIS provides a function to filter the candidate circuits based on some heuristics (read our paper for the exact steps and the supporting rationale).

circuits = ae.core.filter_implausible_circuits(circuits_unfiltered)
circuits

	circuitstring	Parameters
0	[P1,P2]-R3	{'P1w': 1.6594864716689052e-06, 'P1n': 0.97132...
1	P1-R2-P3-[P4,[P5,R6]]	{'P1w': 0.0001709832771605446, 'P1n': 0.999990...
2	R1-L2-[P3,P4]	{'R1': 143.130151058375, 'L2': 8.9946264876557...
3	[P1,P2-L3]-P4-P5-R6	{'P1w': 1.7731744433018003e-06, 'P1n': 0.96184...
4	R1-[[[R2,P3],R4-L5],P6]	{'R1': 142.07487772062197, 'R2': 999016633.188...
5	[[R1,P2],P3]-R4	{'R1': 966022225.0726572, 'P2w': 1.66878355310...
6	R1-[L2-P3,P4]	{'R1': 142.06941780245384, 'L2': 5.0, 'P3w': 4...

Let’s see how well the generated circuits fit the data. You can either use the parameters’ values at the end of the GEP process (stored in the circuits dataframe), or use fit_circuit_parameters to further refine the parameters (recommended).

Note

Normally, the found circuits are good enough, but since we didn’t run the algorithm for long enough (to not timeout our CI on GitHub), we will use a custom circuit for evaluation. If you’re running this notebook on your own data, try using iters >= 200 together with a more stringet tol to get a good pool of circuits. We’re currently working on changing the evolutionary algorithm backend to speed up the process, so you no longer need to wait for hours to get a good pool of circuits.

use_custom_circuit = False

if not use_custom_circuit:
    circuit = circuits.iloc[0]["circuitstring"]
    p = circuits.iloc[0]["Parameters"]
    # Refine the circuit parameters
    p = ae.utils.fit_circuit_parameters(circuit, freq, Z, p0=p)
else:
    circuit = "R4-[P1,R3-P2]"
    p = ae.utils.fit_circuit_parameters(circuit, freq, Z)

# Simulate Z using the circuit and the fitted parameters
circuit_fn = ae.utils.generate_circuit_fn(circuit)
Z_sim = circuit_fn(freq, list(p.values()))

# Plot against ground truth
fig, ax = plt.subplots(figsize=(5.5, 4))
ae.visualization.plot_nyquist(Z_sim, fmt="-", ax=ax, label="simulated")
ae.visualization.plot_nyquist(Z, fmt=".", ax=ax, label="data");
ax.set_title(circuit)

Text(0.5, 1.0, '[P1,P2]-R3')