Response Surface Methods

Response Surface Methodology (RSM) is a collection of statistical and mathematical techniques useful for developing, improving, and optimizing processes where the response of interest is influenced by several factors. RSM quantifies the relationship between input variables and one or more responses, enabling systematic exploration of the factor space to identify optimal operating conditions. The methodology integrates experimental design, regression modeling, and optimization algorithms within a unified mathematical framework.

Quadratic Response Surface Model

Stationary Point and Canonical Analysis

The stationary point of the response surface is found by setting the gradient to zero.

Steepest Ascent Method

When the stationary point lies outside the current experimental region, the method of steepest ascent provides a systematic path toward the optimum.

Ridge Analysis

Ridge analysis finds the stationary point constrained to a specified radius from the design center, particularly useful when the unconstrained optimum lies far from the experimental region.

Desirability Functions

For multi-response optimization, desirability functions transform each response into a dimensionless quality metric.

Python Implementation

import numpy as np
import pandas as pd
from scipy import optimize, linalg
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Fit second-order response surface model
def fit_response_surface(design_matrix, response):
    """Fit second-order RSM model and return coefficients."""
    # Add quadratic and interaction terms
    k = design_matrix.shape[1]
    X_quad = design_matrix.copy()
    
    # Quadratic terms
    for j in range(k):
        X_quad = np.column_stack([X_quad, design_matrix[:, j]**2])
    
    # Interaction terms
    for i in range(k):
        for j in range(i+1, k):
            X_quad = np.column_stack([X_quad, design_matrix[:, i] * design_matrix[:, j]])
    
    # Add intercept
    X_quad = np.column_stack([np.ones(X_quad.shape[0]), X_quad])
    
    # Least squares fit
    coeffs = np.linalg.lstsq(X_quad, response, rcond=None)[0]
    
    return coeffs, X_quad

# Canonical analysis
def canonical_analysis(coeffs, k):
    """Perform canonical analysis of response surface."""
    # Extract coefficient vectors
    b0 = coeffs[0]
    b = coeffs[1:k+1]  # Linear coefficients
    
    # Construct B matrix (quadratic + interaction)
    B = np.zeros((k, k))
    idx = k + 1  # Start of quadratic terms
    
    # Quadratic terms
    for j in range(k):
        B[j, j] = coeffs[idx]
        idx += 1
    
    # Interaction terms (off-diagonal)
    for i in range(k):
        for j in range(i+1, k):
            B[i, j] = coeffs[idx] / 2
            B[j, i] = coeffs[idx] / 2
            idx += 1
    
    # Stationary point
    try:
        x_s = -0.5 * np.linalg.solve(B, b)
    except np.linalg.LinAlgError:
        x_s = None
    
    # Eigenvalues and eigenvectors
    eigenvalues, eigenvectors = np.linalg.eigh(B)
    
    # Predicted response at stationary point
    if x_s is not None:
        Y_s = b0 + np.dot(b, x_s) + 0.5 * np.dot(x_s, B @ x_s)
    else:
        Y_s = None
    
    return {
        'stationary_point': x_s,
        'eigenvalues': eigenvalues,
        'eigenvectors': eigenvectors,
        'Y_stationary': Y_s,
        'B_matrix': B,
        'b_vector': b
    }

# Ridge analysis
def ridge_analysis(coeffs, k, R_values):
    """Compute ridge of stationary points."""
    b0 = coeffs[0]
    b = coeffs[1:k+1]
    
    B = np.zeros((k, k))
    idx = k + 1
    for j in range(k):
        B[j, j] = coeffs[idx]
        idx += 1
    for i in range(k):
        for j in range(i+1, k):
            B[i, j] = coeffs[idx] / 2
            B[j, i] = coeffs[idx] / 2
            idx += 1
    
    # Eigendecomposition
    eigenvalues, P = np.linalg.eigh(B)
    delta = P.T @ b / 2
    
    ridge_points = []
    ridge_responses = []
    
    for R in R_values:
        # Find mu satisfying constraint
        def constraint(mu):
            w = delta / (eigenvalues - mu)
            return np.sum(w**2) - R**2
        
        # Search for mu
        mu_min = max(eigenvalues) + 0.01
        try:
            mu_solution = optimize.brentq(constraint, mu_min, mu_min + 100)
        except ValueError:
            mu_solution = mu_min + 50
        
        # Compute ridge point in original coordinates
        w = delta / (eigenvalues - mu_solution)
        x_r = P @ w
        
        # Predicted response
        Y_r = b0 + np.dot(b, x_r) + 0.5 * np.dot(x_r, B @ x_r)
        
        ridge_points.append(x_r)
        ridge_responses.append(Y_r)
    
    return np.array(ridge_points), np.array(ridge_responses)

# Desirability function
def desirability(y, target, lower, upper, objective='target', r=1):
    """Compute individual desirability."""
    if objective == 'maximize':
        if y <= lower:
            return 0
        elif y >= target:
            return 1
        else:
            return ((y - lower) / (target - lower))**r
    
    elif objective == 'minimize':
        if y >= upper:
            return 0
        elif y <= target:
            return 1
        else:
            return ((upper - y) / (upper - target))**r
    
    else:  # target
        if y <= lower or y >= upper:
            return 0
        elif y <= target:
            return ((y - lower) / (target - lower))**r
        else:
            return ((upper - y) / (upper - target))**r

