Chapter 3: Closed-Loop Optimization

Study time: 25-30 minutes

Introduction

The true value of experimental automation goes beyond mere throughput improvement. By integrating it with Bayesian optimization and active learning, we can build systems that autonomously search for optimal materials without human intervention.

In this chapter, we learn the technology to automate the closed loop of experiment → measurement → analysis → prediction → next-experiment proposal. We will acquire Python implementations that can be applied to real materials-exploration problems, such as optimizing the emission wavelength of quantum dots and maximizing catalyst activity.

Learning Objectives

By studying this chapter, you will be able to master the following:

1. The closed-loop concept: The integrated architecture of experiment and machine learning
2. Integration with Bayesian optimization: Automatic proposal of next experiment candidates via Gaussian processes
3. Active learning: Efficient data-collection strategies
4. Python implementation: Practical code using scikit-optimize and BoTorch
5. Simulation environment: Virtual robotic experiments with PyBullet
6. Demonstration with real data: Automatic optimization of quantum dot emission wavelength

3.1 The Concept of Closed-Loop Optimization

3.1.1 Conventional Manual Optimization vs Closed Loop

Conventional manual optimization: 1. The researcher decides experimental conditions by intuition 2. Runs the experiment (1-2 days) 3. Analyzes the data 4. Devises the next experimental conditions (1-2 days) 5. Repeats steps 1-4

Problems: - Human cognitive bias (tends to fall into local optima) - Difficulty exploring multi-dimensional parameter spaces - Time lag between experiment and analysis - Experiments stop at night and on weekends

Closed-loop optimization:

Advantages: - Operates 24 hours a day, 365 days a year - No cognitive bias (data-driven) - Strong in multi-dimensional exploration - Experiment → analysis → next experiment completes in minutes

3.1.2 Mathematical Formulation of the Closed Loop

Optimization problem: $$ \mathbf{x}^* = \arg\max_{\mathbf{x} \in \mathcal{X}} f(\mathbf{x}) $$

Here: - $\mathbf{x}$: experimental conditions (temperature, time, composition, etc.) - $f(\mathbf{x})$: objective function (catalyst activity, emission wavelength, etc.) - $\mathcal{X}$: search space (the range of allowable experimental conditions)

Challenge: $f(\mathbf{x})$ is an unknown function that can only be evaluated by experiment (evaluation cost is high)

Solution: Bayesian optimization 1. Surrogate model: Approximate $f(\mathbf{x})$ with a Gaussian process 2. Acquisition function: Propose where to experiment next (balancing exploration and exploitation) 3. Iteration: Repeat experiment → model update → next candidate proposal

3.2 Integration with Bayesian Optimization

3.2.1 Fundamentals of Gaussian Processes

A Gaussian Process (GP) is a probability distribution over a function space.

Mathematical definition: $$ f(\mathbf{x}) \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}')) $$

• $m(\mathbf{x})$: mean function (usually 0)
• $k(\mathbf{x}, \mathbf{x}')$: kernel function (similarity between two points)

Prediction: Given observed data $\mathcal{D} = {(\mathbf{x}_i, y_i)}_{i=1}^n$, the prediction at a new point $\mathbf{x}_*$:

$$ \begin{aligned} \mu(\mathbf{x}_*) &= \mathbf{k}_*^\top (\mathbf{K} + \sigma^2 \mathbf{I})^{-1} \mathbf{y} \ \sigma^2(\mathbf{x}_*) &= k(\mathbf{x}_*, \mathbf{x}_*) - \mathbf{k}_*^\top (\mathbf{K} + \sigma^2 \mathbf{I})^{-1} \mathbf{k}_* \end{aligned} $$

• $\mu(\mathbf{x}_*)$: predictive mean (expected value)
• $\sigma^2(\mathbf{x}_*)$: predictive variance (uncertainty)
• $\mathbf{K}$: kernel matrix, $[\mathbf{K}]_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)$

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

# True function (assumed unknown)
def true_function(x):
    """
    The true function to be optimized (e.g., catalyst activity)
    Unknown in a real experiment; only evaluable by experiment
    """
    return np.sin(3*x) + 0.3*np.cos(10*x) + 0.5*x

