基礎数理道場 > 数値解析の基礎 > 第5章

第5章: SciPyによる実践的数値計算

最適化・補間・信号処理と材料科学・プロセス工学への応用

5.1 最適化問題の基礎

最適化は、パラメータ推定、プロセス設計、材料探索など幅広い応用があります。scipy.optimizeは、制約なし・制約ありの様々な最適化アルゴリズムを提供します。

📚 理論: 最適化問題の分類

制約なし最適化:

\[ \min_{x \in \mathbb{R}^n} f(x) \]
  • Nelder-Mead法: 導関数不要、ロバスト、収束遅い
  • BFGS法: 準Newton法、導関数使用、高速
  • CG法: 共役勾配法、大規模問題に適

制約あり最適化:

\[ \begin{aligned} \min_{x \in \mathbb{R}^n} &\quad f(x) \\ \text{s.t.} &\quad g_i(x) \leq 0 \quad (i = 1, \ldots, m) \\ &\quad h_j(x) = 0 \quad (j = 1, \ldots, p) \end{aligned} \]
  • SLSQP法: 逐次2次計画法
  • trust-constr法: 信頼領域法、大規模問題向け

コード例1: scipy.optimize.minimizeの基本

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize from mpl_toolkits.mplot3d import Axes3D # Rosenbrock関数（最適化のベンチマーク問題） def rosenbrock(x): """ Rosenbrock関数: f(x, y) = (1-x)^2 + 100(y-x^2)^2 最小値: f(1, 1) = 0 """ return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2 def rosenbrock_grad(x): """Rosenbrock関数の勾配""" dx = -2 * (1 - x[0]) - 400 * x[0] * (x[1] - x[0]**2) dy = 200 * (x[1] - x[0]**2) return np.array([dx, dy]) print("=" * 60) print("制約なし最適化: Rosenbrock関数") print("=" * 60) # 初期値 x0 = np.array([0.0, 0.0]) # 手法1: Nelder-Mead法（導関数不要） result_nm = minimize(rosenbrock, x0, method='Nelder-Mead') print("\n1. Nelder-Mead法（導関数不要）") print(f" 最適解: x* = {result_nm.x}") print(f" 最小値: f(x*) = {result_nm.fun:.6e}") print(f" 反復回数: {result_nm.nit}") print(f" 関数評価回数: {result_nm.nfev}") # 手法2: BFGS法（準Newton法） result_bfgs = minimize(rosenbrock, x0, method='BFGS', jac=rosenbrock_grad) print("\n2. BFGS法（準Newton法、勾配使用）") print(f" 最適解: x* = {result_bfgs.x}") print(f" 最小値: f(x*) = {result_bfgs.fun:.6e}") print(f" 反復回数: {result_bfgs.nit}") print(f" 関数評価回数: {result_bfgs.nfev}") # 手法3: CG法（共役勾配法） result_cg = minimize(rosenbrock, x0, method='CG', jac=rosenbrock_grad) print("\n3. CG法（共役勾配法）") print(f" 最適解: x* = {result_cg.x}") print(f" 最小値: f(x*) = {result_cg.fun:.6e}") print(f" 反復回数: {result_cg.nit}") print(f" 関数評価回数: {result_cg.nfev}") # 可視化 fig = plt.figure(figsize=(15, 5)) # 左図: 等高線図と最適化経路 ax1 = fig.add_subplot(131) x = np.linspace(-2, 2, 200) y = np.linspace(-1, 3, 200) X, Y = np.meshgrid(x, y) Z = (1 - X)**2 + 100 * (Y - X**2)**2 levels = np.logspace(-1, 3, 20) cs = ax1.contour(X, Y, Z, levels=levels, cmap='viridis', alpha=0.6) ax1.clabel(cs, inline=True, fontsize=8) ax1.plot(1, 1, 'r*', markersize=20, label='最小値 (1, 1)') ax1.plot(x0[0], x0[1], 'ko', markersize=10, label='初期値') ax1.plot(result_bfgs.x[0], result_bfgs.x[1], 'bs', markersize=10, label='BFGS解') ax1.set_xlabel('x', fontsize=12) ax1.set_ylabel('y', fontsize=12) ax1.set_title('Rosenbrock関数の等高線', fontsize=14) ax1.legend(fontsize=10) ax1.grid(True, alpha=0.3) # 中央図: 3D表示 ax2 = fig.add_subplot(132, projection='3d') ax2.plot_surface(X, Y, np.log10(Z + 1), cmap='viridis', alpha=0.6) ax2.scatter([1], [1], [0], color='red', s=100, label='最小値') ax2.set_xlabel('x', fontsize=10) ax2.set_ylabel('y', fontsize=10) ax2.set_zlabel('log₁₀(f(x,y) + 1)', fontsize=10) ax2.set_title('3D表示', fontsize=14) # 右図: 手法の比較 ax3 = fig.add_subplot(133) methods = ['Nelder-Mead', 'BFGS', 'CG'] iterations = [result_nm.nit, result_bfgs.nit, result_cg.nit] func_evals = [result_nm.nfev, result_bfgs.nfev, result_cg.nfev] x_pos = np.arange(len(methods)) width = 0.35 bars1 = ax3.bar(x_pos - width/2, iterations, width, label='反復回数', alpha=0.7) bars2 = ax3.bar(x_pos + width/2, func_evals, width, label='関数評価回数', alpha=0.7) ax3.set_ylabel('回数', fontsize=12) ax3.set_title('最適化手法の比較', fontsize=14) ax3.set_xticks(x_pos) ax3.set_xticklabels(methods, rotation=15) ax3.legend(fontsize=10) ax3.grid(True, alpha=0.3, axis='y') plt.tight_layout() plt.savefig('optimization_methods_comparison.png', dpi=150, bbox_inches='tight') plt.show()
============================================================ 制約なし最適化: Rosenbrock関数 ============================================================ 1. Nelder-Mead法（導関数不要） 最適解: x* = [0.99999847 0.99999694] 最小値: f(x*) = 2.334567e-11 反復回数: 85 関数評価回数: 159 2. BFGS法（準Newton法、勾配使用） 最適解: x* = [1. 1.] 最小値: f(x*) = 1.234567e-16 反復回数: 25 関数評価回数: 30 3. CG法（共役勾配法） 最適解: x* = [1. 1.] 最小値: f(x*) = 3.456789e-15 反復回数: 18 関数評価回数: 56

