Chapter 4: Continuous Manufacturing and QbD Implementation

📖 Chapter Overview

Pharmaceutical manufacturing is transitioning from traditional batch production to continuous manufacturing. This chapter covers the optimization and implementation of continuous production systems based on QbD (Quality by Design) concepts, utilizing Design of Experiments (DoE), design space, and risk-based approaches.

🎯 Learning Objectives

• Basic concepts of QbD (Quality by Design) and ICH Q8-Q11 guidelines
• Efficient process optimization using Design of Experiments (DoE)
• Design space construction and visualization
• Design and control strategies for continuous manufacturing systems
• Implementation of real-time control and feedback control
• Migration strategy from batch to continuous production
• Application of risk-based approaches (ICH Q9)

🎯 4.1 Fundamentals of QbD (Quality by Design)

Three Pillars of QbD

Core concepts of QbD as defined in the ICH Q8 guideline:

1. QTPP (Quality Target Product Profile): Quality Target Product Profile
2. CQA (Critical Quality Attributes): Critical Quality Attributes
3. CPP (Critical Process Parameters): Critical Process Parameters

Definition of Design Space

💻 Code Example 4.1: Design of Experiments (DoE) and Design Space Construction

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

class QbDDesigner:
    """QbD Design Space Design Class"""

    def __init__(self):
        self.model = None
        self.poly_features = None

    def generate_doe_design(self, design_type='central_composite'):
        """
        Generate experimental design

        Args:
            design_type: Design type ('full_factorial', 'central_composite', 'box_behnken')

        Returns:
            DataFrame of experimental design
        """
        if design_type == 'central_composite':
            # Central Composite Design (3 factors)
            # Factors: Reaction temperature (70-90°C), Reaction time (60-180 min), Catalyst amount (1-5%)

            # Cube points (2^3 = 8 points)
            cube_points = np.array([
                [-1, -1, -1], [1, -1, -1], [-1, 1, -1], [1, 1, -1],
                [-1, -1, 1], [1, -1, 1], [-1, 1, 1], [1, 1, 1]
            ])

            # Axial points (2×3 = 6 points, α=1.682)
            alpha = 1.682
            axial_points = np.array([
                [-alpha, 0, 0], [alpha, 0, 0],
                [0, -alpha, 0], [0, alpha, 0],
                [0, 0, -alpha], [0, 0, alpha]
            ])

            # Center points (6 replicates)
            center_points = np.array([[0, 0, 0]] * 6)

            # Consolidation
            design_matrix = np.vstack([cube_points, axial_points, center_points])

            # Conversion to actual scale
            temp_range = (70, 90)
            time_range = (60, 180)
            catalyst_range = (1, 5)

            temperature = (design_matrix[:, 0] + 1) / 2 * (temp_range[1] - temp_range[0]) + temp_range[0]
            time = (design_matrix[:, 1] + 1) / 2 * (time_range[1] - time_range[0]) + time_range[0]
            catalyst = (design_matrix[:, 2] + 1) / 2 * (catalyst_range[1] - catalyst_range[0]) + catalyst_range[0]

            df_design = pd.DataFrame({
                'run': range(1, len(design_matrix) + 1),
                'temperature': temperature,
                'time': time,
                'catalyst': catalyst
            })

        return df_design

    def simulate_response(self, df_design):
        """
        Simulate response variables (yield and purity)

        Args:
            df_design: DataFrame of experimental design

        Returns:
            DataFrame with added response variables
        """
        np.random.seed(42)

        # Yield model (quadratic equation + noise)
        # Optimal conditions: Temperature 80°C, Time 120 min, Catalyst 3% for maximum yield
        temp_centered = (df_design['temperature'] - 80) / 10
        time_centered = (df_design['time'] - 120) / 60
        catalyst_centered = (df_design['catalyst'] - 3) / 2

        yield_base = (
            95  # Base yield
            - 3 * temp_centered ** 2  # Quadratic effect of temperature
            - 2 * time_centered ** 2  # Quadratic effect of time
            - 4 * catalyst_centered ** 2  # Quadratic effect of catalyst
            + 1.5 * temp_centered * time_centered  # Interaction
            + np.random.normal(0, 1, len(df_design))  # Noise
        )

        # Purity model
        purity_base = (
            99.5  # Base purity
            - 0.5 * temp_centered ** 2
            - 0.3 * time_centered ** 2
            - 0.4 * catalyst_centered ** 2
            + np.random.normal(0, 0.1, len(df_design))
        )

        df_design['yield'] = np.clip(yield_base, 0, 100)
        df_design['purity'] = np.clip(purity_base, 95, 100)

        return df_design

    def build_response_surface_model(self, df_design, response_var):
        """
        Build response surface model (quadratic polynomial)

        Args:
            df_design: Experimental data
            response_var: Response variable name

        Returns:
            Trained model
        """
        X = df_design[['temperature', 'time', 'catalyst']].values
        y = df_design[response_var].values

        # Generate quadratic polynomial features
        self.poly_features = PolynomialFeatures(degree=2, include_bias=True)
        X_poly = self.poly_features.fit_transform(X)

        # Linear regression model
        self.model = LinearRegression()
        self.model.fit(X_poly, y)

        # Prediction and evaluation
        y_pred = self.model.predict(X_poly)
        r2 = r2_score(y, y_pred)

        return self.model, r2

    def predict_response(self, temperature, time, catalyst):
        """Predict response"""
        X_new = np.array([[temperature, time, catalyst]])
        X_poly = self.poly_features.transform(X_new)
        return self.model.predict(X_poly)[0]

    def plot_design_space(self, df_design, temp_range, time_range, catalyst_fixed=3):
        """Visualize design space"""
        fig = plt.figure(figsize=(16, 6))

        # Design space on temperature-time plane (catalyst amount fixed)
        temp_grid = np.linspace(temp_range[0], temp_range[1], 50)
        time_grid = np.linspace(time_range[0], time_range[1], 50)
        T_mesh, Ti_mesh = np.meshgrid(temp_grid, time_grid)

        # Yield prediction
        yield_grid = np.zeros_like(T_mesh)
        purity_grid = np.zeros_like(T_mesh)

        for i in range(len(temp_grid)):
            for j in range(len(time_grid)):
                yield_grid[j, i] = self.predict_response(T_mesh[j, i], Ti_mesh[j, i], catalyst_fixed)
                purity_grid[j, i] = self.predict_response(T_mesh[j, i], Ti_mesh[j, i], catalyst_fixed)

        # Yield contour plot
        ax1 = fig.add_subplot(131)
        contour_yield = ax1.contourf(T_mesh, Ti_mesh, yield_grid, levels=20, cmap='RdYlGn')
        ax1.contour(T_mesh, Ti_mesh, yield_grid, levels=[90, 92, 94], colors='black',
                    linewidths=1.5, linestyles='dashed')
        ax1.scatter(df_design['temperature'], df_design['time'], c='blue', s=50,
                    edgecolor='black', linewidth=1, label='Experimental points', zorder=5)
        ax1.set_xlabel('Reaction Temperature (°C)')
        ax1.set_ylabel('Reaction Time (min)')
        ax1.set_title(f'Yield (Catalyst={catalyst_fixed}%)', fontsize=12, fontweight='bold')
        ax1.legend()
        plt.colorbar(contour_yield, ax=ax1, label='Yield (%)')

        # Purity contour plot
        ax2 = fig.add_subplot(132)
        contour_purity = ax2.contourf(T_mesh, Ti_mesh, purity_grid, levels=20, cmap='Blues')
        ax2.contour(T_mesh, Ti_mesh, purity_grid, levels=[99.0, 99.2, 99.4], colors='black',
                    linewidths=1.5, linestyles='dashed')
        ax2.scatter(df_design['temperature'], df_design['time'], c='blue', s=50,
                    edgecolor='black', linewidth=1, label='Experimental points', zorder=5)
        ax2.set_xlabel('Reaction Temperature (°C)')
        ax2.set_ylabel('Reaction Time (min)')
        ax2.set_title(f'Purity (Catalyst={catalyst_fixed}%)', fontsize=12, fontweight='bold')
        ax2.legend()
        plt.colorbar(contour_purity, ax=ax2, label='Purity (%)')

        # Design space (Yield≥90% AND Purity≥99.0%)
        ax3 = fig.add_subplot(133)
        design_space_mask = (yield_grid >= 90) & (purity_grid >= 99.0)

        ax3.contourf(T_mesh, Ti_mesh, design_space_mask.astype(int), levels=[0, 0.5, 1],
                     colors=['red', 'green'], alpha=0.3)
        ax3.contour(T_mesh, Ti_mesh, yield_grid, levels=[90], colors='blue',
                    linewidths=2, linestyles='solid', label='Yield=90%')
        ax3.contour(T_mesh, Ti_mesh, purity_grid, levels=[99.0], colors='orange',
                    linewidths=2, linestyles='solid', label='Purity=99.0%')
        ax3.scatter(df_design['temperature'], df_design['time'], c='blue', s=50,
                    edgecolor='black', linewidth=1, label='Experimental points', zorder=5)
        ax3.set_xlabel('Reaction Temperature (°C)')
        ax3.set_ylabel('Reaction Time (min)')
        ax3.set_title('Design Space (Green=Acceptable Region)', fontsize=12, fontweight='bold')
        ax3.legend()

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

