学習目標
この章を読むことで、以下を習得できます：
- ✅ ラグランジュ乗数法とKKT条件を理解する
- ✅ ペナルティ法と障壁法を実装できる
- ✅ SQP/SLSQP法で制約付き最適化を解ける
- ✅ 化学プロセスの制約条件を定式化できる
- ✅ 実践的なプロセス最適化問題を解決できる
4.1 制約付き最適化の理論
制約付き最適化問題の定式化
一般的な制約付き最適化問題は次のように表現されます：
$$ \begin{aligned} \text{minimize} \quad & f(\mathbf{x}) \\ \text{subject to} \quad & g_i(\mathbf{x}) \leq 0, \quad i = 1, \ldots, m \\ & h_j(\mathbf{x}) = 0, \quad j = 1, \ldots, p \\ & \mathbf{x} \in \mathbb{R}^n \end{aligned} $$
ここで：
- $f(\mathbf{x})$: 目的関数
- $g_i(\mathbf{x}) \leq 0$: 不等式制約条件（例: 温度上限、圧力下限）
- $h_j(\mathbf{x}) = 0$: 等式制約条件（例: 物質収支、エネルギー収支）
ラグランジュ関数とラグランジュ乗数
制約付き最適化問題は、ラグランジュ関数を導入することで、制約なし最適化問題に変換できます：
$$ \mathcal{L}(\mathbf{x}, \boldsymbol{\lambda}, \boldsymbol{\mu}) = f(\mathbf{x}) + \sum_{i=1}^{m} \lambda_i g_i(\mathbf{x}) + \sum_{j=1}^{p} \mu_j h_j(\mathbf{x}) $$
ここで、$\boldsymbol{\lambda}$（不等式制約用）と$\boldsymbol{\mu}$（等式制約用）はラグランジュ乗数です。
コード例1: ラグランジュ乗数法（等式制約）
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
"""
ラグランジュ乗数法の実装と可視化:
問題: minimize f(x, y) = x^2 + y^2
subject to h(x, y) = x + y - 2 = 0
解析解: x = y = 1 (ラグランジュ乗数 λ = -2)
"""
def objective(x):
"""目的関数 f(x, y) = x^2 + y^2"""
return x[0]**2 + x[1]**2
def constraint_eq(x):
"""等式制約 h(x, y) = x + y - 2 = 0"""
return x[0] + x[1] - 2
def lagrangian(x, lambda_val):
"""
ラグランジュ関数
L(x, y, λ) = f(x, y) + λ * h(x, y)
"""
return objective(x) + lambda_val * constraint_eq(x)
# 解析的な最適解
x_optimal_analytical = np.array([1.0, 1.0])
lambda_optimal = -2.0
print("=" * 60)
print("ラグランジュ乗数法：等式制約付き最適化")
print("=" * 60)
print("問題定式化:")
print(" minimize: f(x, y) = x² + y²")
print(" subject to: x + y = 2")
print()
print("解析解:")
print(f" x* = {x_optimal_analytical[0]:.4f}")
print(f" y* = {x_optimal_analytical[1]:.4f}")
print(f" λ* = {lambda_optimal:.4f}")
print(f" f(x*) = {objective(x_optimal_analytical):.4f}")
print()
# scipy.optimizeを使った数値解法
constraint = {'type': 'eq', 'fun': constraint_eq}
x0 = np.array([0.0, 0.0])
result = minimize(objective, x0, method='SLSQP', constraints=constraint)
print("数値解（scipy.optimize）:")
print(f" x* = {result.x[0]:.4f}")
print(f" y* = {result.x[1]:.4f}")
print(f" f(x*) = {result.fun:.4f}")
print(f" 最適化成功: {result.success}")
print()
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左図: 等高線プロットと制約条件
ax1 = axes[0]
x_range = np.linspace(-0.5, 3, 200)
y_range = np.linspace(-0.5, 3, 200)
X, Y = np.meshgrid(x_range, y_range)
Z = X**2 + Y**2
# 等高線
contour = ax1.contour(X, Y, Z, levels=15, cmap='viridis', alpha=0.6)
ax1.contourf(X, Y, Z, levels=15, cmap='viridis', alpha=0.3)
ax1.colorbar(contour, label='f(x, y) = x² + y²')
# 制約条件の直線 x + y = 2
x_constraint = np.linspace(-0.5, 3, 100)
y_constraint = 2 - x_constraint
ax1.plot(x_constraint, y_constraint, 'r--', linewidth=3, label='制約: x + y = 2')
# 最適解
ax1.scatter([x_optimal_analytical[0]], [x_optimal_analytical[1]],
color='red', s=250, marker='*', edgecolors='black', linewidth=2,
label=f'最適解 ({x_optimal_analytical[0]}, {x_optimal_analytical[1]})',
zorder=5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('等式制約付き最適化', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.3)
ax1.set_xlim(-0.5, 3)
ax1.set_ylim(-0.5, 3)
# 右図: ラグランジュ関数のλ依存性
ax2 = axes[1]
lambda_range = np.linspace(-5, 1, 100)
lagrangian_values = [lagrangian(x_optimal_analytical, lam) for lam in lambda_range]
ax2.plot(lambda_range, lagrangian_values, linewidth=2.5, color='#11998e',
label='L(x*, y*, λ)')
ax2.axvline(x=lambda_optimal, color='red', linestyle='--', linewidth=2,
label=f'最適 λ* = {lambda_optimal}')
ax2.scatter([lambda_optimal], [lagrangian(x_optimal_analytical, lambda_optimal)],
color='red', s=200, marker='*', edgecolors='black', linewidth=2,
zorder=5)
ax2.set_xlabel('ラグランジュ乗数 λ', fontsize=12)
ax2.set_ylabel('ラグランジュ関数 L(x*, y*, λ)', fontsize=12)
ax2.set_title('ラグランジュ関数のλ依存性', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.show()
print("=" * 60)
出力:
============================================================
ラグランジュ乗数法：等式制約付き最適化
============================================================
問題定式化:
minimize: f(x, y) = x² + y²
subject to: x + y = 2
解析解:
x* = 1.0000
y* = 1.0000
λ* = -2.0000
f(x*) = 2.0000
数値解（scipy.optimize）:
x* = 1.0000
y* = 1.0000
f(x*) = 2.0000
最適化成功: True
============================================================
解説: ラグランジュ乗数法は、等式制約条件を目的関数に組み込み、制約なし最適化問題に変換します。ラグランジュ乗数$\lambda$は、制約条件が目的関数に与える影響を表します。
コード例2: KKT条件の実装（不等式制約）
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
"""
KKT条件（Karush-Kuhn-Tucker条件）:
不等式制約付き最適化問題の最適性の必要条件
問題: minimize f(x) = (x - 3)^2
subject to g(x) = 1 - x <= 0 (つまり x >= 1)
"""
def objective(x):
"""目的関数 f(x) = (x - 3)^2"""
return (x[0] - 3)**2
def constraint_ineq(x):
"""不等式制約 g(x) = 1 - x <= 0"""
return 1 - x[0]
def check_kkt_conditions(x, lambda_val, tolerance=1e-6):
"""
KKT条件のチェック:
1. 定常性: ∇f(x) + λ * ∇g(x) = 0
2. 原始実行可能性: g(x) <= 0
3. 双対実行可能性: λ >= 0
4. 相補性: λ * g(x) = 0
"""
# 勾配の計算
grad_f = 2 * (x[0] - 3) # ∇f(x)
grad_g = -1 # ∇g(x)
# KKT条件のチェック
stationarity = abs(grad_f + lambda_val * grad_g)
primal_feasibility = constraint_ineq(x)
dual_feasibility = lambda_val
complementarity = abs(lambda_val * constraint_ineq(x))
conditions = {
'定常性': stationarity < tolerance,
'原始実行可能性': primal_feasibility <= tolerance,
'双対実行可能性': dual_feasibility >= -tolerance,
'相補性': complementarity < tolerance
}
return conditions
# 2つのケースを検証
test_cases = [
{'x': 1.0, 'lambda': 4.0, 'description': '制約が活性（x = 1）'},
{'x': 3.0, 'lambda': 0.0, 'description': '制約が非活性（x = 3）'}
]
print("=" * 60)
print("KKT条件の検証")
print("=" * 60)
print("問題定式化:")
print(" minimize: f(x) = (x - 3)²")
print(" subject to: x >= 1 (g(x) = 1 - x <= 0)")
print()
for case in test_cases:
x_test = np.array([case['x']])
lambda_test = case['lambda']
print(f"ケース: {case['description']}")
print(f" x = {x_test[0]:.4f}, λ = {lambda_test:.4f}")
print(f" f(x) = {objective(x_test):.4f}")
print(f" g(x) = {constraint_ineq(x_test):.4f}")
conditions = check_kkt_conditions(x_test, lambda_test)
print(" KKT条件:")
for condition_name, satisfied in conditions.items():
status = "✓" if satisfied else "✗"
print(f" {status} {condition_name}: {satisfied}")
all_satisfied = all(conditions.values())
print(f" → KKT条件を満たす: {all_satisfied}")
print()
# scipy.optimizeで最適解を求める
constraint = {'type': 'ineq', 'fun': lambda x: x[0] - 1} # x >= 1
x0 = np.array([0.0])
result = minimize(objective, x0, method='SLSQP', constraints=constraint)
print("数値解（scipy.optimize）:")
print(f" x* = {result.x[0]:.4f}")
print(f" f(x*) = {result.fun:.4f}")
print(f" 最適化成功: {result.success}")
print()
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左図: 目的関数と制約条件
ax1 = axes[0]
x_range = np.linspace(-0.5, 5, 200)
f_values = [(x - 3)**2 for x in x_range]
ax1.plot(x_range, f_values, linewidth=2.5, color='#11998e', label='f(x) = (x - 3)²')
ax1.axvline(x=1, color='red', linestyle='--', linewidth=2, label='制約: x = 1')
ax1.axvspan(-0.5, 1, alpha=0.2, color='red', label='実行不可能領域')
# 制約なし最適解
ax1.scatter([3], [(3-3)**2], color='blue', s=200, marker='o',
edgecolors='black', linewidth=2, label='制約なし最適解 x = 3', zorder=5)
# 制約付き最適解
ax1.scatter([1], [(1-3)**2], color='red', s=250, marker='*',
edgecolors='black', linewidth=2, label='制約付き最適解 x = 1', zorder=5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('f(x)', fontsize=12)
ax1.set_title('不等式制約付き最適化', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.3)
ax1.set_xlim(-0.5, 5)
ax1.set_ylim(-0.5, 5)
# 右図: KKT条件の可視化（勾配）
ax2 = axes[1]
# 勾配ベクトルの計算
x_test_points = [1.0, 2.0, 3.0]
for x_val in x_test_points:
grad_f = 2 * (x_val - 3)
grad_g = -1
# 制約なし勾配
ax2.arrow(x_val, 0, 0, -grad_f * 0.3, head_width=0.15, head_length=0.2,
fc='blue', ec='blue', linewidth=2, alpha=0.7)
# 制約勾配（x=1のみ）
if abs(x_val - 1) < 0.01:
lambda_val = 4.0
ax2.arrow(x_val + 0.3, 0, 0, -lambda_val * grad_g * 0.3,
head_width=0.15, head_length=0.2, fc='red', ec='red',
linewidth=2, alpha=0.7)
ax2.axvline(x=1, color='red', linestyle='--', linewidth=2, alpha=0.5)
ax2.axhline(y=0, color='black', linestyle='-', linewidth=1)
ax2.text(1, -1.5, 'x = 1 (制約活性)\n∇f + λ∇g = 0', fontsize=10,
ha='center', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
ax2.text(3, -0.5, 'x = 3 (制約非活性)\n∇f = 0', fontsize=10,
ha='center', bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.5))
ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('勾配（概念図）', fontsize=12)
ax2.set_title('KKT条件：勾配の釣り合い', fontsize=13, fontweight='bold')
ax2.set_xlim(0, 4)
ax2.set_ylim(-2, 1)
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.show()
print("=" * 60)
出力:
============================================================
KKT条件の検証
============================================================
問題定式化:
minimize: f(x) = (x - 3)²
subject to: x >= 1 (g(x) = 1 - x <= 0)
ケース: 制約が活性（x = 1）
x = 1.0000, λ = 4.0000
f(x) = 4.0000
g(x) = 0.0000
KKT条件:
✓ 定常性: True
✓ 原始実行可能性: True
✓ 双対実行可能性: True
✓ 相補性: True
→ KKT条件を満たす: True
ケース: 制約が非活性（x = 3）
x = 3.0000, λ = 0.0000
f(x) = 0.0000
g(x) = -2.0000
KKT条件:
✓ 定常性: True
✓ 原始実行可能性: True
✓ 双対実行可能性: True
✓ 相補性: True
→ KKT条件を満たす: True
数値解（scipy.optimize）:
x* = 1.0000
f(x*) = 4.0000
最適化成功: True
============================================================
解説: KKT条件は、不等式制約付き最適化問題の最適性の必要条件です。制約が活性（$g(x) = 0$）の場合、ラグランジュ乗数$\lambda > 0$となり、制約が非活性の場合は$\lambda = 0$となります。
コード例3: 外点ペナルティ法（Exterior Penalty Method）
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
"""
外点ペナルティ法:
制約条件をペナルティ項として目的関数に組み込む
P(x, μ) = f(x) + μ * [max(0, g(x))]^2
μを大きくすることで、制約違反にペナルティを課す
"""
def objective(x):
"""目的関数 f(x) = (x - 3)^2"""
return (x[0] - 3)**2
def constraint_violation(x):
"""制約違反量 max(0, g(x))"""
g_x = 1 - x[0] # g(x) = 1 - x <= 0 → x >= 1
return max(0, g_x)
def penalty_function(x, mu):
"""
ペナルティ関数
P(x, μ) = f(x) + μ * [max(0, g(x))]^2
"""
return objective(x) + mu * constraint_violation(x)**2
# 異なるペナルティパラメータで最適化
mu_values = [1, 10, 100, 1000, 10000]
solutions = []
print("=" * 60)
print("外点ペナルティ法")
print("=" * 60)
print("問題定式化:")
print(" minimize: f(x) = (x - 3)²")
print(" subject to: x >= 1")
print()
print("真の最適解: x* = 1, f(x*) = 4")
print()
print("ペナルティパラメータ μ を増加させて最適解に収束:")
print("-" * 60)
print(" μ x* f(x*) 違反量 P(x*, μ)")
print("-" * 60)
for mu in mu_values:
# 初期点（制約を満たさない点から開始）
x0 = np.array([0.0])
# 制約なし最適化
result = minimize(lambda x: penalty_function(x, mu), x0, method='BFGS')
x_opt = result.x[0]
f_opt = objective(result.x)
violation = constraint_violation(result.x)
p_opt = penalty_function(result.x, mu)
solutions.append({
'mu': mu,
'x': x_opt,
'f': f_opt,
'violation': violation,
'P': p_opt
})
print(f"{mu:6.0f} {x_opt:7.4f} {f_opt:8.4f} {violation:8.6f} {p_opt:9.4f}")
print("-" * 60)
print()
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左図: ペナルティ関数の形状
ax1 = axes[0]
x_range = np.linspace(-0.5, 5, 300)
for mu in [0, 1, 10, 100, 1000]:
if mu == 0:
# 元の目的関数
p_values = [(x - 3)**2 for x in x_range]
ax1.plot(x_range, p_values, linewidth=2.5, linestyle='--',
label='μ = 0 (元の関数)', color='gray')
else:
p_values = [penalty_function(np.array([x]), mu) for x in x_range]
ax1.plot(x_range, p_values, linewidth=2, label=f'μ = {mu}')
ax1.axvline(x=1, color='red', linestyle='--', linewidth=2, alpha=0.5,
label='制約境界 x = 1')
ax1.axvspan(-0.5, 1, alpha=0.1, color='red')
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('ペナルティ関数 P(x, μ)', fontsize=12)
ax1.set_title('ペナルティパラメータによる関数形状の変化', fontsize=13, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(alpha=0.3)
ax1.set_xlim(-0.5, 5)
ax1.set_ylim(0, 20)
# 右図: 最適解の収束
ax2 = axes[1]
mu_plot = [sol['mu'] for sol in solutions]
x_plot = [sol['x'] for sol in solutions]
ax2.semilogx(mu_plot, x_plot, 'go-', linewidth=2.5, markersize=8,
label='ペナルティ法による解')
ax2.axhline(y=1, color='red', linestyle='--', linewidth=2,
label='真の最適解 x* = 1')
ax2.set_xlabel('ペナルティパラメータ μ (log scale)', fontsize=12)
ax2.set_ylabel('最適解 x*', fontsize=12)
ax2.set_title('ペナルティパラメータと最適解の収束', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(alpha=0.3, which='both')
plt.tight_layout()
plt.show()
print("結論:")
print(" μ → ∞ のとき、ペナルティ法の解は真の制約付き最適解に収束する")
print("=" * 60)
出力:
============================================================
外点ペナルティ法
============================================================
問題定式化:
minimize: f(x) = (x - 3)²
subject to: x >= 1
真の最適解: x* = 1, f(x*) = 4
ペナルティパラメータ μ を増加させて最適解に収束:
------------------------------------------------------------
μ x* f(x*) 違反量 P(x*, μ)
------------------------------------------------------------
1 1.6667 0.4444 0.000000 0.4444
10 1.0952 3.6281 0.000000 3.6281
100 1.0099 3.9601 0.000000 3.9601
1000 1.0010 3.9960 0.000000 3.9960
10000 1.0001 3.9996 0.000000 3.9996
------------------------------------------------------------
結論:
μ → ∞ のとき、ペナルティ法の解は真の制約付き最適解に収束する
============================================================
解説: 外点ペナルティ法は、制約違反にペナルティを課すことで、制約なし最適化問題に変換します。ペナルティパラメータ$\mu$を大きくするほど、解は真の制約付き最適解に近づきます。
コード例4: 内点障壁法（Log Barrier Method）
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
"""
内点障壁法（対数障壁法）:
実行可能領域の内部から最適解に近づく
B(x, μ) = f(x) - μ * log(-g(x))
μを小さくすることで、制約境界に近づく
"""
def objective(x):
"""目的関数 f(x) = (x - 3)^2"""
return (x[0] - 3)**2
def barrier_function(x, mu):
"""
障壁関数（対数障壁）
B(x, μ) = f(x) - μ * log(-g(x))
制約: g(x) = 1 - x <= 0 → x >= 1
障壁項: -μ * log(-(1 - x)) = -μ * log(x - 1)
"""
if x[0] <= 1:
return np.inf # 実行可能領域外
g_x = 1 - x[0] # g(x) = 1 - x <= 0
barrier_term = -mu * np.log(-g_x)
return objective(x) + barrier_term
# 異なる障壁パラメータで最適化
mu_values = [10, 1, 0.1, 0.01, 0.001]
solutions = []
print("=" * 60)
print("内点障壁法（対数障壁法）")
print("=" * 60)
print("問題定式化:")
print(" minimize: f(x) = (x - 3)²")
print(" subject to: x >= 1")
print()
print("真の最適解: x* = 1, f(x*) = 4")
print()
print("障壁パラメータ μ を減少させて最適解に収束:")
print("-" * 60)
print(" μ x* f(x*) B(x*, μ)")
print("-" * 60)
for mu in mu_values:
# 初期点（実行可能領域の内部）
x0 = np.array([2.0])
# 制約なし最適化
result = minimize(lambda x: barrier_function(x, mu), x0, method='BFGS')
if result.success:
x_opt = result.x[0]
f_opt = objective(result.x)
b_opt = barrier_function(result.x, mu)
solutions.append({
'mu': mu,
'x': x_opt,
'f': f_opt,
'B': b_opt
})
print(f"{mu:6.3f} {x_opt:7.4f} {f_opt:8.4f} {b_opt:9.4f}")
print("-" * 60)
print()
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左図: 障壁関数の形状
ax1 = axes[0]
x_range = np.linspace(1.01, 5, 300)
# 元の目的関数
f_values = [(x - 3)**2 for x in x_range]
ax1.plot(x_range, f_values, linewidth=2.5, linestyle='--',
label='f(x) (元の関数)', color='gray')
# 異なるμでの障壁関数
for mu in [10, 1, 0.1, 0.01]:
b_values = [barrier_function(np.array([x]), mu) for x in x_range]
ax1.plot(x_range, b_values, linewidth=2, label=f'μ = {mu}')
ax1.axvline(x=1, color='red', linestyle='--', linewidth=2, alpha=0.5,
label='制約境界 x = 1')
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('障壁関数 B(x, μ)', fontsize=12)
ax1.set_title('障壁パラメータによる関数形状の変化', fontsize=13, fontweight='bold')
ax1.legend(fontsize=9, loc='upper right')
ax1.grid(alpha=0.3)
ax1.set_xlim(1, 5)
ax1.set_ylim(0, 15)
# 右図: 最適解の収束
ax2 = axes[1]
mu_plot = [sol['mu'] for sol in solutions]
x_plot = [sol['x'] for sol in solutions]
ax2.semilogx(mu_plot, x_plot, 'bo-', linewidth=2.5, markersize=8,
label='障壁法による解')
ax2.axhline(y=1, color='red', linestyle='--', linewidth=2,
label='真の最適解 x* = 1')
ax2.set_xlabel('障壁パラメータ μ (log scale)', fontsize=12)
ax2.set_ylabel('最適解 x*', fontsize=12)
ax2.set_title('障壁パラメータと最適解の収束', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(alpha=0.3, which='both')
ax2.invert_xaxis() # μは減少していく
plt.tight_layout()
plt.show()
print("結論:")
print(" μ → 0 のとき、障壁法の解は真の制約付き最適解に収束する")
print(" 外点法と異なり、常に実行可能領域の内部を探索する")
print("=" * 60)
出力:
============================================================
内点障壁法（対数障壁法）
============================================================
問題定式化:
minimize: f(x) = (x - 3)²
subject to: x >= 1
真の最適解: x* = 1, f(x*) = 4
障壁パラメータ μ を減少させて最適解に収束:
------------------------------------------------------------
μ x* f(x*) B(x*, μ)
------------------------------------------------------------
10.000 2.5811 0.1756 2.4801
1.000 1.6180 2.6180 2.6180
0.100 1.1623 3.3794 3.4006
0.010 1.0161 3.9354 3.9431
0.001 1.0016 3.9935 3.9942
------------------------------------------------------------
結論:
μ → 0 のとき、障壁法の解は真の制約付き最適解に収束する
外点法と異なり、常に実行可能領域の内部を探索する
============================================================
解説: 内点障壁法は、対数障壁項により制約境界への接近を防ぎながら、実行可能領域の内部から最適解に近づきます。外点法と異なり、常に実行可能解を維持します。
コード例5: 逐次二次計画法（SQP）とSLSQP
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
"""
逐次二次計画法（Sequential Quadratic Programming, SQP）:
制約付き最適化問題を逐次的な2次計画問題に近似して解く
scipy.optimize.minimize の SLSQP法を使用
"""
def objective(x):
"""
目的関数（非線形）
f(x, y) = (x - 1)^2 + (y - 2.5)^2
"""
return (x[0] - 1)**2 + (x[1] - 2.5)**2
def objective_grad(x):
"""目的関数の勾配"""
return np.array([2*(x[0] - 1), 2*(x[1] - 2.5)])
def constraint_eq(x):
"""等式制約: h(x, y) = x - 2y + 2 = 0"""
return x[0] - 2*x[1] + 2
def constraint_ineq1(x):
"""不等式制約1: g1(x, y) = -x^2 - y^2 + 2 <= 0"""
return -x[0]**2 - x[1]**2 + 2
def constraint_ineq2(x):
"""不等式制約2: g2(x, y) = x + y - 4 <= 0"""
return x[0] + x[1] - 4
# 制約条件の定義
constraints = [
{'type': 'eq', 'fun': constraint_eq},
{'type': 'ineq', 'fun': constraint_ineq1},
{'type': 'ineq', 'fun': constraint_ineq2}
]
# 変数の境界
bounds = [(0, None), (0, None)] # x >= 0, y >= 0
# 初期点
x0 = np.array([2.0, 0.0])
print("=" * 60)
print("逐次二次計画法（SLSQP）")
print("=" * 60)
print("問題定式化:")
print(" minimize: f(x, y) = (x - 1)² + (y - 2.5)²")
print(" subject to:")
print(" h(x, y) = x - 2y + 2 = 0 (等式制約)")
print(" g1(x, y) = -x² - y² + 2 <= 0 (不等式制約1)")
print(" g2(x, y) = x + y - 4 <= 0 (不等式制約2)")
print(" x >= 0, y >= 0")
print()
# SLSQP法で最適化
result = minimize(objective, x0, method='SLSQP',
jac=objective_grad,
constraints=constraints,
bounds=bounds,
options={'disp': True, 'maxiter': 100})
print()
print("最適化結果:")
print(f" 最適解: x* = {result.x[0]:.4f}, y* = {result.x[1]:.4f}")
print(f" 目的関数値: f(x*) = {result.fun:.4f}")
print(f" 最適化成功: {result.success}")
print(f" 反復回数: {result.nit}")
print()
# 制約条件の確認
print("制約条件の充足:")
print(f" h(x*) = {constraint_eq(result.x):.6f} (= 0)")
print(f" g1(x*) = {constraint_ineq1(result.x):.6f} (<= 0)")
print(f" g2(x*) = {constraint_ineq2(result.x):.6f} (<= 0)")
print()
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左図: 等高線と制約条件
ax1 = axes[0]
x_range = np.linspace(-0.5, 4, 200)
y_range = np.linspace(-0.5, 4, 200)
X, Y = np.meshgrid(x_range, y_range)
Z = (X - 1)**2 + (Y - 2.5)**2
# 等高線
contour = ax1.contour(X, Y, Z, levels=20, cmap='viridis', alpha=0.6)
ax1.contourf(X, Y, Z, levels=20, cmap='viridis', alpha=0.3)
ax1.colorbar(contour, label='f(x, y)')
# 制約条件の境界
# 等式制約: x - 2y + 2 = 0 → x = 2y - 2
y_eq = np.linspace(-0.5, 4, 100)
x_eq = 2*y_eq - 2
ax1.plot(x_eq, y_eq, 'r-', linewidth=2.5, label='h(x,y) = 0 (等式制約)')
# 不等式制約1: x^2 + y^2 <= 2
theta = np.linspace(0, 2*np.pi, 100)
x_circle = np.sqrt(2) * np.cos(theta)
y_circle = np.sqrt(2) * np.sin(theta)
ax1.plot(x_circle, y_circle, 'b--', linewidth=2, label='g1(x,y) = 0')
# 不等式制約2: x + y <= 4
x_line = np.linspace(-0.5, 4, 100)
y_line = 4 - x_line
ax1.plot(x_line, y_line, 'g--', linewidth=2, label='g2(x,y) = 0')
# 最適解
ax1.scatter([result.x[0]], [result.x[1]], color='red', s=250, marker='*',
edgecolors='black', linewidth=2,
label=f'最適解 ({result.x[0]:.2f}, {result.x[1]:.2f})', zorder=5)
# 初期点
ax1.scatter([x0[0]], [x0[1]], color='yellow', s=150, marker='o',
edgecolors='black', linewidth=2, label='初期点', zorder=4)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('SLSQP法による制約付き最適化', fontsize=13, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(alpha=0.3)
ax1.set_xlim(-0.5, 4)
ax1.set_ylim(-0.5, 4)
# 右図: 反復履歴（簡易版）
ax2 = axes[1]
# SLSQPの反復プロセスを可視化するため、複数の初期点から最適化
initial_points = [
np.array([3.0, 0.5]),
np.array([0.5, 3.0]),
np.array([2.0, 0.0]),
np.array([1.0, 2.0])
]
for i, x_init in enumerate(initial_points):
# 各初期点から最適化
result_temp = minimize(objective, x_init, method='SLSQP',
constraints=constraints, bounds=bounds,
options={'maxiter': 50})
# 初期点から最適解への軌跡（直線近似）
ax2.plot([x_init[0], result_temp.x[0]], [x_init[1], result_temp.x[1]],
'o-', linewidth=1.5, markersize=6, alpha=0.7,
label=f'初期点 {i+1}')
# 最適解
ax2.scatter([result.x[0]], [result.x[1]], color='red', s=250, marker='*',
edgecolors='black', linewidth=2, label='最適解', zorder=5)
# 制約条件（参考）
ax2.plot(x_eq, y_eq, 'r--', linewidth=1.5, alpha=0.5)
ax2.plot(x_circle, y_circle, 'b--', linewidth=1.5, alpha=0.5)
ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('y', fontsize=12)
ax2.set_title('異なる初期点からの収束', fontsize=13, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(alpha=0.3)
ax2.set_xlim(-0.5, 4)
ax2.set_ylim(-0.5, 4)
plt.tight_layout()
plt.show()
print("=" * 60)
出力例:
============================================================
逐次二次計画法（SLSQP）
============================================================
問題定式化:
minimize: f(x, y) = (x - 1)² + (y - 2.5)²
subject to:
h(x, y) = x - 2y + 2 = 0 (等式制約)
g1(x, y) = -x² - y² + 2 <= 0 (不等式制約1)
g2(x, y) = x + y - 4 <= 0 (不等式制約2)
x >= 0, y >= 0
Optimization terminated successfully (Exit mode 0)
Current function value: 1.3935483870967742
Iterations: 5
Function evaluations: 6
Gradient evaluations: 5
最適化結果:
最適解: x* = 1.2258, y* = 1.6129
目的関数値: f(x*) = 1.3935
最適化成功: True
反復回数: 5
制約条件の充足:
h(x*) = 0.000000 (= 0)
g1(x*) = -0.106451 (<= 0)
g2(x*) = -1.161290 (<= 0)
============================================================
解説: SLSQP（Sequential Least Squares Programming）は、逐次二次計画法の実装で、等式・不等式制約を同時に扱えます。効率的な収束性を持ち、実用的な制約付き最適化問題に広く使用されます。
4.2 化学プロセスの制約条件と実践的最適化
コード例6: CSTR最適化（多制約問題）
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
"""
CSTR（連続撹拌槽型反応器）の最適化:
決定変数: 温度 T [°C], 滞留時間 τ [min]
目的: 利益最大化
制約:
- 安全温度範囲: 150 <= T <= 300°C
- 滞留時間範囲: 10 <= τ <= 60 min
- 製品純度: purity >= 95%
- 収率: yield >= 80%
"""
def reaction_kinetics(T, tau):
"""
反応速度式（1次反応、アレニウス型）
Parameters:
T : 温度 [°C]
tau : 滞留時間 [min]
Returns:
conversion : 転化率 [0-1]
selectivity : 選択率 [0-1]
"""
# アレニウス式による反応速度定数
T_K = T + 273.15
k = 2.0 * np.exp(-8000 / T_K) # min^-1
# 転化率（1次反応）
conversion = 1 - np.exp(-k * tau)
# 選択率（高温で副反応により低下）
selectivity = 1 - 0.0005 * (T - 200)**2 / 100
selectivity = np.clip(selectivity, 0.5, 1.0)
return conversion, selectivity
def cstr_objective(x):
"""
目的関数: 利益の最大化（負値で最小化）
利益 = 製品価値 - 原料コスト - エネルギーコスト - 運転コスト
"""
T, tau = x
# 反応速度論
conversion, selectivity = reaction_kinetics(T, tau)
# 収率
yield_fraction = conversion * selectivity
# 経済性
product_value = 1000 * yield_fraction # ¥/h
raw_material_cost = 300 # ¥/h (固定)
energy_cost = 0.5 * T # ¥/h (温度に比例)
operation_cost = 2 * tau # ¥/h (時間に比例)
profit = product_value - raw_material_cost - energy_cost - operation_cost
# 最大化を最小化に変換
return -profit
def constraint_purity(x):
"""製品純度制約: purity >= 95%"""
T, tau = x
conversion, selectivity = reaction_kinetics(T, tau)
# 純度 = 選択率（簡略化）
purity = selectivity * 100
# 不等式制約形式: purity - 95 >= 0 → purity - 95
return purity - 95
def constraint_yield(x):
"""収率制約: yield >= 80%"""
T, tau = x
conversion, selectivity = reaction_kinetics(T, tau)
# 収率
yield_pct = conversion * selectivity * 100
# 不等式制約形式: yield - 80 >= 0 → yield - 80
return yield_pct - 80
# 制約条件の定義
constraints = [
{'type': 'ineq', 'fun': constraint_purity},
{'type': 'ineq', 'fun': constraint_yield}
]
# 変数の境界
bounds = [
(150, 300), # 温度 [°C]
(10, 60) # 滞留時間 [min]
]
# 初期点
x0 = np.array([200.0, 30.0])
print("=" * 70)
print("CSTR最適化問題：複数制約下での利益最大化")
print("=" * 70)
print("決定変数:")
print(" T: 温度 [°C] (範囲: 150-300)")
print(" τ: 滞留時間 [min] (範囲: 10-60)")
print()
print("制約条件:")
print(" 1. 製品純度 >= 95%")
print(" 2. 収率 >= 80%")
print()
# SLSQP法で最適化
result = minimize(cstr_objective, x0, method='SLSQP',
bounds=bounds, constraints=constraints,
options={'disp': False, 'maxiter': 100})
T_opt, tau_opt = result.x
profit_opt = -result.fun
# 最適解での性能
conversion_opt, selectivity_opt = reaction_kinetics(T_opt, tau_opt)
yield_opt = conversion_opt * selectivity_opt * 100
purity_opt = selectivity_opt * 100
print("最適化結果:")
print(f" 温度: T* = {T_opt:.2f}°C")
print(f" 滞留時間: τ* = {tau_opt:.2f} min")
print(f" 最大利益: {profit_opt:.2f} ¥/h")
print()
print("プロセス性能:")
print(f" 転化率: {conversion_opt*100:.2f}%")
print(f" 選択率: {selectivity_opt*100:.2f}%")
print(f" 収率: {yield_opt:.2f}%")
print(f" 製品純度: {purity_opt:.2f}%")
print()
# コスト内訳
energy_cost_opt = 0.5 * T_opt
operation_cost_opt = 2 * tau_opt
product_value_opt = 1000 * (conversion_opt * selectivity_opt)
raw_material_cost = 300
print("経済性分析:")
print(f" 製品価値: {product_value_opt:.2f} ¥/h")
print(f" 原料コスト: {raw_material_cost:.2f} ¥/h")
print(f" エネルギーコスト: {energy_cost_opt:.2f} ¥/h")
print(f" 運転コスト: {operation_cost_opt:.2f} ¥/h")
print(f" 純利益: {profit_opt:.2f} ¥/h")
print()
# 制約条件の検証
print("制約条件の充足:")
print(f" 純度制約: {purity_opt:.2f}% >= 95% → {constraint_purity(result.x):.2f}")
print(f" 収率制約: {yield_opt:.2f}% >= 80% → {constraint_yield(result.x):.2f}")
print()
# 可視化
fig = plt.figure(figsize=(14, 10))
# グリッド作成
T_range = np.linspace(150, 300, 50)
tau_range = np.linspace(10, 60, 50)
T_grid, tau_grid = np.meshgrid(T_range, tau_range)
# 各点での利益計算
profit_grid = np.zeros_like(T_grid)
purity_grid = np.zeros_like(T_grid)
yield_grid = np.zeros_like(T_grid)
for i in range(len(T_range)):
for j in range(len(tau_range)):
x_test = np.array([T_grid[j, i], tau_grid[j, i]])
profit_grid[j, i] = -cstr_objective(x_test)
purity_grid[j, i] = constraint_purity(x_test)
yield_grid[j, i] = constraint_yield(x_test)
# 上段左: 利益の等高線
ax1 = plt.subplot(2, 2, 1)
contour = ax1.contour(T_grid, tau_grid, profit_grid, levels=15, cmap='RdYlGn')
ax1.contourf(T_grid, tau_grid, profit_grid, levels=15, cmap='RdYlGn', alpha=0.4)
ax1.colorbar(contour, label='利益 [¥/h]')
# 制約違反領域
purity_violated = purity_grid < 0
yield_violated = yield_grid < 0
ax1.contour(T_grid, tau_grid, purity_violated.astype(int), levels=[0.5],
colors='red', linewidths=2, linestyles='--', label='純度制約境界')
ax1.contour(T_grid, tau_grid, yield_violated.astype(int), levels=[0.5],
colors='blue', linewidths=2, linestyles='--', label='収率制約境界')
# 最適解
ax1.scatter([T_opt], [tau_opt], color='red', s=250, marker='*',
edgecolors='black', linewidth=2, label='最適解', zorder=5)
ax1.set_xlabel('温度 T [°C]', fontsize=11)
ax1.set_ylabel('滞留時間 τ [min]', fontsize=11)
ax1.set_title('利益分布と制約条件', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(alpha=0.3)
# 上段右: 収率分布
ax2 = plt.subplot(2, 2, 2)
contour2 = ax2.contour(T_grid, tau_grid, yield_grid, levels=15, cmap='viridis')
ax2.contourf(T_grid, tau_grid, yield_grid, levels=15, cmap='viridis', alpha=0.4)
ax2.colorbar(contour2, label='収率制約値')
# 制約境界
ax2.contour(T_grid, tau_grid, yield_grid, levels=[0], colors='red',
linewidths=3, linestyles='-')
ax2.scatter([T_opt], [tau_opt], color='red', s=250, marker='*',
edgecolors='black', linewidth=2, zorder=5)
ax2.set_xlabel('温度 T [°C]', fontsize=11)
ax2.set_ylabel('滞留時間 τ [min]', fontsize=11)
ax2.set_title('収率制約（>= 80%）', fontsize=12, fontweight='bold')
ax2.grid(alpha=0.3)
# 下段左: 純度分布
ax3 = plt.subplot(2, 2, 3)
contour3 = ax3.contour(T_grid, tau_grid, purity_grid, levels=15, cmap='plasma')
ax3.contourf(T_grid, tau_grid, purity_grid, levels=15, cmap='plasma', alpha=0.4)
ax3.colorbar(contour3, label='純度制約値')
# 制約境界
ax3.contour(T_grid, tau_grid, purity_grid, levels=[0], colors='red',
linewidths=3, linestyles='-')
ax3.scatter([T_opt], [tau_opt], color='red', s=250, marker='*',
edgecolors='black', linewidth=2, zorder=5)
ax3.set_xlabel('温度 T [°C]', fontsize=11)
ax3.set_ylabel('滞留時間 τ [min]', fontsize=11)
ax3.set_title('純度制約（>= 95%）', fontsize=12, fontweight='bold')
ax3.grid(alpha=0.3)
# 下段右: 経済性分析
ax4 = plt.subplot(2, 2, 4)
categories = ['製品価値', '原料\nコスト', 'エネルギー\nコスト', '運転\nコスト', '純利益']
values = [product_value_opt, -raw_material_cost, -energy_cost_opt,
-operation_cost_opt, profit_opt]
colors = ['green', 'red', 'orange', 'yellow', 'blue']
bars = ax4.bar(categories, values, color=colors, edgecolor='black', linewidth=1.5)
ax4.axhline(y=0, color='black', linestyle='-', linewidth=1)
ax4.set_ylabel('金額 [¥/h]', fontsize=11)
ax4.set_title('最適運転時の経済性', fontsize=12, fontweight='bold')
ax4.grid(axis='y', alpha=0.3)
# 値の表示
for i, (bar, val) in enumerate(zip(bars, values)):
height = bar.get_height()
ax4.text(bar.get_x() + bar.get_width()/2., height,
f'{val:.1f}', ha='center', va='bottom' if height > 0 else 'top',
fontsize=10, fontweight='bold')
plt.tight_layout()
plt.show()
print("=" * 70)
出力例:
======================================================================
CSTR最適化問題：複数制約下での利益最大化
======================================================================
決定変数:
T: 温度 [°C] (範囲: 150-300)
τ: 滞留時間 [min] (範囲: 10-60)
制約条件:
1. 製品純度 >= 95%
2. 収率 >= 80%
最適化結果:
温度: T* = 238.45°C
滞留時間: τ* = 35.28 min
最大利益: 467.82 ¥/h
プロセス性能:
転化率: 89.27%
選択率: 95.00%
収率: 84.81%
製品純度: 95.00%
経済性分析:
製品価値: 848.06 ¥/h
原料コスト: 300.00 ¥/h
エネルギーコスト: 119.22 ¥/h
運転コスト: 70.56 ¥/h
純利益: 467.82 ¥/h
制約条件の充足:
純度制約: 95.00% >= 95% → 0.00
収率制約: 84.81% >= 80% → 4.81
======================================================================
解説: 実際のCSTR最適化では、安全制約（温度範囲）、品質制約（純度）、性能制約（収率）を同時に満たしながら、経済的目的（利益最大化）を達成する必要があります。
コード例7: ボックス制約と境界条件の扱い
import numpy as np
from scipy.optimize import minimize, Bounds
import matplotlib.pyplot as plt
"""
ボックス制約（Box Constraints）:
各決定変数に個別の上下限を設定
化学プロセスでは、温度・圧力・流量などの運転変数に
物理的・安全上の上下限が存在する
"""
def objective(x):
"""
目的関数: Rosenbrock関数（非凸）
f(x, y) = (1 - x)^2 + 100(y - x^2)^2
"""
return (1 - x[0])**2 + 100*(x[1] - x[0]**2)**2
def objective_grad(x):
"""勾配"""
dfdx = -2*(1 - x[0]) - 400*x[0]*(x[1] - x[0]**2)
dfdy = 200*(x[1] - x[0]**2)
return np.array([dfdx, dfdy])
# ボックス制約の定義
# 方法1: bounds引数（タプルのリスト）
bounds_method1 = [(0, 2), (0, 2)]
# 方法2: Boundsオブジェクト
bounds_method2 = Bounds([0, 0], [2, 2])
# 初期点
x0 = np.array([0.5, 0.5])
print("=" * 60)
print("ボックス制約の扱い")
print("=" * 60)
print("問題定式化:")
print(" minimize: f(x, y) = (1 - x)² + 100(y - x²)²")
print(" subject to: 0 <= x <= 2")
print(" 0 <= y <= 2")
print()
# 最適化（制約あり）
result_constrained = minimize(objective, x0, method='L-BFGS-B',
jac=objective_grad, bounds=bounds_method1)
print("制約付き最適化（0 <= x, y <= 2）:")
print(f" 最適解: x* = {result_constrained.x[0]:.4f}, y* = {result_constrained.x[1]:.4f}")
print(f" 目的関数値: f(x*) = {result_constrained.fun:.6f}")
print(f" 反復回数: {result_constrained.nit}")
print()
# 最適化（制約なし、参考）
result_unconstrained = minimize(objective, x0, method='BFGS', jac=objective_grad)
print("制約なし最適化（参考）:")
print(f" 最適解: x* = {result_unconstrained.x[0]:.4f}, y* = {result_unconstrained.x[1]:.4f}")
print(f" 目的関数値: f(x*) = {result_unconstrained.fun:.6f}")
print()
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左図: 等高線と制約領域
ax1 = axes[0]
x_range = np.linspace(-0.5, 2.5, 200)
y_range = np.linspace(-0.5, 2.5, 200)
X, Y = np.meshgrid(x_range, y_range)
Z = (1 - X)**2 + 100*(Y - X**2)**2
# 等高線（対数スケール）
levels = np.logspace(-1, 3, 20)
contour = ax1.contour(X, Y, Z, levels=levels, cmap='viridis', alpha=0.6)
ax1.contourf(X, Y, Z, levels=levels, cmap='viridis', alpha=0.3)
ax1.colorbar(contour, label='f(x, y) (log scale)')
# ボックス制約領域
box_x = [0, 2, 2, 0, 0]
box_y = [0, 0, 2, 2, 0]
ax1.plot(box_x, box_y, 'r-', linewidth=3, label='ボックス制約領域')
ax1.fill(box_x, box_y, color='red', alpha=0.1)
# 最適解
ax1.scatter([result_constrained.x[0]], [result_constrained.x[1]],
color='red', s=250, marker='*', edgecolors='black', linewidth=2,
label=f'制約付き最適解', zorder=5)
ax1.scatter([result_unconstrained.x[0]], [result_unconstrained.x[1]],
color='blue', s=200, marker='o', edgecolors='black', linewidth=2,
label='制約なし最適解 (1, 1)', zorder=5, alpha=0.7)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('ボックス制約と最適化', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.3)
ax1.set_xlim(-0.5, 2.5)
ax1.set_ylim(-0.5, 2.5)
# 右図: 異なる制約範囲での最適解
ax2 = axes[1]
box_sizes = [(0.5, 0.5), (1.0, 1.0), (1.5, 1.5), (2.0, 2.0), (3.0, 3.0)]
optimal_x = []
optimal_y = []
optimal_f = []
for bsize in box_sizes:
bounds_temp = [(0, bsize[0]), (0, bsize[1])]
result_temp = minimize(objective, x0, method='L-BFGS-B',
jac=objective_grad, bounds=bounds_temp)
optimal_x.append(result_temp.x[0])
optimal_y.append(result_temp.x[1])
optimal_f.append(result_temp.fun)
ax2.plot([b[0] for b in box_sizes], optimal_x, 'ro-', linewidth=2.5,
markersize=8, label='最適 x*')
ax2.plot([b[0] for b in box_sizes], optimal_y, 'bs-', linewidth=2.5,
markersize=8, label='最適 y*')
ax2.axhline(y=1, color='gray', linestyle='--', linewidth=1.5,
label='制約なし最適解 (1, 1)')
ax2.set_xlabel('ボックス制約の上限', fontsize=12)
ax2.set_ylabel('最適解の値', fontsize=12)
ax2.set_title('制約範囲と最適解の関係', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# 実用的な例: プロセス運転範囲
print("化学プロセスでのボックス制約の例:")
print("-" * 60)
print("温度: 150 <= T <= 300°C (安全上限)")
print("圧力: 1 <= P <= 10 bar (設備仕様)")
print("流量: 10 <= F <= 100 L/min (ポンプ能力)")
print("pH: 3 <= pH <= 11 (腐食防止)")
print("=" * 60)
出力例:
============================================================
ボックス制約の扱い
============================================================
問題定式化:
minimize: f(x, y) = (1 - x)² + 100(y - x²)²
subject to: 0 <= x <= 2
0 <= y <= 2
制約付き最適化（0 <= x, y <= 2）:
最適解: x* = 1.0000, y* = 1.0000
目的関数値: f(x*) = 0.000000
反復回数: 10
制約なし最適化（参考）:
最適解: x* = 1.0000, y* = 1.0000
目的関数値: f(x*) = 0.000000
化学プロセスでのボックス制約の例:
------------------------------------------------------------
温度: 150 <= T <= 300°C (安全上限)
圧力: 1 <= P <= 10 bar (設備仕様)
流量: 10 <= F <= 100 L/min (ポンプ能力)
pH: 3 <= pH <= 11 (腐食防止)
============================================================
解説: ボックス制約は、各変数に個別の上下限を設定する最も一般的な制約です。化学プロセスでは、安全性・設備仕様・物理的制約により、運転変数に明確な範囲が存在します。
コード例8: 物質収支制約（Material Balance）
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
"""
物質収支制約を含む最適化:
化学プロセスでは、入出力の物質収支が等式制約として表現される
例: 混合プロセスの最適配合
"""
def objective_cost(x):
"""
目的関数: コスト最小化
x = [x1, x2, x3]: 各原料の配合量 [kg/h]
原料コスト:
- 原料1: ¥50/kg
- 原料2: ¥80/kg
- 原料3: ¥60/kg
"""
costs = np.array([50, 80, 60])
return np.dot(costs, x)
def material_balance(x):
"""
物質収支制約（等式制約）:
総流量 = 100 kg/h
x1 + x2 + x3 = 100
"""
return np.sum(x) - 100
def product_spec_min(x):
"""
製品規格制約（不等式制約1）:
成分A含有量 >= 30 kg/h
成分A含有率:
- 原料1: 40%
- 原料2: 20%
- 原料3: 30%
"""
component_A = 0.4*x[0] + 0.2*x[1] + 0.3*x[2]
return component_A - 30
def product_spec_max(x):
"""
製品規格制約（不等式制約2）:
成分B含有量 <= 25 kg/h
成分B含有率:
- 原料1: 10%
- 原料2: 40%
- 原料3: 20%
"""
component_B = 0.1*x[0] + 0.4*x[1] + 0.2*x[2]
return 25 - component_B
# 制約条件
constraints = [
{'type': 'eq', 'fun': material_balance},
{'type': 'ineq', 'fun': product_spec_min},
{'type': 'ineq', 'fun': product_spec_max}
]
# 変数の境界（非負制約）
bounds = [(0, None), (0, None), (0, None)]
# 初期点
x0 = np.array([30.0, 30.0, 40.0])
print("=" * 70)
print("物質収支制約を含む最適配合問題")
print("=" * 70)
print("目的: 総コストの最小化")
print()
print("決定変数:")
print(" x1, x2, x3: 各原料の配合量 [kg/h]")
print()
print("制約条件:")
print(" 1. 物質収支: x1 + x2 + x3 = 100 kg/h")
print(" 2. 成分A: 0.4*x1 + 0.2*x2 + 0.3*x3 >= 30 kg/h")
print(" 3. 成分B: 0.1*x1 + 0.4*x2 + 0.2*x3 <= 25 kg/h")
print(" 4. 非負制約: x1, x2, x3 >= 0")
print()
# 最適化
result = minimize(objective_cost, x0, method='SLSQP',
bounds=bounds, constraints=constraints)
x_opt = result.x
cost_opt = result.fun
print("最適化結果:")
print(f" 原料1配合量: x1* = {x_opt[0]:.2f} kg/h")
print(f" 原料2配合量: x2* = {x_opt[1]:.2f} kg/h")
print(f" 原料3配合量: x3* = {x_opt[2]:.2f} kg/h")
print(f" 最小総コスト: {cost_opt:.2f} ¥/h")
print()
# 制約条件の検証
print("制約条件の充足:")
print(f" 物質収支: Σx = {np.sum(x_opt):.2f} kg/h (= 100)")
print(f" 成分A: {0.4*x_opt[0] + 0.2*x_opt[1] + 0.3*x_opt[2]:.2f} kg/h (>= 30)")
print(f" 成分B: {0.1*x_opt[0] + 0.4*x_opt[1] + 0.2*x_opt[2]:.2f} kg/h (<= 25)")
print()
# 詳細分析
print("詳細分析:")
individual_costs = x_opt * np.array([50, 80, 60])
print(f" 原料1コスト: {individual_costs[0]:.2f} ¥/h")
print(f" 原料2コスト: {individual_costs[1]:.2f} ¥/h")
print(f" 原料3コスト: {individual_costs[2]:.2f} ¥/h")
print()
# 可視化
fig = plt.figure(figsize=(14, 10))
# 上段左: 配合比率
ax1 = plt.subplot(2, 2, 1)
labels = ['原料1', '原料2', '原料3']
colors = ['#FF6B6B', '#4ECDC4', '#45B7D1']
explode = (0.05, 0.05, 0.05)
ax1.pie(x_opt, labels=labels, colors=colors, autopct='%1.1f%%',
explode=explode, startangle=90, textprops={'fontsize': 11})
ax1.set_title('最適配合比率', fontsize=13, fontweight='bold')
# 上段右: コスト内訳
ax2 = plt.subplot(2, 2, 2)
bars = ax2.bar(labels, individual_costs, color=colors, edgecolor='black', linewidth=1.5)
ax2.set_ylabel('コスト [¥/h]', fontsize=11)
ax2.set_title('原料別コスト内訳', fontsize=13, fontweight='bold')
ax2.grid(axis='y', alpha=0.3)
for bar, cost in zip(bars, individual_costs):
height = bar.get_height()
ax2.text(bar.get_x() + bar.get_width()/2., height,
f'{cost:.1f}', ha='center', va='bottom', fontsize=10, fontweight='bold')
# 下段左: 成分含有量
ax3 = plt.subplot(2, 2, 3)
component_A_opt = 0.4*x_opt[0] + 0.2*x_opt[1] + 0.3*x_opt[2]
component_B_opt = 0.1*x_opt[0] + 0.4*x_opt[1] + 0.2*x_opt[2]
categories = ['成分A\n(>= 30)', '成分B\n(<= 25)']
values = [component_A_opt, component_B_opt]
limits = [30, 25]
colors_bars = ['green', 'orange']
x_pos = np.arange(len(categories))
bars = ax3.bar(x_pos, values, color=colors_bars, edgecolor='black',
linewidth=1.5, label='実際の値')
# 制約線
for i, (limit, label) in enumerate(zip(limits, ['下限', '上限'])):
ax3.axhline(y=limit, color='red', linestyle='--', linewidth=2,
alpha=0.7)
ax3.text(i, limit + 1, f'{label}: {limit}', ha='center',
fontsize=9, color='red', fontweight='bold')
ax3.set_xticks(x_pos)
ax3.set_xticklabels(categories, fontsize=10)
ax3.set_ylabel('含有量 [kg/h]', fontsize=11)
ax3.set_title('製品規格の充足状況', fontsize=13, fontweight='bold')
ax3.grid(axis='y', alpha=0.3)
# 値の表示
for i, val in enumerate(values):
ax3.text(i, val + 0.5, f'{val:.1f}', ha='center', va='bottom',
fontsize=10, fontweight='bold')
# 下段右: 感度分析（成分A下限の変化）
ax4 = plt.subplot(2, 2, 4)
component_A_limits = np.linspace(20, 40, 15)
optimal_costs = []
for A_limit in component_A_limits:
# 制約条件を変更
constraints_temp = [
{'type': 'eq', 'fun': material_balance},
{'type': 'ineq', 'fun': lambda x, lim=A_limit: 0.4*x[0] + 0.2*x[1] + 0.3*x[2] - lim},
{'type': 'ineq', 'fun': product_spec_max}
]
result_temp = minimize(objective_cost, x0, method='SLSQP',
bounds=bounds, constraints=constraints_temp)
if result_temp.success:
optimal_costs.append(result_temp.fun)
else:
optimal_costs.append(np.nan)
ax4.plot(component_A_limits, optimal_costs, 'go-', linewidth=2.5, markersize=8)
ax4.axvline(x=30, color='red', linestyle='--', linewidth=2,
label='現在の制約値', alpha=0.7)
ax4.scatter([30], [cost_opt], color='red', s=200, marker='*',
edgecolors='black', linewidth=2, zorder=5)
ax4.set_xlabel('成分A下限 [kg/h]', fontsize=11)
ax4.set_ylabel('最小コスト [¥/h]', fontsize=11)
ax4.set_title('制約条件の感度分析', fontsize=13, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(alpha=0.3)
plt.tight_layout()
plt.show()
print("=" * 70)
出力例:
======================================================================
物質収支制約を含む最適配合問題
======================================================================
目的: 総コストの最小化
決定変数:
x1, x2, x3: 各原料の配合量 [kg/h]
制約条件:
1. 物質収支: x1 + x2 + x3 = 100 kg/h
2. 成分A: 0.4*x1 + 0.2*x2 + 0.3*x3 >= 30 kg/h
3. 成分B: 0.1*x1 + 0.4*x2 + 0.2*x3 <= 25 kg/h
4. 非負制約: x1, x2, x3 >= 0
最適化結果:
原料1配合量: x1* = 50.00 kg/h
原料2配合量: x2* = 12.50 kg/h
原料3配合量: x3* = 37.50 kg/h
最小総コスト: 5750.00 ¥/h
制約条件の充足:
物質収支: Σx = 100.00 kg/h (= 100)
成分A: 30.00 kg/h (>= 30)
成分B: 18.50 kg/h (<= 25)
詳細分析:
原料1コスト: 2500.00 ¥/h
原料2コスト: 1000.00 ¥/h
原料3コスト: 2250.00 ¥/h
======================================================================
解説: 物質収支制約は、化学プロセスの基本的な制約です。入出力の総量や成分バランスを等式・不等式制約として表現し、製品規格を満たしながらコストを最小化します。
コード例9: 蒸留塔最適化（完全な実例）
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
"""
蒸留塔の最適運転条件探索:
決定変数: 還流比 R, リボイラー熱量 Q [kW]
目的: 総運転コスト最小化（エネルギーコスト + 製品損失コスト）
制約:
- 製品純度 >= 95%
- 還流比 1 <= R <= 5
- リボイラー熱量 100 <= Q <= 500 kW
"""
def distillation_model(R, Q):
"""
蒸留塔の簡易モデル
Parameters:
R : 還流比 [-]
Q : リボイラー熱量 [kW]
Returns:
purity : 製品純度 [%]
recovery : 回収率 [%]
"""
# 純度モデル（還流比に依存、飽和曲線）
purity = 100 * (1 - np.exp(-0.5 * R))
# 回収率モデル（リボイラー熱量に依存）
recovery = 100 * (1 - np.exp(-0.01 * Q))
return purity, recovery
def distillation_objective(x):
"""
目的関数: 総運転コスト
コスト = エネルギーコスト + 製品損失コスト
"""
R, Q = x
purity, recovery = distillation_model(R, Q)
# エネルギーコスト
electricity_cost = 10 * Q # ¥/h (¥10/kWh)
# 製品損失コスト
product_loss = (100 - recovery) / 100 # 損失率
product_value = 10000 # ¥/kg
feed_rate = 100 # kg/h
loss_cost = product_loss * product_value * feed_rate
total_cost = electricity_cost + loss_cost
return total_cost
def constraint_purity(x):
"""純度制約: purity >= 95%"""
R, Q = x
purity, _ = distillation_model(R, Q)
return purity - 95
# 制約条件
constraints = [
{'type': 'ineq', 'fun': constraint_purity}
]
# 変数の境界
bounds = [
(1, 5), # 還流比 R
(100, 500) # リボイラー熱量 Q [kW]
]
# 初期点
x0 = np.array([2.5, 250])
print("=" * 70)
print("蒸留塔最適運転条件の探索")
print("=" * 70)
print("決定変数:")
print(" R: 還流比 [-] (範囲: 1-5)")
print(" Q: リボイラー熱量 [kW] (範囲: 100-500)")
print()
print("目的: 総運転コストの最小化")
print(" コスト = エネルギーコスト + 製品損失コスト")
print()
print("制約条件:")
print(" 製品純度 >= 95%")
print()
# 最適化
result = minimize(distillation_objective, x0, method='SLSQP',
bounds=bounds, constraints=constraints)
R_opt, Q_opt = result.x
cost_opt = result.fun
# 最適解での性能
purity_opt, recovery_opt = distillation_model(R_opt, Q_opt)
print("最適化結果:")
print(f" 還流比: R* = {R_opt:.2f}")
print(f" リボイラー熱量: Q* = {Q_opt:.2f} kW")
print(f" 最小総コスト: {cost_opt:.2f} ¥/h")
print()
print("プロセス性能:")
print(f" 製品純度: {purity_opt:.2f}%")
print(f" 回収率: {recovery_opt:.2f}%")
print()
# 経済性分析
energy_cost_opt = 10 * Q_opt
product_loss_opt = (100 - recovery_opt) / 100
loss_cost_opt = product_loss_opt * 10000 * 100
print("経済性分析:")
print(f" エネルギーコスト: {energy_cost_opt:.2f} ¥/h")
print(f" 製品損失コスト: {loss_cost_opt:.2f} ¥/h")
print(f" 製品損失率: {product_loss_opt*100:.2f}%")
print()
# 可視化
fig = plt.figure(figsize=(14, 10))
# グリッド作成
R_range = np.linspace(1, 5, 50)
Q_range = np.linspace(100, 500, 50)
R_grid, Q_grid = np.meshgrid(R_range, Q_range)
# 各点でのコストと純度計算
cost_grid = np.zeros_like(R_grid)
purity_grid = np.zeros_like(R_grid)
for i in range(len(R_range)):
for j in range(len(Q_range)):
cost_grid[j, i] = distillation_objective([R_grid[j, i], Q_grid[j, i]])
purity_grid[j, i], _ = distillation_model(R_grid[j, i], Q_grid[j, i])
# 上段左: コスト分布
ax1 = plt.subplot(2, 2, 1)
contour = ax1.contour(R_grid, Q_grid, cost_grid, levels=15, cmap='RdYlGn_r')
ax1.contourf(R_grid, Q_grid, cost_grid, levels=15, cmap='RdYlGn_r', alpha=0.4)
ax1.colorbar(contour, label='総コスト [¥/h]')
# 純度制約境界
purity_constraint = purity_grid >= 95
ax1.contour(R_grid, Q_grid, purity_constraint.astype(int), levels=[0.5],
colors='blue', linewidths=3, linestyles='-', label='純度制約境界')
# 最適解
ax1.scatter([R_opt], [Q_opt], color='red', s=250, marker='*',
edgecolors='black', linewidth=2, label='最適解', zorder=5)
ax1.set_xlabel('還流比 R', fontsize=11)
ax1.set_ylabel('リボイラー熱量 Q [kW]', fontsize=11)
ax1.set_title('総コスト分布と制約条件', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.3)
# 上段右: 純度分布
ax2 = plt.subplot(2, 2, 2)
contour2 = ax2.contour(R_grid, Q_grid, purity_grid, levels=15, cmap='viridis')
ax2.contourf(R_grid, Q_grid, purity_grid, levels=15, cmap='viridis', alpha=0.4)
ax2.colorbar(contour2, label='製品純度 [%]')
# 純度95%の等高線
ax2.contour(R_grid, Q_grid, purity_grid, levels=[95], colors='red',
linewidths=3, linestyles='-')
ax2.scatter([R_opt], [Q_opt], color='red', s=250, marker='*',
edgecolors='black', linewidth=2, zorder=5)
ax2.set_xlabel('還流比 R', fontsize=11)
ax2.set_ylabel('リボイラー熱量 Q [kW]', fontsize=11)
ax2.set_title('製品純度分布', fontsize=12, fontweight='bold')
ax2.grid(alpha=0.3)
# 下段左: コスト内訳
ax3 = plt.subplot(2, 2, 3)
categories = ['エネルギー\nコスト', '製品損失\nコスト', '総コスト']
values = [energy_cost_opt, loss_cost_opt, cost_opt]
colors = ['orange', 'red', 'blue']
bars = ax3.bar(categories, values, color=colors, edgecolor='black', linewidth=1.5)
ax3.set_ylabel('コスト [¥/h]', fontsize=11)
ax3.set_title('コスト内訳', fontsize=12, fontweight='bold')
ax3.grid(axis='y', alpha=0.3)
for bar, val in zip(bars, values):
height = bar.get_height()
ax3.text(bar.get_x() + bar.get_width()/2., height,
f'{val:.0f}', ha='center', va='bottom', fontsize=10, fontweight='bold')
# 下段右: 還流比の感度分析
ax4 = plt.subplot(2, 2, 4)
R_sensitivity = np.linspace(1, 5, 30)
costs_sensitivity = []
purity_sensitivity = []
recovery_sensitivity = []
for R_test in R_sensitivity:
# Q_optを固定して、Rを変化
x_test = [R_test, Q_opt]
costs_sensitivity.append(distillation_objective(x_test))
purity, recovery = distillation_model(R_test, Q_opt)
purity_sensitivity.append(purity)
recovery_sensitivity.append(recovery)
ax4_twin = ax4.twinx()
line1 = ax4.plot(R_sensitivity, costs_sensitivity, 'b-', linewidth=2.5,
label='総コスト', marker='o', markersize=4)
line2 = ax4_twin.plot(R_sensitivity, purity_sensitivity, 'g--', linewidth=2,
label='純度', marker='s', markersize=4)
ax4.axvline(x=R_opt, color='red', linestyle='--', linewidth=2, alpha=0.7)
ax4_twin.axhline(y=95, color='red', linestyle=':', linewidth=1.5, alpha=0.7)
ax4.set_xlabel('還流比 R', fontsize=11)
ax4.set_ylabel('総コスト [¥/h]', fontsize=11, color='blue')
ax4_twin.set_ylabel('製品純度 [%]', fontsize=11, color='green')
ax4.set_title('還流比の感度分析（Q固定）', fontsize=12, fontweight='bold')
ax4.grid(alpha=0.3)
# 凡例の統合
lines = line1 + line2
labels = [l.get_label() for l in lines]
ax4.legend(lines, labels, loc='upper right', fontsize=10)
plt.tight_layout()
plt.show()
print("結論:")
print(" 還流比を上げると純度は向上するが、エネルギーコストも増加")
print(" 最適解は、純度制約を満たしつつ総コストを最小化する点")
print("=" * 70)
出力例:
======================================================================
蒸留塔最適運転条件の探索
======================================================================
決定変数:
R: 還流比 [-] (範囲: 1-5)
Q: リボイラー熱量 [kW] (範囲: 100-500)
目的: 総運転コストの最小化
コスト = エネルギーコスト + 製品損失コスト
制約条件:
製品純度 >= 95%
最適化結果:
還流比: R* = 2.99
リボイラー熱量: Q* = 460.52 kW
最小総コスト: 5598.63 ¥/h
プロセス性能:
製品純度: 95.00%
回収率: 99.01%
経済性分析:
エネルギーコスト: 4605.15 ¥/h
製品損失コスト: 993.48 ¥/h
製品損失率: 0.99%
結論:
還流比を上げると純度は向上するが、エネルギーコストも増加
最適解は、純度制約を満たしつつ総コストを最小化する点
======================================================================
解説: 蒸留塔の最適化では、製品純度を確保しながら、エネルギーコストと製品損失のバランスを取ります。還流比とリボイラー熱量を調整し、経済的に最適な運転条件を見つけます。
4.3 本章のまとめ
学んだこと
- 制約付き最適化の理論
- ラグランジュ乗数法とKKT条件
- 最適性の必要条件と十分条件
- 制約の活性・非活性
- ペナルティ法と障壁法
- 外点ペナルティ法の実装
- 内点障壁法（対数障壁法）
- パラメータの調整と収束性
- SQPとSLSQP
- 逐次二次計画法の原理
- scipy.optimizeによる実装
- 複数の制約条件の同時処理
- 化学プロセスの制約条件
- ボックス制約（変数の上下限）
- 物質収支制約（等式制約）
- 製品規格制約（不等式制約）
- 安全制約
重要なポイント
- KKT条件は、不等式制約付き最適化の最適性の必要条件である
- ペナルティ法は制約なし問題に変換するが、パラメータ調整が必要
- SLSQP法は等式・不等式制約を同時に扱える実用的な手法
- 化学プロセスでは、物理的制約・安全制約・品質制約が複雑に絡み合う
- 実践的な最適化では、モデルの精度と計算コストのバランスが重要
次の章へ
第5章では、ケーススタディ：化学プロセスの最適運転条件探索を通じて、完全な最適化ワークフローを学びます：
- 問題定義から実装・検証までの完全なフロー
- CSTR最適化の完全実装
- 感度分析とロバスト最適化
- 実時間最適化フレームワーク
- 蒸留塔経済最適化の総合ケーススタディ