scipy-optimization

エンジニアリング設計のためのSciPy最適化

概要

scipy.optimize モジュールは、エンジニアリング設計の問題を解くための包括的な最適化アルゴリズム群を提供します。このスキルでは、ポンプ設計の最適化、システムパラメーターの調整、性能分析にこれらの手法を適用する方法に焦点を当てます。

主な機能：

• 制約なし最適化：制限なしで最適なパラメーターを見つけます
• 制約付き最適化：設計上の制約や境界条件を考慮して最適化します
• 曲線フィッティング：実験データにモデルを当てはめます
• 大域最適化：多峰性問題で大域的最小値を見つけます
• 多目的最適化：相反する設計目標のバランスを取ります

最適化手法

1. minimize() - 制約なし・制約付き最適化

ほとんどの最適化タスクで使用する中核的な関数です。

from scipy.optimize import minimize
import numpy as np

# Basic usage
def objective(x):
    """Objective function to minimize"""
    return x[0]**2 + x[1]**2

result = minimize(objective, x0=[1.0, 1.0])
print(f"Optimal point: {result.x}")
print(f"Optimal value: {result.fun}")

一般的なアルゴリズム：

• 'Nelder-Mead'：導関数不要で頑健ですが、処理は遅めです
• 'BFGS'：準ニュートン法で、滑らかな関数に対して高速です
• 'L-BFGS-B'：境界条件付きのBFGSです
• 'SLSQP'：制約付きの逐次最小二乗法です
• 'trust-constr'：一般的な制約付き最適化です（推奨）

2. least_squares() - 曲線フィッティングと残差の最小化

次の形式の問題に特化しています：残差の二乗和を最小化する

from scipy.optimize import least_squares

def residuals(params, x_data, y_data):
    """Calculate residuals between model and data"""
    a, b, c = params
    y_model = a * x_data**2 + b * x_data + c
    return y_model - y_data

# Fit quadratic to data
x_data = np.array([1, 2, 3, 4, 5])
y_data = np.array([2.1, 3.9, 9.2, 15.8, 25.1])
result = least_squares(residuals, x0=[1, 1, 1], args=(x_data, y_data))

主な機能：

• 変数の境界条件に対応します
• ロバスト損失関数を使用できます（外れ値に対応）
• 残差が大規模な問題でも効率的です
• ヤコビアン情報を返します

3. differential_evolution() - 大域最適化

進化的アルゴリズムを使用して大域的最小値を見つけます。多峰性問題で不可欠です。

from scipy.optimize import differential_evolution

def complex_function(x):
    """Function with multiple local minima"""
    return np.sin(x[0]) * np.cos(x[1]) + (x[0] - 1)**2 + (x[1] + 2)**2

# Search over bounds
bounds = [(-5, 5), (-5, 5)]
result = differential_evolution(complex_function, bounds)

使用する場面：

• 複数の局所最小値が予想される場合
• 勾配情報が得られない場合
• 頑健な大域解が必要な場合
• 計算時間が長くなっても問題ない場合

4. 制約付き最適化

等式制約と不等式制約のもとで最適化します。

from scipy.optimize import minimize, NonlinearConstraint, LinearConstraint

def objective(x):
    return x[0]**2 + x[1]**2

# Equality constraint: x[0] + x[1] = 1
def eq_constraint(x):
    return x[0] + x[1] - 1

# Inequality constraint: x[0]^2 + x[1]^2 >= 0.5
def ineq_constraint(x):
    return x[0]**2 + x[1]**2 - 0.5

constraints = [
    {'type': 'eq', 'fun': eq_constraint},
    {'type': 'ineq', 'fun': ineq_constraint}
]

result = minimize(objective, x0=[0.5, 0.5],
                 method='SLSQP', constraints=constraints)

エンジニアリングへの応用

1. ポンプ効率の最大化

設計点で効率を最大化するよう、ポンプの形状（インペラー径、羽根角度、幅）を最適化します。

問題の定式化：

• 目的関数：効率 η(D, β, b) を最大化します
• 制約条件：
  • 最小流量の要件
  • 最大消費電力
  • 幾何学的制約（例：D_min ≤ D ≤ D_max）
  • NPSHの要件

アプローチ：

def pump_efficiency(params):
    """Calculate pump efficiency (negative for minimization)"""
    D, beta, b = params  # diameter, blade angle, width

    # Simplified efficiency model
    Q = calculate_flow_rate(D, beta, b)
    H = calculate_head(D, beta, b)
    P = calculate_power(D, beta, b)

    eta = (rho * g * Q * H) / P
    return -eta  # Negative for maximization

bounds = [(0.1, 0.5), (20, 45), (0.02, 0.1)]  # D, beta, b
result = minimize(pump_efficiency, x0=[0.3, 30, 0.05],
                 method='L-BFGS-B', bounds=bounds)

2. インペラー形状の最適化

特定の運転点に合わせて、羽根形状、シュラウド輪郭、ハブ形状を最適化します。

主な考慮事項：

• 入口・出口の速度三角形
• 羽根負荷の分布
• キャビテーション性能（NPSH）
• 構造健全性（応力限界）

3. システムコストの最小化

総ライフサイクルコストを最小化する、最適なポンプサイズと運転速度を求めます。

コスト項目：

• 設備費（ポンプ＋モーター＋VFD）
• エネルギー費（運転時間 × 効率）
• メンテナンス費
• 交換費（NPV）

定式化：

def total_cost(params):
    """Total lifecycle cost"""
    pump_size, speed_rpm = params

    # Capital cost
    capital = cost_pump(pump_size) + cost_motor(pump_size, speed_rpm)

    # Operating cost (20 year life)
    annual_energy = operating_hours * power(pump_size, speed_rpm) * electricity_rate
    operating = annual_energy * 20

    # Maintenance
    maintenance = maintenance_rate * capital * 20

    return capital + operating + maintenance