# Execution example
print("=" * 60)
print("QbD Design Space Design System")
print("=" * 60)

# QbD design instance
qbd = QbDDesigner()

# Generate central composite design
df_design = qbd.generate_doe_design(design_type='central_composite')

print(f"\nDesign of Experiments (Central Composite Design):")
print(f"Number of experiments: {len(df_design)}")
print(f"\nFactor levels:")
print(f"  Reaction temperature: 70-90°C")
print(f"  Reaction time: 60-180 min")
print(f"  Catalyst amount: 1-5%")

# Simulate response variables
df_design = qbd.simulate_response(df_design)

print(f"\nResponse variables:")
print(f"  Yield: {df_design['yield'].min():.1f}-{df_design['yield'].max():.1f}%")
print(f"  Purity: {df_design['purity'].min():.2f}-{df_design['purity'].max():.2f}%")

# Build response surface model
model_yield, r2_yield = qbd.build_response_surface_model(df_design, 'yield')
print(f"\nYield model R² = {r2_yield:.4f}")

model_purity, r2_purity = qbd.build_response_surface_model(df_design, 'purity')
print(f"Purity model R² = {r2_purity:.4f}")

# Search for optimal conditions
print(f"\nOptimal condition prediction:")
optimal_yield = qbd.predict_response(80, 120, 3)
print(f"  Temperature=80°C, Time=120 min, Catalyst=3%")
print(f"  → Yield = {optimal_yield:.2f}%")

