Chapter 3: Experiencing PI with Python - Process Optimization Practice

Learn the workflow of combining predictive models with multi-objective optimization to handle the yield-versus-cost trade-off. Build a practical sense of reducing the number of experiments with Bayesian optimization.

💡 Supplement: "The optimal solution is not necessarily unique." Lay out your options along the Pareto front and align the criteria with the decision-makers.

Learning Objectives

By reading this article, you will acquire the following: - Set up a Python environment and install PI libraries - Preprocess and visualize process data - Predict process characteristics with five or more regression models - Optimize process conditions using Bayesian optimization - Analyze trade-offs among multiple objectives with multi-objective optimization - Troubleshoot errors on your own

1. Environment Setup: Three Options

There are three ways to build a Python environment for chemical process optimization, depending on your situation.

1.1 Option 1: Anaconda (Recommended for Beginners)

Features: - Scientific computing libraries are bundled from the start - Environment management is easy (a GUI is available) - Supports Windows/Mac/Linux

Installation steps:

# 1. Download Anaconda
# Official site: https://www.anaconda.com/download
# Select Python 3.11 or later

# 2. After installation, launch Anaconda Prompt

# 3. Create a virtual environment (a dedicated PI environment)
conda create -n pi-env python=3.11 numpy pandas matplotlib scikit-learn jupyter scipy

# 4. Activate the environment
conda activate pi-env

# 5. Install additional libraries
conda install -c conda-forge lightgbm scikit-optimize pymoo

# 6. Verify the installation
python --version
# Output: Python 3.11.x

Advantages of Anaconda: - ✅ NumPy, SciPy, etc. are included from the start - ✅ Fewer dependency issues - ✅ Can be managed visually with Anaconda Navigator - ❌ Large file size (3GB or more)

1.2 Option 2: venv (Python Standard)

Features: - Standard Python tool (no extra installation required) - Lightweight (install only what you need) - Isolate environments per project

Installation steps:

# 1. Check that Python 3.11 or later is installed
python3 --version
# Output: Python 3.11.x or later is required

# 2. Create a virtual environment
python3 -m venv pi-env

# 3. Activate the environment
# macOS/Linux:
source pi-env/bin/activate

# Windows (PowerShell):
pi-env\Scripts\Activate.ps1

# Windows (Command Prompt):
pi-env\Scripts\activate.bat

# 4. Upgrade pip
pip install --upgrade pip

# 5. Install the required libraries
pip install numpy pandas matplotlib scikit-learn scipy jupyter
pip install lightgbm scikit-optimize pymoo

# 6. Verify the installation
pip list

Advantages of venv: - ✅ Lightweight (tens of MB) - ✅ Standard Python tool (no extra installation required) - ✅ Independent per project - ❌ You must resolve dependencies manually

1.3 Option 3: Google Colab (No Installation)

Features: - Runs entirely in the browser - No installation required (runs in the cloud) - Free GPU/TPU available

How to use:

1. Go to Google Colab: https://colab.research.google.com
2. Create a new notebook
3. Run the following code (the required libraries are already installed)

# On Google Colab, the following are installed from the start
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor

# Libraries that require additional installation
!pip install scikit-optimize pymoo lightgbm

print("Libraries imported successfully!")
print(f"NumPy version: {np.__version__}")
print(f"Pandas version: {pd.__version__}")

Advantages of Google Colab: - ✅ No installation required (start immediately) - ✅ Free GPU available - ✅ Integrates with Google Drive (easy data storage) - ❌ Requires an internet connection - ❌ Sessions reset after 12 hours

1.4 Environment Selection Guide

1.5 Installation Verification and Troubleshooting

Verification command:

# Runnable in all environments
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sklearn
import scipy

print("===== Environment Check =====")
print(f"Python version: {sys.version}")
print(f"NumPy version: {np.__version__}")
print(f"Pandas version: {pd.__version__}")
print(f"Matplotlib version: {plt.matplotlib.__version__}")
print(f"scikit-learn version: {sklearn.__version__}")
print(f"SciPy version: {scipy.__version__}")

# Check PI-specific libraries
try:
    import skopt
    print(f"scikit-optimize version: {skopt.__version__}")
except ImportError:
    print("⚠️ scikit-optimize not installed (pip install scikit-optimize)")

try:
    import pymoo
    print(f"pymoo version: {pymoo.__version__}")
except ImportError:
    print("⚠️ pymoo not installed (pip install pymoo)")

print("\n✅ The basic libraries are installed correctly!")

Common errors and solutions:

2. Preparing and Visualizing Process Data

We simulate real chemical process data and perform preprocessing and visualization.

2.1 Example 1: Generating and Loading Process Data

Overview: Generate simulated data for a chemical reaction process (temperature, pressure, catalyst amount → yield).

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
import time

# Generate simulated data for a chemical reaction process
np.random.seed(42)
n_samples = 200

# Process conditions (input variables)
temperature = np.random.uniform(300, 500, n_samples)  # Temperature [K]
pressure = np.random.uniform(1, 10, n_samples)  # Pressure [bar]
catalyst = np.random.uniform(0.1, 5.0, n_samples)  # Catalyst amount [wt%]

# Yield model (nonlinear relationship + noise)
# Yield = f(temperature, pressure, catalyst amount) + noise
yield_percentage = (
    20  # Base yield
    + 0.15 * temperature  # Temperature effect (positive correlation)
    - 0.0002 * temperature**2  # Quadratic term for temperature (an optimum exists)
    + 5.0 * pressure  # Pressure effect (positive correlation)
    + 3.0 * catalyst  # Catalyst effect (positive correlation)
    - 0.3 * catalyst**2  # Quadratic term for catalyst (excess addition reduces effect)
    + 0.01 * temperature * pressure  # Temperature-pressure interaction
    + np.random.normal(0, 3, n_samples)  # Noise (measurement error)
)

# Store data in a DataFrame
process_data = pd.DataFrame({
    'temperature_K': temperature,
    'pressure_bar': pressure,
    'catalyst_wt%': catalyst,
    'yield_%': yield_percentage
})

print("===== Check Process Data =====")
print(process_data.head(10))
print(f"\nNumber of entries: {len(process_data)}")
print(f"\nBasic statistics:")
print(process_data.describe())

# Save as CSV (real processes typically provide data in this format)
process_data.to_csv('process_data.csv', index=False)
print("\n✅ Saved data to process_data.csv")

Code explanation: 1. Process conditions: temperature (300-500 K), pressure (1-10 bar), catalyst amount (0.1-5.0 wt%) 2. Yield model: nonlinear relationship (quadratic terms, interaction term) + noise 3. Real-data emulation: typical behavior of chemical reactions (an optimum exists, excess addition reduces effect)

2.2 Example 2: Data Visualization (Scatter Plot Matrix)

import seaborn as sns

# Check relationships between variables with a scatter plot matrix
fig = plt.figure(figsize=(12, 10))
sns.pairplot(
    process_data,
    diag_kind='hist',  # Histogram on the diagonal
    plot_kws={'alpha': 0.6, 's': 50},  # Scatter plot settings
    diag_kws={'bins': 20, 'edgecolor': 'black'}  # Histogram settings
)
plt.suptitle('Scatter Plot Matrix of Process Data', y=1.01, fontsize=16)
plt.tight_layout()
plt.show()

print("===== Correlation Matrix =====")
correlation_matrix = process_data.corr()
print(correlation_matrix)