4. 多目的最適化

相反する目標（例：効率とコスト、流量と揚程）のバランスを取ります。

パレートフロントのアプローチ：

def weighted_objectives(params, weight):
    """Combine objectives with weighting"""
    eta = pump_efficiency(params)
    cost = pump_cost(params)

    # Normalize and combine
    return -weight * eta + (1 - weight) * cost

# Sweep weights to find Pareto front
pareto_solutions = []
for w in np.linspace(0, 1, 21):
    result = minimize(lambda x: weighted_objectives(x, w),
                     x0=[0.3, 30, 0.05], bounds=bounds)
    pareto_solutions.append(result.x)

5. 実験データへの曲線フィッティング

ポンプ特性曲線（H-Q、η-Q、P-Q）を試験データに当てはめます。

用途：

• 試験データから性能モデルを作成します
• 設計点外の条件に対する性能を外挿します
• 解析モデルの経験係数を特定します

例：

from scipy.optimize import curve_fit

def pump_curve(Q, a, b, c):
    """Quadratic head-flow relationship"""
    return a - b * Q - c * Q**2

# Fit to test data
Q_test = np.array([0, 50, 100, 150, 200])  # m³/h
H_test = np.array([120, 115, 105, 90, 70])  # m

params, covariance = curve_fit(pump_curve, Q_test, H_test)
print(f"H = {params[0]:.2f} - {params[1]:.4f}*Q - {params[2]:.6f}*Q²")

ベストプラクティス

1. 問題の定式化

• 変数のスケーリング：同程度の大きさ（0～1または-1～1）に正規化します
• 適切な初期推定値の選択：予想される解に近い値から開始します
• 明確な境界条件の定義：物理的に不可能な解を防ぎます
• 滑らかな目的関数の定式化：不連続性を避けます

2. アルゴリズムの選択

• 滑らかで勾配が利用可能な場合：'BFGS' または 'L-BFGS-B' を使用します
• 制約がある場合：'SLSQP' または 'trust-constr' を使用します
• ノイズが多い／不連続な場合：'Nelder-Mead' または導関数不要の手法を使用します
• 複数の局所最小値がある場合：まず differential_evolution を使用し、その後、解を精緻化します
• 曲線フィッティングの場合：curve_fit または least_squares を使用します

3. 収束のヒント

• 複数の初期推定値を試してください
• 適切な許容誤差を設定してください（厳しすぎないようにします）
• 可能であれば、解析的な勾配を指定してください
• 最適化の成功フラグを確認してください
• 結果が物理的に実現可能か検証してください

4. 制約の定式化

• 不等式制約は g(x) ≥ 0 の形式で表してください
• 等式制約は h(x) = 0 の形式で表してください
• 制約は滑らかで微分可能なものにしてください
• 冗長な制約は避けてください

よくある落とし穴

1. スケーリングが不適切：変数の大きさが大きく異なる
2. 初期推定値が不適切：実行可能領域から離れた位置から開始する
3. アルゴリズムの選択ミス：滑らかでない関数に勾配ベースの手法を使用する
4. 許容誤差が厳しすぎる：達成可能な精度を超える精度を求める
5. 境界の無視：物理的な限界を制約に含めない
6. 検証を行わない：解が工学的に妥当か確認しない

パフォーマンス最適化のヒント

1. 演算をベクトル化する：ループの代わりに NumPy 配列を使用する
2. ヤコビアンを指定する：解析的な勾配はより高速で精度も高くなります
3. 適切な手法を使用する：局所最適化で十分な場合は、大域最適化を使用しないでください
4. コストの高い計算をキャッシュする：中間結果を保存する
5. シンプルなモデルから始める：段階的に複雑さを高める

ポンプ設計ワークフローとの統合

1. Define design variables (D, β, b, ω, etc.)
2. Establish objective function (maximize η, minimize cost)
3. Define constraints (Q_min, NPSH, stress limits)
4. Set bounds (physical limits)
5. Choose optimization algorithm
6. Run optimization with multiple initial guesses
7. Validate results (CFD, analytical checks)
8. Perform sensitivity analysis
9. Generate documentation and plots

ワークフローの例

import numpy as np
from scipy.optimize import minimize, differential_evolution
import matplotlib.pyplot as plt

# Step 1: Define problem
def objective(x):
    D, beta = x
    return -calculate_efficiency(D, beta)

def constraint_flow(x):
    D, beta = x
    return calculate_flow_rate(D, beta) - Q_min

def constraint_power(x):
    D, beta = x
    return P_max - calculate_power(D, beta)

# Step 2: Set up optimization
bounds = [(0.2, 0.6), (20, 50)]
constraints = [
    {'type': 'ineq', 'fun': constraint_flow},
    {'type': 'ineq', 'fun': constraint_power}
]

# Step 3: Solve with global optimizer first
result_global = differential_evolution(
    objective, bounds, constraints=constraints
)

# Step 4: Refine with local optimizer
result_local = minimize(
    objective, x0=result_global.x,
    method='SLSQP', bounds=bounds,
    constraints=constraints
)

# Step 5: Report results
print(f"Optimal diameter: {result_local.x[0]:.3f} m")
print(f"Optimal blade angle: {result_local.x[1]:.1f}°")
print(f"Maximum efficiency: {-result_local.fun:.2%}")

参考資料

• SciPy 最適化ドキュメント
• 最適化手法のチュートリアル
• 問題の定式化に関する工学最適化の教科書
• 現実的な制約やモデルに関するポンプ設計ハンドブック