Skip to content

Reaction process optimization

Electrochemical flow hydrogenation

In this tutorial, we’ll use GraphFLA to explore how process settings affect an electrochemical reaction. To work with ordered variables, we’ll compare concentration, temperature and residence time in a flow reactor, then analyse the landscape formed by their combinations.

Open In Colab

Download notebook Notebook + data

Run the notebook in Google Colab, or download it to run locally. It downloads its data when no data/ folder is beside it.

1. Setting up

On Google Colab, the cell below installs GraphFLA. In any environment, it also downloads the dataset into a data/ folder next to the notebook if the files are not already there.

We then import GraphFLA for landscape construction and analysis, and pandas for working with the data.

import sys
from pathlib import Path
from urllib.request import urlretrieve

if "google.colab" in sys.modules:
    %pip install -q graphfla==0.4.0

DATA_URL = "https://raw.githubusercontent.com/COLA-Laboratory/GraphFLA/v0.4.0/tutorials/datasets/data/"
Path("data").mkdir(exist_ok=True)
for name in ["flow_semihydrogenation.csv"]:
    if not Path("data", name).exists():
        urlretrieve(DATA_URL + name, Path("data", name))

from graphfla import analysis
from graphfla.landscape import OrdinalLandscape
import pandas as pd

2. Loading the dataset

Selective hydrogenation aims to add hydrogen to a molecule while favouring a particular product. To explore how operating conditions affect this process, we’ll use the study by Ashikari et al. (2022). The researchers used an electrochemical flow reactor to convert diphenylacetylene to cis-stilbene, testing 54 combinations of concentration, temperature and residence time.

We’ll maximise their integrated current-efficiency/selectivity index, which rewards both efficient use of current and selectivity towards cis-stilbene. It measures this combined objective rather than isolated product yield.

Column Sampled levels
concentration_M 0.25, 0.5 and 1 M.
temperature_C 30, 40, 50, 60, 70 and 80 °C.
residence_time_s 6, 12 and 24 s.
integrated_index_pct Reported index on the percentage scale.

Let’s load Table S5. We retain its reported index of 64 at 0.25 M, 80 °C and 12 s, although that entry is inconsistent with its displayed components.

df = pd.read_csv("data/flow_semihydrogenation.csv")
df.head()

Output

concentration_M temperature_C flow_rate_mL_min residence_time_s current_efficiency_pct selectivity_cis_stilbene_pct integrated_index_pct
0 1.0 30 1.00 6 95 92 87
1 1.0 30 0.50 12 97 92 89
2 1.0 30 0.25 24 >99 92 92
3 1.0 40 1.00 6 86 95 82
4 1.0 40 0.50 12 90 94 85

3. Preparing the inputs

To construct a landscape, GraphFLA needs two aligned inputs:

  • X: one row per configuration and one column per variable.
  • f: one measured or calculated outcome for each row of X.

The three independent process settings form X, and the reported index forms f. Their complete Cartesian product contains 3 × 6 × 3 = 54 points. Flow rate determines residence time in this setup, so we do not include both as independent variables.

X = df[["concentration_M", "temperature_C", "residence_time_s"]]
f = df["integrated_index_pct"]
X.head()

Output

concentration_M temperature_C residence_time_s
0 1.0 30 6
1 1.0 30 12
2 1.0 30 24
3 1.0 40 6
4 1.0 40 12

4. Constructing the landscape

We can use OrdinalLandscape because each process variable has ordered levels. Neighbours move to the next sampled level on one axis while retaining the other settings.

Unequal physical spacings still count as one rank step: 0.25→0.5 M and 0.5→1 M are both neighbouring changes. We set maximize=True to favour larger integrated indices.

landscape = OrdinalLandscape(maximize=True)
landscape.build_from_data(X, f, verbose=False)

Output

OrdinalLandscape(maximize=True)

5. Inspecting the landscape

Let’s first look at what we’ve built. Printing the landscape shows its variables, configurations, improving edges and local optima:

print(landscape)

Output

OrdinalLandscape(kind='ordinal'): 3 variables, 54 configurations, 111 edges, 2 local optima

For a closer look at individual process conditions, we can call get_data(). The outcome appears as fitness; out_degree counts directly improving moves, and is_lo indicates local-optimum membership.

landscape.get_data().head()

Output

concentration_M temperature_C residence_time_s fitness plateau_id plateau_size in_degree out_degree is_lo
0 1.0 30 6 87.0 -1 1 2 1 False
1 1.0 30 12 89.0 -1 1 3 1 False
2 1.0 30 24 92.0 0 2 2 0 False
3 1.0 40 6 82.0 1 2 1 2 False
4 1.0 40 12 85.0 -1 1 2 3 False

We can also summarise the integrated indices with pandas:

landscape.get_data()["fitness"].describe()

Output

count    54.000000
mean     76.203704
std      14.694858
min      44.000000
25%      64.000000
50%      79.500000
75%      87.750000
max      96.000000
Name: fitness, dtype: float64

