pyomo

Pyomo - 数理最適化モデリング

Pyomoを使うと、自然な数学的構文（集合、パラメーター、変数、制約）で最適化問題を定義できます。モデルとソルバーを分離できるため、コードを変更せずに同じモデルを異なるエンジンで解くことができます。

最初に：前提条件を確認する

pip install pyomo
# Also install a solver (e.g., GLPK for linear/integer problems)
# Conda: conda install -c conda-forge glpk ipopt

使用する場面

• 戦略計画：長期的なリソース配分や投資計画。
• プロセスエンジニアリング：化学プラントや製油所の操業の最適化（非線形）。
• エネルギーシステム：電力網のディスパッチやユニットコミットメント問題。
• サプライチェーン最適化：複数期間・複数品目のフロー問題。
• 非線形計画法（NLP）：制約や目的関数に滑らかな曲線（例：x²、log(x)）が含まれる場合。
• 確率計画法：最適化における不確実性のモデリング。
• カスタムソルバーの統合：IPOPT、SCIP、Baronなどの特定のソルバーを使用する必要がある場合。

リファレンスドキュメント

公式ドキュメント：http://www.pyomo.org/
GitHub：https://github.com/Pyomo/pyomo
検索パターン：pyo.ConcreteModel、pyo.Constraint、pyo.Objective、pyo.SolverFactory

基本原則

Concrete モデルと Abstract モデル

• ConcreteModel：モデルの構築時にデータを定義します（Python／データサイエンスでは最も一般的です）。
• AbstractModel：最初に構造を定義し、その後でデータを読み込みます（大規模な産業用モデルで標準的です）。

コンポーネント

• Var：ソルバーが求める未知数です。
• Set/Param：問題のインスタンスを定義するデータです。
• Objective：最小化または最大化する関数です。
• Constraint：変数が従わなければならない規則です。

ソルバー

Pyomoには独自のソルバーはありません。システムにインストールされた外部ソルバー（LP／MIP用のglpkやNLP用のipoptなど）が必要です。

クイックリファレンス

インストール

pip install pyomo
# Also install a solver (e.g., GLPK for linear/integer problems)
# Conda: conda install -c conda-forge glpk ipopt

標準的なインポート

import pyomo.environ as pyo
from pyomo.opt import SolverFactory

基本パターン - シンプルな線形モデル（Concrete）

import pyomo.environ as pyo

# 1. Create Model
model = pyo.ConcreteModel()

# 2. Define Variables
model.x = pyo.Var(within=pyo.NonNegativeReals)
model.y = pyo.Var(within=pyo.NonNegativeReals)

# 3. Define Objective (Minimize x + 2*y)
model.obj = pyo.Objective(expr=model.x + 2*model.y, sense=pyo.minimize)

# 4. Define Constraints
model.con1 = pyo.Constraint(expr=3*model.x + 4*model.y >= 12)
model.con2 = pyo.Constraint(expr=2*model.x + 5*model.y >= 10)

# 5. Solve
solver = pyo.SolverFactory('glpk')
results = solver.solve(model)

# 6. Access Results
print(f"x = {pyo.value(model.x)}, y = {pyo.value(model.y)}")

重要なルール

✅ 推奨事項

• 問題をスケーリングする - あらゆる最適化と同様に、変数と制約の桁数が近い範囲（例：0.1～100）にあると、Pyomoは最も効果的に動作します。
• pyo.value()を使う - 求解後に数値結果を得るには、変数を必ずpyo.value(model.x)でラップします。
• ソルバーのステータスを確認する - results.solver.termination_conditionを使って、ソルバーが最適解に到達したことを確認します。
• pyo.Setとpyo.Paramを使う - 複雑なモデルでは、通常のPythonのリストや辞書ではなく、これらのコンポーネントを使ってデータを整理します。
• ソルバーのオプションを指定する - 実行時間を管理するには、solver.options['tm_limit'] = 60（またはソルバーに応じた同様の設定）を使います。
• model.dualを使う - ソルバーが対応している場合は、双対値を取得できるように双対サフィックスを有効にします（経済分析で重要です）。

❌ 避けるべき事項

• Pythonの数学関数を使わない - 式の中ではmath.log()やnp.log()ではなく、pyo.log()、pyo.exp()、pyo.sqrt()を使います。
• 制約内でif文を使わない - 制約は代数的でなければなりません。条件ロジックには「Big-M」表記や論理的な再定式化を使います。
• 「扱いにくい」非線形性を避ける - abs(x)やmax(x, y)などの関数は非滑らかで、多くのソルバーで問題を引き起こす可能性があります。再定式化を使います。
• 切片／定数項を忘れない - x + y = 10のような制約は問題ありませんが、単位に一貫性があることを確認します。

アンチパターン（絶対に避ける）

import pyomo.environ as pyo
import numpy as np

# ❌ BAD: Using NumPy/Math functions in constraints
# model.con = pyo.Constraint(expr=np.sin(model.x) <= 0.5)

# ✅ GOOD: Use Pyomo-compatible functions
model.con = pyo.Constraint(expr=pyo.sin(model.x) <= 0.5)

# ❌ BAD: Using Python IF for conditional constraints
# if model.x > 10:
#     model.con = pyo.Constraint(expr=model.y <= 5)

# ✅ GOOD: Using Big-M notation (for binary variable z)
# y <= 5 + M * (1 - z)
# x >= 10 - M * (1 - z)

# ❌ BAD: Printing a variable directly
# print(model.x) # Returns a reference object, not a number!

# ✅ GOOD: Use value()
print(pyo.value(model.x))

