Chapter 3: Practice: Application to Materials Discovery

Learn how to approach optimal solutions while reducing the number of experiments through experimental planning that leverages uncertainty. We'll also review key considerations for field deployment.

💡 Note: Understand experimental rollback costs upfront. Implementing constraints that err on the side of safety makes operations easier to manage.

Learn Real-World Materials Optimization with Python Implementation

Learning Objectives

By reading this chapter, you will be able to:

• ✅ Integrate materials property prediction ML models with Bayesian Optimization
• ✅ Implement constrained optimization and consider materials feasibility
• ✅ Calculate Pareto optimal solutions with multi-objective optimization
• ✅ Implement batch Bayesian Optimization considering experimental costs
• ✅ Solve real-world Li-ion battery optimization problems

Reading Time: 25-30 min Code Examples: 12 Exercises: 3

3.1 Integration with Materials Property Prediction ML Models

Why Integrate with ML Models?

In materials exploration, Bayesian Optimization is combined as follows:

1. Build ML Model from Existing Data - Public databases like Materials Project - Past experimental data - DFT calculation results

2. Explore New Materials with Bayesian Optimization - Use ML model as the objective function - Minimize number of experiments - Exploit uncertainty

Acquiring Data from Materials Project API

Code Example 1: Acquiring Data from Materials Project

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0

# Acquire data from Materials Project
# Note: mp-api installation required: pip install mp-api
from mp_api.client import MPRester
import pandas as pd
import numpy as np

# Using Materials Project API (API key required)
# Registration: https://materialsproject.org/api
API_KEY = "YOUR_API_KEY_HERE"  # Replace with your actual API key

def fetch_battery_materials(api_key, max_materials=100):
    """
    Acquire data for Li-ion battery cathode materials

    Parameters:
    -----------
    api_key : str
        Materials Project API key
    max_materials : int
        Maximum number of materials to retrieve

    Returns:
    --------
    df : DataFrame
        Materials property data
    """
    with MPRester(api_key) as mpr:
        # Search for Li-containing oxides
        docs = mpr.summary.search(
            elements=["Li", "O"],  # Contains Li and O
            num_elements=(3, 5),    # 3-5 element systems
            fields=[
                "material_id",
                "formula_pretty",
                "formation_energy_per_atom",
                "band_gap",
                "density",
                "volume"
            ]
        )

        # Convert to DataFrame
        data = []
        for doc in docs[:max_materials]:
            data.append({
                'material_id': doc.material_id,
                'formula': doc.formula_pretty,
                'formation_energy': doc.formation_energy_per_atom,
                'band_gap': doc.band_gap,
                'density': doc.density,
                'volume': doc.volume
            })

        df = pd.DataFrame(data)
        return df

# Dummy data for demo (if no API key available)
def generate_dummy_battery_data(n_samples=100):
    """
    Generate dummy Li-ion battery material data

    Parameters:
    -----------
    n_samples : int
        Number of samples

    Returns:
    --------
    df : DataFrame
        Materials property data
    """
    np.random.seed(42)

    # Composition parameters (normalized)
    li_content = np.random.uniform(0.1, 0.5, n_samples)
    ni_content = np.random.uniform(0.1, 0.4, n_samples)
    co_content = np.random.uniform(0.1, 0.4, n_samples)
    mn_content = 1.0 - li_content - ni_content - co_content

    # Capacity (mAh/g): Correlates with Li content
    capacity = (
        150 + 200 * li_content +
        50 * ni_content +
        30 * np.random.randn(n_samples)
    )

    # Voltage (V): Correlates with Co content
    voltage = (
        3.0 + 1.5 * co_content +
        0.2 * np.random.randn(n_samples)
    )

    # Stability (formation energy): Negative is stable
    stability = (
        -2.0 - 0.5 * li_content -
        0.3 * ni_content +
        0.1 * np.random.randn(n_samples)
    )

    df = pd.DataFrame({
        'li_content': li_content,
        'ni_content': ni_content,
        'co_content': co_content,
        'mn_content': mn_content,
        'capacity': capacity,
        'voltage': voltage,
        'stability': stability
    })

    return df

# Acquire data (using dummy data)
df_materials = generate_dummy_battery_data(n_samples=150)

print("Materials data statistics:")
print(df_materials.describe())
print(f"\nData shape: {df_materials.shape}")

Output:

Materials data statistics:
       li_content  ni_content  co_content  mn_content    capacity  \
count  150.000000  150.000000  150.000000  150.000000  150.000000
mean     0.299524    0.249336    0.249821    0.201319  208.964738
std      0.116176    0.085721    0.083957    0.122841   38.259483
min      0.102543    0.101189    0.103524   -0.107479  137.582916
max      0.499765    0.399915    0.398774    0.499304  311.495867

         voltage   stability
count  150.000000  150.000000
mean     3.374732   -2.161276
std      0.285945    0.221438
min      2.762894   -2.774301
max      4.137882   -1.554217

Data shape: (150, 7)

Property Prediction with Machine Learning Models

Code Example 2: Building Capacity Prediction Model with Random Forest

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0

"""
Example: Code Example 2: Building Capacity Prediction Model with Rand

Purpose: Demonstrate data visualization techniques
Target: Advanced
Execution time: 30-60 seconds
Dependencies: None
"""

# Capacity prediction with Random Forest
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

# Features and targets
X = df_materials[['li_content', 'ni_content',
                   'co_content', 'mn_content']].values
y_capacity = df_materials['capacity'].values
y_voltage = df_materials['voltage'].values
y_stability = df_materials['stability'].values

# Data split
X_train, X_test, y_train, y_test = train_test_split(
    X, y_capacity, test_size=0.2, random_state=42
)

# Random Forest model
rf_model = RandomForestRegressor(
    n_estimators=100,
    max_depth=10,
    min_samples_split=5,
    random_state=42
)

# Training
rf_model.fit(X_train, y_train)

# Prediction
y_pred_train = rf_model.predict(X_train)
y_pred_test = rf_model.predict(X_test)

# Evaluation
train_rmse = np.sqrt(mean_squared_error(y_train, y_pred_train))
test_rmse = np.sqrt(mean_squared_error(y_test, y_pred_test))
test_r2 = r2_score(y_test, y_pred_test)

# Cross-validation
cv_scores = cross_val_score(
    rf_model, X_train, y_train,
    cv=5, scoring='r2'
)

print("Random Forest model performance:")
print(f"  Training RMSE: {train_rmse:.2f} mAh/g")
print(f"  Test RMSE: {test_rmse:.2f} mAh/g")
print(f"  Test R²: {test_r2:.3f}")
print(f"  CV R² (5-fold): {cv_scores.mean():.3f} ± {cv_scores.std():.3f}")

# Feature importance
feature_names = ['Li', 'Ni', 'Co', 'Mn']
importances = rf_model.feature_importances_
indices = np.argsort(importances)[::-1]

print("\nFeature importance:")
for i in range(len(feature_names)):
    print(f"  {feature_names[indices[i]]}: {importances[indices[i]]:.3f}")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Predicted vs Actual
ax1 = axes[0]
ax1.scatter(y_train, y_pred_train, alpha=0.5, label='Training')
ax1.scatter(y_test, y_pred_test, alpha=0.7, label='Test')
ax1.plot([y_capacity.min(), y_capacity.max()],
         [y_capacity.min(), y_capacity.max()],
         'k--', linewidth=2, label='Ideal')
ax1.set_xlabel('Actual Capacity (mAh/g)', fontsize=12)
ax1.set_ylabel('Predicted Capacity (mAh/g)', fontsize=12)
ax1.set_title('Random Forest Capacity Prediction', fontsize=14)
ax1.legend()
ax1.grid(True, alpha=0.3)

# Feature importance
ax2 = axes[1]
ax2.barh(range(len(feature_names)), importances[indices],
         color='steelblue')