To inspect the setting with the highest reported index, we can use the global-optimum node index:

landscape[landscape.go_index]

Output

{'concentration_M': np.float64(1.0),
 'temperature_C': np.int64(70),
 'residence_time_s': np.int64(12),
 'fitness': np.float64(96.0),
 'plateau_id': 2,
 'plateau_size': 2,
 'in_degree': 4,
 'out_degree': 0,
 'is_lo': True}

6. Analysing the landscape

Next, we’ll count local optima, examine local fitness similarity and interaction orders, and assess the trend towards the best observed configuration.

6.1 Number of local optima

We can start with landscape.n_lo. Each connected neutral optimum plateau counts once. A local optimum has no improving move out of it, including through its neutral plateau. Multiple optima mean successive improvements can end at different configurations.

print(f"Number of local optima: {landscape.n_lo}")

Output

Number of local optima: 2

6.2 Fitness autocorrelation

How similar are responses at neighbouring steps? autocorrelation() measures persistence along sampled walks. Higher values indicate more persistence under the same walk settings.

We use 200 walks of up to 20 visited configurations, lag one and a fixed seed. Walks traverse stored edges in either direction; separately stored neutral pairs are excluded.

walk_autocorrelation = analysis.autocorrelation(
    landscape, walk_length=20, walk_times=200, lag=1, seed=42,
)
print(f"Lag-1 fitness autocorrelation: {walk_autocorrelation:.4f}")

Output

Lag-1 fitness autocorrelation: 0.6308

6.3 Walsh–Hadamard decomposition

Does changing temperature have the same effect at every concentration? walsh_hadamard() fits separate contributions from each variable and interactions between pairs. It treats the sampled levels as discrete states, rather than assuming a linear response per degree or second.

Compare first- and second-order training r2 using delta_r2. Three-way interactions remain outside this fit. model_variance_fraction describes the fitted model under uniform combinations of the sampled levels; it does not measure a separate experimental variance component.

wh = analysis.walsh_hadamard(landscape, max_order=2)
wh["order_summary"]

Output

order r2 delta_r2 rmse n_terms rank alpha model_variance_fraction n_nonzero
0 0 0.000000 0.000000 14.558159 1 1 NaN 0.000000 <NA>
1 1 0.877437 0.877437 5.096665 10 10 NaN 0.889408 <NA>
2 2 0.986541 0.109104 1.688946 34 34 NaN 0.110592 <NA>

We can also inspect the largest fitted nonconstant coefficients. Their units follow the response, and their signs depend on the encoded contrasts:

wh["coefficients"].query("order > 0").sort_values(
    "coefficient", key=abs, ascending=False,
).head(8)

Output

order positions term coefficient
1 1 (1,) 1.0_1_0.25 -28.555556
9 1 (3,) 6_3_24 15.777778
13 2 (1, 2) 1.0_1_0.25-30_2_70 -15.666667
12 2 (1, 2) 1.0_1_0.25-30_2_60 -15.333333
21 2 (1, 2) 1.0_1_0.5-30_2_80 15.000000
23 2 (1, 3) 1.0_1_0.5-6_3_24 12.833333
2 1 (1,) 1.0_1_0.5 -11.333333
11 2 (1, 2) 1.0_1_0.25-30_2_50 -10.666667

6.4 Fitness-distance correlation

Distance sums the sampled rank steps in concentration, temperature and residence time. It does not add their physical units together.

We use Spearman fdc to compare ranks. A negative value means higher responses tend to occur nearer the optimum. A value near zero indicates little monotonic association. This trend does not guarantee an improving path from every configuration.

fitness_distance_r = analysis.fdc(landscape, method="spearman")
print(f"Fitness-distance correlation: {fitness_distance_r:.4f}")

Output

Fitness-distance correlation: -0.5820

The two local optima coexist with positive autocorrelation (0.631) and negative FDC (−0.582). Nearby process settings often have similar responses, and the overall trend favours settings closer to the best observed point.

7. Running several analyses together

We can collect autocorrelation and FDC with analysis.profile(), keeping the same walk settings. landscape.n_lo supplies the count directly; the W–H result keeps its separate coefficient and order-summary tables.

analysis.profile(
    landscape,
    metrics=["autocorrelation", "fdc"],
    params={"autocorrelation": {"walk_length": 20, "walk_times": 200, "lag": 1}},
    seed=42, progress=False,
)

Output

autocorrelation    0.630833
fdc               -0.582001
dtype: float64

8. Analysis reference

For further exploration, analysis.list_metrics() lists the metrics available to profile(). The worked W–H result also exposes coefficients and fit_info; its order summary separates cumulative training fit from the fitted model’s variance spectrum.

Data sources

Ashikari Y. et al. (2022). Investigation of Parameter Control for Electrocatalytic Semihydrogenation in a Proton-Exchange Membrane Reactor Utilizing Bayesian Optimization. We use the complete 54-point grid in supplementary Table S5, preserving the authors’ reported index rather than recomputing it from rounded component values.