# Initial experimental data (3 points)
np.random.seed(42)
X_init = np.array([0.2, 0.5, 0.8]).reshape(-1, 1)
y_init = true_function(X_init).ravel() + np.random.normal(0, 0.05, 3)  # noise

# Build the Gaussian process model
kernel = C(1.0, (1e-3, 1e3)) * RBF(length_scale=0.2, length_scale_bounds=(1e-2, 1e2))
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, alpha=0.05**2)

# Train the model
gp.fit(X_init, y_init)

# Predict (over the entire search space)
X_pred = np.linspace(0, 1, 100).reshape(-1, 1)
y_pred, sigma = gp.predict(X_pred, return_std=True)

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

# True function
plt.plot(X_pred, true_function(X_pred), 'k--', label='True function (unknown)', linewidth=2)

# Gaussian process prediction
plt.plot(X_pred, y_pred, 'b-', label='GP prediction (mean)', linewidth=2)
plt.fill_between(X_pred.ravel(),
                 y_pred - 1.96*sigma,  # 95% confidence interval
                 y_pred + 1.96*sigma,
                 alpha=0.2, color='blue', label='95% confidence interval')

# Observed points
plt.plot(X_init, y_init, 'ro', markersize=12, label='Initial observations', zorder=10)

plt.xlabel('Experimental condition x', fontsize=12)
plt.ylabel('Objective function f(x)', fontsize=12)
plt.title('Function approximation by a Gaussian process', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('gaussian_process_approximation.png', dpi=300, bbox_inches='tight')
plt.show()

print("Gaussian process training complete")
print(f"Optimized kernel parameters: {gp.kernel_}")

Code explanation: 1. Kernel function: Uses the RBF (Radial Basis Function) kernel - length_scale: controls the smoothness of the function 2. Noise: The alpha parameter accounts for measurement noise 3. Prediction: Computes the mean $\mu(\mathbf{x})$ and standard deviation $\sigma(\mathbf{x})$ simultaneously 4. Confidence interval: 95% confidence interval as $\mu \pm 1.96\sigma$

3.2.2 Acquisition Function

A function that decides where to experiment next. It automatically adjusts the trade-off between exploration and exploitation.

Main acquisition functions:

1. Expected Improvement (EI): $$ \text{EI}(\mathbf{x}) = \mathbb{E}[\max(f(\mathbf{x}) - f^+, 0)] $$ The expected improvement over the best observed value $f^+$

2. Upper Confidence Bound (UCB): $$ \text{UCB}(\mathbf{x}) = \mu(\mathbf{x}) + \kappa \sigma(\mathbf{x}) $$ Predictive mean + uncertainty (tuned by $\kappa$)

3. Probability of Improvement (PI): $$ \text{PI}(\mathbf{x}) = P(f(\mathbf{x}) > f^+) $$

from scipy.stats import norm

def expected_improvement(X, gp, f_best, xi=0.01):
    """
    Expected Improvement acquisition function

    Args:
        X: candidate points
        gp: Gaussian process model
        f_best: current best value
        xi: exploitation-exploration trade-off (smaller favors exploitation)

    Returns:
        EI values
    """
    mu, sigma = gp.predict(X, return_std=True)

    # Avoid sigma=0 (already observed)
    sigma = np.maximum(sigma, 1e-9)

    # Compute EI
    z = (mu - f_best - xi) / sigma
    ei = (mu - f_best - xi) * norm.cdf(z) + sigma * norm.pdf(z)

    return ei

def upper_confidence_bound(X, gp, kappa=2.0):
    """
    Upper Confidence Bound acquisition function

    Args:
        X: candidate points
        gp: Gaussian process model
        kappa: exploration weight (larger favors exploration)

    Returns:
        UCB values
    """
    mu, sigma = gp.predict(X, return_std=True)
    return mu + kappa * sigma


# Visualize the acquisition functions
X_pred = np.linspace(0, 1, 100).reshape(-1, 1)
y_pred, sigma = gp.predict(X_pred, return_std=True)

# Current best value
f_best = y_init.max()

# Compute the acquisition functions
ei_values = expected_improvement(X_pred, gp, f_best)
ucb_values = upper_confidence_bound(X_pred, gp, kappa=2.0)

# Visualization
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 12))