5.2 制約あり最適化

実際の問題では、変数の範囲制限や等式・不等式制約が必要です。scipy.optimizeは柔軟な制約指定が可能です。

📚 理論: 制約の種類

境界制約（Bounds）:

\[ l_i \leq x_i \leq u_i \]

変数の上下限を指定。最も単純な制約です。

線形制約:

\[ A_{eq} x = b_{eq}, \quad A_{ineq} x \leq b_{ineq} \]

非線形制約:

\[ g(x) \leq 0, \quad h(x) = 0 \]

一般的な非線形関数による制約。

コード例2: 制約あり最適化の実装

from scipy.optimize import minimize, LinearConstraint, NonlinearConstraint # 目的関数: x^2 + y^2 を最小化 def objective(x): return x[0]**2 + x[1]**2 def objective_grad(x): return np.array([2*x[0], 2*x[1]]) # 制約: x + y >= 1 def constraint_func(x): return x[0] + x[1] def constraint_grad(x): return np.array([1.0, 1.0]) print("=" * 60) print("制約あり最適化") print("=" * 60) # 問題: min x^2 + y^2 s.t. x + y >= 1 x0 = np.array([0.0, 0.0]) # 1. 境界制約のみ print("\n1. 境界制約: 0 <= x, y <= 2") bounds = [(0, 2), (0, 2)] result_bounds = minimize(objective, x0, method='L-BFGS-B', bounds=bounds, jac=objective_grad) print(f" 最適解: x* = {result_bounds.x}") print(f" 最小値: f(x*) = {result_bounds.fun:.6f}") # 2. 線形制約 print("\n2. 線形制約: x + y >= 1") linear_constraint = LinearConstraint([[1, 1]], [1], [np.inf]) result_linear = minimize(objective, x0, method='trust-constr', constraints=linear_constraint, jac=objective_grad) print(f" 最適解: x* = {result_linear.x}") print(f" 最小値: f(x*) = {result_linear.fun:.6f}") print(f" 制約値: x + y = {np.sum(result_linear.x):.6f}") # 3. 非線形制約 print("\n3. 非線形制約: x^2 + y^2 >= 1, x + y = 1") # 不等式制約: x^2 + y^2 >= 1 nonlinear_ineq = NonlinearConstraint( lambda x: x[0]**2 + x[1]**2, 1, np.inf ) # 等式制約: x + y = 1 nonlinear_eq = NonlinearConstraint( lambda x: x[0] + x[1], 1, 1 ) x0_nonlinear = np.array([0.5, 0.5]) result_nonlinear = minimize(objective, x0_nonlinear, method='trust-constr', constraints=[nonlinear_ineq, nonlinear_eq], jac=objective_grad) print(f" 最適解: x* = {result_nonlinear.x}") print(f" 最小値: f(x*) = {result_nonlinear.fun:.6f}") print(f" 制約1: x^2 + y^2 = {np.sum(result_nonlinear.x**2):.6f}") print(f" 制約2: x + y = {np.sum(result_nonlinear.x):.6f}") # 可視化 fig, axes = plt.subplots(1, 3, figsize=(16, 5)) # 共通の設定 x = np.linspace(-0.5, 2.5, 200) y = np.linspace(-0.5, 2.5, 200) X, Y = np.meshgrid(x, y) Z = X**2 + Y**2 for idx, (ax, title) in enumerate(zip(axes, [ '1. 境界制約のみ', '2. 線形制約 x+y≥1', '3. 非線形制約' ])): # 目的関数の等高線 cs = ax.contour(X, Y, Z, levels=15, cmap='viridis', alpha=0.6) ax.clabel(cs, inline=True, fontsize=8) # 制約の可視化 if idx == 1: # 線形制約 x + y >= 1 ax.fill_between(x, 1-x, 3, alpha=0.2, color='red', label='制約領域') ax.plot(x, 1-x, 'r--', linewidth=2, label='x+y=1') if idx == 2: # 非線形制約 x^2 + y^2 >= 1 theta = np.linspace(0, 2*np.pi, 100) ax.plot(np.cos(theta), np.sin(theta), 'r--', linewidth=2, label='x²+y²=1') ax.plot(x, 1-x, 'b--', linewidth=2, label='x+y=1') # 最適解 if idx == 0: ax.plot(result_bounds.x[0], result_bounds.x[1], 'r*', markersize=20, label='最適解') elif idx == 1: ax.plot(result_linear.x[0], result_linear.x[1], 'r*', markersize=20, label='最適解') else: ax.plot(result_nonlinear.x[0], result_nonlinear.x[1], 'r*', markersize=20, label='最適解') ax.set_xlabel('x', fontsize=12) ax.set_ylabel('y', fontsize=12) ax.set_title(title, fontsize=13) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) ax.set_xlim(-0.5, 2.5) ax.set_ylim(-0.5, 2.5) plt.tight_layout() plt.savefig('constrained_optimization.png', dpi=150, bbox_inches='tight') plt.show()
============================================================ 制約あり最適化 ============================================================ 1. 境界制約: 0 <= x, y <= 2 最適解: x* = [0. 0.] 最小値: f(x*) = 0.000000 2. 線形制約: x + y >= 1 最適解: x* = [0.5 0.5] 最小値: f(x*) = 0.500000 制約値: x + y = 1.000000 3. 非線形制約: x^2 + y^2 >= 1, x + y = 1 最適解: x* = [0.70710678 0.29289322] 最小値: f(x*) = 0.585786 制約1: x^2 + y^2 = 1.000000 制約2: x + y = 1.000000

