Chapter 3 ProcessOptimization

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


            <div class="code-header">📊 Code Examples1: Bayesian Optimizationによる加熱条件のOptimization</div>
            <pre><code class="language-python">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
from scipy.stats import norm
import warnings
warnings.filterwarnings('ignore')

# Food加熱Processの目的関数（実際のProcessをシミュレート）
def food_processing_objective(temperature, time):
    """
    目的関数: 品質スコアを最大化（温度と hoursの関数）
    最適点付近: temperature=90°C, time=20 minutes
    """
    # 複雑な非線形関数（実際のFoodProcessを模擬）
    quality = (
        100 * np.exp(-0.5 * ((temperature - 90)/10)**2) * 
        np.exp(-0.5 * ((time - 20)/5)**2) +
        10 * np.sin(temperature / 10) * np.cos(time / 5) +
        np.random.normal(0, 2)  # 測定ノイズ
    )
    return quality

# Bayesian Optimizationの実装
class BayesianOptimizer:
    def __init__(self, bounds, n_init=5):
        self.bounds = np.array(bounds)
        self.n_init = n_init
        self.X_sample = []
        self.y_sample = []
        
        # ガウス過程回帰モデル
        kernel = C(1.0, (1e-3, 1e3)) * RBF([10, 10], (1e-2, 1e2))
        self.gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10,
                                           alpha=1e-6, normalize_y=True)
    
    def acquisition_function(self, X, xi=0.01):
        """
        獲得関数: Expected Improvement (EI)
        """
        mu, sigma = self.gp.predict(X, return_std=True)
        mu_sample_opt = np.max(self.y_sample)
        
        with np.errstate(divide='warn'):
            imp = mu - mu_sample_opt - xi
            Z = imp / sigma
            ei = imp * norm.cdf(Z) + sigma * norm.pdf(Z)
            ei[sigma == 0.0] = 0.0
        
        return ei
    
    def propose_location(self):
        """
        次の評価点を提案（EIを最大化）
        """
        dim = self.bounds.shape[0]
        min_val = float('inf')
        min_x = None
        
        # ランダムサンプリング + Optimization
        n_restarts = 25
        for _ in range(n_restarts):
            x0 = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1], size=dim)
            res = self._minimize_acquisition(x0)
            if res < min_val:
                min_val = res
                min_x = x0
        
        return min_x
    
    def _minimize_acquisition(self, x0):
        """
        獲得関数の最小化（負のEIを最小化）
        """
        from scipy.optimize import minimize
        
        def objective(x):
            return -self.acquisition_function(x.reshape(1, -1))
        
        bounds_list = [(self.bounds[i, 0], self.bounds[i, 1]) 
                       for i in range(self.bounds.shape[0])]
        
        res = minimize(objective, x0, bounds=bounds_list, method='L-BFGS-B')
        return res.fun
    
    def optimize(self, objective_func, n_iter=20):
        """
        Bayesian Optimizationの実行
        """
        # 初期ランダムサンプリング
        for _ in range(self.n_init):
            x = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1], 
                                  size=self.bounds.shape[0])
            y = objective_func(x[0], x[1])
            self.X_sample.append(x)
            self.y_sample.append(y)
        
        # Bayesian Optimizationループ
        for i in range(n_iter):
            # GPモデル更新
            self.gp.fit(np.array(self.X_sample), np.array(self.y_sample))
            
            # 次の評価点を提案
            x_next = self.propose_location()
            y_next = objective_func(x_next[0], x_next[1])
            
            self.X_sample.append(x_next)
            self.y_sample.append(y_next)
            
            if (i + 1) % 5 == 0:
                print(f"Iteration {i+1}: Best quality = {np.max(self.y_sample):.2f} "
                      f"at T={self.X_sample[np.argmax(self.y_sample)][0]:.1f}°C, "
                      f"t={self.X_sample[np.argmax(self.y_sample)][1]:.1f}min")
        
        return np.array(self.X_sample), np.array(self.y_sample)