# (1) Gaussian process prediction
ax1.plot(X_pred, true_function(X_pred), 'k--', label='True function', linewidth=2)
ax1.plot(X_pred, y_pred, 'b-', label='GP prediction', linewidth=2)
ax1.fill_between(X_pred.ravel(), y_pred - 1.96*sigma, y_pred + 1.96*sigma,
                 alpha=0.2, color='blue')
ax1.plot(X_init, y_init, 'ro', markersize=12, label='Observations')
ax1.axhline(y=f_best, color='red', linestyle=':', label=f'Current best ({f_best:.2f})')
ax1.set_ylabel('Objective function f(x)', fontsize=12)
ax1.set_title('(1) Gaussian process prediction', fontsize=13, fontweight='bold')
ax1.legend()
ax1.grid(alpha=0.3)

# (2) Expected Improvement
next_x_ei = X_pred[np.argmax(ei_values)]
ax2.plot(X_pred, ei_values, 'g-', linewidth=2)
ax2.axvline(x=next_x_ei, color='red', linestyle='--', label=f'Next candidate (x={next_x_ei[0]:.3f})')
ax2.fill_between(X_pred.ravel(), 0, ei_values.ravel(), alpha=0.3, color='green')
ax2.set_ylabel('EI(x)', fontsize=12)
ax2.set_title('(2) Expected Improvement acquisition function', fontsize=13, fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)