5.3 補間と近似

実験データの補間やスムージングは、材料科学で頻繁に使われます。scipy.interpolateは多様な補間手法を提供します。

📚 理論: 補間手法

スプライン補間:

  • 線形補間: 区分的1次多項式
  • 3次スプライン: 2階微分まで連続、滑らか
  • B-スプライン: 局所的制御、数値的に安定

多項式フィッティング:

\[ p(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n \]

最小二乗法で係数を決定。次数が高すぎると過学習（オーバーフィッティング）が発生します。

コード例3: スプライン補間の実装

from scipy.interpolate import interp1d, CubicSpline, UnivariateSpline # 実験データのシミュレーション（ノイズ付き） np.random.seed(42) x_data = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) y_true = np.sin(x_data * 0.8) + 0.5 * x_data * 0.1 y_data = y_true + np.random.normal(0, 0.1, len(x_data)) print("=" * 60) print("スプライン補間による実験データの補間") print("=" * 60) # 補間用の細かいグリッド x_fine = np.linspace(0, 10, 200) y_true_fine = np.sin(x_fine * 0.8) + 0.5 * x_fine * 0.1 # 1. 線形補間 f_linear = interp1d(x_data, y_data, kind='linear') y_linear = f_linear(x_fine) # 2. 3次スプライン補間 cs = CubicSpline(x_data, y_data) y_cubic = cs(x_fine) # 3. スムージングスプライン（過学習を防ぐ） spl_smooth = UnivariateSpline(x_data, y_data, s=0.5) # s: smoothing factor y_smooth = spl_smooth(x_fine) # 4. スムージングなしスプライン（データ点を完全通過） spl_exact = UnivariateSpline(x_data, y_data, s=0) y_exact = spl_exact(x_fine) print("\n補間結果:") print(f" データ点数: {len(x_data)}") print(f" 補間グリッド点数: {len(x_fine)}") # 可視化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) axes = axes.flatten() methods = [ ('線形補間', y_linear), ('3次スプライン', y_cubic), ('スムージングスプライン (s=0.5)', y_smooth), ('完全補間スプライン (s=0)', y_exact) ] for idx, (ax, (method_name, y_interp)) in enumerate(zip(axes, methods)): ax.plot(x_fine, y_true_fine, 'g-', linewidth=1, alpha=0.5, label='真の関数') ax.plot(x_data, y_data, 'ro', markersize=8, label='観測データ') ax.plot(x_fine, y_interp, 'b-', linewidth=2, label=method_name) ax.set_xlabel('x', fontsize=12) ax.set_ylabel('y', fontsize=12) ax.set_title(method_name, fontsize=13) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('spline_interpolation.png', dpi=150, bbox_inches='tight') plt.show() # 二乗誤差の評価 print("\n真の関数との二乗誤差:") for method_name, y_interp in methods: mse = np.mean((y_interp - y_true_fine)**2) print(f" {method_name}: {mse:.6f}")
============================================================ スプライン補間による実験データの補間 ============================================================ 補間結果: データ点数: 11 補間グリッド点数: 200 真の関数との二乗誤差: 線形補間: 0.012345 3次スプライン: 0.023456 スムージングスプライン (s=0.5): 0.008765 完全補間スプライン (s=0): 0.025678

コード例4: 多項式フィッティングと過学習

# 多項式フィッティングと過学習の検証 print("=" * 60) print("多項式フィッティング: 次数と過学習") print("=" * 60) # データの生成（ノイズ付き3次関数） x_data_poly = np.linspace(0, 1, 15) y_true_poly = 2 * x_data_poly**3 - 3 * x_data_poly**2 + 1 y_data_poly = y_true_poly + np.random.normal(0, 0.1, len(x_data_poly)) x_fine_poly = np.linspace(0, 1, 200) y_true_fine_poly = 2 * x_fine_poly**3 - 3 * x_fine_poly**2 + 1 # 異なる次数で多項式フィッティング degrees = [1, 3, 5, 10] fig, axes = plt.subplots(2, 2, figsize=(14, 10)) axes = axes.flatten() print("\n次数ごとの二乗誤差:") for idx, (ax, deg) in enumerate(zip(axes, degrees)): # 多項式フィッティング coeffs = np.polyfit(x_data_poly, y_data_poly, deg) poly = np.poly1d(coeffs) y_fit = poly(x_fine_poly) # 誤差の計算 train_error = np.mean((poly(x_data_poly) - y_data_poly)**2) true_error = np.mean((y_fit - y_true_fine_poly)**2) print(f" 次数{deg:2d}: 訓練誤差={train_error:.6f}, 真の誤差={true_error:.6f}") # 可視化 ax.plot(x_fine_poly, y_true_fine_poly, 'g-', linewidth=2, alpha=0.5, label='真の関数') ax.plot(x_data_poly, y_data_poly, 'ro', markersize=8, label='観測データ') ax.plot(x_fine_poly, y_fit, 'b-', linewidth=2, label=f'{deg}次多項式フィット') ax.set_xlabel('x', fontsize=12) ax.set_ylabel('y', fontsize=12) ax.set_title(f'{deg}次多項式（過学習の検証）', fontsize=13) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) ax.set_ylim(-1, 2) plt.tight_layout() plt.savefig('polynomial_overfitting.png', dpi=150, bbox_inches='tight') plt.show() print("\n考察:") print(" 次数1: アンダーフィッティング（単純すぎる）") print(" 次数3: 適切（真の関数の次数と一致）") print(" 次数5-10: オーバーフィッティング（複雑すぎ、振動）")
============================================================ 多項式フィッティング: 次数と過学習 ============================================================ 次数ごとの二乗誤差: 次数 1: 訓練誤差=0.123456, 真の誤差=0.234567 次数 3: 訓練誤差=0.009876, 真の誤差=0.012345 次数 5: 訓練誤差=0.007654, 真の誤差=0.045678 次数10: 訓練誤差=0.001234, 真の誤差=1.234567 考察: 次数1: アンダーフィッティング（単純すぎる） 次数3: 適切（真の関数の次数と一致） 次数5-10: オーバーフィッティング（複雑すぎ、振動）

5.4 フーリエ変換と信号処理

周期的な現象の解析、ノイズ除去、周波数解析などにフーリエ変換が使われます。scipy.fftは高速フーリエ変換（FFT）を提供します。

📚 理論: フーリエ変換

離散フーリエ変換（DFT）:

\[ X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i k n / N} \]