def composite_desirability(responses, targets, lowers, uppers, 
                          objectives, weights=None):
    """Compute composite desirability across multiple responses."""
    m = len(responses)
    if weights is None:
        weights = np.ones(m)
    
    d_values = []
    for i in range(m):
        d = desirability(responses[i], targets[i], lowers[i], uppers[i],
                        objectives[i])
        d_values.append(d)
    
    D = np.prod(np.array(d_values)**weights)**(1/np.sum(weights))
    return D, d_values

# Example: Fit response surface to simulated data
np.random.seed(42)

# Generate design points (CCD)
factorial = np.array([[-1, -1], [-1, 1], [1, -1], [1, 1]])
axial = np.array([[-1.414, 0], [1.414, 0], [0, -1.414], [0, 1.414]])
center = np.array([[0, 0], [0, 0], [0, 0]])
design = np.vstack([factorial, axial, center])

# True model: Y = 85 + 2.1*X1 + 1.8*X2 - 3.2*X1^2 - 2.8*X2^2 + 1.5*X1*X2
def true_model(x1, x2):
    return 85 + 2.1*x1 + 1.8*x2 - 3.2*x1**2 - 2.8*x2**2 + 1.5*x1*x2

response = np.array([true_model(x[0], x[1]) + np.random.normal(0, 1) 
                     for x in design])

# Fit model
k = 2
coeffs, X_quad = fit_response_surface(design, response)
print("Fitted coefficients:", np.round(coeffs, 3))

# Canonical analysis
result = canonical_analysis(coeffs, k)
print(f"\nStationary point: {np.round(result['stationary_point'], 3)}")
print(f"Eigenvalues: {np.round(result['eigenvalues'], 3)}")
print(f"Predicted response at stationary point: {result['Y_stationary']:.3f}")

# Classify stationary point
if all(result['eigenvalues'] < 0):
    print("Classification: Maximum")
elif all(result['eigenvalues'] > 0):
    print("Classification: Minimum")
else:
    print("Classification: Saddle point")

# Ridge analysis
R_values = np.linspace(0, 2, 20)
ridge_pts, ridge_resp = ridge_analysis(coeffs, k, R_values)

# Plot response surface and ridge
fig = plt.figure(figsize=(12, 5))

# 3D response surface
ax1 = fig.add_subplot(121, projection='3d')
x1_range = np.linspace(-1.5, 1.5, 50)
x2_range = np.linspace(-1.5, 1.5, 50)
X1, X2 = np.meshgrid(x1_range, x2_range)
Y_surface = 85 + 2.1*X1 + 1.8*X2 - 3.2*X1**2 - 2.8*X2**2 + 1.5*X1*X2

ax1.plot_surface(X1, X2, Y_surface, alpha=0.7, cmap='viridis')
ax1.scatter(ridge_pts[:, 0], ridge_pts[:, 1], ridge_resp, 
           c='red', s=50, label='Ridge path')
ax1.scatter(*result['stationary_point'], result['Y_stationary'], 
           c='red', s=200, marker='*', label='Stationary point')
ax1.set_xlabel('X1')
ax1.set_ylabel('X2')
ax1.set_zlabel('Y')
ax1.set_title('Response Surface with Ridge Path')
ax1.legend()

# Ridge path in factor space
ax2 = fig.add_subplot(122)
ax2.plot(ridge_pts[:, 0], ridge_pts[:, 1], 'r-o', markersize=3, label='Ridge path')
ax2.scatter(*result['stationary_point'], c='red', s=200, marker='*', 
           label='Stationary point')
ax2.contour(X1, X2, Y_surface, levels=20, alpha=0.6)
ax2.set_xlabel('X1')
ax2.set_ylabel('X2')
ax2.set_title('Contour Plot with Ridge Path')
ax2.legend()
ax2.set_aspect('equal')

plt.tight_layout()
plt.savefig('rsm_analysis.png', dpi=150)
plt.show()

# Multi-response desirability optimization
# Example: Optimize for max strength and min cost
def multi_response_optimization():
    """Find optimal conditions using desirability."""
    # Simulated response functions
    def strength(x1, x2):
        return 85 + 2.1*x1 + 1.8*x2 - 3.2*x1**2 - 2.8*x2**2
    
    def cost(x1, x2):
        return 50 - 3*x1 + 2*x2 + x1**2 + 0.5*x2**2
    
    # Search grid
    x1_grid = np.linspace(-1, 1, 20)
    x2_grid = np.linspace(-1, 1, 20)
    
    best_D = 0
    best_x = None
    
    for x1 in x1_grid:
        for x2 in x2_grid:
            s = strength(x1, x2)
            c = cost(x1, x2)
            
            D, _ = composite_desirability(
                responses=[s, c],
                targets=[90, 45],  # max strength, min cost
                lowers=[80, 40],
                uppers=[100, 60],
                objectives=['maximize', 'minimize'],
                weights=[0.6, 0.4]
            )
            
            if D > best_D:
                best_D = D
                best_x = (x1, x2)
    
    return best_x, best_D

optimal, D_opt = multi_response_optimization()
print(f"\nOptimal conditions: X1={optimal[0]:.3f}, X2={optimal[1]:.3f}")
print(f"Composite desirability: {D_opt:.4f}")