# (3) Upper Confidence Bound
next_x_ucb = X_pred[np.argmax(ucb_values)]
ax3.plot(X_pred, ucb_values, 'm-', linewidth=2, label='UCB')
ax3.plot(X_pred, y_pred, 'b--', linewidth=1, alpha=0.5, label='GP mean')
ax3.axvline(x=next_x_ucb, color='red', linestyle='--', label=f'Next candidate (x={next_x_ucb[0]:.3f})')
ax3.set_xlabel('Experimental condition x', fontsize=12)
ax3.set_ylabel('UCB(x)', fontsize=12)
ax3.set_title('(3) Upper Confidence Bound acquisition function', fontsize=13, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('acquisition_functions.png', dpi=300, bbox_inches='tight')
plt.show()

print(f"Next candidate by EI: x = {next_x_ei[0]:.3f}")
print(f"Next candidate by UCB: x = {next_x_ucb[0]:.3f}")

Characteristics of the acquisition functions: - EI: balanced, the most widely used - UCB: exploration-focused, tunable via $\kappa$ - PI: conservative, prioritizes improvement in known regions

3.3 Implementing the Closed Loop

3.3.1 Implementation with scikit-optimize

from skopt import gp_minimize
from skopt.space import Real
from skopt.utils import use_named_args
from skopt.plots import plot_convergence, plot_objective

# The function to optimize (simulation of a robotic experiment)
def robot_experiment(x):
    """
    Simulation function for a robotic experiment

    Args:
        x: experimental conditions (e.g., [temperature, time])

    Returns:
        -f(x): converted to a minimization problem (scikit-optimize minimizes)
    """
    # In a real experiment, this part is the robotic synthesis and measurement
    result = true_function(np.array([[x[0]]]))[0]

    # Noise (experimental error)
    result += np.random.normal(0, 0.05)

    print(f"Running experiment: x={x[0]:.3f}, result={result:.3f}")

    # Convert maximization to minimization (scikit-optimize minimizes)
    return -result


# Define the search space
space = [Real(0.0, 1.0, name='x')]

# Run Bayesian optimization
n_calls = 20  # number of experiments
result = gp_minimize(
    robot_experiment,  # objective function (robotic experiment)
    space,  # search space
    n_calls=n_calls,  # total number of experiments
    n_initial_points=5,  # number of initial random experiments
    acq_func='EI',  # acquisition function
    random_state=42
)

print(f"\nOptimization complete!")
print(f"Optimal condition: x = {result.x[0]:.3f}")
print(f"Optimal value: f(x) = {-result.fun:.3f}")
print(f"Number of experiments: {n_calls}")

# Convergence plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# (1) Convergence history
plot_convergence(result, ax=ax1)
ax1.set_title('Optimization convergence', fontsize=14, fontweight='bold')

# (2) Explored points
X_evaluated = np.array([x[0] for x in result.x_iters])
y_evaluated = -np.array(result.func_vals)  # restore the original scale

X_plot = np.linspace(0, 1, 100)
y_true = true_function(X_plot.reshape(-1, 1)).ravel()

ax2.plot(X_plot, y_true, 'k--', label='True function', linewidth=2)
ax2.plot(X_evaluated, y_evaluated, 'ro-', markersize=8, label='Evaluated points', alpha=0.6)
ax2.plot(result.x[0], -result.fun, 'g*', markersize=20, label='Optimal solution')
ax2.set_xlabel('Experimental condition x', fontsize=12)
ax2.set_ylabel('Objective function f(x)', fontsize=12)
ax2.set_title('Search trajectory', fontsize=14, fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('bayesian_optimization_result.png', dpi=300, bbox_inches='tight')
plt.show()

3.3.2 A Complete Closed-Loop System

class ClosedLoopOptimization:
    """
    Closed-loop optimization system
    Experiment -> Measurement -> Analysis -> Prediction -> Next experiment
    """

    def __init__(self, robot_controller, sensor_controller):
        """
        Args:
            robot_controller: robot control interface
            sensor_controller: sensor control interface
        """
        self.robot = robot_controller
        self.sensor = sensor_controller
        self.data = []

    def run_experiment(self, conditions):
        """
        A single experiment cycle

        Args:
            conditions: experimental conditions (dict)

        Returns:
            measurement: the measured value
        """
        print(f"\n--- Experiment {len(self.data)+1} ---")
        print(f"Conditions: {conditions}")

        # (1) Sample preparation
        print("  [1/4] Robotic sample preparation...")
        self.robot.prepare_sample(conditions)

        # (2) Measurement
        print("  [2/4] Sensor measurement...")
        measurement = self.sensor.measure()

        # (3) Data recording
        print("  [3/4] Recording data...")
        self.data.append({'conditions': conditions, 'measurement': measurement})

        # (4) Display result
        print(f"  [4/4] Measurement result: {measurement:.3f}")

        return measurement

    def optimize(self, n_iterations=10):
        """
        Run closed-loop optimization

        Args:
            n_iterations: number of optimization iterations
        """
        from skopt import Optimizer

        # Search space (e.g., temperature and pH)
        space = [
            Real(50, 150, name='temperature'),  # temperature (degrees C)
            Real(4, 10, name='pH')  # pH
        ]

        optimizer = Optimizer(space, base_estimator='GP', acq_func='EI')

        print("=" * 50)
        print("Starting closed-loop optimization")
        print("=" * 50)

        for i in range(n_iterations):
            # Propose the next experimental conditions
            next_conditions = optimizer.ask()
            conditions_dict = {'temperature': next_conditions[0], 'pH': next_conditions[1]}

            # Run the experiment
            result = self.run_experiment(conditions_dict)

            # Report the result to the optimizer
            optimizer.tell(next_conditions, -result)  # maximization -> minimization

        # Optimal conditions
        best_result = optimizer.get_result()
        print("\n" + "=" * 50)
        print("Optimization complete")
        print("=" * 50)
        print(f"Optimal conditions: temperature={best_result.x[0]:.1f} C, pH={best_result.x[1]:.1f}")
        print(f"Optimal value: {-best_result.fun:.3f}")

        return best_result


# Controllers for simulation
class RobotControllerSimulator:
    def prepare_sample(self, conditions):
        """Simulation of sample preparation"""
        time.sleep(0.1)  # actually several minutes

class SensorControllerSimulator:
    def __init__(self):
        self.measurement_count = 0

    def measure(self):
        """Simulation of measurement (two-variable function)"""
        # A hypothetical objective function
        # Optimal value: temperature=100 C, pH=7
        temp = np.random.uniform(50, 150)
        ph = np.random.uniform(4, 10)
        result = -(temp - 100)**2 / 1000 - (ph - 7)**2 + 10
        result += np.random.normal(0, 0.2)  # noise
        self.measurement_count += 1
        time.sleep(0.1)  # actually tens of seconds to several minutes
        return result


# Run the closed loop
robot = RobotControllerSimulator()
sensor = SensorControllerSimulator()
closed_loop = ClosedLoopOptimization(robot, sensor)

result = closed_loop.optimize(n_iterations=15)

# Visualize the explored points (2D)
import pandas as pd

df = pd.DataFrame([
    {'temperature': d['conditions']['temperature'],
     'pH': d['conditions']['pH'],
     'measurement': d['measurement']}
    for d in closed_loop.data
])

fig = plt.figure(figsize=(12, 5))

# Temperature vs measurement
ax1 = fig.add_subplot(121)
ax1.scatter(df['temperature'], df['measurement'], c=range(len(df)), cmap='viridis', s=100, edgecolors='black')
ax1.set_xlabel('Temperature (C)', fontsize=12)
ax1.set_ylabel('Measured value', fontsize=12)
ax1.set_title('Relationship between temperature and measured value', fontsize=13, fontweight='bold')
ax1.grid(alpha=0.3)

# pH vs measurement
ax2 = fig.add_subplot(122)
scatter = ax2.scatter(df['pH'], df['measurement'], c=range(len(df)), cmap='viridis', s=100, edgecolors='black')
ax2.set_xlabel('pH', fontsize=12)
ax2.set_ylabel('Measured value', fontsize=12)
ax2.set_title('Relationship between pH and measured value', fontsize=13, fontweight='bold')
ax2.grid(alpha=0.3)

# Colorbar (experiment order)
cbar = plt.colorbar(scatter, ax=[ax1, ax2])
cbar.set_label('Experiment order', fontsize=11)

plt.tight_layout()
plt.savefig('closed_loop_exploration.png', dpi=300, bbox_inches='tight')
plt.show()

3.4 Automatic Optimization of Quantum Dot Emission Wavelength

A real materials-exploration case study: automatically optimizing the emission wavelength of quantum dots.

3.4.1 Problem Setup

Goal: Optimize the emission wavelength of CdSeS quantum dots to 520 nm (green)

Experimental parameters: - Cd/Se ratio (0.5-2.0) - Reaction temperature (150-300 degrees C) - Reaction time (5-60 minutes)

def quantum_dot_synthesis_simulator(cd_se_ratio, temperature, reaction_time):
    """
    Simulator for quantum dot synthesis and emission wavelength measurement

    Args:
        cd_se_ratio: Cd/Se ratio (0.5-2.0)
        temperature: reaction temperature (150-300 degrees C)
        reaction_time: reaction time (5-60 minutes)

    Returns:
        emission_wavelength: emission wavelength (nm)
    """
    # A simulation model based on empirical rules
    # In a real experiment, this part is robotic synthesis + fluorescence spectroscopy

    # Base wavelength (composition-dependent)
    base_wavelength = 480 + 100 * (cd_se_ratio - 0.5) / 1.5

    # Temperature effect (particle-size control)
    temp_effect = 0.2 * (temperature - 225)

    # Time effect (growth time)
    time_effect = 0.3 * (reaction_time - 32.5)

    # Overall emission wavelength
    emission_wavelength = base_wavelength + temp_effect + time_effect

    # Noise (experimental error)
    emission_wavelength += np.random.normal(0, 3)

    return emission_wavelength


# Closed-loop optimization
from skopt import gp_minimize
from skopt.space import Real

# Target wavelength
target_wavelength = 520  # nm (green)

def objective_function(params):
    """
    Objective function: minimize the difference from the target wavelength

    Args:
        params: [cd_se_ratio, temperature, reaction_time]

    Returns:
        error: difference from the target (to minimize)
    """
    cd_se_ratio, temperature, reaction_time = params

    # Quantum dot synthesis (simulation)
    emission = quantum_dot_synthesis_simulator(cd_se_ratio, temperature, reaction_time)

    # Error (absolute difference from the target wavelength)
    error = abs(emission - target_wavelength)

    print(f"Cd/Se={cd_se_ratio:.2f}, T={temperature:.0f}C, t={reaction_time:.0f}min -> lambda={emission:.1f}nm (error: {error:.1f}nm)")

    return error


# Search space
space = [
    Real(0.5, 2.0, name='cd_se_ratio'),
    Real(150, 300, name='temperature'),
    Real(5, 60, name='reaction_time')
]

print("=" * 70)
print("Automatic optimization of quantum dot emission wavelength")
print(f"Target wavelength: {target_wavelength}nm (green)")
print("=" * 70)

# Run Bayesian optimization
result = gp_minimize(
    objective_function,
    space,
    n_calls=30,  # 30 experiments
    n_initial_points=10,  # first 10 are random
    acq_func='EI',
    random_state=42
)

print("\n" + "=" * 70)
print("Optimization complete")
print("=" * 70)
print(f"Optimal conditions:")
print(f"  Cd/Se ratio: {result.x[0]:.2f}")
print(f"  Temperature: {result.x[1]:.0f} C")
print(f"  Reaction time: {result.x[2]:.0f} min")

# Measured value at the optimal conditions
optimized_wavelength = quantum_dot_synthesis_simulator(result.x[0], result.x[1], result.x[2])
print(f"\nEmission wavelength at optimal conditions: {optimized_wavelength:.1f}nm")
print(f"Error from target: {abs(optimized_wavelength - target_wavelength):.1f}nm")

# Visualize the results
fig = plt.figure(figsize=(14, 10))

# (1) Convergence history
ax1 = fig.add_subplot(221)
errors = result.func_vals
ax1.plot(errors, 'o-', linewidth=2, markersize=8, color='steelblue')
ax1.axhline(y=5, color='red', linestyle='--', label='Target error (<5nm)')
ax1.set_xlabel('Number of experiments', fontsize=12)
ax1.set_ylabel('Error (nm)', fontsize=12)
ax1.set_title('(1) Optimization convergence', fontsize=13, fontweight='bold')
ax1.legend()
ax1.grid(alpha=0.3)

# (2) Cd/Se ratio vs emission wavelength
ax2 = fig.add_subplot(222)
cd_se_ratios = [x[0] for x in result.x_iters]
wavelengths = [target_wavelength - e if i % 2 == 0 else target_wavelength + e
               for i, e in enumerate(errors)]  # simple reconstruction
ax2.scatter(cd_se_ratios, wavelengths, c=range(len(cd_se_ratios)), cmap='viridis',
            s=100, edgecolors='black')
ax2.axhline(y=target_wavelength, color='red', linestyle='--', label='Target wavelength')
ax2.set_xlabel('Cd/Se ratio', fontsize=12)
ax2.set_ylabel('Emission wavelength (nm)', fontsize=12)
ax2.set_title('(2) Cd/Se ratio and emission wavelength', fontsize=13, fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)

# (3) Temperature vs emission wavelength
ax3 = fig.add_subplot(223)
temperatures = [x[1] for x in result.x_iters]
ax3.scatter(temperatures, wavelengths, c=range(len(temperatures)), cmap='plasma',
            s=100, edgecolors='black')
ax3.axhline(y=target_wavelength, color='red', linestyle='--', label='Target wavelength')
ax3.set_xlabel('Temperature (C)', fontsize=12)
ax3.set_ylabel('Emission wavelength (nm)', fontsize=12)
ax3.set_title('(3) Temperature and emission wavelength', fontsize=13, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3)

# (4) Reaction time vs emission wavelength
ax4 = fig.add_subplot(224)
reaction_times = [x[2] for x in result.x_iters]
scatter = ax4.scatter(reaction_times, wavelengths, c=range(len(reaction_times)),
                      cmap='coolwarm', s=100, edgecolors='black')
ax4.axhline(y=target_wavelength, color='red', linestyle='--', label='Target wavelength')
ax4.set_xlabel('Reaction time (min)', fontsize=12)
ax4.set_ylabel('Emission wavelength (nm)', fontsize=12)
ax4.set_title('(4) Reaction time and emission wavelength', fontsize=13, fontweight='bold')
ax4.legend()
ax4.grid(alpha=0.3)

# Colorbar
cbar = plt.colorbar(scatter, ax=[ax2, ax3, ax4])
cbar.set_label('Experiment order', fontsize=11)

plt.tight_layout()
plt.savefig('quantum_dot_optimization.png', dpi=300, bbox_inches='tight')
plt.show()

Interpretation of the results: - First 10 experiments: random exploration sampling a wide range - From the 11th onward: Bayesian optimization concentrates on promising regions - 20th-30th: converges to within the target wavelength +/-5 nm

3.5 Exercises

Exercise 1: Comparing Acquisition Functions (Difficulty: Medium)

Solve the same optimization problem with the two acquisition functions EI (Expected Improvement) and UCB (Upper Confidence Bound), and compare the difference in their exploration behavior.

# Optimize the same problem with EI and UCB
results = {}

for acq_func in ['EI', 'LCB']:  # LCB = -UCB (because it is a minimization problem)
    print(f"\nAcquisition function: {acq_func}")
    result = gp_minimize(
        robot_experiment,
        space,
        n_calls=20,
        n_initial_points=5,
        acq_func=acq_func,
        random_state=42
    )
    results[acq_func] = result
    print(f"Optimal value: {-result.fun:.3f}, optimal condition: {result.x[0]:.3f}")

# Comparison plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

colors = {'EI': 'blue', 'LCB': 'red'}

# Convergence comparison
for acq_func, result in results.items():
    ax1.plot(-np.array(result.func_vals), label=acq_func, color=colors[acq_func], linewidth=2)

ax1.set_xlabel('Number of experiments', fontsize=12)
ax1.set_ylabel('Objective function value', fontsize=12)
ax1.set_title('Difference in convergence by acquisition function', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(alpha=0.3)

# Distribution of explored points
X_true = np.linspace(0, 1, 100)
y_true = true_function(X_true.reshape(-1, 1)).ravel()

ax2.plot(X_true, y_true, 'k--', label='True function', linewidth=2)

for acq_func, result in results.items():
    X_eval = np.array([x[0] for x in result.x_iters])
    y_eval = -np.array(result.func_vals)
    ax2.scatter(X_eval, y_eval, label=f'{acq_func}', s=80, alpha=0.6, color=colors[acq_func])

ax2.set_xlabel('Experimental condition x', fontsize=12)
ax2.set_ylabel('Objective function f(x)', fontsize=12)
ax2.set_title('Distribution of explored points', fontsize=14, fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('acquisition_function_comparison.png', dpi=300, bbox_inches='tight')
plt.show()

print("\nDiscussion:")
print("EI: balanced, efficiently explores a wide range")
print("UCB: exploration-focused, actively investigates high-uncertainty regions")

Exercise 2: Constrained Optimization (Difficulty: Hard)

In quantum dot synthesis, add the following constraints: - Temperature ≤ 250 degrees C (safety) - Reaction time ≤ 30 minutes (throughput)

def constrained_objective(params):
    """
    Constrained objective function

    Args:
        params: [cd_se_ratio, temperature, reaction_time]

    Returns:
        error: error + constraint-violation penalty
    """
    cd_se_ratio, temperature, reaction_time = params

    # Base objective function
    emission = quantum_dot_synthesis_simulator(cd_se_ratio, temperature, reaction_time)
    error = abs(emission - target_wavelength)

    # Constraint check
    penalty = 0
    constraints_violated = []

    if temperature > 250:
        penalty += 1000 * (temperature - 250)
        constraints_violated.append(f"temperature constraint violated({temperature:.0f}C > 250C)")

    if reaction_time > 30:
        penalty += 1000 * (reaction_time - 30)
        constraints_violated.append(f"time constraint violated({reaction_time:.0f}min > 30min)")

    if constraints_violated:
        print(f"  Constraint violation: {', '.join(constraints_violated)} -> penalty={penalty:.0f}")
    else:
        print(f"Cd/Se={cd_se_ratio:.2f}, T={temperature:.0f}C, t={reaction_time:.0f}min -> lambda={emission:.1f}nm")

    return error + penalty


# Constrained optimization
space_constrained = [
    Real(0.5, 2.0, name='cd_se_ratio'),
    Real(150, 300, name='temperature'),  # do not change the search space
    Real(5, 60, name='reaction_time')
]

print("Constrained optimization:")
print("  Temperature <= 250 C")
print("  Reaction time <= 30 min\n")

result_constrained = gp_minimize(
    constrained_objective,
    space_constrained,
    n_calls=30,
    n_initial_points=10,
    acq_func='EI',
    random_state=42
)

print(f"\nOptimal conditions:")
print(f"  Cd/Se ratio: {result_constrained.x[0]:.2f}")
print(f"  Temperature: {result_constrained.x[1]:.0f} C")
print(f"  Reaction time: {result_constrained.x[2]:.0f} min")

# Confirm constraint satisfaction
if result_constrained.x[1] <= 250 and result_constrained.x[2] <= 30:
    print("✅ All constraints are satisfied")
else:
    print("❌ There are constraint violations")

Exercise 3: Multi-Objective Optimization (Difficulty: Hard)

Simultaneously optimize the emission wavelength and the quantum yield (emission efficiency) of quantum dots.

Goals: - Emission wavelength: 520 nm - Quantum yield: maximize (0-100%)

def multi\_objective\_synthesis(cd\_se\_ratio, temperature, reaction\_time):
    """
    Quantum dot synthesis (emission wavelength + quantum yield)

    Returns:
        emission\_wavelength, quantum\_yield
    """
    # Emission wavelength
    emission = quantum\_dot\_synthesis\_simulator(cd\_se\_ratio, temperature, reaction\_time)

    # Quantum yield (hypothetical model)
    # High QY at low temperature and short time (fewer surface defects)
    qy = 80 - 0.2 * (temperature - 200) - 0.5 * (reaction\_time - 15)
    qy = np.clip(qy + np.random.normal(0, 2), 0, 100)

    return emission, qy


def multi\_objective\_function(params, w1=1.0, w2=0.1):
    """
    Multi-objective function (weighted sum)

    Args:
        params: [cd\_se\_ratio, temperature, reaction\_time]
        w1: weight for the wavelength error
        w2: weight for the quantum yield

    Returns:
        weighted\_error: weighted error
    """
    cd\_se\_ratio, temperature, reaction\_time = params

    emission, qy = multi\_objective\_synthesis(cd\_se\_ratio, temperature, reaction\_time)

    # Two objectives
    wavelength\_error = abs(emission - 520)
    qy\_loss = 100 - qy  # maximization -> minimization

    # Weighted sum
    weighted\_error = w1 * wavelength\_error + w2 * qy\_loss

    print(f"T={temperature:.0f}C, t={reaction\_time:.0f}min -> lambda={emission:.1f}nm, QY={qy:.1f}%")

    return weighted\_error


# Multi-objective optimization
result\_multi = gp\_minimize(
    multi\_objective\_function,
    space,
    n\_calls=30,
    n\_initial\_points=10,
    acq\_func='EI',
    random\_state=42
)

print(f"\nOptimal conditions:")
print(f"  Cd/Se ratio: {result\_multi.x[0]:.2f}")
print(f"  Temperature: {result\_multi.x[1]:.0f} C")
print(f"  Reaction time: {result\_multi.x[2]:.0f} min")

# Confirm performance at the optimal conditions
opt\_emission, opt\_qy = multi\_objective\_synthesis(result\_multi.x[0], result\_multi.x[1], result\_multi.x[2])
print(f"\nPerformance:")
print(f"  Emission wavelength: {opt\_emission:.1f}nm (target 520nm, error {abs(opt\_emission-520):.1f}nm)")
print(f"  Quantum yield: {opt\_qy:.1f}%")

Chapter Summary

In this chapter, we learned the theory and implementation of closed-loop optimization.

Key Points

1. The closed-loop concept: - An automated cycle of experiment → measurement → analysis → prediction → next experiment - Acceleration through 24/7/365 operation

2. Bayesian optimization: - Function approximation via Gaussian processes - Next-candidate selection via acquisition functions (EI, UCB) - Balancing exploration and exploitation

3. Implementation: - scikit-optimize: simple optimization - A custom closed-loop system - Robot and sensor integration

4. Real applications: - Quantum dot emission wavelength optimization - Target achieved in 30 experiments - Conventional trial and error (several weeks) → automation (1 day)

5. Advanced methods: - Constrained optimization - Multi-objective optimization

Preview of the Next Chapter

In Chapter 4, we will learn about remote experiments using a cloud lab (Emerald Cloud Lab). Without owning any instruments, you will experience the cutting-edge experimental environment where experiments are requested via an API and data is acquired automatically.

References

1. Shahriari, B. et al. (2016). "Taking the Human Out of the Loop: A Review of Bayesian Optimization." Proceedings of the IEEE, 104(1), 148-175.
2. Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.
3. MacLeod, B. P. et al. (2020). "Self-driving laboratory for accelerated discovery of thin-film materials." Science Advances, 6(20), eaaz8867.
4. Häse, F. et al. (2018). "Next-Generation Experimentation with Self-Driving Laboratories." Trends in Chemistry, 1(3), 282-291.

To the next chapter: Chapter 4: Cloud Labs and Remote Experiments

Back to contents