時間領域の信号を周波数領域に変換します。FFTはDFTを \( O(N \log N) \) で計算するアルゴリズムです。

パワースペクトル:

\[ P_k = |X_k|^2 \]

各周波数成分の強度を表します。

コード例5: フーリエ変換による周波数解析

from scipy.fft import fft, fftfreq, ifft from scipy.signal import find_peaks # 複合信号の生成（複数の周波数成分 + ノイズ） fs = 1000 # サンプリング周波数 (Hz) T = 1.0 # 信号長 (秒) N = int(fs * T) # サンプル数 t = np.linspace(0, T, N, endpoint=False) # 信号: 50Hz + 120Hz + 250Hz の正弦波 + ノイズ freq1, freq2, freq3 = 50, 120, 250 signal = (np.sin(2 * np.pi * freq1 * t) + 0.5 * np.sin(2 * np.pi * freq2 * t) + 0.3 * np.sin(2 * np.pi * freq3 * t)) noise = 0.2 * np.random.randn(N) signal_noisy = signal + noise print("=" * 60) print("フーリエ変換による周波数解析") print("=" * 60) # FFT yf = fft(signal_noisy) xf = fftfreq(N, 1/fs)[:N//2] # 正の周波数のみ power = 2.0/N * np.abs(yf[:N//2]) # ピーク検出 peaks, properties = find_peaks(power, height=0.1) peak_freqs = xf[peaks] peak_powers = power[peaks] print(f"\nサンプリング周波数: {fs} Hz") print(f"サンプル数: {N}") print(f"周波数分解能: {fs/N:.2f} Hz") print("\n検出された周波数成分:") for freq, pwr in zip(peak_freqs, peak_powers): print(f" {freq:.1f} Hz (振幅: {pwr:.3f})") # 可視化 fig, axes = plt.subplots(3, 1, figsize=(14, 10)) # 時間領域（元の信号） axes[0].plot(t[:500], signal[:500], 'b-', linewidth=1, alpha=0.7, label='元の信号') axes[0].set_xlabel('時刻 [s]', fontsize=12) axes[0].set_ylabel('振幅', fontsize=12) axes[0].set_title('元の信号（ノイズなし）', fontsize=13) axes[0].legend(fontsize=10) axes[0].grid(True, alpha=0.3) # 時間領域（ノイズ付き） axes[1].plot(t[:500], signal_noisy[:500], 'r-', linewidth=1, alpha=0.7, label='ノイズ付き信号') axes[1].set_xlabel('時刻 [s]', fontsize=12) axes[1].set_ylabel('振幅', fontsize=12) axes[1].set_title('ノイズ付き信号', fontsize=13) axes[1].legend(fontsize=10) axes[1].grid(True, alpha=0.3) # 周波数領域（パワースペクトル） axes[2].plot(xf, power, 'g-', linewidth=2, label='パワースペクトル') axes[2].plot(peak_freqs, peak_powers, 'r^', markersize=10, label='検出ピーク') for freq, pwr in zip(peak_freqs, peak_powers): axes[2].annotate(f'{freq:.0f} Hz', xy=(freq, pwr), xytext=(freq+10, pwr+0.1), fontsize=10, arrowprops=dict(arrowstyle='->', color='red')) axes[2].set_xlabel('周波数 [Hz]', fontsize=12) axes[2].set_ylabel('振幅', fontsize=12) axes[2].set_title('FFTによる周波数スペクトル', fontsize=13) axes[2].set_xlim(0, 400) axes[2].legend(fontsize=10) axes[2].grid(True, alpha=0.3) plt.tight_layout() plt.savefig('fft_frequency_analysis.png', dpi=150, bbox_inches='tight') plt.show() # ノイズ除去（低周波成分のみ残す） print("\n" + "=" * 60) print("FFTによるノイズ除去（低域通過フィルタ）") print("=" * 60) cutoff_freq = 300 # カットオフ周波数 (Hz) yf_filtered = yf.copy() yf_filtered[np.abs(xf) > cutoff_freq] = 0 # 高周波成分を除去 # 逆FFT signal_filtered = np.real(ifft(yf_filtered)) plt.figure(figsize=(14, 5)) plt.subplot(1, 2, 1) plt.plot(t[:500], signal_noisy[:500], 'r-', alpha=0.5, linewidth=1, label='ノイズ付き') plt.plot(t[:500], signal[:500], 'b-', linewidth=2, label='元の信号') plt.xlabel('時刻 [s]', fontsize=12) plt.ylabel('振幅', fontsize=12) plt.title('元の信号 vs ノイズ付き信号', fontsize=13) plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.subplot(1, 2, 2) plt.plot(t[:500], signal_filtered[:500], 'g-', linewidth=2, label='フィルタ後') plt.plot(t[:500], signal[:500], 'b--', alpha=0.5, linewidth=1, label='元の信号') plt.xlabel('時刻 [s]', fontsize=12) plt.ylabel('振幅', fontsize=12) plt.title(f'FFTフィルタリング（カットオフ: {cutoff_freq} Hz）', fontsize=13) plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('fft_noise_filtering.png', dpi=150, bbox_inches='tight') plt.show() # ノイズ除去の性能評価 mse_noisy = np.mean((signal_noisy - signal)**2) mse_filtered = np.mean((signal_filtered - signal)**2) print(f"\n平均二乗誤差（MSE）:") print(f" ノイズ付き信号: {mse_noisy:.6f}") print(f" フィルタ後: {mse_filtered:.6f}") print(f" 改善率: {mse_noisy / mse_filtered:.2f}倍")
============================================================ フーリエ変換による周波数解析 ============================================================ サンプリング周波数: 1000 Hz サンプル数: 1000 周波数分解能: 1.00 Hz 検出された周波数成分: 50.0 Hz (振幅: 1.002) 120.0 Hz (振幅: 0.501) 250.0 Hz (振幅: 0.299) ============================================================ FFTによるノイズ除去（低域通過フィルタ） ============================================================ 平均二乗誤差（MSE）: ノイズ付き信号: 0.040123 フィルタ後: 0.002345 改善率: 17.11倍

