Bayesian Optimization of UV Curable Elastomer

Introduction

Let's use Bayesian optimization to efficiently design a UV curable elastomer. Note for this post I'll use a mix of real and synthetic data for illustrative purposes.

If you're not familiar with the concept of Bayesian optimization and how it compares to OFAT and DOE, check out: https://en.wikipedia.org/wiki/Bayesian_optimization, https://distill.pub/2020/bayesian-optimization/

High Level Comparison

There are several open and closed source options for performing Bayesian optimization, for this post I'll be using BayBE (https://emdgroup.github.io/baybe/stable/index.html)

The maintainers have a good intro video here:
https://www.youtube.com/watch?v=EcVvYhh4wI0&t=1258s&pp=ygUMYmF5YmUgcHl0aG9u

Installation

I recommend installing BayBE with chem and insights modules (can leverage molecular descriptors and fingerprints in model during optimization from chem module, Shapley values from insights module, save for another post)

pip install baybe[chem,insights]

I'll also be using RDKit for some early screening, but not needed.

pip install rdkit

Targets

As an example, we will design a UV-curable elastomer targeting the properties outlined below. Note this is intended for demonstrative purposes only and does not represent a complete formulation or an exhaustive list of physical and mechanical property requirements.

Table 1: Key Properties

Raw materials:

• 9 monomers (Tg, SMILES known)
• Oligomer_1 (Polyurethane methacrylate elastomer, high viscosity)
• Oligomer_2 (Polyurethane methacrylate elastomer, medium viscosity)
• Oligomer_3 (Polyurethane methacrylate elastomer, low viscosity)

Prescreening

Operating under the following assumptions:

• High Tg monomers needed to impart strength but will worsen resilience / energy return
• Higher viscosity oligomers potentially higher molecular weight, possibly higher resilience, but lower crosslink density and strength

I'll start with a theory that while bulky steric groups impart strength they will also reduce resilience of the cured elastomer, loss modulus will increase (think rubber-band snapping back), thus need to find good trade-off between Tg and bulkiness in monomers. For this I'll use RDKit to calculate van der waals volume for each monomer based on its SMILES string then divide its literature Tg by calculated volume.

from rdkit import Chem  
from rdkit.Chem import AllChem  
import pandas as pd  
  
def get_vdw_volume(smiles: str) -> float:  
    mol = Chem.MolFromSmiles(smiles)  
    mol = Chem.AddHs(mol)  
    AllChem.EmbedMolecule(mol)  
    AllChem.MMFFOptimizeMolecule(mol)  
    vdw_volume = AllChem.ComputeMolVolume(mol) # A^3  
    return vdw_volume  
  
high_tg_monomers = pd.DataFrame(  
    {  
        "monomers": [  
            "IBOMA",  
            "Monomer D",    
            "Monomer C",  
            "Monomer B",  
            "Monomer A",  
            "HEMA",  
            "tert-Butyl methacrylate",  
            "Neopentyl methacrylate",  
        ],  
        "smiles": [  
            "CC(=C)C(=O)OC1CC2CCC1(C)C2(C)C",  
            "CC1(C)[C@H]2CC[C@]1(C)[C@H](C2)OC(=O)C=C",  
            "CC1(C)[C@H]2CC[C@]1(C)[C@H](C2)OC(=O)C=C",
            "CC1(C)[C@H]2CC[C@]1(C)[C@H](C2)OC(=O)C=C",
            "CC1(C)[C@H]2CC[C@]1(C)[C@H](C2)OC(=O)C=C",
            "CC1(C)[C@H]2CC[C@]1(C)[C@H](C2)OC(=O)C=C",  
            "CC(=C)C(=O)OCCO",  
            "CC(=C)C(=O)OC(C)(C)C",  
            "CC(=C)C(=O)OCC(C)(C)C",  
        ],  
        "Tg": [115, 125, 145, 142, 147, 175, 55, 107, 110],  
    }  
)  
  
vdw = []  
for i in high_tg_monomers["smiles"]:  
    vdw.append(get_vdw_volume(i))  
high_tg_monomers["vdw_volume"] = vdw  
high_tg_monomers["Tg/vdw_volume"] = high_tg_monomers["Tg"] / high_tg_monomers["vdw_volume"]  
print(high_tg_monomers.sort_values("Tg/vdw_volume", ascending=False).to_string())

Which gives us:

We'll filter out monomers with lower Tg/vdw ratios to reduce search space with the exception of IBOMA as we expect it to perform poorly (will our theory hold water...).

Bayesian Optimization

The Setup

In order to perform Bayesian optimization we must set the following:

• Targets / objectives
  • Selecting multi-objective for: UTS, resilience, viscosity
• Objective type
  • Selecting Pareto Objective (This approach allows one to select the optimal trade-offs in target properties after optimization has been completed. An alternative is to combine objectives into a single scoring function using user-defined weights. However, the later approach (DesirabilityObjective) requires those weights to be specified before optimization, meaning that if the relative importance of the objectives or target properties changes, the optimization must be rerun)