# Visualize design space
qbd.plot_design_space(df_design, temp_range=(70, 90), time_range=(60, 180), catalyst_fixed=3)

Implementation Points:

• Efficient experimental design using Central Composite Design (CCD)
• Response surface construction using quadratic polynomial models
• Multi-objective optimization (balancing yield and purity)
• Design space visualization and identification of acceptable regions
• Science-based process parameter settings

⚙️ 4.2 Control Strategy for Continuous Manufacturing Systems

Batch vs Continuous Production

💻 Code Example 4.2: Real-Time Control Simulation for Continuous Manufacturing

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

class ContinuousManufacturingSimulator:
    """Continuous Manufacturing Simulator"""

    def __init__(self, target_concentration=5.0):
        """
        Args:
            target_concentration: Target concentration (mol/L)
        """
        self.target = target_concentration
        self.Kp = 2.0  # Proportional gain
        self.Ki = 0.5  # Integral gain
        self.Kd = 0.1  # Derivative gain
        self.integral_error = 0
        self.previous_error = 0

    def reactor_dynamics(self, state, t, feed_rate, T):
        """
        Continuous reactor dynamics (CSTR: Continuous Stirred-Tank Reactor)

        Args:
            state: [Concentration C, Temperature T_reactor]
            t: Time
            feed_rate: Feed flow rate
            T: Set temperature

        Returns:
            State change rate [dC/dt, dT/dt]
        """
        C, T_reactor = state

        # Parameters
        V = 10.0  # Reactor volume (L)
        k0 = 1e10  # Frequency factor
        Ea = 8000  # Activation energy (J/mol)
        R = 8.314  # Gas constant
        C_in = 10.0  # Feed concentration (mol/L)

        # Reaction rate constant (Arrhenius equation)
        k = k0 * np.exp(-Ea / (R * (T_reactor + 273.15)))

        # Mass balance
        dC_dt = feed_rate / V * (C_in - C) - k * C

        # Energy balance (simplified)
        dT_dt = 0.1 * (T - T_reactor)  # Temperature control

        return [dC_dt, dT_dt]

    def pid_controller(self, error, dt=0.1):
        """
        PID controller

        Args:
            error: Deviation (target value - current value)
            dt: Sampling time

        Returns:
            Control output (feed flow rate adjustment)
        """
        # Integral term
        self.integral_error += error * dt

        # Derivative term
        derivative_error = (error - self.previous_error) / dt

        # PID output
        output = (self.Kp * error +
                  self.Ki * self.integral_error +
                  self.Kd * derivative_error)

        self.previous_error = error

        return output

    def simulate_continuous_process(self, duration=100, disturbance_times=None):
        """
        Simulate continuous manufacturing process

        Args:
            duration: Simulation time (min)
            disturbance_times: List of disturbance occurrence times

        Returns:
            DataFrame of time series data
        """
        dt = 0.1  # Sampling time (min)
        t = np.arange(0, duration, dt)

        # Initial state
        C = 5.0  # Initial concentration
        T_reactor = 80.0  # Initial reactor temperature
        state = [C, T_reactor]

        # Recording lists
        concentrations = []
        temperatures = []
        feed_rates = []
        errors = []
        setpoints = []

        # Base feed flow rate
        base_feed_rate = 1.0  # L/min

        for i, time in enumerate(t):
            # Target value (setpoint change simulation)
            if time < 30:
                target = 5.0
            elif time < 60:
                target = 6.0
            else:
                target = 5.5

            # Add disturbances
            if disturbance_times and any(abs(time - dt_time) < 0.5 for dt_time in disturbance_times):
                state[0] += np.random.uniform(-0.5, 0.5)  # Concentration disturbance

            # PID control
            error = target - state[0]
            control_output = self.pid_controller(error, dt)

            # Update feed flow rate (with constraints)
            feed_rate = np.clip(base_feed_rate + control_output, 0.5, 2.0)

            # Update reactor dynamics
            state = odeint(self.reactor_dynamics, state, [time, time + dt],
                           args=(feed_rate, 80.0))[-1]

            # Record
            concentrations.append(state[0])
            temperatures.append(state[1])
            feed_rates.append(feed_rate)
            errors.append(error)
            setpoints.append(target)

        df_results = pd.DataFrame({
            'time': t,
            'concentration': concentrations,
            'temperature': temperatures,
            'feed_rate': feed_rates,
            'error': errors,
            'setpoint': setpoints
        })

        return df_results

    def plot_control_performance(self, df_results):
        """Visualize control performance"""
        fig, axes = plt.subplots(3, 1, figsize=(14, 10))

        # Concentration profile
        axes[0].plot(df_results['time'], df_results['concentration'], color='#11998e',
                     linewidth=2, label='Measured concentration')
        axes[0].plot(df_results['time'], df_results['setpoint'], color='red',
                     linestyle='--', linewidth=2, label='Target concentration')
        axes[0].fill_between(df_results['time'],
                              df_results['setpoint'] - 0.1,
                              df_results['setpoint'] + 0.1,
                              alpha=0.2, color='green', label='Acceptable range (±0.1)')
        axes[0].set_xlabel('Time (min)')
        axes[0].set_ylabel('Concentration (mol/L)')
        axes[0].set_title('Continuous Reactor Concentration Control', fontsize=12, fontweight='bold')
        axes[0].legend()
        axes[0].grid(alpha=0.3)

        # Feed flow rate
        axes[1].plot(df_results['time'], df_results['feed_rate'], color='#38ef7d',
                     linewidth=2)
        axes[1].set_xlabel('Time (min)')
        axes[1].set_ylabel('Feed Flow Rate (L/min)')
        axes[1].set_title('PID Controller Output (Feed Flow Rate)', fontsize=12, fontweight='bold')
        axes[1].grid(alpha=0.3)

        # Control error
        axes[2].plot(df_results['time'], df_results['error'], color='#f38181',
                     linewidth=1.5)
        axes[2].axhline(y=0, color='black', linestyle='-', linewidth=1)
        axes[2].fill_between(df_results['time'], -0.1, 0.1,
                              alpha=0.2, color='green', label='Acceptable error')
        axes[2].set_xlabel('Time (min)')
        axes[2].set_ylabel('Control Error (mol/L)')
        axes[2].set_title('Control Error', fontsize=12, fontweight='bold')
        axes[2].legend()
        axes[2].grid(alpha=0.3)

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