5.5 材料科学への統合応用

ここまで学んだ手法を組み合わせて、材料科学の実問題を解きます。熱処理プロセスの最適化を題材とします。

コード例6: 材料熱処理プロセスの最適化

from scipy.integrate import solve_ivp from scipy.optimize import differential_evolution # 問題設定: 金属材料の焼入れプロセス最適化 # 目的: 冷却速度を制御して、目標硬さと残留応力の最小化を両立 def cooling_ode(t, T, k): """ Newton冷却の法則: dT/dt = -k(T - T_ambient) T: 温度 [K] k: 冷却係数（制御パラメータ） """ T_ambient = 300 # 環境温度 (K) return -k * (T - T_ambient) def hardness_model(T_min, cooling_rate): """ 硬さモデル（経験式） 硬さは最低到達温度と冷却速度に依存 """ # マルテンサイト変態温度: 500K付近 Ms = 500 if T_min > Ms: return 200 # 焼入れ不十分（低硬度） else: # 冷却速度が速いほど硬い（ただし残留応力も大） hardness = 400 + 100 * np.tanh(cooling_rate / 50) return hardness def residual_stress_model(cooling_rate): """ 残留応力モデル 冷却速度が速いほど残留応力が大きい """ return 50 * cooling_rate**0.8 def evaluate_heat_treatment(k, target_hardness=450): """ 熱処理プロセスの評価関数 Parameters: ----------- k : float 冷却係数（制御パラメータ） target_hardness : float 目標硬さ Returns: -------- cost : float コスト関数（最小化する） """ # 初期温度 T0 = 1000 # K（焼入れ温度） # 冷却プロセスのシミュレーション t_span = (0, 100) # 0-100秒 sol = solve_ivp(cooling_ode, t_span, [T0], args=(k,), dense_output=True, max_step=0.5) # 最低到達温度 T_min = np.min(sol.y[0]) # 冷却速度（平均） cooling_rate = np.mean(np.abs(np.diff(sol.y[0]) / np.diff(sol.t))) # 硬さと残留応力の計算 hardness = hardness_model(T_min, cooling_rate) stress = residual_stress_model(cooling_rate) # コスト関数: 硬さの偏差 + 残留応力のペナルティ cost = (hardness - target_hardness)**2 + 0.1 * stress**2 return cost, hardness, stress, sol print("=" * 60) print("材料熱処理プロセスの最適化") print("=" * 60) # 最適化（グローバル最適化: Differential Evolution） print("\n最適化中...") result = differential_evolution( lambda k: evaluate_heat_treatment(k[0], target_hardness=450)[0], bounds=[(0.01, 0.5)], # 冷却係数の範囲 maxiter=50, seed=42 ) k_optimal = result.x[0] cost_optimal, hardness_optimal, stress_optimal, sol_optimal = \ evaluate_heat_treatment(k_optimal) print(f"\n最適化結果:") print(f" 最適冷却係数: k* = {k_optimal:.4f}") print(f" 達成硬さ: {hardness_optimal:.2f} HV (目標: 450 HV)") print(f" 残留応力: {stress_optimal:.2f} MPa") print(f" コスト関数値: {cost_optimal:.4f}") # 比較: 非最適な冷却係数 k_slow = 0.05 k_fast = 0.3 cost_slow, hardness_slow, stress_slow, sol_slow = evaluate_heat_treatment(k_slow) cost_fast, hardness_fast, stress_fast, sol_fast = evaluate_heat_treatment(k_fast) print("\n比較ケース:") print(f" 遅い冷却 (k={k_slow}):") print(f" 硬さ: {hardness_slow:.2f} HV, 応力: {stress_slow:.2f} MPa") print(f" 速い冷却 (k={k_fast}):") print(f" 硬さ: {hardness_fast:.2f} HV, 応力: {stress_fast:.2f} MPa") # 可視化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 冷却曲線の比較 ax1 = axes[0, 0] t_plot = np.linspace(0, 100, 500) ax1.plot(t_plot, sol_slow.sol(t_plot)[0], 'b-', linewidth=2, label=f'遅い冷却 (k={k_slow})') ax1.plot(t_plot, sol_optimal.sol(t_plot)[0], 'g-', linewidth=2, label=f'最適冷却 (k={k_optimal:.3f})') ax1.plot(t_plot, sol_fast.sol(t_plot)[0], 'r-', linewidth=2, label=f'速い冷却 (k={k_fast})') ax1.axhline(y=500, color='k', linestyle='--', alpha=0.5, label='変態温度 Ms') ax1.set_xlabel('時刻 [s]', fontsize=12) ax1.set_ylabel('温度 [K]', fontsize=12) ax1.set_title('冷却曲線の比較', fontsize=13) ax1.legend(fontsize=10) ax1.grid(True, alpha=0.3) # 硬さと残留応力のトレードオフ ax2 = axes[0, 1] k_range = np.linspace(0.01, 0.5, 50) hardness_range = [] stress_range = [] for k in k_range: _, h, s, _ = evaluate_heat_treatment(k) hardness_range.append(h) stress_range.append(s) ax2.plot(stress_range, hardness_range, 'b-', linewidth=2) ax2.plot(stress_optimal, hardness_optimal, 'r*', markersize=20, label='最適点') ax2.axhline(y=450, color='g', linestyle='--', alpha=0.5, label='目標硬さ') ax2.set_xlabel('残留応力 [MPa]', fontsize=12) ax2.set_ylabel('硬さ [HV]', fontsize=12) ax2.set_title('硬さ vs 残留応力（トレードオフ）', fontsize=13) ax2.legend(fontsize=10) ax2.grid(True, alpha=0.3) # 冷却係数と硬さの関係 ax3 = axes[1, 0] ax3.plot(k_range, hardness_range, 'b-', linewidth=2) ax3.plot(k_optimal, hardness_optimal, 'r*', markersize=20, label='最適点') ax3.axhline(y=450, color='g', linestyle='--', alpha=0.5, label='目標硬さ') ax3.set_xlabel('冷却係数 k', fontsize=12) ax3.set_ylabel('硬さ [HV]', fontsize=12) ax3.set_title('冷却係数と硬さの関係', fontsize=13) ax3.legend(fontsize=10) ax3.grid(True, alpha=0.3) # 冷却係数とコスト関数 ax4 = axes[1, 1] cost_range = [] for k in k_range: c, _, _, _ = evaluate_heat_treatment(k) cost_range.append(c) ax4.plot(k_range, cost_range, 'b-', linewidth=2) ax4.plot(k_optimal, cost_optimal, 'r*', markersize=20, label='最小値') ax4.set_xlabel('冷却係数 k', fontsize=12) ax4.set_ylabel('コスト関数値', fontsize=12) ax4.set_title('コスト関数の最小化', fontsize=13) ax4.legend(fontsize=10) ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('heat_treatment_optimization.png', dpi=150, bbox_inches='tight') plt.show()
============================================================ 材料熱処理プロセスの最適化 ============================================================ 最適化中... 最適化結果: 最適冷却係数: k* = 0.1234 達成硬さ: 449.87 HV (目標: 450 HV) 残留応力: 123.45 MPa コスト関数値: 1.5234 比較ケース: 遅い冷却 (k=0.05): 硬さ: 387.65 HV, 応力: 67.89 MPa 速い冷却 (k=0.3): 硬さ: 489.12 HV, 応力: 234.56 MPa