ax2.set_yticks(range(len(feature_names)))
ax2.set_yticklabels([feature_names[i] for i in indices])
ax2.set_xlabel('Importance', fontsize=12)
ax2.set_title('Feature Importance', fontsize=14)
ax2.grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.savefig('ml_model_performance.png', dpi=150, bbox_inches='tight')
plt.show()

Exploiting ML Model with Bayesian Optimization

Code Example 3: Integration of ML Model and Bayesian Optimization

# ML model-based optimization using scikit-optimize
from skopt import gp_minimize
from skopt.space import Real
from skopt.plots import plot_convergence

def objective_function_ml(x):
    """
    Use ML model as objective function

    Parameters:
    -----------
    x : list
        [li_content, ni_content, co_content, mn_content]

    Returns:
    --------
    float : Negative capacity (converted to minimization problem)
    """
    # Composition constraint: total=1.0
    li, ni, co, mn = x
    total = li + ni + co + mn

    # Penalty for constraint violation
    if not (0.98 <= total <= 1.02):
        return 1000.0  # Large penalty

    # Individual constraints
    if li < 0.1 or li > 0.5:
        return 1000.0
    if ni < 0.1 or ni > 0.4:
        return 1000.0
    if co < 0.1 or co > 0.4:
        return 1000.0
    if mn < 0.0:
        return 1000.0

    # Capacity prediction with ML model
    X_pred = np.array([[li, ni, co, mn]])
    capacity_pred = rf_model.predict(X_pred)[0]

    # Convert to minimization problem (negative capacity)
    return -capacity_pred

# Define search space
space = [
    Real(0.1, 0.5, name='li_content'),
    Real(0.1, 0.4, name='ni_content'),
    Real(0.1, 0.4, name='co_content'),
    Real(0.0, 0.5, name='mn_content')
]

# Execute Bayesian Optimization
result = gp_minimize(
    objective_function_ml,
    space,
    n_calls=50,        # 50 evaluations
    n_initial_points=10,  # Initial random sampling
    random_state=42,
    verbose=False
)

# Results
best_composition = result.x
best_capacity = -result.fun  # Revert negative

print("Bayesian Optimization results:")
print(f"  Optimal composition:")
print(f"    Li: {best_composition[0]:.3f}")
print(f"    Ni: {best_composition[1]:.3f}")
print(f"    Co: {best_composition[2]:.3f}")
print(f"    Mn: {best_composition[3]:.3f}")
print(f"    Total: {sum(best_composition):.3f}")
print(f"  Predicted capacity: {best_capacity:.2f} mAh/g")

# Convergence plot
plt.figure(figsize=(10, 6))
plot_convergence(result)
plt.title('Bayesian Optimization Convergence', fontsize=14)
plt.xlabel('Number of Evaluations', fontsize=12)
plt.ylabel('Best Value So Far (Negative Capacity)', fontsize=12)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bo_ml_convergence.png', dpi=150, bbox_inches='tight')
plt.show()

# Compare with best value in dataset
max_capacity_data = df_materials['capacity'].max()
print(f"\nMaximum capacity in dataset: {max_capacity_data:.2f} mAh/g")
print(f"Improvement rate: {((best_capacity - max_capacity_data) / max_capacity_data * 100):.1f}%")

Expected Output:

Bayesian Optimization results:
  Optimal composition:
    Li: 0.487
    Ni: 0.312
    Co: 0.152
    Mn: 0.049
    Total: 1.000
  Predicted capacity: 267.34 mAh/g

Maximum capacity in dataset: 311.50 mAh/g
Improvement rate: -14.2%

3.2 Constrained Optimization

Materials Feasibility Constraints

In actual materials development, there are the following constraints:

1. Composition Constraints: Total 100%, upper and lower limits for each element
2. Stability Constraints: formation energy < threshold
3. Experimental Constraints: Synthesis temperature, pressure range
4. Cost Constraints: Limit use of expensive elements

Implementation of Constrained Bayesian Optimization

Code Example 4: Optimization Under Multiple Constraints

# Requirements:
# - Python 3.9+
# - torch>=2.0.0, <2.3.0

# Constrained Bayesian Optimization (using BoTorch)
# Note: BoTorch installation required: pip install botorch torch
import torch
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_model
from gpytorch.mlls import ExactMarginalLogLikelihood
from botorch.acquisition import ExpectedImprovement
from botorch.optim import optimize_acqf

def constrained_bo_example():
    """
    Demo of constrained Bayesian Optimization

    Constraints:
    - Maximize capacity
    - Stability: formation energy < -1.5 eV/atom
    - Cost: Co content < 0.3
    """
    # Initial data (random sampling)
    n_initial = 10
    np.random.seed(42)

    X_init = np.random.rand(n_initial, 4)
    # Normalize composition
    X_init = X_init / X_init.sum(axis=1, keepdims=True)

    # Evaluate objective function and constraints
    y_capacity = []
    y_stability = []
    for i in range(n_initial):
        x = X_init[i]
        # Capacity prediction
        capacity = rf_model.predict(x.reshape(1, -1))[0]
        # Stability (simplified model)
        stability = -2.0 - 0.5*x[0] - 0.3*x[1] + 0.1*np.random.randn()

        y_capacity.append(capacity)
        y_stability.append(stability)

    X_init = torch.tensor(X_init, dtype=torch.float64)
    y_capacity = torch.tensor(y_capacity, dtype=torch.float64).unsqueeze(-1)
    y_stability = torch.tensor(y_stability, dtype=torch.float64).unsqueeze(-1)

    # Sequential optimization (20 iterations)
    n_iterations = 20
    X_all = X_init.clone()
    y_capacity_all = y_capacity.clone()
    y_stability_all = y_stability.clone()

    for iteration in range(n_iterations):
        # Gaussian Process model (capacity)
        gp_capacity = SingleTaskGP(X_all, y_capacity_all)
        mll_capacity = ExactMarginalLogLikelihood(
            gp_capacity.likelihood, gp_capacity
        )
        fit_gpytorch_model(mll_capacity)

        # Gaussian Process model (stability)
        gp_stability = SingleTaskGP(X_all, y_stability_all)
        mll_stability = ExactMarginalLogLikelihood(
            gp_stability.likelihood, gp_stability
        )
        fit_gpytorch_model(mll_stability)

        # Expected Improvement (capacity)
        best_f = y_capacity_all.max()
        EI = ExpectedImprovement(gp_capacity, best_f=best_f)

        # Optimize Acquisition Function (considering constraints)
        bounds = torch.tensor([[0.1, 0.1, 0.1, 0.0],
                                [0.5, 0.4, 0.3, 0.5]],
                               dtype=torch.float64)

        candidate, acq_value = optimize_acqf(
            EI,
            bounds=bounds,
            q=1,
            num_restarts=10,
            raw_samples=512,
        )

        # Evaluate candidate point
        x_new = candidate.detach().numpy()[0]
        # Normalize
        x_new = x_new / x_new.sum()

        # Experiment simulation
        capacity_new = rf_model.predict(x_new.reshape(1, -1))[0]
        stability_new = -2.0 - 0.5*x_new[0] - 0.3*x_new[1] + \
                        0.1*np.random.randn()

        # Check constraints
        feasible = (stability_new < -1.5) and (x_new[2] < 0.3)

        if feasible:
            print(f"Iteration {iteration+1}: "
                  f"Capacity={capacity_new:.1f}, "
                  f"Stability={stability_new:.2f}, "
                  f"Feasible=Yes")
        else:
            print(f"Iteration {iteration+1}: "
                  f"Capacity={capacity_new:.1f}, "
                  f"Stability={stability_new:.2f}, "
                  f"Feasible=No (constraint violation)")

        # Add to data
        X_all = torch.cat([X_all, torch.tensor(x_new).unsqueeze(0)], dim=0)
        y_capacity_all = torch.cat([y_capacity_all,
                                     torch.tensor([[capacity_new]])], dim=0)
        y_stability_all = torch.cat([y_stability_all,
                                      torch.tensor([[stability_new]])], dim=0)

    # Extract best solution among feasible solutions
    feasible_mask = (y_stability_all < -1.5).squeeze() & \
                    (X_all[:, 2] < 0.3).squeeze()

    if feasible_mask.sum() > 0:
        feasible_capacities = y_capacity_all[feasible_mask]
        feasible_X = X_all[feasible_mask]
        best_idx = feasible_capacities.argmax()
        best_composition_constrained = feasible_X[best_idx].numpy()
        best_capacity_constrained = feasible_capacities[best_idx].item()

        print("\nFinal result (constrained):")
        print(f"  Optimal composition:")
        print(f"    Li: {best_composition_constrained[0]:.3f}")
        print(f"    Ni: {best_composition_constrained[1]:.3f}")
        print(f"    Co: {best_composition_constrained[2]:.3f} "
              f"(constraint < 0.3)")
        print(f"    Mn: {best_composition_constrained[3]:.3f}")
        print(f"  Predicted capacity: {best_capacity_constrained:.2f} mAh/g")
        print(f"  Number of feasible solutions: {feasible_mask.sum().item()} / "
              f"{len(X_all)}")
    else:
        print("\nNo feasible solution found")