• Search space
  • Individual monomers may be included between 0-50% by wt.
  • Individual oligomers may be included between 0-60% by wt.
  • Total monomer content must be between 20-50% by wt.
  • The above bounds and constraints are based on prior knowledge / screening, monomer loadings below threshold are not likely to exhibit adequate strength, monomer loadings above threshold are likely to restrict too much movement
  • Combined monomer and oligomer content must sum to 100% by wt.
  • Recommended formulations must have raw material cost below $26/kg
• Recommender
  • We initially start with randomly sampling the search space to get first set of recommendations then switch to Bayesian optimization via BotorchRecommender

We instantiate the Campaign class with the specified parameters and generate the initial set of recommended formulations. Because optimization is an iterative process, requiring formulations to be physically tested, resulting data fed back into the system, and new recommendations generated, the campaign state and recommendations must persist throughout the experiment cycle. Thus we serialize the campaign object to a JSON file while the recommended formulations are exported to a CSV file. This allows experiment data to be collected and subsequently reintroduced into the optimization process across iterations.

create_campaign.py

import json  
from baybe.targets import NumericalTarget  
from baybe.objectives import ParetoObjective  
from baybe.parameters import NumericalContinuousParameter  
from baybe.searchspace import SearchSpace  
from baybe.constraints import ContinuousLinearConstraint  
from baybe.recommenders import (  
    BotorchRecommender,  
    TwoPhaseMetaRecommender,  
    RandomRecommender,  
)  
from baybe import Campaign  
  
try:  
    with open("campaign.json", "r", encoding="utf-8") as file:  
        campaign = Campaign.from_json(str(json.load(file)))  
    print("Campaign already exists")  
except:  
    target_1 = NumericalTarget(name="UTS")  
	target_2 = NumericalTarget(name="resilience") 
    target_3 = NumericalTarget.match_absolute(  
        name="viscosity", match_value=17500, match_mode="<="  
    )  
    objective = ParetoObjective(  
        targets=[target_1, target_2, target_3]  
    )  
  
    parameters = [  
        NumericalContinuousParameter(name="Monomer_A", bounds=(0, 0.5)),  
        NumericalContinuousParameter(name="Monomer_B", bounds=(0, 0.5)),  
        NumericalContinuousParameter(name="Monomer_C", bounds=(0, 0.5)),  
        NumericalContinuousParameter(name="Monomer_D", bounds=(0, 0.5)),  
        NumericalContinuousParameter(name="IBOMA", bounds=(0, 0.5)),  
        NumericalContinuousParameter(name="Oligomer_1", bounds=(0, 0.6)),  
        NumericalContinuousParameter(name="Oligomer_2", bounds=(0, 0.6)),  
        NumericalContinuousParameter(name="Oligomer_3", bounds=(0, 0.6)),  
    ]  
  
    constraints = [  
        ContinuousLinearConstraint(  
            parameters=[  
                "Monomer_A",  
                "Monomer_B",  
                "Monomer_C",  
                "Monomer_D",  
                "IBOMA",  
                "Oligomer_1",  
                "Oligomer_2",  
                "Oligomer_3",  
            ],  
            operator="=",  
            coefficients=(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0),  
            rhs=1.0,  # All components must sum to 100%  
        ),  
        ContinuousLinearConstraint(  
            parameters=["Monomer_A", "Monomer_B", "Monomer_C", "Monomer_D", "IBOMA"],  
            operator="<=",  
            coefficients=(1.0, 1.0, 1.0, 1.0, 1.0),  
            rhs=0.5,  # Total monomer loading must be less than or equal to 50%  
        ),  
        ContinuousLinearConstraint(  
            parameters=["Monomer_A", "Monomer_B", "Monomer_C", "Monomer_D", "IBOMA"],  
            operator=">=",  
            coefficients=(1.0, 1.0, 1.0, 1.0, 1.0),  
            rhs=0.2,  # Total monomer loading must be greater than or equal to 20%  
        ),  
        ContinuousLinearConstraint(  
            parameters=[  
                "Monomer_A",  
                "Monomer_B",  
                "Monomer_C",  
                "Monomer_D",  
                "IBOMA",  
                "Oligomer_1",  
                "Oligomer_2",  
                "Oligomer_3",  
            ],  
            operator="<=",  
            coefficients=(  
                25,  
                14,  
                32,  
                17,  
                23,  
                24,  
                20,  
                23,  
            ),  # Each coefficient represents raw material cost $/kg  
            rhs=2600,  # Total raw cost must be less than or equal to $26/kg  
        ),  
    ]  
  
    searchspace = SearchSpace.from_product(  
        parameters=parameters, constraints=constraints  
    )  
  
    recommender = TwoPhaseMetaRecommender(  
        initial_recommender=RandomRecommender(),  
        recommender=BotorchRecommender(),  
    )  
  
    campaign = Campaign(searchspace, objective, recommender)  
    recommended_experiments = campaign.recommend(batch_size=3)  
  
    with open("campaign.json", "w") as file:  
        json.dump(campaign.to_json(), file)  
  
    for i in range(0, len(campaign.targets)):  
        recommended_experiments[campaign.targets[i].to_dict()["name"]] = ""  
    recommended_experiments.to_csv("recommended_experiments.csv")  
  
    print("Created campaign, wrote to campaign.json")  
    print("Created initial recommendations, wrote to recommended_experiments.csv")