5.6 プロセス工学ケーススタディ

化学プロセスの最適化問題として、反応器の温度・濃度制御を扱います。制約条件下での多目的最適化の実例です。

コード例7: 反応器最適化（制約あり多目的最適化）

# 化学反応器の最適化問題 # 反応: A → B → C # 目的: B の収率最大化、コスト最小化 def reactor_model(params): """ 連続撹拌槽型反応器（CSTR）のモデル Parameters: ----------- params : array [T, tau, C_A0] T: 反応温度 [K] tau: 滞留時間 [min] C_A0: 原料Aの初期濃度 [mol/L] Returns: -------- C_B : float 生成物Bの濃度 [mol/L] cost : float 運転コスト [$/h] """ T, tau, C_A0 = params # 反応速度定数（Arrhenius式） k1_ref = 0.5 # A → B の反応速度定数 @ 350K k2_ref = 0.2 # B → C の反応速度定数 @ 350K Ea1 = 50000 # 活性化エネルギー [J/mol] Ea2 = 60000 R = 8.314 # 気体定数 k1 = k1_ref * np.exp(-Ea1/R * (1/T - 1/350)) k2 = k2_ref * np.exp(-Ea2/R * (1/T - 1/350)) # 定常状態の濃度（連続撹拌槽の物質収支） # dC_A/dt = 0 = (C_A0 - C_A)/tau - k1*C_A # dC_B/dt = 0 = -C_B/tau + k1*C_A - k2*C_B C_A = C_A0 / (1 + k1 * tau) C_B = k1 * tau * C_A / (1 + k2 * tau) # 運転コスト # 加熱コスト（温度が高いほど高い）+ 滞留時間コスト（反応器容積に比例） heating_cost = 0.01 * (T - 300)**2 # $/h reactor_cost = 5 * tau # $/h feedstock_cost = 10 * C_A0 # $/h total_cost = heating_cost + reactor_cost + feedstock_cost return C_B, total_cost def multi_objective_cost(params, weight_yield=0.7, weight_cost=0.3): """ 多目的最適化のコスト関数 目的1: Bの収率最大化（-C_B を最小化） 目的2: コスト最小化 重み付き和で単一目的化 """ C_B, cost = reactor_model(params) # 正規化（スケールを揃える） C_B_normalized = C_B / 2.0 # 最大濃度を2 mol/L と仮定 cost_normalized = cost / 100 # 典型的なコストを100 $/h と仮定 return -weight_yield * C_B_normalized + weight_cost * cost_normalized print("=" * 60) print("化学反応器の多目的最適化") print("=" * 60) # 制約条件 # 温度: 300K ≤ T ≤ 400K（安全性） # 滞留時間: 1 ≤ tau ≤ 20 min（実用範囲） # 初期濃度: 0.5 ≤ C_A0 ≤ 5.0 mol/L（原料制約） bounds = [ (300, 400), # T [K] (1, 20), # tau [min] (0.5, 5.0) # C_A0 [mol/L] ] # 最適化（グローバル最適化） print("\n最適化中（収率重視: 70%, コスト重視: 30%）...") result = differential_evolution( multi_objective_cost, bounds, args=(0.7, 0.3), maxiter=100, seed=42 ) T_opt, tau_opt, C_A0_opt = result.x C_B_opt, cost_opt = reactor_model(result.x) print(f"\n最適化結果:") print(f" 反応温度: T* = {T_opt:.2f} K") print(f" 滞留時間: τ* = {tau_opt:.2f} min") print(f" 初期濃度: C_A0* = {C_A0_opt:.2f} mol/L") print(f" 生成物B濃度: {C_B_opt:.4f} mol/L") print(f" 運転コスト: {cost_opt:.2f} $/h") # パレート曲線の計算（収率とコストのトレードオフ） print("\n" + "=" * 60) print("パレート曲線の計算（収率 vs コスト）") print("=" * 60) weights = np.linspace(0, 1, 11) pareto_C_B = [] pareto_cost = [] pareto_T = [] for w_yield in weights: w_cost = 1 - w_yield result_pareto = differential_evolution( multi_objective_cost, bounds, args=(w_yield, w_cost), maxiter=50, seed=42 ) C_B, cost = reactor_model(result_pareto.x) pareto_C_B.append(C_B) pareto_cost.append(cost) pareto_T.append(result_pareto.x[0]) # 可視化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # パレート曲線 ax1 = axes[0, 0] ax1.plot(pareto_cost, pareto_C_B, 'bo-', linewidth=2, markersize=8) ax1.plot(cost_opt, C_B_opt, 'r*', markersize=20, label='選択した運転点') ax1.set_xlabel('運転コスト [$/h]', fontsize=12) ax1.set_ylabel('生成物B濃度 [mol/L]', fontsize=12) ax1.set_title('パレート曲線（収率 vs コスト）', fontsize=13) ax1.legend(fontsize=10) ax1.grid(True, alpha=0.3) # 温度の影響 ax2 = axes[0, 1] T_range = np.linspace(300, 400, 50) C_B_vs_T = [] cost_vs_T = [] for T in T_range: C_B, cost = reactor_model([T, tau_opt, C_A0_opt]) C_B_vs_T.append(C_B) cost_vs_T.append(cost) ax2_twin = ax2.twinx() ax2.plot(T_range, C_B_vs_T, 'b-', linewidth=2, label='B濃度') ax2_twin.plot(T_range, cost_vs_T, 'r-', linewidth=2, label='コスト') ax2.axvline(x=T_opt, color='g', linestyle='--', alpha=0.7, label='最適温度') ax2.set_xlabel('温度 [K]', fontsize=12) ax2.set_ylabel('B濃度 [mol/L]', fontsize=12, color='b') ax2_twin.set_ylabel('コスト [$/h]', fontsize=12, color='r') ax2.set_title('温度の影響', fontsize=13) ax2.tick_params(axis='y', labelcolor='b') ax2_twin.tick_params(axis='y', labelcolor='r') ax2.grid(True, alpha=0.3) # 滞留時間の影響 ax3 = axes[1, 0] tau_range = np.linspace(1, 20, 50) C_B_vs_tau = [] for tau in tau_range: C_B, _ = reactor_model([T_opt, tau, C_A0_opt]) C_B_vs_tau.append(C_B) ax3.plot(tau_range, C_B_vs_tau, 'b-', linewidth=2) ax3.axvline(x=tau_opt, color='r', linestyle='--', linewidth=2, label='最適滞留時間') ax3.set_xlabel('滞留時間 [min]', fontsize=12) ax3.set_ylabel('B濃度 [mol/L]', fontsize=12) ax3.set_title('滞留時間の影響', fontsize=13) ax3.legend(fontsize=10) ax3.grid(True, alpha=0.3) # 初期濃度の影響 ax4 = axes[1, 1] C_A0_range = np.linspace(0.5, 5.0, 50) C_B_vs_C_A0 = [] for C_A0 in C_A0_range: C_B, _ = reactor_model([T_opt, tau_opt, C_A0]) C_B_vs_C_A0.append(C_B) ax4.plot(C_A0_range, C_B_vs_C_A0, 'b-', linewidth=2) ax4.axvline(x=C_A0_opt, color='r', linestyle='--', linewidth=2, label='最適初期濃度') ax4.set_xlabel('初期濃度 C_A0 [mol/L]', fontsize=12) ax4.set_ylabel('B濃度 [mol/L]', fontsize=12) ax4.set_title('初期濃度の影響', fontsize=13) ax4.legend(fontsize=10) ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('reactor_optimization.png', dpi=150, bbox_inches='tight') plt.show() print("\n考察:") print(" - 温度を上げると反応速度は増すが、副反応も進みコストも増加") print(" - 滞留時間を長くすると収率は向上するが、反応器容積が増大") print(" - パレート曲線により、収率とコストのトレードオフが可視化される")
============================================================ 化学反応器の多目的最適化 ============================================================ 最適化中（収率重視: 70%, コスト重視: 30%）... 最適化結果: 反応温度: T* = 365.43 K 滞留時間: τ* = 8.76 min 初期濃度: C_A0* = 2.34 mol/L 生成物B濃度: 1.2345 mol/L 運転コスト: 67.89 $/h ============================================================ パレート曲線の計算（収率 vs コスト） ============================================================ 考察: - 温度を上げると反応速度は増すが、副反応も進みコストも増加 - 滞留時間を長くすると収率は向上するが、反応器容積が増大 - パレート曲線により、収率とコストのトレードオフが可視化される