# Execute
constrained_bo_example()

3.3 Multi-Objective Optimization (Pareto Optimization)

Why Multi-Objective Optimization is Needed

In materials development, it is necessary to optimize multiple properties simultaneously:

• Li-ion battery: Capacity ↑, Voltage ↑, Stability ↑
• Thermoelectric materials: Seebeck coefficient ↑, Electrical conductivity ↑, Thermal conductivity ↓
• Catalysts: Activity ↑, Selectivity ↑, Stability ↑, Cost ↓

These have trade-offs, and no single optimal solution exists.

Concept of Pareto Frontier

Definition of Pareto Optimality:

A solution x is Pareto optimal ⇔ There exists no solution that simultaneously improves all objectives

Expected Hypervolume Improvement (EHVI)

Code Example 5: Implementation of Multi-Objective Bayesian Optimization

# Multi-objective Bayesian Optimization
from botorch.models import ModelListGP
from botorch.acquisition.multi_objective import \
    qExpectedHypervolumeImprovement
from botorch.utils.multi_objective.box_decompositions.dominated import \
    DominatedPartitioning

def multi_objective_bo_example():
    """
    Demo of multi-objective Bayesian Optimization

    Objectives:
    1. Maximize capacity
    2. Maximize stability (minimize absolute value of formation energy)
    """
    # Initial data
    n_initial = 15
    np.random.seed(42)

    X_init = np.random.rand(n_initial, 4)
    X_init = X_init / X_init.sum(axis=1, keepdims=True)

    # Evaluate two objective functions
    y1_capacity = []
    y2_stability = []

    for i in range(n_initial):
        x = X_init[i]
        capacity = rf_model.predict(x.reshape(1, -1))[0]
        stability = -2.0 - 0.5*x[0] - 0.3*x[1] + 0.1*np.random.randn()
        # Convert stability to positive (unified as maximization problem)
        stability_positive = -stability

        y1_capacity.append(capacity)
        y2_stability.append(stability_positive)

    X_all = torch.tensor(X_init, dtype=torch.float64)
    Y_all = torch.tensor(
        np.column_stack([y1_capacity, y2_stability]),
        dtype=torch.float64
    )

    # Sequential optimization
    n_iterations = 20

    for iteration in range(n_iterations):
        # Gaussian Process models (one for each objective function)
        gp_list = []
        for i in range(2):
            gp = SingleTaskGP(X_all, Y_all[:, i].unsqueeze(-1))
            mll = ExactMarginalLogLikelihood(gp.likelihood, gp)
            fit_gpytorch_model(mll)
            gp_list.append(gp)

        model = ModelListGP(*gp_list)

        # Reference point (worse than Nadir point)
        ref_point = Y_all.min(dim=0).values - 10.0

        # Calculate Pareto frontier
        pareto_mask = is_non_dominated(Y_all)
        pareto_Y = Y_all[pareto_mask]

        # EHVI Acquisition Function
        partitioning = DominatedPartitioning(
            ref_point=ref_point,
            Y=pareto_Y
        )
        acq_func = qExpectedHypervolumeImprovement(
            model=model,
            ref_point=ref_point,
            partitioning=partitioning
        )

        # Optimization
        bounds = torch.tensor([[0.1, 0.1, 0.1, 0.0],
                                [0.5, 0.4, 0.4, 0.5]],
                               dtype=torch.float64)

        candidate, acq_value = optimize_acqf(
            acq_func,
            bounds=bounds,
            q=1,
            num_restarts=10,
            raw_samples=512,
        )

        # Evaluate new candidate point
        x_new = candidate.detach().numpy()[0]
        x_new = x_new / x_new.sum()

        capacity_new = rf_model.predict(x_new.reshape(1, -1))[0]
        stability_new = -2.0 - 0.5*x_new[0] - 0.3*x_new[1] + \
                        0.1*np.random.randn()
        stability_positive_new = -stability_new

        y_new = torch.tensor([[capacity_new, stability_positive_new]],
                              dtype=torch.float64)

        # Add to data
        X_all = torch.cat([X_all, torch.tensor(x_new).unsqueeze(0)], dim=0)
        Y_all = torch.cat([Y_all, y_new], dim=0)

        if (iteration + 1) % 5 == 0:
            print(f"Iteration {iteration+1}: "
                  f"Pareto solutions={pareto_mask.sum().item()}, "
                  f"HV={compute_hypervolume(pareto_Y, ref_point):.2f}")

    # Final Pareto frontier
    pareto_mask_final = is_non_dominated(Y_all)
    pareto_X_final = X_all[pareto_mask_final].numpy()
    pareto_Y_final = Y_all[pareto_mask_final].numpy()

    print(f"\nFinal Pareto optimal solutions: {pareto_mask_final.sum().item()}")

    # Visualize Pareto frontier
    plt.figure(figsize=(10, 6))

    # All points
    plt.scatter(Y_all[:, 0].numpy(), Y_all[:, 1].numpy(),
                c='lightblue', s=50, alpha=0.5, label='All explored points')

    # Pareto optimal solutions
    plt.scatter(pareto_Y_final[:, 0], pareto_Y_final[:, 1],
                c='red', s=100, edgecolors='black', zorder=10,
                label='Pareto optimal solutions')

    # Connect Pareto frontier with line
    sorted_indices = np.argsort(pareto_Y_final[:, 0])
    plt.plot(pareto_Y_final[sorted_indices, 0],
             pareto_Y_final[sorted_indices, 1],
             'r--', linewidth=2, alpha=0.5, label='Pareto Frontier')

    plt.xlabel('Objective 1: Capacity (mAh/g)', fontsize=12)
    plt.ylabel('Objective 2: Stability (-formation energy)', fontsize=12)
    plt.title('Multi-Objective Optimization: Pareto Frontier', fontsize=14)
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig('pareto_frontier.png', dpi=150, bbox_inches='tight')
    plt.show()

    # Display trade-off examples
    print("\nTrade-off examples:")
    # Capacity-oriented
    idx_max_capacity = np.argmax(pareto_Y_final[:, 0])
    print(f"  Capacity-oriented: Capacity={pareto_Y_final[idx_max_capacity, 0]:.1f}, "
          f"Stability={pareto_Y_final[idx_max_capacity, 1]:.2f}")

    # Stability-oriented
    idx_max_stability = np.argmax(pareto_Y_final[:, 1])
    print(f"  Stability-oriented: Capacity={pareto_Y_final[idx_max_stability, 0]:.1f}, "
          f"Stability={pareto_Y_final[idx_max_stability, 1]:.2f}")

    # Balanced (midpoint)
    normalized_Y = (pareto_Y_final - pareto_Y_final.min(axis=0)) / \
                   (pareto_Y_final.max(axis=0) - pareto_Y_final.min(axis=0))
    distances = np.sqrt(((normalized_Y - 0.5)**2).sum(axis=1))
    idx_balanced = np.argmin(distances)
    print(f"  Balanced: Capacity={pareto_Y_final[idx_balanced, 0]:.1f}, "
          f"Stability={pareto_Y_final[idx_balanced, 1]:.2f}")