# Optimization実行
bounds = [[75, 105], [10, 30]]  # [温度範囲,  hours範囲]
optimizer = BayesianOptimizer(bounds, n_init=5)
X_sample, y_sample = optimizer.optimize(food_processing_objective, n_iter=25)

# 最適解
best_idx = np.argmax(y_sample)
best_temp, best_time = X_sample[best_idx]
best_quality = y_sample[best_idx]

print(f"\n=== Optimization結果 ===")
print(f"最適温度: {best_temp:.2f}°C")
print(f"最適 hours: {best_time:.2f} minutes")
print(f"最大品質スコア: {best_quality:.2f}")
print(f"総評価回数: {len(X_sample)}")

# 可視化
fig = plt.figure(figsize=(16, 10))
gs = fig.add_gridspec(2, 2, hspace=0.3, wspace=0.3)

# 1. 真の目的関数（参考）
ax1 = fig.add_subplot(gs[0, 0])
T_grid = np.linspace(75, 105, 100)
t_grid = np.linspace(10, 30, 100)
T_mesh, t_mesh = np.meshgrid(T_grid, t_grid)
quality_grid = np.zeros_like(T_mesh)
for i in range(T_mesh.shape[0]):
    for j in range(T_mesh.shape[1]):
        quality_grid[i, j] = food_processing_objective(T_mesh[i, j], t_mesh[i, j])

contour = ax1.contourf(T_mesh, t_mesh, quality_grid, levels=20, cmap='viridis', alpha=0.8)
ax1.scatter(X_sample[:5, 0], X_sample[:5, 1], c='red', s=100, 
            marker='o', edgecolors='white', linewidth=2, label='初期サンプル', zorder=5)
ax1.scatter(X_sample[5:, 0], X_sample[5:, 1], c='white', s=80, 
            marker='s', edgecolors='black', linewidth=1.5, label='BO提案点', zorder=5)
ax1.scatter(best_temp, best_time, c='gold', s=300, marker='*', 
            edgecolors='red', linewidth=2, label='最適解', zorder=10)