# Visualize correlations with a heatmap
plt.figure(figsize=(8, 6))
sns.heatmap(
    correlation_matrix,
    annot=True,  # Show values
    fmt='.3f',  # 3 decimal places
    cmap='coolwarm',  # Colormap
    center=0,  # Center at 0
    square=True,  # Square cells
    linewidths=1,  # Cell borders
    cbar_kws={'label': 'Correlation Coefficient'}
)
plt.title('Correlation Heatmap of Process Variables', fontsize=14)
plt.tight_layout()
plt.show()

Interpretation points: - Temperature vs yield: curved (an optimal temperature exists) - Pressure vs yield: positive correlation (pressure↑ → yield↑) - Catalyst vs yield: curved (excess addition reduces effect)

2.3 Example 3: Data Preprocessing (Missing Values and Outliers)

# Check and handle missing values
print("===== Check Missing Values =====")
print(process_data.isnull().sum())

# Artificially add missing values (common in real data)
process_data_with_missing = process_data.copy()
missing_indices = np.random.choice(process_data.index, size=10, replace=False)
process_data_with_missing.loc[missing_indices, 'catalyst_wt%'] = np.nan

print(f"\nAfter adding missing values: {process_data_with_missing.isnull().sum()['catalyst_wt%']} entries")

# Impute missing values (with the mean)
process_data_filled = process_data_with_missing.copy()
process_data_filled['catalyst_wt%'].fillna(
    process_data_filled['catalyst_wt%'].mean(),
    inplace=True
)

print("✅ Imputed missing values with the mean")

# Outlier detection (Z-score method)
from scipy import stats

z_scores = np.abs(stats.zscore(process_data[['yield_%']]))
outliers = (z_scores > 3).any(axis=1)

print(f"\n===== Outlier Detection =====")
print(f"Number of outliers: {outliers.sum()} entries")
print(f"Outlier ratio: {outliers.sum() / len(process_data) * 100:.1f}%")

# Remove outliers
process_data_clean = process_data[~outliers].copy()
print(f"Data count after removal: {len(process_data_clean)} entries")

2.4 Example 4: Feature Engineering (Interaction Terms)

# Create process-specific features
process_data_enhanced = process_data_clean.copy()

# Interaction terms (temperature × pressure, temperature × catalyst)
process_data_enhanced['temp_pressure'] = (
    process_data_enhanced['temperature_K'] * process_data_enhanced['pressure_bar']
)
process_data_enhanced['temp_catalyst'] = (
    process_data_enhanced['temperature_K'] * process_data_enhanced['catalyst_wt%']
)

# Quadratic terms (temperature², catalyst²)
process_data_enhanced['temp_squared'] = process_data_enhanced['temperature_K'] ** 2
process_data_enhanced['catalyst_squared'] = process_data_enhanced['catalyst_wt%'] ** 2

# Ratio (catalyst / pressure)
process_data_enhanced['catalyst_pressure_ratio'] = (
    process_data_enhanced['catalyst_wt%'] / (process_data_enhanced['pressure_bar'] + 1e-10)
)

print("===== Data After Feature Engineering =====")
print(process_data_enhanced.head())
print(f"\nNumber of features: {len(process_data_enhanced.columns) - 1} (original 3 → 8)")

2.5 Example 5: Data Splitting (Training/Test)

# Separate features and target
X = process_data_enhanced.drop('yield_%', axis=1)  # Input: process conditions
y = process_data_enhanced['yield_%']  # Output: yield

# Split into training and test data (80% vs 20%)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

print("===== Data Split =====")
print(f"Training data: {len(X_train)} entries")
print(f"Test data: {len(X_test)} entries")
print(f"\nFeature column names:")
print(list(X.columns))

2.6 Example 6: Data Standardization

from sklearn.preprocessing import StandardScaler

# Create scaler (transform to mean 0, std 1)
scaler = StandardScaler()

# Learn scaling parameters on training data
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)  # Transform test data with the same parameters

print("===== Standardization =====")
print("Training data (before standardization):")
print(f"  Mean temperature: {X_train['temperature_K'].mean():.1f} K")
print(f"  Std temperature: {X_train['temperature_K'].std():.1f} K")

print("\nTraining data (after standardization):")
print(f"  Mean temperature: {X_train_scaled[:, 0].mean():.3f}")
print(f"  Std temperature: {X_train_scaled[:, 0].std():.3f}")

print("\n✅ Standardization made all features have mean 0 and std 1")

2.7 Example 7: Visualizing Time-Series Process Data

# Generate time-series process data (batch process example)
np.random.seed(42)
time_hours = np.arange(0, 24, 0.5)  # 24 hours, 0.5-hour intervals
n_points = len(time_hours)

# Time-series variation of process variables
temp_time = 350 + 50 * np.sin(2 * np.pi * time_hours / 24) + np.random.normal(0, 2, n_points)
pressure_time = 5 + 2 * np.sin(2 * np.pi * time_hours / 12 + np.pi/4) + np.random.normal(0, 0.3, n_points)
yield_time = 60 + 10 * np.sin(2 * np.pi * time_hours / 24 - np.pi/2) + np.random.normal(0, 1.5, n_points)

# Visualization
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

axes[0].plot(time_hours, temp_time, 'r-', linewidth=2, label='Temperature')
axes[0].set_ylabel('Temperature [K]', fontsize=12)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(time_hours, pressure_time, 'b-', linewidth=2, label='Pressure')
axes[1].set_ylabel('Pressure [bar]', fontsize=12)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].plot(time_hours, yield_time, 'g-', linewidth=2, label='Yield')
axes[2].set_xlabel('Time [h]', fontsize=12)
axes[2].set_ylabel('Yield [%]', fontsize=12)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.suptitle('Time-Series Data of a Batch Process', fontsize=16)
plt.tight_layout()
plt.show()

3. Yield Prediction with Regression Models

We implement five machine learning models to predict yield from process conditions.

3.1 Example 8: Linear Regression (Baseline)

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

# Build linear regression model
start_time = time.time()
model_lr = LinearRegression()
model_lr.fit(X_train, y_train)
training_time_lr = time.time() - start_time

# Prediction
y_pred_lr = model_lr.predict(X_test)

# Evaluation
mae_lr = mean_absolute_error(y_test, y_pred_lr)
rmse_lr = np.sqrt(mean_squared_error(y_test, y_pred_lr))
r2_lr = r2_score(y_test, y_pred_lr)

print("===== Linear Regression Model Performance =====")
print(f"Training time: {training_time_lr:.4f} s")
print(f"Mean Absolute Error (MAE): {mae_lr:.2f} %")
print(f"Root Mean Squared Error (RMSE): {rmse_lr:.2f} %")
print(f"Coefficient of Determination (R²): {r2_lr:.4f}")

# Display learned coefficients
print("\n===== Learned Coefficients (Top 3) =====")
coefficients = pd.DataFrame({
    'Feature': X.columns,
    'Coefficient': model_lr.coef_
}).sort_values('Coefficient', key=abs, ascending=False)
print(coefficients.head(3))

3.2 Example 9: Random Forest Regression

from sklearn.ensemble import RandomForestRegressor

# Build random forest model
start_time = time.time()
model_rf = RandomForestRegressor(
    n_estimators=100,
    max_depth=15,
    min_samples_split=5,
    min_samples_leaf=2,
    random_state=42,
    n_jobs=-1
)
model_rf.fit(X_train, y_train)
training_time_rf = time.time() - start_time