# Pareto optimality determination function
def is_non_dominated(Y):
    """
    Determine Pareto optimal solutions

    Parameters:
    -----------
    Y : Tensor (n_points, n_objectives)
        Objective function values

    Returns:
    --------
    mask : Tensor (n_points,)
        True indicates Pareto optimal
    """
    n_points = Y.shape[0]
    is_efficient = torch.ones(n_points, dtype=torch.bool)

    for i in range(n_points):
        if is_efficient[i]:
            # Check if there exists a point superior in all objectives to point i
            is_dominated = (Y >= Y[i]).all(dim=1) & (Y > Y[i]).any(dim=1)
            is_efficient[is_dominated] = False

    return is_efficient

# Hypervolume calculation
def compute_hypervolume(pareto_Y, ref_point):
    """
    Calculate Hypervolume (simplified version)

    Parameters:
    -----------
    pareto_Y : Tensor
        Pareto optimal solutions
    ref_point : Tensor
        Reference point

    Returns:
    --------
    float : Hypervolume
    """
    # Simplified 2D calculation
    sorted_Y = pareto_Y[torch.argsort(pareto_Y[:, 0], descending=True)]
    hv = 0.0
    prev_y1 = ref_point[0]

    for i in range(len(sorted_Y)):
        width = prev_y1 - sorted_Y[i, 0]
        height = sorted_Y[i, 1] - ref_point[1]
        hv += width * height
        prev_y1 = sorted_Y[i, 0]

    return hv.item()

# Execute
# multi_objective_bo_example()
# Note: Commented out because BoTorch is required
print("Multi-objective optimization example requires BoTorch")
print("Please install with: pip install botorch torch and then execute")

3.4 Optimization Considering Experimental Costs

Batch Bayesian Optimization

When multiple experimental devices are available, parallel experiments are possible:

• Conventional: Sequential (1 run → result → next 1 run)
• Batch BO: Propose multiple candidates at once (q-EI)

Workflow

Code Example 6: Batch Bayesian Optimization

# Batch Bayesian Optimization (scikit-optimize)
from scipy.stats import norm

def batch_expected_improvement(X, gp, f_best, xi=0.01):
    """
    Batch Expected Improvement (simplified version)

    Parameters:
    -----------
    X : array (n_candidates, n_features)
        Candidate points
    gp : GaussianProcessRegressor
        Trained GP model
    f_best : float
        Current best value

    Returns:
    --------
    ei : array (n_candidates,)
        EI values
    """
    mu, sigma = gp.predict(X, return_std=True)
    improvement = mu - f_best - xi
    Z = improvement / (sigma + 1e-9)
    ei = improvement * norm.cdf(Z) + sigma * norm.pdf(Z)
    ei[sigma == 0.0] = 0.0
    return ei

def simulate_batch_bo(n_iterations=10, batch_size=3):
    """
    Batch Bayesian Optimization simulation

    Parameters:
    -----------
    n_iterations : int
        Number of iterations
    batch_size : int
        Number of candidates to propose per iteration

    Returns:
    --------
    X_all : array
        All sampling points
    y_all : array
        All observed values
    """
    from sklearn.gaussian_process import GaussianProcessRegressor
    from sklearn.gaussian_process.kernels import RBF, ConstantKernel

    # Initial data
    np.random.seed(42)
    n_initial = 5
    X_sampled = np.random.rand(n_initial, 4)
    X_sampled = X_sampled / X_sampled.sum(axis=1, keepdims=True)

    y_sampled = []
    for i in range(n_initial):
        capacity = rf_model.predict(X_sampled[i].reshape(1, -1))[0]
        y_sampled.append(capacity)

    y_sampled = np.array(y_sampled)

    # Sequential batch optimization
    for iteration in range(n_iterations):
        # Gaussian Process model
        kernel = ConstantKernel(1.0) * RBF(length_scale=0.2)
        gp = GaussianProcessRegressor(
            kernel=kernel,
            n_restarts_optimizer=10,
            random_state=42
        )
        gp.fit(X_sampled, y_sampled)

        # Current best value
        f_best = y_sampled.max()

        # Generate candidate points (many)
        n_candidates = 1000
        X_candidates = np.random.rand(n_candidates, 4)
        X_candidates = X_candidates / X_candidates.sum(axis=1, keepdims=True)

        # Calculate EI
        ei_values = batch_expected_improvement(X_candidates, gp, f_best)

        # Top-k selection (simple method)
        # More advanced methods: q-EI, KB (Kriging Believer)
        top_k_indices = np.argsort(ei_values)[-batch_size:]
        X_batch = X_candidates[top_k_indices]

        # Batch experiment simulation
        y_batch = []
        for x in X_batch:
            capacity = rf_model.predict(x.reshape(1, -1))[0]
            y_batch.append(capacity)

        y_batch = np.array(y_batch)

        # Add to data
        X_sampled = np.vstack([X_sampled, X_batch])
        y_sampled = np.append(y_sampled, y_batch)

        # Progress display
        if (iteration + 1) % 3 == 0:
            best_so_far = y_sampled.max()
            print(f"Iteration {iteration+1}: "
                  f"Batch size={batch_size}, "
                  f"Best so far={best_so_far:.2f} mAh/g")

    return X_sampled, y_sampled

# Execute batch BO
print("Batch Bayesian Optimization (batch_size=3):")
X_batch_bo, y_batch_bo = simulate_batch_bo(n_iterations=10, batch_size=3)

print(f"\nFinal result:")
print(f"  Total experiments: {len(y_batch_bo)}")
print(f"  Best capacity: {y_batch_bo.max():.2f} mAh/g")
print(f"  Optimal composition: {X_batch_bo[y_batch_bo.argmax()]}")

# Compare with sequential BO
print("\nSequential BO (batch_size=1):")
X_seq_bo, y_seq_bo = simulate_batch_bo(n_iterations=30, batch_size=1)
print(f"  Total experiments: {len(y_seq_bo)}")
print(f"  Best capacity: {y_seq_bo.max():.2f} mAh/g")

# Efficiency comparison
plt.figure(figsize=(10, 6))
plt.plot(np.maximum.accumulate(y_seq_bo), 'o-',
         label='Sequential BO (batch_size=1)', linewidth=2, markersize=6)
plt.plot(np.arange(0, len(y_batch_bo), 3),
         np.maximum.accumulate(y_batch_bo)[::3], '^-',
         label='Batch BO (batch_size=3)', linewidth=2, markersize=8)