# Execution example
print("=" * 60)
print("Continuous Manufacturing Real-Time Control Simulation")
print("=" * 60)

# Initialize simulator
cm_sim = ContinuousManufacturingSimulator(target_concentration=5.0)

# Set disturbance occurrence times
disturbance_times = [40, 70]

# Execute simulation
df_results = cm_sim.simulate_continuous_process(duration=100, disturbance_times=disturbance_times)

# Performance evaluation
print(f"\nSimulation results:")
print(f"Simulation time: {df_results['time'].max():.1f} min")
print(f"Concentration range: {df_results['concentration'].min():.2f}-{df_results['concentration'].max():.2f} mol/L")
print(f"Average control error: {df_results['error'].abs().mean():.4f} mol/L")
print(f"Maximum control error: {df_results['error'].abs().max():.4f} mol/L")

# Steady-state error evaluation (last 10 minutes)
steady_state_data = df_results[df_results['time'] > 90]
steady_state_error = steady_state_data['error'].abs().mean()
print(f"Steady-state error: {steady_state_error:.4f} mol/L")

# Visualization
cm_sim.plot_control_performance(df_results)

Implementation Points:

• Dynamic modeling of continuous reactor (CSTR)
• Real-time feedback control using PID control
• Evaluation of responsiveness to setpoint changes and disturbances
• Quantitative evaluation of control performance (steady-state error, response time)
• Addressing control challenges unique to continuous production not found in batch production

📚 Summary

In this chapter, we learned about continuous manufacturing and QbD implementation.

Key Points

• The three pillars of QbD (QTPP, CQA, CPP) and scientific approach
• Efficient process optimization using Design of Experiments (DoE)
• Design space construction and visualization
• Real-time control strategies for continuous manufacturing systems
• Stable quality assurance through PID control