# Prediction and evaluation
y_pred_rf = model_rf.predict(X_test)
mae_rf = mean_absolute_error(y_test, y_pred_rf)
rmse_rf = np.sqrt(mean_squared_error(y_test, y_pred_rf))
r2_rf = r2_score(y_test, y_pred_rf)

print("\n===== Random Forest Performance =====")
print(f"Training time: {training_time_rf:.4f} s")
print(f"Mean Absolute Error (MAE): {mae_rf:.2f} %")
print(f"Root Mean Squared Error (RMSE): {rmse_rf:.2f} %")
print(f"Coefficient of Determination (R²): {r2_rf:.4f}")

# Feature importance
feature_importance = pd.DataFrame({
    'Feature': X.columns,
    'Importance': model_rf.feature_importances_
}).sort_values('Importance', ascending=False)

print("\n===== Feature Importance (Top 3) =====")
print(feature_importance.head(3))

3.3 Example 10: LightGBM Regression

import lightgbm as lgb

# Build LightGBM model
start_time = time.time()
model_lgb = lgb.LGBMRegressor(
    n_estimators=100,
    learning_rate=0.1,
    max_depth=10,
    num_leaves=31,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42,
    verbose=-1
)
model_lgb.fit(
    X_train, y_train,
    eval_set=[(X_test, y_test)],
    eval_metric='rmse',
    callbacks=[lgb.early_stopping(stopping_rounds=10, verbose=False)]
)
training_time_lgb = time.time() - start_time

# Prediction and evaluation
y_pred_lgb = model_lgb.predict(X_test)
mae_lgb = mean_absolute_error(y_test, y_pred_lgb)
rmse_lgb = np.sqrt(mean_squared_error(y_test, y_pred_lgb))
r2_lgb = r2_score(y_test, y_pred_lgb)

print("\n===== LightGBM Performance =====")
print(f"Training time: {training_time_lgb:.4f} s")
print(f"Mean Absolute Error (MAE): {mae_lgb:.2f} %")
print(f"Root Mean Squared Error (RMSE): {rmse_lgb:.2f} %")
print(f"Coefficient of Determination (R²): {r2_lgb:.4f}")

3.4 Example 11: Support Vector Regression (SVR)

from sklearn.svm import SVR

# Build SVR model (using standardized data)
start_time = time.time()
model_svr = SVR(
    kernel='rbf',
    C=100,
    gamma='scale',
    epsilon=0.1
)
model_svr.fit(X_train_scaled, y_train)
training_time_svr = time.time() - start_time

# Prediction and evaluation
y_pred_svr = model_svr.predict(X_test_scaled)
mae_svr = mean_absolute_error(y_test, y_pred_svr)
rmse_svr = np.sqrt(mean_squared_error(y_test, y_pred_svr))
r2_svr = r2_score(y_test, y_pred_svr)

print("\n===== SVR Performance =====")
print(f"Training time: {training_time_svr:.4f} s")
print(f"Mean Absolute Error (MAE): {mae_svr:.2f} %")
print(f"Root Mean Squared Error (RMSE): {rmse_svr:.2f} %")
print(f"Coefficient of Determination (R²): {r2_svr:.4f}")

3.5 Example 12: Neural Network (MLP)

from sklearn.neural_network import MLPRegressor

# Build MLP model
start_time = time.time()
model_mlp = MLPRegressor(
    hidden_layer_sizes=(64, 32, 16),
    activation='relu',
    solver='adam',
    alpha=0.001,
    learning_rate_init=0.01,
    max_iter=500,
    random_state=42,
    early_stopping=True,
    validation_fraction=0.2,
    verbose=False
)
model_mlp.fit(X_train_scaled, y_train)
training_time_mlp = time.time() - start_time

# Prediction and evaluation
y_pred_mlp = model_mlp.predict(X_test_scaled)
mae_mlp = mean_absolute_error(y_test, y_pred_mlp)
rmse_mlp = np.sqrt(mean_squared_error(y_test, y_pred_mlp))
r2_mlp = r2_score(y_test, y_pred_mlp)

print("\n===== MLP Performance =====")
print(f"Training time: {training_time_mlp:.4f} s")
print(f"Mean Absolute Error (MAE): {mae_mlp:.2f} %")
print(f"Root Mean Squared Error (RMSE): {rmse_mlp:.2f} %")
print(f"Coefficient of Determination (R²): {r2_mlp:.4f}")
print(f"Number of iterations: {model_mlp.n_iter_}")

3.6 Example 13: Comparing Model Performance

# Model performance comparison table
comparison = pd.DataFrame({
    'Model': ['Linear Regression', 'Random Forest', 'LightGBM', 'SVR', 'MLP'],
    'MAE (%)': [mae_lr, mae_rf, mae_lgb, mae_svr, mae_mlp],
    'RMSE (%)': [rmse_lr, rmse_rf, rmse_lgb, rmse_svr, rmse_mlp],
    'R²': [r2_lr, r2_rf, r2_lgb, r2_svr, r2_mlp],
    'Training Time (s)': [training_time_lr, training_time_rf, training_time_lgb,
                  training_time_svr, training_time_mlp]
})

print("\n===== Overall Model Performance Comparison =====")
print(comparison.to_string(index=False))

# Visualization
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# MAE comparison
axes[0].bar(comparison['Model'], comparison['MAE (%)'],
            color=['blue', 'green', 'orange', 'purple', 'red'])
axes[0].set_ylabel('MAE (%)', fontsize=12)
axes[0].set_title('Mean Absolute Error (lower is better)', fontsize=14)
axes[0].tick_params(axis='x', rotation=45)
axes[0].grid(True, alpha=0.3, axis='y')

# R² comparison
axes[1].bar(comparison['Model'], comparison['R²'],
            color=['blue', 'green', 'orange', 'purple', 'red'])
axes[1].set_ylabel('R²', fontsize=12)
axes[1].set_title('Coefficient of Determination (closer to 1 is better)', fontsize=14)
axes[1].tick_params(axis='x', rotation=45)
axes[1].grid(True, alpha=0.3, axis='y')

# Training time comparison
axes[2].bar(comparison['Model'], comparison['Training Time (s)'],
            color=['blue', 'green', 'orange', 'purple', 'red'])
axes[2].set_ylabel('Training Time (s)', fontsize=12)
axes[2].set_title('Training Time (shorter is better)', fontsize=14)
axes[2].tick_params(axis='x', rotation=45)
axes[2].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

3.7 Example 14: Predicted vs Actual Plot

# Visualize prediction results of the best model (LightGBM)
plt.figure(figsize=(10, 8))
plt.scatter(y_test, y_pred_lgb, alpha=0.6, s=100, c='green', edgecolors='k', linewidth=0.5)
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()],
         'r--', lw=2, label='Perfect Prediction')
plt.xlabel('Actual Yield (%)', fontsize=14)
plt.ylabel('Predicted Yield (%)', fontsize=14)
plt.title('LightGBM: Yield Prediction Accuracy', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)