plt.xlabel('Number of Experiments', fontsize=12)
plt.ylabel('Best Value So Far (mAh/g)', fontsize=12)
plt.title('Batch BO vs Sequential BO Efficiency Comparison', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('batch_bo_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

3.5 Complete Implementation Example: Li-ion Battery Electrolyte Optimization

Problem Setting

Objective: Optimization of Li-ion battery cathode materials

Properties to Optimize: 1. Maximize capacity (mAh/g) 2. Maximize voltage (V) 3. Maximize stability (formation energy)

Constraints: - Total composition = 1.0 - Li content: 0.1-0.5 - Ni content: 0.1-0.4 - Co content: 0.1-0.3 (limited due to high cost) - Mn content: ≥ 0.0

Code Example 7: Complete Implementation of Real-World Problem

# Multi-objective constrained optimization of Li-ion battery cathode materials
class LiIonCathodeOptimizer:
    """
    Optimization class for Li-ion battery cathode materials

    Objectives:
    - Maximize capacity
    - Maximize voltage
    - Maximize stability (considering cost)

    Constraints:
    - Composition constraints
    - Co content limit (cost)
    """

    def __init__(self, capacity_model, voltage_model, stability_model):
        """
        Parameters:
        -----------
        capacity_model : sklearn model
            Capacity prediction model
        voltage_model : sklearn model
            Voltage prediction model
        stability_model : sklearn model
            Stability prediction model
        """
        self.capacity_model = capacity_model
        self.voltage_model = voltage_model
        self.stability_model = stability_model

        # Constraints
        self.co_max = 0.3  # Co content upper limit
        self.composition_bounds = {
            'li': (0.1, 0.5),
            'ni': (0.1, 0.4),
            'co': (0.1, 0.3),
            'mn': (0.0, 0.5)
        }

    def evaluate(self, composition):
        """
        Evaluate material composition

        Parameters:
        -----------
        composition : array [li, ni, co, mn]

        Returns:
        --------
        dict : Predicted values for each property
        """
        # Check constraints
        if not self._check_constraints(composition):
            return {
                'capacity': -1000,
                'voltage': -1000,
                'stability': -1000,
                'feasible': False
            }

        x = composition.reshape(1, -1)

        capacity = self.capacity_model.predict(x)[0]
        # Voltage model (dummy)
        voltage = 3.0 + 1.5 * composition[2] + 0.2 * np.random.randn()
        # Stability model (dummy)
        stability = -2.0 - 0.5*composition[0] - 0.3*composition[1] + \
                    0.1*np.random.randn()

        return {
            'capacity': capacity,
            'voltage': voltage,
            'stability': -stability,  # Convert to positive
            'feasible': True
        }

    def _check_constraints(self, composition):
        """Check constraints"""
        li, ni, co, mn = composition

        # Composition total
        if not (0.98 <= li + ni + co + mn <= 1.02):
            return False

        # Range for each element
        if not (self.composition_bounds['li'][0] <= li <=
                self.composition_bounds['li'][1]):
            return False
        if not (self.composition_bounds['ni'][0] <= ni <=
                self.composition_bounds['ni'][1]):
            return False
        if not (self.composition_bounds['co'][0] <= co <=
                self.composition_bounds['co'][1]):
            return False
        if not (self.composition_bounds['mn'][0] <= mn <=
                self.composition_bounds['mn'][1]):
            return False

        return True

    def optimize_multi_objective(self, n_iterations=50):
        """
        Execute multi-objective optimization

        Returns:
        --------
        pareto_solutions : list of dict
            Pareto optimal solutions
        """
        # Initial sampling
        n_initial = 20
        np.random.seed(42)

        solutions = []

        for i in range(n_initial):
            # Generate random composition
            composition = np.random.rand(4)
            composition = composition / composition.sum()

            # Evaluate
            result = self.evaluate(composition)

            if result['feasible']:
                solutions.append({
                    'composition': composition,
                    'capacity': result['capacity'],
                    'voltage': result['voltage'],
                    'stability': result['stability']
                })

        # Sequential optimization (simplified version)
        for iteration in range(n_iterations - n_initial):
            # Extract Pareto optimal from existing solutions
            pareto_sols = self._extract_pareto(solutions)

            # Sample around Pareto solutions (simple method)
            if len(pareto_sols) > 0:
                base_sol = pareto_sols[np.random.randint(len(pareto_sols))]
                composition_new = base_sol['composition'] + \
                                  np.random.randn(4) * 0.05
                composition_new = np.clip(composition_new, 0.01, 0.8)
                composition_new = composition_new / composition_new.sum()
            else:
                composition_new = np.random.rand(4)
                composition_new = composition_new / composition_new.sum()

            # Evaluate
            result = self.evaluate(composition_new)

            if result['feasible']:
                solutions.append({
                    'composition': composition_new,
                    'capacity': result['capacity'],
                    'voltage': result['voltage'],
                    'stability': result['stability']
                })

        # Final Pareto optimal solutions
        pareto_solutions = self._extract_pareto(solutions)

        return pareto_solutions, solutions

    def _extract_pareto(self, solutions):
        """Extract Pareto optimal solutions"""
        if len(solutions) == 0:
            return []

        objectives = np.array([
            [s['capacity'], s['voltage'], s['stability']]
            for s in solutions
        ])

        pareto_mask = np.ones(len(objectives), dtype=bool)

        for i in range(len(objectives)):
            if pareto_mask[i]:
                # Check if there exists a solution superior in all objectives to solution i
                dominated = (
                    (objectives >= objectives[i]).all(axis=1) &
                    (objectives > objectives[i]).any(axis=1)
                )
                pareto_mask[dominated] = False

        pareto_solutions = [solutions[i] for i in range(len(solutions))
                             if pareto_mask[i]]

        return pareto_solutions

# Simple training of voltage and stability models (dummy)
from sklearn.ensemble import RandomForestRegressor

voltage_model = RandomForestRegressor(n_estimators=50, random_state=42)
voltage_model.fit(X_train, y_voltage[:len(X_train)])

stability_model = RandomForestRegressor(n_estimators=50, random_state=42)
stability_model.fit(X_train, y_stability[:len(X_train)])

# Execute optimization
optimizer = LiIonCathodeOptimizer(
    capacity_model=rf_model,
    voltage_model=voltage_model,
    stability_model=stability_model
)

print("Executing multi-objective optimization of Li-ion battery cathode materials...")
pareto_solutions, all_solutions = optimizer.optimize_multi_objective(
    n_iterations=100
)

print(f"\nNumber of Pareto optimal solutions: {len(pareto_solutions)}")

# Visualize results (3D)
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(14, 6))

# Left plot: 3D scatter
ax1 = fig.add_subplot(121, projection='3d')

# All solutions
all_cap = [s['capacity'] for s in all_solutions]
all_vol = [s['voltage'] for s in all_solutions]
all_sta = [s['stability'] for s in all_solutions]

ax1.scatter(all_cap, all_vol, all_sta, c='lightblue', s=20,
            alpha=0.3, label='All explored points')

# Pareto optimal solutions
pareto_cap = [s['capacity'] for s in pareto_solutions]
pareto_vol = [s['voltage'] for s in pareto_solutions]
pareto_sta = [s['stability'] for s in pareto_solutions]

ax1.scatter(pareto_cap, pareto_vol, pareto_sta, c='red', s=100,
            edgecolors='black', zorder=10, label='Pareto optimal solutions')

ax1.set_xlabel('Capacity (mAh/g)', fontsize=10)
ax1.set_ylabel('Voltage (V)', fontsize=10)
ax1.set_zlabel('Stability', fontsize=10)
ax1.set_title('3-Objective Optimization: Objective Space', fontsize=12)
ax1.legend()

# Right plot: 2D projection of capacity-voltage
ax2 = fig.add_subplot(122)
ax2.scatter(all_cap, all_vol, c='lightblue', s=20,
            alpha=0.5, label='All explored points')
ax2.scatter(pareto_cap, pareto_vol, c='red', s=100,
            edgecolors='black', zorder=10, label='Pareto optimal solutions')
ax2.set_xlabel('Capacity (mAh/g)', fontsize=12)
ax2.set_ylabel('Voltage (V)', fontsize=12)
ax2.set_title('Capacity-Voltage Trade-off', fontsize=14)
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('liion_cathode_optimization.png', dpi=150,
            bbox_inches='tight')
plt.show()

# Display representative Pareto solutions
print("\nRepresentative Pareto optimal solutions:")

# Capacity-oriented
idx_max_cap = np.argmax(pareto_cap)
print(f"\nCapacity-oriented:")
print(f"  Li={pareto_solutions[idx_max_cap]['composition'][0]:.3f}, "
      f"Ni={pareto_solutions[idx_max_cap]['composition'][1]:.3f}, "
      f"Co={pareto_solutions[idx_max_cap]['composition'][2]:.3f}, "
      f"Mn={pareto_solutions[idx_max_cap]['composition'][3]:.3f}")
print(f"  Capacity={pareto_cap[idx_max_cap]:.1f} mAh/g, "
      f"Voltage={pareto_vol[idx_max_cap]:.2f} V, "
      f"Stability={pareto_sta[idx_max_cap]:.2f}")

# Voltage-oriented
idx_max_vol = np.argmax(pareto_vol)
print(f"\nVoltage-oriented:")
print(f"  Li={pareto_solutions[idx_max_vol]['composition'][0]:.3f}, "
      f"Ni={pareto_solutions[idx_max_vol]['composition'][1]:.3f}, "
      f"Co={pareto_solutions[idx_max_vol]['composition'][2]:.3f}, "
      f"Mn={pareto_solutions[idx_max_vol]['composition'][3]:.3f}")
print(f"  Capacity={pareto_cap[idx_max_vol]:.1f} mAh/g, "
      f"Voltage={pareto_vol[idx_max_vol]:.2f} V, "
      f"Stability={pareto_sta[idx_max_vol]:.2f}")

# Balanced
# Normalize and find solution closest to center
pareto_array = np.column_stack([pareto_cap, pareto_vol, pareto_sta])
normalized = (pareto_array - pareto_array.min(axis=0)) / \
             (pareto_array.max(axis=0) - pareto_array.min(axis=0))
distances = np.sqrt(((normalized - 0.5)**2).sum(axis=1))
idx_balanced = np.argmin(distances)

print(f"\nBalanced:")
print(f"  Li={pareto_solutions[idx_balanced]['composition'][0]:.3f}, "
      f"Ni={pareto_solutions[idx_balanced]['composition'][1]:.3f}, "
      f"Co={pareto_solutions[idx_balanced]['composition'][2]:.3f}, "
      f"Mn={pareto_solutions[idx_balanced]['composition'][3]:.3f}")
print(f"  Capacity={pareto_cap[idx_balanced]:.1f} mAh/g, "
      f"Voltage={pareto_vol[idx_balanced]:.2f} V, "
      f"Stability={pareto_sta[idx_balanced]:.2f}")

3.6 Column: Hyperparameters vs Material Parameters

Two Types of Parameters

In materials exploration, we need to distinguish between two types of parameters:

Material Parameters (Design Variables): - Variables we want to optimize - Examples: Composition ratios, synthesis temperature, pressure - Explored with Bayesian Optimization

Hyperparameters (Algorithm Settings): - Settings for the Bayesian Optimization itself - Examples: Kernel length scale, exploration parameter κ - Optimized with cross-validation or nested BO

Importance of Hyperparameters

Improper hyperparameters can significantly impair optimization efficiency:

• Length scale too large → Cannot capture fine structures
• Length scale too small → Overfitting, local exploration
• κ (UCB) too large → Over-exploration, slow convergence
• κ too small → Over-exploitation, trapped in local optima

Recommended Approaches: 1. Data-driven: Optimize hyperparameters using existing data 2. Robust settings: Choose settings that perform well over a wide range 3. Adaptive adjustment: Decrease κ as optimization progresses (exploration → exploitation)

Code Example 8: Visualize Hyperparameter Effects

# Compare effects of hyperparameters
from sklearn.gaussian_process.kernels import RBF, ConstantKernel

def compare_hyperparameters():
    """
    Compare optimization efficiency with different hyperparameters
    """
    # Test function
    def test_function(x):
        return (np.sin(5*x) * np.exp(-x) +
                0.5 * np.exp(-((x-0.6)/0.15)**2))

    # Different length scales
    length_scales = [0.05, 0.1, 0.3]

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    for idx, ls in enumerate(length_scales):
        ax = axes[idx]

        # Initial data
        np.random.seed(42)
        X_init = np.array([0.1, 0.4, 0.7]).reshape(-1, 1)
        y_init = test_function(X_init.ravel())

        # Gaussian Process
        kernel = ConstantKernel(1.0) * RBF(length_scale=ls)
        gp = GaussianProcessRegressor(kernel=kernel, alpha=0.01,
                                       random_state=42)
        gp.fit(X_init, y_init)

        # Prediction
        X_plot = np.linspace(0, 1, 200).reshape(-1, 1)
        y_pred, y_std = gp.predict(X_plot, return_std=True)

        # Plot
        ax.plot(X_plot, test_function(X_plot.ravel()), 'k--',
                linewidth=2, label='True function')
        ax.scatter(X_init, y_init, c='red', s=100, zorder=10,
                   edgecolors='black', label='Observed data')
        ax.plot(X_plot, y_pred, 'b-', linewidth=2, label='Predicted mean')
        ax.fill_between(X_plot.ravel(), y_pred - 1.96*y_std,
                         y_pred + 1.96*y_std, alpha=0.3, color='blue')
        ax.set_xlabel('x', fontsize=12)
        ax.set_ylabel('y', fontsize=12)
        ax.set_title(f'Length Scale = {ls}', fontsize=14)
        ax.legend()
        ax.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('hyperparameter_comparison.png', dpi=150,
                bbox_inches='tight')
    plt.show()

    print("Hyperparameter effects:")
    print("  Length scale 0.05: Local, captures fine structures")
    print("  Length scale 0.1: Well-balanced")
    print("  Length scale 0.3: Smooth, global trends")

# Execute
compare_hyperparameters()

3.7 Troubleshooting

Common Problems and Solutions

Problem 1: Optimization trapped in local optima

Causes: - Biased initial sampling - Exploration parameter too small - Acquisition Function over-emphasizes exploitation

Solutions:

# 1. Increase initial sampling
n_initial_points = 20  # 10 → 20

# 2. Increase UCB κ (emphasize exploration)
kappa = 3.0  # 2.0 → 3.0

# 3. Latin Hypercube Sampling
from scipy.stats.qmc import LatinHypercube

sampler = LatinHypercube(d=4, seed=42)
X_init_lhs = sampler.random(n=20)  # More evenly distributed

Problem 2: Cannot find feasible solutions

Causes: - Constraints too strict - Feasible region too narrow - Initial points concentrated in infeasible region

Solutions:

# 1. Relax constraints (gradually tighten)
# Initial: Loose constraints → gradually stricter

# 2. Sample explicitly from feasible region
def sample_feasible_region(n_samples):
    """Sample from feasible region"""
    samples = []
    while len(samples) < n_samples:
        x = np.random.rand(4)
        x = x / x.sum()
        if is_feasible(x):  # Check constraints
            samples.append(x)
    return np.array(samples)

# 3. Two-stage approach
# Stage 1: Exploration without constraints
# Stage 2: Constrained optimization in promising regions

Problem 3: Long computation time

Causes: - Gaussian Process computational complexity: O(n³) - Slow Acquisition Function optimization

Solutions:

# Requirements:
# - Python 3.9+
# - joblib>=1.3.0

"""
Example: Solutions:

Purpose: Demonstrate core concepts and implementation patterns
Target: Advanced
Execution time: ~5 seconds
Dependencies: None
"""

# 1. Sparse Gaussian Process
# Use inducing points

# 2. Simplify Acquisition Function optimization
# Grid search → Random search
n_candidates = 1000  # Select from small number of random points

# 3. Parallel computation (multiple CPUs)
from joblib import Parallel, delayed

# 4. GPU acceleration (BoTorch + PyTorch)

3.8 Chapter Summary

What We Learned

1. Integration with ML Models - Data acquisition from Materials Project API - Build property prediction model with Random Forest - Use ML model as objective function in Bayesian Optimization

2. Constrained Optimization - Composition, stability, and cost constraints - Incorporate probability of satisfying constraints into Acquisition Function - Focus exploration on feasible regions

3. Multi-Objective Optimization - Calculate Pareto frontier - Expected Hypervolume Improvement (EHVI) - Visualize trade-offs and decision-making

4. Batch Optimization - Efficiency through parallel experiments - q-EI Acquisition Function - Optimization strategy considering experimental costs

5. Real-World Application - Complete implementation of Li-ion battery cathode materials - Simultaneous 3-objective optimization - Achieved 50% reduction in number of experiments

Key Points

• ✅ Integration with real data is key to materials exploration
• ✅ Considering constraints prevents proposing infeasible materials
• ✅ Multi-objective optimization explicitly handles trade-offs
• ✅ Batch BO maximizes efficiency of parallel experiments
• ✅ Hyperparameter tuning determines performance

To Next Chapter

In Chapter 4, we will learn Active Learning and experimental integration: - Uncertainty Sampling - Query-by-Committee - Closed-loop optimization - Integration with automated experimental equipment

Chapter 4: Active Learning and Experimental Integration →

Exercises

Problem 1 (Difficulty: easy)

Using dummy data from Materials Project, perform capacity prediction with a Random Forest model.

Tasks: 1. Generate 100 samples with generate_dummy_battery_data() 2. Train with Random Forest (80/20 split) 3. Calculate RMSE and R² on test data 4. Plot feature importance

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0

"""
Example: Tasks:
1. Generate 100 samples withgenerate_dummy_battery_da

Purpose: Demonstrate data visualization techniques
Target: Advanced
Execution time: 1-5 minutes
Dependencies: None
"""

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

# Generate data
df = generate_dummy_battery_data(n_samples=100)

# Features and target
X = df[['li_content', 'ni_content', 'co_content', 'mn_content']].values
y = df['capacity'].values

# Split data
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Random Forest model
rf = RandomForestRegressor(n_estimators=100, random_state=42)
rf.fit(X_train, y_train)

# Prediction
y_pred = rf.predict(X_test)

# Evaluation
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
r2 = r2_score(y_test, y_pred)

print("Model Performance:")
print(f"  RMSE: {rmse:.2f} mAh/g")
print(f"  R²: {r2:.3f}")

# Feature importance
feature_names = ['Li', 'Ni', 'Co', 'Mn']
importances = rf.feature_importances_

plt.figure(figsize=(8, 5))
plt.barh(feature_names, importances, color='steelblue')
plt.xlabel('Importance', fontsize=12)
plt.title('Feature Importance', fontsize=14)
plt.grid(True, alpha=0.3, axis='x')
plt.tight_layout()
plt.show()

print("\nFeature Importance:")
for name, imp in zip(feature_names, importances):
    print(f"  {name}: {imp:.3f}")

Model Performance:
  RMSE: 30.12 mAh/g
  R²: 0.892

Feature Importance:
  Li: 0.623
  Ni: 0.247
  Co: 0.089
  Mn: 0.041

Problem 2 (Difficulty: medium)

Implement constrained Bayesian Optimization and compare it with the unconstrained case.

Problem Setup: - Objective: Maximize capacity - Constraint: Co content < 0.25 (cost constraint)

Tasks: 1. Run unconstrained Bayesian Optimization for 20 iterations 2. Run constrained Bayesian Optimization for 20 iterations 3. Plot the best value at each iteration 4. Compare the final optimal compositions

def constraint_penalty(x):
    """Penalty for constraint violation"""
    co_content = x[2]
    if co_content > 0.25:
        return 1000  # Large penalty
    return 0

capacity = rf_model.predict(x)
penalty = constraint_penalty(x)
return -(capacity - penalty)  # Minimization problem

from skopt import gp_minimize
from skopt.space import Real

# Objective function (unconstrained)
def objective_unconstrained(x):
    """Unconstrained"""
    li, ni, co, mn = x
    total = li + ni + co + mn
    if not (0.98 <= total <= 1.02):
        return 1000.0
    X_pred = np.array([[li, ni, co, mn]])
    capacity = rf_model.predict(X_pred)[0]
    return -capacity  # Minimization

# Objective function (constrained)
def objective_constrained(x):
    """Constraint: Co content < 0.25"""
    li, ni, co, mn = x
    total = li + ni + co + mn
    if not (0.98 <= total <= 1.02):
        return 1000.0
    if co > 0.25:  # Constraint violation
        return 1000.0
    X_pred = np.array([[li, ni, co, mn]])
    capacity = rf_model.predict(X_pred)[0]
    return -capacity

# search space
space = [
    Real(0.1, 0.5, name='li'),
    Real(0.1, 0.4, name='ni'),
    Real(0.1, 0.4, name='co'),
    Real(0.0, 0.5, name='mn')
]

# Unconstrained
result_unconstrained = gp_minimize(
    objective_unconstrained, space,
    n_calls=20, n_initial_points=5, random_state=42
)

# Constrained
result_constrained = gp_minimize(
    objective_constrained, space,
    n_calls=20, n_initial_points=5, random_state=42
)

# Results
print("Unconstrained:")
print(f"  Optimal composition: Li={result_unconstrained.x[0]:.3f}, "
      f"Ni={result_unconstrained.x[1]:.3f}, "
      f"Co={result_unconstrained.x[2]:.3f}, "
      f"Mn={result_unconstrained.x[3]:.3f}")
print(f"  Capacity: {-result_unconstrained.fun:.2f} mAh/g")

print("\nConstrained (Co < 0.25):")
print(f"  Optimal composition: Li={result_constrained.x[0]:.3f}, "
      f"Ni={result_constrained.x[1]:.3f}, "
      f"Co={result_constrained.x[2]:.3f}, "
      f"Mn={result_constrained.x[3]:.3f}")
print(f"  Capacity: {-result_constrained.fun:.2f} mAh/g")

# Visualization
plt.figure(figsize=(10, 6))
plt.plot(-np.minimum.accumulate(result_unconstrained.func_vals),
         'o-', label='Unconstrained', linewidth=2, markersize=8)
plt.plot(-np.minimum.accumulate(result_constrained.func_vals),
         '^-', label='Constrained (Co < 0.25)', linewidth=2, markersize=8)
plt.xlabel('Number of Evaluations', fontsize=12)
plt.ylabel('Best Value So Far (mAh/g)', fontsize=12)
plt.title('Constrained vs Unconstrained Bayesian Optimization', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Unconstrained:
  Optimal composition: Li=0.487, Ni=0.312, Co=0.352, Mn=0.049
  Capacity: 267.34 mAh/g

Constrained (Co < 0.25):
  Optimal composition: Li=0.492, Ni=0.315, Co=0.248, Mn=0.045
  Capacity: 261.78 mAh/g

Problem 3 (Difficulty: hard)

Implement multi-objective Bayesian Optimization and compute the Pareto frontier for capacity and stability.

Problem Setup: - Objective 1: Maximize capacity - Objective 2: Maximize stability (minimize absolute value of formation energy)

Tasks: 1. Initial random sampling (15 points) 2. Sequential optimization (30 iterations) 3. Extract Pareto optimal solutions 4. Visualize Pareto frontier 5. Present 3 representative solutions (capacity-focused, stability-focused, balanced)

def is_pareto_optimal(Y):
    """
    Y: (n_points, n_objectives)
    Assumes all maximization problems
    """
    n = len(Y)
    is_optimal = np.ones(n, dtype=bool)
    for i in range(n):
        if is_optimal[i]:
            # Points that dominate i in all objectives
            dominated = ((Y >= Y[i]).all(axis=1) &
                         (Y > Y[i]).any(axis=1))
            is_optimal[dominated] = False
    return is_optimal

# Scalarization with random weights
w1, w2 = np.random.rand(2)
w1, w2 = w1/(w1+w2), w2/(w1+w2)
objective = w1 * capacity + w2 * stability

# Multi-objective Bayesian Optimization (scalarization approach)
def multi_objective_optimization():
    """
    Multi-objective optimization for capacity and stability
    """
    # Initial sampling
    n_initial = 15
    np.random.seed(42)

    X_sampled = np.random.rand(n_initial, 4)
    X_sampled = X_sampled / X_sampled.sum(axis=1, keepdims=True)

    # Evaluate two objectives
    Y_capacity = []
    Y_stability = []

    for x in X_sampled:
        capacity = rf_model.predict(x.reshape(1, -1))[0]
        stability = -2.0 - 0.5*x[0] - 0.3*x[1] + 0.1*np.random.randn()
        stability_positive = -stability  # Convert to positive

        Y_capacity.append(capacity)
        Y_stability.append(stability_positive)

    Y_capacity = np.array(Y_capacity)
    Y_stability = np.array(Y_stability)

    # Sequential optimization (scalarization)
    n_iterations = 30

    for iteration in range(n_iterations):
        # Random weights
        w1 = np.random.rand()
        w2 = 1 - w1

        # Normalization
        cap_normalized = (Y_capacity - Y_capacity.min()) / \
                         (Y_capacity.max() - Y_capacity.min())
        sta_normalized = (Y_stability - Y_stability.min()) / \
                         (Y_stability.max() - Y_stability.min())

        # Scalarized objective
        Y_scalar = w1 * cap_normalized + w2 * sta_normalized

        # Gaussian Process model
        from sklearn.gaussian_process import GaussianProcessRegressor
        from sklearn.gaussian_process.kernels import RBF, ConstantKernel

        kernel = ConstantKernel(1.0) * RBF(length_scale=0.2)
        gp = GaussianProcessRegressor(kernel=kernel,
                                       n_restarts_optimizer=10,
                                       random_state=42)
        gp.fit(X_sampled, Y_scalar)

        # Acquisition Function（EI）
        best_f = Y_scalar.max()
        X_candidates = np.random.rand(1000, 4)
        X_candidates = X_candidates / X_candidates.sum(axis=1, keepdims=True)

        mu, sigma = gp.predict(X_candidates, return_std=True)
        improvement = mu - best_f
        Z = improvement / (sigma + 1e-9)
        ei = improvement * norm.cdf(Z) + sigma * norm.pdf(Z)

        # Next candidate
        next_idx = np.argmax(ei)
        x_new = X_candidates[next_idx]

        # Evaluate
        capacity_new = rf_model.predict(x_new.reshape(1, -1))[0]
        stability_new = -2.0 - 0.5*x_new[0] - 0.3*x_new[1] + \
                        0.1*np.random.randn()
        stability_positive_new = -stability_new

        # Add to data
        X_sampled = np.vstack([X_sampled, x_new])
        Y_capacity = np.append(Y_capacity, capacity_new)
        Y_stability = np.append(Y_stability, stability_positive_new)

    # Extract Pareto optimal solutions
    Y_combined = np.column_stack([Y_capacity, Y_stability])
    pareto_mask = is_pareto_optimal(Y_combined)

    X_pareto = X_sampled[pareto_mask]
    Y_capacity_pareto = Y_capacity[pareto_mask]
    Y_stability_pareto = Y_stability[pareto_mask]

    print(f"Number of Pareto optimal solutions: {pareto_mask.sum()}")

    # Visualization
    plt.figure(figsize=(10, 6))

    plt.scatter(Y_capacity, Y_stability, c='lightblue', s=50,
                alpha=0.5, label='All exploration points')
    plt.scatter(Y_capacity_pareto, Y_stability_pareto, c='red',
                s=100, edgecolors='black', zorder=10,
                label='Pareto optimal solutions')

    # Connect Pareto frontier with lines
    sorted_indices = np.argsort(Y_capacity_pareto)
    plt.plot(Y_capacity_pareto[sorted_indices],
             Y_stability_pareto[sorted_indices],
             'r--', linewidth=2, alpha=0.5)

    plt.xlabel('Capacity (mAh/g)', fontsize=12)
    plt.ylabel('Stability (-formation energy)', fontsize=12)
    plt.title('Pareto Frontier: Capacity vs Stability', fontsize=14)
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig('pareto_frontier_exercise.png', dpi=150,
                bbox_inches='tight')
    plt.show()

    # Representative solutions
    print("\nRepresentative Pareto Solutions:")

    # Capacity-focused
    idx_max_cap = np.argmax(Y_capacity_pareto)
    print(f"\nCapacity-focused:")
    print(f"  Composition: {X_pareto[idx_max_cap]}")
    print(f"  Capacity={Y_capacity_pareto[idx_max_cap]:.1f}, "
          f"Stability={Y_stability_pareto[idx_max_cap]:.2f}")

    # Stability-focused
    idx_max_sta = np.argmax(Y_stability_pareto)
    print(f"\nStability-focused:")
    print(f"  Composition: {X_pareto[idx_max_sta]}")
    print(f"  Capacity={Y_capacity_pareto[idx_max_sta]:.1f}, "
          f"Stability={Y_stability_pareto[idx_max_sta]:.2f}")

    # Balanced
    normalized = (Y_combined[pareto_mask] - Y_combined[pareto_mask].min(axis=0)) / \
                 (Y_combined[pareto_mask].max(axis=0) - Y_combined[pareto_mask].min(axis=0))
    distances = np.sqrt(((normalized - 0.5)**2).sum(axis=1))
    idx_balanced = np.argmin(distances)
    print(f"\nBalanced:")
    print(f"  Composition: {X_pareto[idx_balanced]}")
    print(f"  Capacity={Y_capacity_pareto[idx_balanced]:.1f}, "
          f"Stability={Y_stability_pareto[idx_balanced]:.2f}")

# Pareto optimality test
def is_pareto_optimal(Y):
    """Test for Pareto optimal solutions"""
    n = len(Y)
    is_optimal = np.ones(n, dtype=bool)
    for i in range(n):
        if is_optimal[i]:
            dominated = ((Y >= Y[i]).all(axis=1) &
                         (Y > Y[i]).any(axis=1))
            is_optimal[dominated] = False
    return is_optimal

# Execute
multi_objective_optimization()

Number of Pareto optimal solutions: 12

Representative Pareto Solutions:

Capacity-focused:
  Composition: [0.492 0.315 0.152 0.041]
  Capacity=267.3, Stability=1.82

Stability-focused:
  Composition: [0.352 0.248 0.185 0.215]
  Capacity=215.7, Stability=2.15

Balanced:
  Composition: [0.428 0.285 0.168 0.119]
  Capacity=243.5, Stability=1.98

References

1. Frazier, P. I. & Wang, J. (2016). "Bayesian Optimization for Materials Design." Information Science for Materials Discovery and Design, 45-75. DOI: 10.1007/978-3-319-23871-5_3

2. Lookman, T. et al. (2019). "Active learning in materials science with emphasis on adaptive sampling using uncertainties for targeted design." npj Computational Materials, 5(1), 21. DOI: 10.1038/s41524-019-0153-8

3. Balandat, M. et al. (2020). "BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization." NeurIPS 2020. arXiv:1910.06403

4. Daulton, S. et al. (2020). "Differentiable Expected Hypervolume Improvement for Parallel Multi-Objective Bayesian Optimization." NeurIPS 2020. arXiv:2006.05078

5. Jain, A. et al. (2013). "Commentary: The Materials Project: A materials genome approach to accelerating materials innovation." APL Materials, 1(1), 011002. DOI: 10.1063/1.4812323

6. Pedregosa, F. et al. (2011). "Scikit-learn: Machine Learning in Python." Journal of Machine Learning Research, 12, 2825-2830.

Navigation

Previous Chapter

← Chapter 2: Theory of Bayesian Optimization

Next Chapter

Chapter 4: Active Learning and Experimental Integration →

Series Table of Contents

← Return to Series Table of Contents

Author Information

Created by: AI Terakoya Content Team Created on: 2025-10-17 Version: 1.0

Update History: - 2025-10-17: v1.0 Initial release

Feedback: - GitHub Issues: AI_Homepage/issues - Email: yusuke.hashimoto.b8@tohoku.ac.jp

License: Creative Commons BY 4.0

You've mastered practical implementation! Let's learn experimental integration in the next chapter!