recommended_experiments.csv

Go in the lab, test formulations, update CSV with results:

Now we'd like to update our campaign with test results:

update_campaign.py

import json  
from baybe import Campaign  
import numpy as np  
import pandas as pd  
  
try:  
    with open("campaign.json", "r", encoding="utf-8") as file:  
        campaign = Campaign.from_json(str(json.load(file)))  
    print("Loaded campaign...")  
except:  
    raise ValueError("Campaign does not exist")  
  
try:  
    recommended_experiments = pd.read_csv("recommended_experiments.csv")  
    recommended_experiments = recommended_experiments.reset_index(drop=True)  
    recommended_experiments = recommended_experiments.replace("", np.nan).dropna()  
  
    filter = campaign.measurements.filter(regex=r"^(?!.*BatchNr)").filter(  
        regex=r"^(?!.*FitNr)"  
    )  
    new_measurements_only = pd.concat(  
        [filter, recommended_experiments]  
    ).drop_duplicates(keep=False)  
  
    if not new_measurements_only.empty:  
        print(  
            "New measurements imported from recommended_experiments.csv and added to campaign"  
        )  
        campaign.add_measurements(new_measurements_only)  
        with open("campaign.json", "w") as file:  
            json.dump(campaign.to_json(), file)  
        print("Updated campaign.json")  
    else:  
        print("No new measurements in recommended_experiments.csv, doing nothing")  
  
except:  
    print("No recommended_experiments.csv to update from")

• This will only import test results for formulations that contain data for every target, missing data will cause entire formulation row to be skipped

Now we perform Bayesian optimization with initial results, recommend next batch of formulations to test, and update both campaign.json and recommended_experiments.csv files:

recommend_formulations.py

import json  
from baybe import Campaign  
  
try:  
    with open("campaign.json", "r", encoding="utf-8") as file:  
        campaign = Campaign.from_json(str(json.load(file)))  
    print("Loaded campaign...")  
except:  
    raise ValueError("Campaign does not exist")  
  
recommended_experiments = campaign.recommend(batch_size=3)  
  
for i in range(0, len(campaign.targets)):  
    recommended_experiments[campaign.targets[i].to_dict()["name"]] = ""  
  
recommended_experiments.to_csv("recommended_experiments.csv", mode="a", header=False)  
print("Added recommendations to recommended_experiments.csv")  
  
with open("campaign.json", "w") as file:  
    json.dump(campaign.to_json(), file)  
print("Updated campaign.json")

Once the next batch of experimental data has been collected only update_campaign.py and recommend_formulations.py are required to continue the optimization process.

After multiple iterations of optimization we have the following formulations and collected properties:

ID11 tested with 36% IBOMA loading resulted in ~5% resilience and causes IBOMA to primarily be avoided in future tests.

Since we selected Pareto Objective we end up with a set of optimal compromises where we can seek out a formulation matching our original targets, eg: ID22, but we also have access to other formulations that may be of interest, eg: if product requirements shift mid project perhaps ID20 with a lower UTS and higher resilience value addresses our updated needs:

Pareto Front of UTS vs Resilience for each optimization cycle:

I'll end with a few points for consideration:

Bayesian optimization becomes increasingly challenging to leverage as the number of objectives grows, eg: 4, 5, 6... Each additional objective generally requires its own surrogate model, increasing computational cost and model complexity. At the same time, the size of the Pareto front typically grows with the number of objectives, creating a larger set of potentially optimal trade-offs. This makes it more difficult to efficiently explore the design space, distinguish meaningful improvements between candidates, and achieve convergence.

BayBE offers a broad range of capabilities, and its documentation is well worth reviewing. Two features that are particularly noteworthy:

• Transfer learning allows existing experimental data and prior learnings from previous campaigns to be used in new campaigns and ultimately accelerate optimization
• Molecular descriptors can also be incorporated as parameter representations, either by supplying SMILES strings and selecting Scikit-fingerprint fingerprint type, or through custom descriptor dictionaries, eg: including Tg, MW, etc. Incorporating these descriptors provides the model with chemically meaningful information and can improve predictions compared with relying solely on categorical material labels--such incorporation does however increase dimensionality and require more compute.