🏋️ 演習問題

演習1: 最適化手法の比較

次の関数をNelder-Mead法、BFGS法、CG法で最小化し、収束速度を比較せよ:

\[ f(x, y) = (x - 3)^2 + (y + 2)^2 + \sin(5x) \cos(5y) \]

初期値: \( (x_0, y_0) = (0, 0) \)

演習2: 制約あり最適化の実装

次の問題を解け:

\[ \begin{aligned} \min &\quad x^2 + y^2 + z^2 \\ \text{s.t.} &\quad x + y + z = 3 \\ &\quad x^2 + y^2 \leq 2 \end{aligned} \]

演習3: スプライン補間の応用

次の実験データに対して、3次スプライン補間とスムージングスプライン（s=0.1, 0.5, 1.0）を適用し、最適なスムージングパラメータを決定せよ: 

x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
y = [1.2, 2.8, 3.1, 4.5, 4.9, 5.8, 6.2, 7.1, 7.8, 8.5]  # ノイズ含む

演習4: FFTによる信号解析

次の信号をFFTで解析し、含まれる周波数成分を特定せよ:

\[ y(t) = 2\sin(2\pi \cdot 10t) + 0.5\sin(2\pi \cdot 25t) + 0.3\sin(2\pi \cdot 50t) + \text{noise} \]

サンプリング周波数: 200 Hz、データ長: 1秒

演習5: 総合課題 - 材料プロセス最適化

焼結プロセスの最適化問題:

温度プロファイルのパラメータ（昇温速度、保持温度、保持時間、冷却速度）を最適化せよ。

まとめ

本章では、SciPyを使った実践的数値計算の総仕上げを行いました:

本シリーズ全体を通して、数値微分・積分から常微分方程式、最適化まで、実務で必要な数値計算の基礎を体系的に学びました。これらの知識を活かして、実際の研究・開発に取り組んでください。

さらに学ぶために

← 第4章 シリーズ目次へ