ax1.set_xlabel('温度 (°C)', fontsize=12)
ax1.set_ylabel(' hours ( minutes)', fontsize=12)
ax1.set_title('Bayesian Optimizationの探索過程', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
plt.colorbar(contour, ax=ax1, label='品質スコア')

# 2. 品質スコアの推移
ax2 = fig.add_subplot(gs[0, 1])
iterations = np.arange(1, len(y_sample) + 1)
best_so_far = np.maximum.accumulate(y_sample)
ax2.plot(iterations, y_sample, 'o-', color='#11998e', alpha=0.6, 
         linewidth=2, markersize=6, label='各評価点の品質')
ax2.plot(iterations, best_so_far, 'r-', linewidth=3, label='最良値の推移')
ax2.axvline(x=5, color='gray', linestyle='--', alpha=0.5, label='初期サンプル終了')
ax2.set_xlabel('評価回数', fontsize=12)
ax2.set_ylabel('品質スコア', fontsize=12)
ax2.set_title('Optimizationの収束過程', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 3. ガウス過程によるPrediction（温度方向のスライス、 hours=best_time）
ax3 = fig.add_subplot(gs[1, 0])
T_test = np.linspace(75, 105, 200).reshape(-1, 1)
t_test = np.full_like(T_test, best_time)
X_test = np.hstack([T_test, t_test])

mu, sigma = optimizer.gp.predict(X_test, return_std=True)

ax3.plot(T_test, mu, 'b-', linewidth=2, label='GPPrediction平均')
ax3.fill_between(T_test.ravel(), mu - 1.96*sigma, mu + 1.96*sigma,
                 alpha=0.3, color='blue', label='95%信頼区間')
ax3.scatter(X_sample[:, 0], y_sample, c='red', s=80, marker='o', 
            edgecolors='white', linewidth=1.5, label='観測点', zorder=5)
ax3.axvline(x=best_temp, color='green', linestyle='--', linewidth=2, label='最適温度')
ax3.set_xlabel('温度 (°C)', fontsize=12)
ax3.set_ylabel('品質スコア', fontsize=12)
ax3.set_title(f'ガウス過程Prediction（ hours={best_time:.1f} minutes固定）', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 4. 獲得関数（Expected Improvement）
ax4 = fig.add_subplot(gs[1, 1])
ei = optimizer.acquisition_function(X_test)
ax4.plot(T_test, ei, 'g-', linewidth=2)
ax4.fill_between(T_test.ravel(), 0, ei, alpha=0.3, color='green')
ax4.set_xlabel('温度 (°C)', fontsize=12)
ax4.set_ylabel('Expected Improvement', fontsize=12)
ax4.set_title('獲得関数（次の評価点の選択基準）', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)

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

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

"""
Example: FoodProcessでは、品質、コスト、エネルギー効率など複数の目的を同時にOptimizationする必要があります

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

            <div class="code-header">📊 Code Examples2: 多目的Optimization（品質 vs コスト）</div>
            <pre><code class="language-python">import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution

# 多目的Optimization問題の定義
class FoodProcessMultiObjective:
    def __init__(self):
        # パラメータ範囲
        self.bounds = [
            (75, 105),   # 温度 (°C)
            (10, 30),    #  hours ( minutes)
            (3.0, 5.0),  # pH
            (50, 150),   # 流量 (L/min)
        ]
    
    def quality_objective(self, params):
        """
        目的関数1: 品質スコアを最大化（最小化問題なので負値）
        """
        temp, time, ph, flow = params
        
        # 品質関数（非線形）
        quality = (
            100 * np.exp(-0.5 * ((temp - 90)/10)**2) *
            np.exp(-0.5 * ((time - 20)/5)**2) *
            np.exp(-0.5 * ((ph - 4.0)/0.5)**2) *
            (1 + 0.1 * np.sin(flow / 20))
        )
        
        return -quality  # 最小化問題に変換
    
    def cost_objective(self, params):
        """
        目的関数2: コストを最小化
        """
        temp, time, ph, flow = params
        
        # エネルギーコスト（温度と hoursに依存）
        energy_cost = 0.1 * (temp - 70) * time
        
        # pH調整剤コスト
        ph_cost = 5 * abs(ph - 4.0)
        
        # ポンプ運転コスト
        pump_cost = 0.05 * flow * time
        
        total_cost = energy_cost + ph_cost + pump_cost
        
        return total_cost
    
    def constraints_check(self, params):
        """
        制約条件のチェック（HACCP基準など）
        """
        temp, time, ph, flow = params
        
        # F値（殺菌効果）の計算: F0 = time * 10^((temp-121)/10)
        F0 = time * (10 ** ((temp - 121) / 10))
        
        # 最小F値要求（例: F0 >= 3）
        if F0 < 3:
            return False
        
        # pH範囲制約
        if not (3.0 <= ph <= 5.0):
            return False
        
        return True

# 簡易的なパレートOptimization（グリッドサーチ + パレート判定）
problem = FoodProcessMultiObjective()

# ランダムサンプリングでパレートフロント近似
np.random.seed(42)
n_samples = 1000

# ランダムサンプル生成
samples = []
qualities = []
costs = []

for _ in range(n_samples):
    params = [
        np.random.uniform(low, high) 
        for low, high in problem.bounds
    ]
    
    if problem.constraints_check(params):
        quality = -problem.quality_objective(params)  # 元に戻す（最大化）
        cost = problem.cost_objective(params)
        
        samples.append(params)
        qualities.append(quality)
        costs.append(cost)

samples = np.array(samples)
qualities = np.array(qualities)
costs = np.array(costs)

# パレート最適解の判定
def is_pareto_efficient(costs_qualities):
    """
    パレート最適解を判定（2目的の場合）
    cost: 最小化, quality: 最大化
    """
    is_efficient = np.ones(costs_qualities.shape[0], dtype=bool)
    for i, c in enumerate(costs_qualities):
        if is_efficient[i]:
            # 他の解がこの解を支配しているかチェック
            # 支配: cost <= c[0] かつ quality >= c[1] でどちらかが真に優れている
            dominated = np.logical_and(
                costs_qualities[:, 0] <= c[0],  # コストが同じか少ない
                costs_qualities[:, 1] >= c[1]   # 品質が同じか高い
            )
            # 真に優れている（等しくない）
            strictly_better = np.logical_or(
                costs_qualities[:, 0] < c[0],
                costs_qualities[:, 1] > c[1]
            )
            dominated = np.logical_and(dominated, strictly_better)
            
            if np.any(dominated):
                is_efficient[i] = False
    
    return is_efficient

# パレート最適解の抽出
costs_qualities = np.column_stack([costs, qualities])
pareto_mask = is_pareto_efficient(costs_qualities)
pareto_samples = samples[pareto_mask]
pareto_costs = costs[pareto_mask]
pareto_qualities = qualities[pareto_mask]

# パレートフロントをコストでソート
sorted_indices = np.argsort(pareto_costs)
pareto_costs_sorted = pareto_costs[sorted_indices]
pareto_qualities_sorted = pareto_qualities[sorted_indices]
pareto_samples_sorted = pareto_samples[sorted_indices]

# 代表的な3つの解を選択
n_pareto = len(pareto_costs_sorted)
representative_indices = [0, n_pareto // 2, n_pareto - 1]
representative_solutions = []

for idx in representative_indices:
    sol = {
        'params': pareto_samples_sorted[idx],
        'quality': pareto_qualities_sorted[idx],
        'cost': pareto_costs_sorted[idx],
        'label': ['コスト重視', 'バランス型', '品質重視'][representative_indices.index(idx)]
    }
    representative_solutions.append(sol)

# 可視化
fig = plt.figure(figsize=(16, 10))
gs = fig.add_gridspec(2, 2, hspace=0.3, wspace=0.3)

# 1. パレートフロント
ax1 = fig.add_subplot(gs[0, :])
ax1.scatter(costs, qualities, c='lightgray', alpha=0.4, s=30, label='全サンプル')
ax1.plot(pareto_costs_sorted, pareto_qualities_sorted, 'r-', linewidth=2.5, 
         label='パレートフロント', zorder=10)
ax1.scatter(pareto_costs_sorted, pareto_qualities_sorted, c='red', s=100, 
            marker='o', edgecolors='white', linewidth=2, label='パレート最適解', zorder=11)

# 代表解をハイライト
colors_rep = ['#1f77b4', '#ff7f0e', '#2ca02c']
for i, sol in enumerate(representative_solutions):
    ax1.scatter(sol['cost'], sol['quality'], c=colors_rep[i], s=300, 
                marker='*', edgecolors='black', linewidth=2, 
                label=sol['label'], zorder=12)
    ax1.annotate(sol['label'], 
                xy=(sol['cost'], sol['quality']),
                xytext=(15, 15), textcoords='offset points',
                fontsize=11, fontweight='bold',
                bbox=dict(boxstyle='round,pad=0.5', fc=colors_rep[i], alpha=0.7),
                arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0'))

ax1.set_xlabel('コスト（円/バッチ）', fontsize=13)
ax1.set_ylabel('品質スコア', fontsize=13)
ax1.set_title('多目的Optimization：パレートフロント', fontsize=15, fontweight='bold')
ax1.legend(fontsize=11, loc='upper right')
ax1.grid(True, alpha=0.3)

# 2-4. 代表解のパラメータ比較
param_names = ['温度 (°C)', ' hours ( minutes)', 'pH', '流量 (L/min)']
axes = [fig.add_subplot(gs[1, 0]), fig.add_subplot(gs[1, 1])]

# 温度と hours
ax2 = axes[0]
x_pos = np.arange(len(representative_solutions))
width = 0.35

temp_values = [sol['params'][0] for sol in representative_solutions]
time_values = [sol['params'][1] for sol in representative_solutions]

ax2_twin = ax2.twinx()
bars1 = ax2.bar(x_pos - width/2, temp_values, width, label='温度', color='#ff6b6b', alpha=0.8)
bars2 = ax2_twin.bar(x_pos + width/2, time_values, width, label=' hours', color='#4ecdc4', alpha=0.8)

ax2.set_xlabel('解のタイプ', fontsize=12)
ax2.set_ylabel('温度 (°C)', fontsize=12, color='#ff6b6b')
ax2_twin.set_ylabel(' hours ( minutes)', fontsize=12, color='#4ecdc4')
ax2.set_title('代表解のパラメータ比較（温度・ hours）', fontsize=13, fontweight='bold')
ax2.set_xticks(x_pos)
ax2.set_xticklabels([sol['label'] for sol in representative_solutions])
ax2.tick_params(axis='y', labelcolor='#ff6b6b')
ax2_twin.tick_params(axis='y', labelcolor='#4ecdc4')
ax2.grid(True, alpha=0.3, axis='y')

# pH と流量
ax3 = axes[1]
ph_values = [sol['params'][2] for sol in representative_solutions]
flow_values = [sol['params'][3] for sol in representative_solutions]

ax3_twin = ax3.twinx()
bars3 = ax3.bar(x_pos - width/2, ph_values, width, label='pH', color='#95e1d3', alpha=0.8)
bars4 = ax3_twin.bar(x_pos + width/2, flow_values, width, label='流量', color='#f38181', alpha=0.8)

ax3.set_xlabel('解のタイプ', fontsize=12)
ax3.set_ylabel('pH', fontsize=12, color='#95e1d3')
ax3_twin.set_ylabel('流量 (L/min)', fontsize=12, color='#f38181')
ax3.set_title('代表解のパラメータ比較（pH・流量）', fontsize=13, fontweight='bold')
ax3.set_xticks(x_pos)
ax3.set_xticklabels([sol['label'] for sol in representative_solutions])
ax3.tick_params(axis='y', labelcolor='#95e1d3')
ax3_twin.tick_params(axis='y', labelcolor='#f38181')
ax3.grid(True, alpha=0.3, axis='y')

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

# レポート出力
print("=== 多目的Optimization結果 ===")
print(f"総サンプル数: {n_samples}")
print(f"制約を満たすサンプル数: {len(samples)}")
print(f"パレート最適解数: {len(pareto_samples)}")

print("\n=== 代表的な3つの解 ===")
for i, sol in enumerate(representative_solutions):
    print(f"\n{i+1}. {sol['label']}")
    print(f"   温度: {sol['params'][0]:.2f}°C")
    print(f"    hours: {sol['params'][1]:.2f} minutes")
    print(f"   pH: {sol['params'][2]:.2f}")
    print(f"   流量: {sol['params'][3]:.2f} L/min")
    print(f"   品質スコア: {sol['quality']:.2f}")
    print(f"   コスト: {sol['cost']:.2f} 円/バッチ")
    
    # F値計算
    F0 = sol['params'][1] * (10 ** ((sol['params'][0] - 121) / 10))
    print(f"   殺菌効果（F0値）: {F0:.2f}")

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


            <div class="code-header">📊 Code Examples3: 応答曲面法（Box-Behnken計画）</div>
            <pre><code class="language-python">import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from mpl_toolkits.mplot3d import Axes3D

# Box-Behnken実験計画（3因子）
def box_behnken_design(factors):
    """
    Box-Behnken実験計画を生成（3因子の場合）
    """
    # 符号化Level: -1, 0, +1
    design = np.array([
        [-1, -1,  0], [1, -1,  0], [-1,  1,  0], [1,  1,  0],
        [-1,  0, -1], [1,  0, -1], [-1,  0,  1], [1,  0,  1],
        [ 0, -1, -1], [0,  1, -1], [ 0, -1,  1], [0,  1,  1],
        [ 0,  0,  0], [0,  0,  0], [ 0,  0,  0]  # 中心点（3回反復）
    ])
    
    # 実際の値に変換
    actual_design = []
    for row in design:
        actual_row = []
        for i, level in enumerate(row):
            low, center, high = factors[i]
            if level == -1:
                actual_row.append(low)
            elif level == 0:
                actual_row.append(center)
            else:  # level == 1
                actual_row.append(high)
        actual_design.append(actual_row)
    
    return np.array(actual_design)

# 因子Level設定（低・中・高）
factors = [
    (80, 90, 100),   # 温度 (°C)
    (15, 20, 25),    #  hours ( minutes)
    (3.5, 4.0, 4.5), # pH
]
factor_names = ['温度', ' hours', 'pH']

# 実験計画生成
X_design = box_behnken_design(factors)

# 応答変数（品質スコア）のシミュレート
def response_function(temp, time, ph):
    """
    真の応答関数（2次多項式 + 交互作用 + ノイズ）
    """
    response = (
        50 +
        2 * (temp - 90) + 0.5 * (time - 20) + 10 * (ph - 4.0) +
        -0.05 * (temp - 90)**2 - 0.1 * (time - 20)**2 - 15 * (ph - 4.0)**2 +
        0.02 * (temp - 90) * (time - 20) +
        0.5 * (temp - 90) * (ph - 4.0) +
        np.random.normal(0, 0.5)  # 実験誤差
    )
    return response

y_response = np.array([response_function(x[0], x[1], x[2]) for x in X_design])

# データフレーム作成
df_experiment = pd.DataFrame(X_design, columns=factor_names)
df_experiment['品質スコア'] = y_response

print("=== Box-Behnken実験計画 ===")
print(df_experiment)

# 2次応答曲面モデルの構築
poly = PolynomialFeatures(degree=2, include_bias=True)
X_poly = poly.fit_transform(X_design)
model = LinearRegression()
model.fit(X_poly, y_response)

# モデル性能
y_pred = model.predict(X_poly)
r2_score = model.score(X_poly, y_response)
print(f"\n=== 応答曲面モデル ===")
print(f"R² スコア: {r2_score:.4f}")

# 係数の表示
feature_names = poly.get_feature_names_out(factor_names)
coefficients = pd.DataFrame({
    '項': feature_names,
    '係数': model.coef_
})
print("\n=== モデル係数 ===")
print(coefficients)

# Optimization：グリッドサーチで最大値を探索
temp_range = np.linspace(80, 100, 50)
time_range = np.linspace(15, 25, 50)
ph_fixed = 4.0  # pHを固定

T_mesh, t_mesh = np.meshgrid(temp_range, time_range)
response_mesh = np.zeros_like(T_mesh)

for i in range(T_mesh.shape[0]):
    for j in range(T_mesh.shape[1]):
        X_test = np.array([[T_mesh[i, j], t_mesh[i, j], ph_fixed]])
        X_test_poly = poly.transform(X_test)
        response_mesh[i, j] = model.predict(X_test_poly)[0]

# 最適点
max_idx = np.unravel_index(np.argmax(response_mesh), response_mesh.shape)
optimal_temp = T_mesh[max_idx]
optimal_time = t_mesh[max_idx]
optimal_response = response_mesh[max_idx]

print(f"\n=== 最適条件（pH={ph_fixed}固定） ===")
print(f"最適温度: {optimal_temp:.2f}°C")
print(f"最適 hours: {optimal_time:.2f} minutes")
print(f"Prediction品質スコア: {optimal_response:.2f}")

# 可視化
fig = plt.figure(figsize=(16, 12))

# 1. 3D応答曲面
ax1 = fig.add_subplot(2, 2, 1, projection='3d')
surf = ax1.plot_surface(T_mesh, t_mesh, response_mesh, cmap='viridis', alpha=0.8)
ax1.scatter(X_design[:, 0], X_design[:, 1], y_response, c='red', s=100, 
            marker='o', edgecolors='white', linewidth=2, label='実験点')
ax1.scatter(optimal_temp, optimal_time, optimal_response, c='gold', s=300, 
            marker='*', edgecolors='red', linewidth=2, label='最適点')
ax1.set_xlabel('温度 (°C)', fontsize=11)
ax1.set_ylabel(' hours ( minutes)', fontsize=11)
ax1.set_zlabel('品質スコア', fontsize=11)
ax1.set_title(f'応答曲面（pH={ph_fixed}固定）', fontsize=13, fontweight='bold')
plt.colorbar(surf, ax=ax1, shrink=0.5, label='品質スコア')

# 2. 等高線図
ax2 = fig.add_subplot(2, 2, 2)
contour = ax2.contourf(T_mesh, t_mesh, response_mesh, levels=15, cmap='viridis', alpha=0.9)
contour_lines = ax2.contour(T_mesh, t_mesh, response_mesh, levels=10, colors='white', 
                             linewidths=0.5, alpha=0.7)
ax2.clabel(contour_lines, inline=True, fontsize=8, fmt='%.1f')
ax2.scatter(X_design[:, 0], X_design[:, 1], c='red', s=100, 
            marker='o', edgecolors='white', linewidth=2, label='実験点')
ax2.scatter(optimal_temp, optimal_time, c='gold', s=300, 
            marker='*', edgecolors='red', linewidth=2, label='最適点')
ax2.set_xlabel('温度 (°C)', fontsize=12)
ax2.set_ylabel(' hours ( minutes)', fontsize=12)
ax2.set_title(f'等高線図（pH={ph_fixed}固定）', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
plt.colorbar(contour, ax=ax2, label='品質スコア')

# 3. Predictionvs実測値
ax3 = fig.add_subplot(2, 2, 3)
ax3.scatter(y_response, y_pred, c='#11998e', s=100, alpha=0.7, edgecolors='white', linewidth=1.5)
ax3.plot([y_response.min(), y_response.max()], [y_response.min(), y_response.max()],
         'r--', linewidth=2, label='理想ライン')
ax3.set_xlabel('実測値', fontsize=12)
ax3.set_ylabel('Prediction値', fontsize=12)
ax3.set_title(f'Prediction精度（R²={r2_score:.4f}）', fontsize=13, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 4. 主効果プロット（温度の効果）
ax4 = fig.add_subplot(2, 2, 4)
temp_test = np.linspace(80, 100, 100)
time_fixed = 20
ph_test_fixed = 4.0

response_temp_effect = []
for t in temp_test:
    X_test = np.array([[t, time_fixed, ph_test_fixed]])
    X_test_poly = poly.transform(X_test)
    response_temp_effect.append(model.predict(X_test_poly)[0])

ax4.plot(temp_test, response_temp_effect, linewidth=2.5, color='#11998e', label='Prediction応答')
ax4.scatter(X_design[:, 0], y_response, c='red', s=80, alpha=0.6, 
            edgecolors='white', linewidth=1.5, label='実験点')
ax4.axvline(x=optimal_temp, color='green', linestyle='--', linewidth=2, label='最適温度')
ax4.set_xlabel('温度 (°C)', fontsize=12)
ax4.set_ylabel('品質スコア', fontsize=12)
ax4.set_title(f'主効果プロット（ hours={time_fixed} minutes, pH={ph_test_fixed}固定）', 
              fontsize=13, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

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