集合とインデックスを使ったモデリング

大規模モデルの構築

model = pyo.ConcreteModel()

# Data
products = ['A', 'B', 'C']
profit = {'A': 10, 'B': 20, 'C': 15}
limit = 100

# Components
model.P = pyo.Set(initialize=products)
model.x = pyo.Var(model.P, within=pyo.NonNegativeReals)

# Indexed Objective
def obj_rule(model):
    return sum(profit[p] * model.x[p] for p in model.P)
model.obj = pyo.Objective(rule=obj_rule, sense=pyo.maximize)

# Indexed Constraint
def limit_rule(model):
    return sum(model.x[p] for p in model.P) <= limit
model.con = pyo.Constraint(rule=limit_rule)

非線形計画法（NLP）

IPOPTソルバーの使用

# Solving: Minimize (x-2)^2 + (y-2)^2
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=0) # Initialization is CRUCIAL for NLP
model.y = pyo.Var(initialize=0)

model.obj = pyo.Objective(expr=(model.x - 2)**2 + (model.y - 2)**2)

# Constraint: x^2 + y <= 1
model.con = pyo.Constraint(expr=model.x**2 + model.y <= 1)

# Solve with IPOPT
solver = pyo.SolverFactory('ipopt')
solver.solve(model)

混合整数線形計画法（MILP）

離散的な意思決定

# Binary variable: 1 if we open a warehouse, 0 otherwise
model.use_warehouse = pyo.Var(within=pyo.Binary)

# Integer variable: Number of trucks to buy
model.num_trucks = pyo.Var(within=pyo.NonNegativeIntegers)

# Conditional logic: If warehouse is not used, trucks must be 0
# trucks <= Capacity * use_warehouse
model.cap_con = pyo.Constraint(expr=model.num_trucks <= 100 * model.use_warehouse)

実践的なワークフロー

1. 食事問題（古典的な線形計画法）

def solve_diet(foods, nutrients, costs, requirements):
    model = pyo.ConcreteModel()
    model.F = pyo.Set(initialize=foods)
    model.N = pyo.Set(initialize=nutrients)
    
    model.x = pyo.Var(model.F, within=pyo.NonNegativeReals)
    
    # Minimize cost
    model.obj = pyo.Objective(expr=sum(costs[f] * model.x[f] for f in model.F))
    
    # Meet nutrient requirements
    def nutrient_rule(model, n):
        return sum(nutrients[f][n] * model.x[f] for f in model.F) >= requirements[n]
    model.con = pyo.Constraint(model.N, rule=nutrient_rule)
    
    pyo.SolverFactory('glpk').solve(model)
    return {f: pyo.value(model.x[f]) for f in model.F}

2. ブレンド問題（化学／プロセスエンジニアリング）

# Balancing component fractions in a mixture
# Note: Often becomes non-linear (NLP) if both flow and fraction are variables
def blend_optimization(inputs, target_purity):
    model = pyo.ConcreteModel()
    # ... model setup ...
    # con: sum(flow[i] * purity[i]) / sum(flow[i]) == target_purity
    # becomes: sum(flow[i] * purity[i]) == target_purity * sum(flow[i]) (Linearized)

パフォーマンスの最適化

ソルバーの選択

• GLPK/CBC：無償／オープンソースのLP／MIPに適しています。
• Gurobi/CPLEX：産業標準です（大規模なMIPを非常に高速に解きます）。
• IPOPT：滑らかな非線形問題（NLP）に最適です。

ウォームスタート

反復最適化では、前回の解を初期値として使用してください。

# For NLP solvers like IPOPT
model.x.set_value(prev_x_value)
solver.solve(model)

よくある落とし穴と解決策

終了条件

ソルバーが停止した理由を必ず確認してください。

from pyomo.opt import TerminationCondition

results = solver.solve(model)
if results.solver.termination_condition == TerminationCondition.optimal:
    print("Success")
elif results.solver.termination_condition == TerminationCondition.infeasible:
    print("Check your constraints!")

NLP における非凸性

NLP モデルに複数の局所最小値がある場合、IPOPT は局所最小値から抜け出せなくなることがあります。

# ✅ Solution: 
# 1. Provide multiple different initial guesses (multistart).
# 2. Use a global solver like BARON or SCIP.

変数によるインデックス指定

model.x[model.y] のように、y が Var の場合は無効です。

# ✅ Solution: Use model.AddElement or binary variable reformulations.

ベストプラクティス

1. 問題は必ずスケーリングしてください。変数と制約の桁数は、同程度にする必要があります
2. 求解後に変数から数値結果を取り出すには、pyo.value() を使用してください
3. 最適性の確認や問題の診断のために、ソルバーの終了条件を確認してください
4. 複雑なモデルのデータ整理には、pyo.Set と pyo.Param を使用してください
5. NLP 問題では変数を初期化してください。適切な初期値が重要です
6. 式では NumPy や math の関数ではなく、Pyomo に対応した関数（pyo.log()、pyo.exp() など）を使用してください
7. 条件分岐のロジックは、Big-M 記法またはバイナリ変数を使って再定式化してください
8. abs() や max() などの非滑らかな関数は避け、再定式化してください
9. 経済分析でシャドープライスが必要な場合は、dual suffix を有効にしてください
10. 問題の規模に応じて、ソルバーの制限時間とオプションを適切に設定してください

Pyomo は、高水準の数学的抽象化を、解が得られるビジネスや科学の問題へと変換するための究極のツールです。代数モデリングと高性能ソルバーの間をつなぐ機能により、現代の規範的分析の基盤となっています。