# Add performance metrics as text
textstr = f'R² = {r2_lgb:.3f}\nMAE = {mae_lgb:.2f} %\nRMSE = {rmse_lgb:.2f} %'
plt.text(0.05, 0.95, textstr, transform=plt.gca().transAxes,
         fontsize=12, verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.tight_layout()
plt.show()

4. Process Optimization

4.1 Example 15: Searching for Optimal Conditions via Grid Search

from scipy.optimize import minimize

# Objective function: maximize yield (minimize negative yield)
def objective_yield(params):
    """
    params = [temperature, pressure, catalyst]
    """
    # Unpack parameters
    temp = params[0]
    press = params[1]
    cat = params[2]

    # Build features (same order as during training)
    features = np.array([[
        temp, press, cat,
        temp * press,  # temp_pressure
        temp * cat,    # temp_catalyst
        temp**2,       # temp_squared
        cat**2,        # catalyst_squared
        cat / (press + 1e-10)  # catalyst_pressure_ratio
    ]])

    # Predict with the model (using LightGBM)
    predicted_yield = model_lgb.predict(features)[0]

    # Return negative value for maximization
    return -predicted_yield

# Constraints (process operating ranges)
bounds = [
    (300, 500),  # Temperature [K]
    (1, 10),     # Pressure [bar]
    (0.1, 5.0)   # Catalyst amount [wt%]
]

# Initial guess
x0 = [400, 5, 2.5]

# Run optimization
result = minimize(
    objective_yield,
    x0,
    method='L-BFGS-B',
    bounds=bounds
)

print("===== Grid Search Optimization Result =====")
print(f"Optimal conditions:")
print(f"  Temperature: {result.x[0]:.1f} K")
print(f"  Pressure: {result.x[1]:.2f} bar")
print(f"  Catalyst: {result.x[2]:.2f} wt%")
print(f"\nMaximum predicted yield: {-result.fun:.2f} %")
print(f"Optimization success: {result.success}")
print(f"Number of iterations: {result.nit}")

4.2 Example 16: Bayesian Optimization (Efficient Search)

from skopt import gp_minimize
from skopt.space import Real
from skopt.utils import use_named_args

# Define search space
space = [
    Real(300, 500, name='temperature'),
    Real(1, 10, name='pressure'),
    Real(0.1, 5.0, name='catalyst')
]

# Objective function (for Bayesian optimization)
@use_named_args(space)
def objective_bayes(**params):
    temp = params['temperature']
    press = params['pressure']
    cat = params['catalyst']

    # Build features
    features = np.array([[
        temp, press, cat,
        temp * press,
        temp * cat,
        temp**2,
        cat**2,
        cat / (press + 1e-10)
    ]])

    # Predicted yield (negative for maximization)
    predicted_yield = model_lgb.predict(features)[0]
    return -predicted_yield

# Run Bayesian optimization
result_bayes = gp_minimize(
    objective_bayes,
    space,
    n_calls=30,  # 30 evaluations
    random_state=42,
    verbose=False
)

print("\n===== Bayesian Optimization Result =====")
print(f"Optimal conditions:")
print(f"  Temperature: {result_bayes.x[0]:.1f} K")
print(f"  Pressure: {result_bayes.x[1]:.2f} bar")
print(f"  Catalyst: {result_bayes.x[2]:.2f} wt%")
print(f"\nMaximum predicted yield: {-result_bayes.fun:.2f} %")

# Optimization convergence history
plt.figure(figsize=(10, 6))
plt.plot(range(1, len(result_bayes.func_vals) + 1),
         -result_bayes.func_vals, 'b-o', linewidth=2, markersize=6)
plt.axhline(y=-result_bayes.fun, color='r', linestyle='--',
            label=f'Best value: {-result_bayes.fun:.2f}%')
plt.xlabel('Number of Evaluations', fontsize=12)
plt.ylabel('Predicted Yield (%)', fontsize=12)
plt.title('Bayesian Optimization Convergence History', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

4.3 Example 17: Design of Experiments (DoE)

from itertools import product

# Two-level design of experiments (2^3 = 8 experiments)
levels = {
    'temperature': [350, 450],  # Low level, high level
    'pressure': [3, 8],
    'catalyst': [1.0, 4.0]
}

# Generate all combinations
experiments = list(product(levels['temperature'], levels['pressure'], levels['catalyst']))

# Compute the predicted yield for each experiment
results_doe = []
for temp, press, cat in experiments:
    features = np.array([[
        temp, press, cat,
        temp * press,
        temp * cat,
        temp**2,
        cat**2,
        cat / (press + 1e-10)
    ]])
    predicted_yield = model_lgb.predict(features)[0]
    results_doe.append({
        'Temperature [K]': temp,
        'Pressure [bar]': press,
        'Catalyst [wt%]': cat,
        'Predicted Yield [%]': predicted_yield
    })

# Convert results to a DataFrame
df_doe = pd.DataFrame(results_doe).sort_values('Predicted Yield [%]', ascending=False)

print("\n===== Design of Experiments (2^3 DoE) Result =====")
print(df_doe.to_string(index=False))

# Extract the best condition
best_condition = df_doe.iloc[0]
print(f"\nBest condition:")
print(f"  Temperature: {best_condition['Temperature [K]']:.0f} K")
print(f"  Pressure: {best_condition['Pressure [bar]']:.0f} bar")
print(f"  Catalyst: {best_condition['Catalyst [wt%]']:.1f} wt%")
print(f"  Predicted yield: {best_condition['Predicted Yield [%]']:.2f} %")

4.4 Example 18: Response Surface Method

from scipy.interpolate import griddata

# Create a grid over temperature and pressure ranges (catalyst fixed)
temp_range = np.linspace(300, 500, 50)
press_range = np.linspace(1, 10, 50)
temp_grid, press_grid = np.meshgrid(temp_range, press_range)

# Predict yield at each grid point (catalyst fixed at optimal value)
catalyst_fixed = result_bayes.x[2]
yield_grid = np.zeros_like(temp_grid)

for i in range(len(temp_range)):
    for j in range(len(press_range)):
        temp = temp_grid[j, i]
        press = press_grid[j, i]
        cat = catalyst_fixed

        features = np.array([[
            temp, press, cat,
            temp * press,
            temp * cat,
            temp**2,
            cat**2,
            cat / (press + 1e-10)
        ]])

        yield_grid[j, i] = model_lgb.predict(features)[0]

# Visualize the response surface
fig = plt.figure(figsize=(14, 6))

# Contour plot
ax1 = fig.add_subplot(1, 2, 1)
contour = ax1.contourf(temp_grid, press_grid, yield_grid, levels=20, cmap='viridis')
ax1.set_xlabel('Temperature [K]', fontsize=12)
ax1.set_ylabel('Pressure [bar]', fontsize=12)
ax1.set_title(f'Response Surface (Catalyst = {catalyst_fixed:.2f} wt%)', fontsize=14)
plt.colorbar(contour, ax=ax1, label='Predicted Yield [%]')

# 3D surface
ax2 = fig.add_subplot(1, 2, 2, projection='3d')
surf = ax2.plot_surface(temp_grid, press_grid, yield_grid,
                        cmap='viridis', alpha=0.8)
ax2.set_xlabel('Temperature [K]', fontsize=10)
ax2.set_ylabel('Pressure [bar]', fontsize=10)
ax2.set_zlabel('Predicted Yield [%]', fontsize=10)
ax2.set_title('3D Response Surface', fontsize=14)
plt.colorbar(surf, ax=ax2, label='Predicted Yield [%]', shrink=0.5)

plt.tight_layout()
plt.show()

4.5 Example 19: Constrained Optimization

from scipy.optimize import NonlinearConstraint

# Objective function (maximize yield)
def objective_constrained(params):
    temp, press, cat = params
    features = np.array([[
        temp, press, cat,
        temp * press,
        temp * cat,
        temp**2,
        cat**2,
        cat / (press + 1e-10)
    ]])
    predicted_yield = model_lgb.predict(features)[0]
    return -predicted_yield

# Constraint function: energy cost < 100 [arbitrary units]
# Cost = 0.1 * temperature + 2.0 * pressure
def energy_cost_constraint(params):
    temp, press, cat = params
    cost = 0.1 * temp + 2.0 * press
    return cost

# Constraint: energy cost <= 100
constraint = NonlinearConstraint(energy_cost_constraint, -np.inf, 100)

# Run optimization
result_constrained = minimize(
    objective_constrained,
    x0=[400, 5, 2.5],
    method='SLSQP',
    bounds=bounds,
    constraints=constraint
)

print("\n===== Constrained Optimization Result =====")
print(f"Optimal conditions (under energy cost constraint):")
print(f"  Temperature: {result_constrained.x[0]:.1f} K")
print(f"  Pressure: {result_constrained.x[1]:.2f} bar")
print(f"  Catalyst: {result_constrained.x[2]:.2f} wt%")
print(f"\nMaximum predicted yield: {-result_constrained.fun:.2f} %")
print(f"Energy cost: {energy_cost_constraint(result_constrained.x):.2f}")
print(f"Constraint satisfied: {energy_cost_constraint(result_constrained.x) <= 100}")

4.6 Example 20: Multi-Objective Optimization (Yield vs Cost)

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.core.problem import Problem
from pymoo.optimize import minimize as pymoo_minimize

# Define the multi-objective optimization problem
class ProcessOptimizationProblem(Problem):
    def __init__(self):
        super().__init__(
            n_var=3,  # Number of variables (temperature, pressure, catalyst)
            n_obj=2,  # Number of objectives (yield, cost)
            xl=np.array([300, 1, 0.1]),  # Lower bounds
            xu=np.array([500, 10, 5.0])  # Upper bounds
        )

    def _evaluate(self, X, out, *args, **kwargs):
        # X: array of shape (n_samples, 3)
        n_samples = X.shape[0]
        f1 = np.zeros(n_samples)  # Objective 1: -yield (minimize)
        f2 = np.zeros(n_samples)  # Objective 2: cost (minimize)

        for i in range(n_samples):
            temp, press, cat = X[i]

            # Predict yield
            features = np.array([[
                temp, press, cat,
                temp * press,
                temp * cat,
                temp**2,
                cat**2,
                cat / (press + 1e-10)
            ]])
            predicted_yield = model_lgb.predict(features)[0]

            # Objective 1: maximize yield → minimize -yield
            f1[i] = -predicted_yield

            # Objective 2: minimize cost
            # Cost = energy cost + catalyst cost
            energy_cost = 0.1 * temp + 2.0 * press
            catalyst_cost = 5.0 * cat
            f2[i] = energy_cost + catalyst_cost

        out["F"] = np.column_stack([f1, f2])

# Optimize with the NSGA-II algorithm
problem = ProcessOptimizationProblem()
algorithm = NSGA2(pop_size=50)

result_nsga2 = pymoo_minimize(
    problem,
    algorithm,
    ('n_gen', 100),  # 100 generations
    verbose=False
)

# Obtain the Pareto-optimal solutions
pareto_front = result_nsga2.F
pareto_solutions = result_nsga2.X

print("\n===== Multi-Objective Optimization (NSGA-II) Result =====")
print(f"Number of Pareto-optimal solutions: {len(pareto_solutions)}")
print(f"\nExamples of Pareto-optimal solutions (first 3):")
for i in range(min(3, len(pareto_solutions))):
    temp, press, cat = pareto_solutions[i]
    yield_val = -pareto_front[i, 0]
    cost_val = pareto_front[i, 1]
    print(f"\nSolution {i+1}:")
    print(f"  Temperature: {temp:.1f} K, Pressure: {press:.2f} bar, Catalyst: {cat:.2f} wt%")
    print(f"  Yield: {yield_val:.2f} %, Cost: {cost_val:.2f}")

# Visualize the Pareto front
plt.figure(figsize=(10, 6))
plt.scatter(-pareto_front[:, 0], pareto_front[:, 1],
            c='blue', s=50, alpha=0.6, edgecolors='k', linewidth=0.5)
plt.xlabel('Yield [%]', fontsize=12)
plt.ylabel('Cost [arbitrary units]', fontsize=12)
plt.title('Pareto Front (Yield vs Cost)', fontsize=14)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

4.7 Example 21: Comparing Optimization Results

# Compare the results of all optimization methods
optimization_results = pd.DataFrame({
    'Method': [
        'Grid Search',
        'Bayesian Optimization',
        'DoE (2^3)',
        'Constrained Optimization'
    ],
    'Temperature [K]': [
        result.x[0],
        result_bayes.x[0],
        best_condition['Temperature [K]'],
        result_constrained.x[0]
    ],
    'Pressure [bar]': [
        result.x[1],
        result_bayes.x[1],
        best_condition['Pressure [bar]'],
        result_constrained.x[1]
    ],
    'Catalyst [wt%]': [
        result.x[2],
        result_bayes.x[2],
        best_condition['Catalyst [wt%]'],
        result_constrained.x[2]
    ],
    'Predicted Yield [%]': [
        -result.fun,
        -result_bayes.fun,
        best_condition['Predicted Yield [%]'],
        -result_constrained.fun
    ]
})

print("\n===== Comparison of Optimization Methods =====")
print(optimization_results.to_string(index=False))

4.8 Example 22: Flowchart of Optimization Methods

# Mermaid flowchart (displayed in Markdown)
print("""
<div class="mermaid">
graph TD
    A[Process Optimization Task] --> B{Number of objectives?}
    B -->|Single objective| C{Constraints?}
    B -->|Multiple objectives| D[NSGA-II / Genetic Algorithm]

    C -->|None| E{Evaluation cost?}
    C -->|Yes| F[Constrained Optimization<br>SLSQP/COBYLA]

    E -->|Low| G[Grid Search<br>or DoE]
    E -->|High| H[Bayesian Optimization<br>Gaussian Process]

    D --> I[Obtain Pareto Front]
    F --> J[Obtain Optimal Conditions]
    G --> J
    H --> J
    I --> K[Trade-off Analysis]
    J --> L[Experimental Validation]
    K --> L

    style A fill:#e3f2fd
    style D fill:#c8e6c9
    style F fill:#fff9c4
    style H fill:#ffccbc
    style I fill:#f3e5f5
</div>
""")

5. Advanced Methods

5.9 Example 23: Hyperparameter Tuning (Grid Search)

from sklearn.model_selection import GridSearchCV

# Hyperparameter candidates for random forest
param_grid = {
    'n_estimators': [50, 100, 200],
    'max_depth': [10, 15, 20],
    'min_samples_split': [2, 5, 10]
}

# Grid Search settings
grid_search = GridSearchCV(
    estimator=RandomForestRegressor(random_state=42),
    param_grid=param_grid,
    cv=5,  # 5-fold cross-validation
    scoring='neg_mean_absolute_error',
    n_jobs=-1,
    verbose=1
)

# Run Grid Search
print("===== Grid Search Start =====")
grid_search.fit(X_train, y_train)

print(f"\n===== Best Hyperparameters =====")
for param, value in grid_search.best_params_.items():
    print(f"{param}: {value}")

print(f"\nCross-validation MAE: {-grid_search.best_score_:.2f} %")

# Evaluate the best model on test data
best_model_gs = grid_search.best_estimator_
y_pred_gs = best_model_gs.predict(X_test)
mae_gs = mean_absolute_error(y_test, y_pred_gs)
r2_gs = r2_score(y_test, y_pred_gs)

print(f"\nPerformance on test data:")
print(f"  MAE: {mae_gs:.2f} %")
print(f"  R²: {r2_gs:.4f}")

5.10 Example 24: Anomaly Detection in Time-Series Processes

from sklearn.ensemble import IsolationForest

# Time-series process data (using the data generated in Example 7)
process_time_series = pd.DataFrame({
    'time_h': time_hours,
    'temperature_K': temp_time,
    'pressure_bar': pressure_time,
    'yield_%': yield_time
})

# Anomaly detection with Isolation Forest
iso_forest = IsolationForest(
    contamination=0.1,  # Assume 10% anomaly ratio
    random_state=42
)

# Features (temperature, pressure, yield)
X_anomaly = process_time_series[['temperature_K', 'pressure_bar', 'yield_%']]

# Compute anomaly scores
anomaly_scores = iso_forest.fit_predict(X_anomaly)
process_time_series['anomaly'] = anomaly_scores

# Extract anomalous data
anomalies = process_time_series[process_time_series['anomaly'] == -1]

print(f"\n===== Anomaly Detection Result =====")
print(f"Number of anomalies: {len(anomalies)} / {len(process_time_series)}")
print(f"Anomaly ratio: {len(anomalies) / len(process_time_series) * 100:.1f}%")

# Visualization
plt.figure(figsize=(14, 5))
plt.plot(process_time_series['time_h'], process_time_series['yield_%'],
         'b-', linewidth=1.5, label='Normal data')
plt.scatter(anomalies['time_h'], anomalies['yield_%'],
            c='red', s=100, marker='x', linewidth=2, label='Anomalous data')
plt.xlabel('Time [h]', fontsize=12)
plt.ylabel('Yield [%]', fontsize=12)
plt.title('Anomaly Detection in Time-Series Processes', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

5.11 Example 25: Interpretability Analysis with SHAP Values

# Analyze feature impact with SHAP (SHapley Additive exPlanations)
try:
    import shap

    # Create a SHAP Explainer (for LightGBM)
    explainer = shap.TreeExplainer(model_lgb)
    shap_values = explainer.shap_values(X_test)

    # SHAP value summary plot
    plt.figure(figsize=(10, 6))
    shap.summary_plot(shap_values, X_test, feature_names=X.columns, show=False)
    plt.title('SHAP Value Summary Plot (Feature Impact)', fontsize=14)
    plt.tight_layout()
    plt.show()

    print("\n✅ SHAP value analysis complete")
    print("Visualized how each feature affects predictions.")

except ImportError:
    print("\n⚠️ SHAP library not installed")
    print("Please install it with pip install shap.")

5.12 Example 26: Process Simulation (PID Control)

# Simulation of a simple PID controller
class PIDController:
    def __init__(self, Kp, Ki, Kd, setpoint):
        self.Kp = Kp  # Proportional gain
        self.Ki = Ki  # Integral gain
        self.Kd = Kd  # Derivative gain
        self.setpoint = setpoint  # Setpoint
        self.integral = 0
        self.prev_error = 0

    def update(self, measured_value, dt):
        # Compute error
        error = self.setpoint - measured_value

        # Integral term
        self.integral += error * dt

        # Derivative term
        derivative = (error - self.prev_error) / dt

        # PID output
        output = (
            self.Kp * error +
            self.Ki * self.integral +
            self.Kd * derivative
        )

        # Save error for the next step
        self.prev_error = error

        return output

# Process model (first-order lag system)
def process_model(input_val, current_temp, tau=5.0, K=1.0, dt=0.1):
    """
    First-order lag process model
    tau: time constant, K: gain
    """
    dT = (K * input_val - current_temp) / tau
    new_temp = current_temp + dT * dt
    return new_temp

# Simulation settings
dt = 0.1  # Time step [s]
t_end = 50  # Simulation duration [s]
time_sim = np.arange(0, t_end, dt)

# Initialize the PID controller (target temperature: 400 K)
pid = PIDController(Kp=2.0, Ki=0.5, Kd=1.0, setpoint=400)

# Run simulation
temperature = 350  # Initial temperature [K]
temperatures = []
inputs = []

for t in time_sim:
    # Compute the PID control input
    control_input = pid.update(temperature, dt)

    # Update temperature with the process model
    temperature = process_model(control_input, temperature, dt=dt)

    # Record
    temperatures.append(temperature)
    inputs.append(control_input)

# Visualization
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# Temperature evolution
axes[0].plot(time_sim, temperatures, 'b-', linewidth=2, label='Process temperature')
axes[0].axhline(y=400, color='r', linestyle='--', linewidth=1.5, label='Target temperature')
axes[0].set_ylabel('Temperature [K]', fontsize=12)
axes[0].set_title('PID Control Simulation', fontsize=14)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Control input evolution
axes[1].plot(time_sim, inputs, 'g-', linewidth=2, label='Control input')
axes[1].set_xlabel('Time [s]', fontsize=12)
axes[1].set_ylabel('Control input', fontsize=12)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\n===== PID Control Simulation Result =====")
print(f"Final temperature: {temperatures[-1]:.2f} K (target: 400 K)")
print(f"Steady-state error: {abs(400 - temperatures[-1]):.2f} K")

6. Troubleshooting Guide

6.1 List of Common Errors (Extended)

6.2 Debugging Checklist (Extended)

Step 1: Check the data

# Basic statistics of the data
print("===== Data Statistics =====")
print(process_data.describe())

# Check missing values
print("\n===== Missing Values =====")
print(process_data.isnull().sum())

# Check data types
print("\n===== Data Types =====")
print(process_data.dtypes)

# Check data shape
print(f"\nData shape: {process_data.shape}")
print(f"Rows: {len(process_data)}, Columns: {len(process_data.columns)}")

# Check memory usage
print(f"\nMemory usage: {process_data.memory_usage(deep=True).sum() / 1024**2:.2f} MB")

Step 2: Simplify the model

# If a complex model fails, first try linear regression
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_absolute_error

model_simple = LinearRegression()
model_simple.fit(X_train, y_train)

y_pred_simple = model_simple.predict(X_test)
r2_simple = r2_score(y_test, y_pred_simple)
mae_simple = mean_absolute_error(y_test, y_pred_simple)

print(f"Linear regression R²: {r2_simple:.4f}")
print(f"Linear regression MAE: {mae_simple:.2f}")

# If R² < 0.5, there may be a data quality issue
if r2_simple < 0.5:
    print("⚠️ Poor performance even with linear regression → check data quality")

Step 3: Check scaling

# Standardization is essential for SVR and MLP
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Compare statistics before and after scaling
print("Before scaling:")
print(f"  Mean: {X_train.mean().values[:3]}")  # First 3 features
print(f"  Std: {X_train.std().values[:3]}")

print("After scaling:")
print(f"  Mean: {X_train_scaled.mean(axis=0)[:3]}")  # ≈ 0
print(f"  Std: {X_train_scaled.std(axis=0)[:3]}")  # ≈ 1

Step 4: Overfitting check

# Check the performance gap between training and test data
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)

r2_train = r2_score(y_train, y_train_pred)
r2_test = r2_score(y_test, y_test_pred)

print(f"Training R²: {r2_train:.4f}")
print(f"Test R²: {r2_test:.4f}")

# If the train-test gap is 0.1 or more, overfitting is possible
if r2_train - r2_test > 0.1:
    print("⚠️ Possible overfitting → add regularization")

Step 5: Check feature importance

# Visualize feature importance with random forest
import matplotlib.pyplot as plt

if hasattr(model, 'feature_importances_'):
    importance = model.feature_importances_
    feature_names = X.columns

    plt.figure(figsize=(10, 6))
    plt.barh(feature_names, importance)
    plt.xlabel('Importance')
    plt.title('Feature Importance')
    plt.tight_layout()
    plt.show()

    # Consider removing low-importance features (< 0.01)
    low_importance = feature_names[importance < 0.01]
    if len(low_importance) > 0:
        print(f"⚠️ Low-importance features: {list(low_importance)}")

6.3 Practical Pitfalls

Pitfall 1: Convergence issues in Bayesian optimization

Problem: Bayesian optimization gets stuck in a local optimum

# ❌ Wrong: too few evaluations
result = gp_minimize(objective, space, n_calls=10)  # insufficient

# ✅ Correct: sufficient evaluations + initial random exploration
result = gp_minimize(
    objective,
    space,
    n_calls=50,  # Increase the number of evaluations
    n_initial_points=20,  # More initial random exploration
    acq_func='EI',  # Expected Improvement (balance of exploration and exploitation)
    random_state=42
)

Pitfall 2: Data leakage in LSTM time-series prediction

Problem: Using future information to predict the past

# ❌ Wrong: standardize on all data
scaler = MinMaxScaler()
data_scaled = scaler.fit_transform(data)  # Uses all data
train, test = data_scaled[:800], data_scaled[800:]

# ✅ Correct: standardize on training data only
train, test = data[:800], data[800:]
scaler = MinMaxScaler()
train_scaled = scaler.fit_transform(train)  # Training only
test_scaled = scaler.transform(test)  # Transform with the same parameters

Pitfall 3: Selecting a Pareto solution in multi-objective optimization

Problem: Too many Pareto-optimal solutions to choose from

# Select three representative solutions from the Pareto-optimal set
# 1. Prioritize objective 1 (maximum yield)
idx_yield = np.argmin(pareto_front[:, 0])  # minimum f1 (maximum -yield)
# 2. Prioritize objective 2 (minimum cost)
idx_cost = np.argmin(pareto_front[:, 1])
# 3. Balanced (minimum Euclidean distance)
utopia = np.min(pareto_front, axis=0)  # Ideal point
distances = np.linalg.norm(pareto_front - utopia, axis=1)
idx_balance = np.argmin(distances)

print(f"Yield-focused solution: {pareto_solutions[idx_yield]}")
print(f"Cost-focused solution: {pareto_solutions[idx_cost]}")
print(f"Balanced solution: {pareto_solutions[idx_balance]}")

End-of-Chapter Checklist (50 Items)

1. Environment Setup (8 items)

• [ ] 1.1 Can install Python 3.11 or later
• [ ] 1.2 Can explain the differences between Anaconda, venv, and Google Colab
• [ ] 1.3 Can create and activate a virtual environment
• [ ] 1.4 Can install required libraries with pip install
• [ ] 1.5 Can run installation verification code
• [ ] 1.6 Understand the environment selection guide (beginners → Anaconda, etc.)
• [ ] 1.7 Can resolve common installation errors (DLL load failed, etc.)
• [ ] 1.8 Can check and record library versions

2. Data Preprocessing (8 items)

• [ ] 2.1 Can generate and load process data
• [ ] 2.2 Can visualize relationships between variables with a scatter plot matrix
• [ ] 2.3 Can create a correlation coefficient heatmap
• [ ] 2.4 Can check and impute missing values (mean, linear interpolation)
• [ ] 2.5 Can detect (Z-score method) and remove outliers
• [ ] 2.6 Can implement feature engineering (interaction terms, quadratic terms)
• [ ] 2.7 Can correctly perform data splitting (train_test_split)
• [ ] 2.8 Can implement standardization (StandardScaler) without leakage

3. Regression Model Building (10 items)

• [ ] 3.1 Can build and evaluate a linear regression model
• [ ] 3.2 Can interpret random forest feature importance
• [ ] 3.3 Can use early_stopping in LightGBM
• [ ] 3.4 Understand the need for standardization with SVR
• [ ] 3.5 Can tune MLP hyperparameters (hidden_layer_sizes, etc.)
• [ ] 3.6 Can compare the performance of five models (MAE, RMSE, R²)
• [ ] 3.7 Can create and evaluate predicted-vs-actual plots
• [ ] 3.8 Understand model performance guidelines (R²>0.8 good, >0.9 excellent)
• [ ] 3.9 Can detect and address overfitting (gap between training and test R²)
• [ ] 3.10 Can remove variables with low feature importance

4. Process Optimization (12 items)

• [ ] 4.1 Can optimize with grid search (L-BFGS-B)
• [ ] 4.2 Can implement Bayesian optimization (gp_minimize)
• [ ] 4.3 Can visualize Bayesian optimization convergence history
• [ ] 4.4 Can cover conditions with design of experiments (2^3 DoE)
• [ ] 4.5 Can create a two-variable contour plot with the response surface method
• [ ] 4.6 Can visualize a 3D response surface
• [ ] 4.7 Can implement constrained optimization (NonlinearConstraint)
• [ ] 4.8 Can run multi-objective optimization (NSGA-II)
• [ ] 4.9 Can visualize and interpret the Pareto front
• [ ] 4.10 Can select three representative solutions from the Pareto set
• [ ] 4.11 Can draw a flowchart of optimization methods
• [ ] 4.12 Can choose methods based on evaluation cost (high/low)

5. Advanced Methods (6 items)

• [ ] 5.1 Can do hyperparameter tuning with GridSearchCV
• [ ] 5.2 Can perform anomaly detection with Isolation Forest
• [ ] 5.3 Can analyze feature impact with SHAP values (optional)
• [ ] 5.4 Can simulate a PID controller
• [ ] 5.5 Can explain the effects of PID parameters (Kp, Ki, Kd)
• [ ] 5.6 Can visualize time-series process data

6. Troubleshooting (6 items)

• [ ] 6.1 Know how to solve 8 types of common errors
• [ ] 6.2 Can execute the 5-step debugging checklist
• [ ] 6.3 Can solve convergence issues in Bayesian optimization
• [ ] 6.4 Can avoid data leakage in LSTM time-series prediction
• [ ] 6.5 Can select a representative solution from the Pareto set in multi-objective optimization
• [ ] 6.6 Can perform overfitting checks and add regularization

7. Project Challenge: Optimizing a Chemical Reactor

Integrate what you have learned and tackle a practical project.

7.1 Project Overview

Goal: Optimize the operating conditions of a chemical reactor to maximize yield

Target performance: - Predictive model: R² > 0.85 - Optimization: yield > 90%

7.2 Step-by-Step Guide

Step 1: Data generation (more realistic data)

# Generate more complex reactor data
np.random.seed(42)
n_reactor = 300

temp_reactor = np.random.uniform(320, 480, n_reactor)
press_reactor = np.random.uniform(2, 12, n_reactor)
cat_reactor = np.random.uniform(0.5, 6.0, n_reactor)
residence_time = np.random.uniform(5, 30, n_reactor)  # Residence time [min]

# A more complex yield model (4 variables, interactions, with optima)
yield_reactor = (
    25
    + 0.18 * temp_reactor
    - 0.00025 * temp_reactor**2
    + 6.0 * press_reactor
    - 0.3 * press_reactor**2
    + 4.0 * cat_reactor
    - 0.4 * cat_reactor**2
    + 1.5 * residence_time
    - 0.03 * residence_time**2
    + 0.015 * temp_reactor * press_reactor
    + 0.008 * cat_reactor * residence_time
    + np.random.normal(0, 2.5, n_reactor)
)

reactor_data = pd.DataFrame({
    'temperature': temp_reactor,
    'pressure': press_reactor,
    'catalyst': cat_reactor,
    'residence_time': residence_time,
    'yield': yield_reactor
})

print("===== Reactor Data =====")
print(reactor_data.describe())

Step 2: Feature engineering

# Add features
reactor_data['temp_press'] = reactor_data['temperature'] * reactor_data['pressure']
reactor_data['cat_time'] = reactor_data['catalyst'] * reactor_data['residence_time']
reactor_data['temp_sq'] = reactor_data['temperature'] ** 2
reactor_data['press_sq'] = reactor_data['pressure'] ** 2

X_reactor = reactor_data.drop('yield', axis=1)
y_reactor = reactor_data['yield']

X_train_r, X_test_r, y_train_r, y_test_r = train_test_split(
    X_reactor, y_reactor, test_size=0.2, random_state=42
)

Step 3: Model training (LightGBM)

model_reactor = lgb.LGBMRegressor(
    n_estimators=200,
    learning_rate=0.05,
    max_depth=15,
    random_state=42,
    verbose=-1
)
model_reactor.fit(X_train_r, y_train_r)

y_pred_r = model_reactor.predict(X_test_r)
r2_reactor = r2_score(y_test_r, y_pred_r)
mae_reactor = mean_absolute_error(y_test_r, y_pred_r)

print(f"\n===== Reactor Model Performance =====")
print(f"R²: {r2_reactor:.3f}")
print(f"MAE: {mae_reactor:.2f}%")

if r2_reactor > 0.85:
    print("🎉 Goal achieved! (R² > 0.85)")

Step 4: Condition search with Bayesian optimization

# Search for optimal conditions
space_reactor = [
    Real(320, 480, name='temperature'),
    Real(2, 12, name='pressure'),
    Real(0.5, 6.0, name='catalyst'),
    Real(5, 30, name='residence_time')
]

@use_named_args(space_reactor)
def objective_reactor(**params):
    temp = params['temperature']
    press = params['pressure']
    cat = params['catalyst']
    res_time = params['residence_time']

    features = np.array([[
        temp, press, cat, res_time,
        temp * press,
        cat * res_time,
        temp**2,
        press**2
    ]])

    predicted_yield = model_reactor.predict(features)[0]
    return -predicted_yield

result_reactor = gp_minimize(
    objective_reactor,
    space_reactor,
    n_calls=50,
    random_state=42,
    verbose=False
)

print(f"\n===== Optimal Conditions =====")
print(f"Temperature: {result_reactor.x[0]:.1f} K")
print(f"Pressure: {result_reactor.x[1]:.2f} bar")
print(f"Catalyst: {result_reactor.x[2]:.2f} wt%")
print(f"Residence time: {result_reactor.x[3]:.1f} min")
print(f"\nMaximum predicted yield: {-result_reactor.fun:.2f}%")

if -result_reactor.fun > 90:
    print("🎉 Goal achieved! (yield > 90%)")

8. Summary

What You Learned in This Chapter

1. Environment setup - Three options: Anaconda, venv, Google Colab - Installing PI-specific libraries (scikit-optimize, pymoo)

2. Process data handling - Data generation and visualization (scatter plot matrix, heatmap) - Preprocessing (missing value imputation, outlier removal, standardization) - Feature engineering (interaction terms, quadratic terms)

3. Five regression models - Linear regression, random forest, LightGBM, SVR, MLP - Comparison of model performance (MAE, RMSE, R²)

4. Process optimization methods - Grid search, Bayesian optimization, DoE, response surface method - Constrained optimization, multi-objective optimization (NSGA-II)

5. Advanced methods - Hyperparameter tuning - Anomaly detection (Isolation Forest) - Interpretability analysis (SHAP values) - PID control simulation

Next Steps

Having completed this tutorial, you can: - ✅ Preprocess and visualize process data - ✅ Use five or more types of regression models appropriately - ✅ Optimize process conditions with Bayesian optimization - ✅ Analyze trade-offs with multi-objective optimization

What to learn next: 1. Application to real processes - Data acquisition from a DCS (distributed control system) - Real-time optimization

2. Applications of deep learning - LSTM (time-series prediction) - Autoencoders (anomaly detection)

3. Autonomous process control - Control with reinforcement learning - Model Predictive Control (MPC)

Exercises

Problem 1 (Difficulty: easy)

In process optimization, give three reasons why Bayesian optimization is superior to grid search.

Problem 2 (Difficulty: medium)

Explain what the Pareto front obtained from multi-objective optimization (NSGA-II) is, and give one application example in chemical processes.

Problem 3 (Difficulty: hard)

Explain how the PID control parameters (Kp, Ki, Kd) affect the behavior of temperature control, and state the advantages and disadvantages of increasing each one.

# Automatically tune PID parameters with Bayesian optimization
space_pid = [
    Real(0.1, 10.0, name='Kp'),
    Real(0.01, 1.0, name='Ki'),
    Real(0.01, 5.0, name='Kd')
]

@use_named_args(space_pid)
def objective_pid(**params):
    # Run the PID simulation
    # Goal: minimum overshoot + shortest settling time
    overshoot, settling_time = simulate_pid(
        Kp=params['Kp'],
        Ki=params['Ki'],
        Kd=params['Kd']
    )
    return overshoot + 0.1 * settling_time

# Run optimization
result_pid = gp_minimize(objective_pid, space_pid, n_calls=50)

References

1. Pedregosa, F., et al. (2011). "Scikit-learn: Machine Learning in Python." Journal of Machine Learning Research, 12, 2825-2830. URL: https://scikit-learn.org

2. Brochu, E., Cora, V. M., & de Freitas, N. (2010). "A Tutorial on Bayesian Optimization of Expensive Cost Functions." arXiv:1012.2599. URL: https://arxiv.org/abs/1012.2599

3. Deb, K., et al. (2002). "A fast and elitist multiobjective genetic algorithm: NSGA-II." IEEE Transactions on Evolutionary Computation, 6(2), 182-197. DOI: 10.1109/4235.996017

4. 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. DOI: 10.1109/JPROC.2015.2494218

5. Lundberg, S. M., & Lee, S. I. (2017). "A Unified Approach to Interpreting Model Predictions." Advances in Neural Information Processing Systems, 30. URL: https://github.com/slundberg/shap

6. Åström, K. J., & Hägglund, T. (2006). Advanced PID Control. ISA-The Instrumentation, Systems, and Automation Society. ISBN: 978-1556179426

7. scikit-optimize Documentation. (2024). "Bayesian Optimization." URL: https://scikit-optimize.github.io/stable/

8. pymoo Documentation. (2024). "Multi-objective Optimization." URL: https://pymoo.org/

Created: 2025-10-16 Version: 1.0 Series: PI Introduction Series v1.0 Author: MI Knowledge Hub Project