Python数理最適化ライブラリPyomo詳解

Pythonで複雑な最適化問題をモデル化し解決するための強力なツール

Pyomoとは？
- 数理最適化とは？
- Pyomoを使うメリット
インストール
基本的な概念とワークフロー
- ConcreteModel vs AbstractModel
- 主要なモデリングコンポーネント
モデリングコンポーネント詳解
簡単な例題: 線形計画問題（生産計画）
- 問題設定
- Pyomoコード (ConcreteModel)
ソルバーとの連携
求解と結果の分析
高度な機能（一部紹介）
応用分野と事例
まとめ

Pyomoとは？

Pyomo (Python Optimization Modeling Objects) は、Pythonベースのオープンソース数理最適化モデリング言語です。線形計画法(LP)、混合整数計画法(MIP)、非線形計画法(NLP)、混合整数非線形計画法(MINLP)など、多様な最適化問題に対応しています。Pyomoを使うことで、研究者やエンジニアは複雑な実世界の応用問題をPythonの構文を用いて直感的かつ数学的に厳密な方法でモデル化し、様々なソルバーを用いて解くことができます。

Pyomoは、以前はCooprソフトウェアライブラリの一部としてリリースされていました。BSDライセンスの下で利用可能です。

数理最適化とは？

数理最適化は、与えられた制約条件の下で、特定の目的関数（コスト、利益、効率など）を最大化または最小化する変数の値を見つけ出す数学的な手法です。サプライチェーン管理、生産計画、スケジューリング、金融ポートフォリオ最適化、エネルギーシステム設計、機械学習モデルのパラメータ調整など、幅広い分野で活用されています。

最適化問題は通常、以下の3つの要素で構成されます。

目的関数 (Objective Function): 最大化または最小化したい対象を表す数式。
決定変数 (Decision Variables): 最適化プロセスで値を決定する変数。
制約条件 (Constraints): 決定変数が満たすべき等式または不等式。

Pyomoを使うメリット

Pythonicな構文: Pythonに慣れているユーザーにとって学習曲線が緩やかです。標準的なPython構文を使用し、Pythonの豊富なライブラリ（データ処理、可視化など）と連携できます。
高い表現力と柔軟性: 複雑な数理モデルを直感的に表現できます。標準的なPythonの構造を使用して、複雑な関係性をモデル化できます。
幅広い問題タイプへの対応: LP, MILP, NLP, MINLP, GDP (Generalized Disjunctive Programming), DAE (Differential Algebraic Equations), MPEC (Mathematical Programming with Equilibrium Constraints) など、多様な問題タイプをサポートしています。
ソルバー連携: GLPK, CBC, IPOPTなどのオープンソースソルバーや、Gurobi, CPLEXなどの商用ソルバーと簡単に連携できます。SolverFactoryを通じて、統一されたインターフェースでソルバーを切り替えることが可能です。
拡張性: Pyomoは高レベルの最適化・分析ツール開発のフレームワークとしても機能します。例えば、確率計画法のための`mpi-sppy`パッケージなど、特定の領域向けの拡張パッケージが存在します。
活発なコミュニティ: 継続的な開発と改善が行われており、ユーザーコミュニティによるサポートや事例共有も活発です。

インストール

前提条件

Pyomoを利用するには、Pythonがインストールされている必要があります。科学技術計算向けのPythonディストリビューション（Anacondaなど）を利用すると、NumPy, SciPy, Matplotlibなどの便利なライブラリも一緒にインストールされるため推奨されます。

Pyomoは現在、CPython 3.9から3.13およびPyPy 3.10でテストされています。

Pyomo本体のインストール

Pyomoはpipを使って簡単にインストールできます。ターミナルまたはコマンドプロンプトで以下のコマンドを実行します。

pip install pyomo

Conda環境を使用している場合は、conda-forgeチャンネルからインストールすることもできます。

conda install -c conda-forge pyomo

ソルバーのインストール

Pyomoはモデリング言語であり、実際に問題を解くためには別途「ソルバー」が必要です。Pyomoはソルバーを同梱していないため、ユーザー自身でインストールする必要があります。どのソルバーを選択するかは、解きたい問題の種類（線形、非線形、整数など）やライセンス（オープンソース、商用）によって異なります。

代表的なオープンソースソルバー:

GLPK (GNU Linear Programming Kit): 線形計画(LP)、混合整数線形計画(MILP)に対応。多くのLinuxディストリビューションやmacOSのパッケージマネージャ (apt, yum, brewなど) でインストールできます。
```
# macOS (Homebrew) の場合
brew install glpk

# Debian/Ubuntu の場合
sudo apt-get update
sudo apt-get install glpk-utils

# Conda の場合
conda install -c conda-forge glpk
```

CBC (COIN-OR Branch and Cut): LP, MILPに対応。GLPKよりも高速な場合があります。

# macOS (Homebrew) の場合
brew install cbc

# Debian/Ubuntu の場合
sudo apt-get install coinor-cbc

# Conda の場合
conda install -c conda-forge coincbc

IPOPT (Interior Point Optimizer): 非線形計画(NLP)に対応。大規模なNLP問題に強力です。

# Conda の場合 (推奨)
conda install -c conda-forge ipopt

# ソースからビルドすることも可能 (上級者向け)

Bonmin / Couenne: 混合整数非線形計画(MINLP)に対応するオープンソースソルバー。
```
# Conda の場合 (推奨)
conda install -c conda-forge bonmin couenne
```
SCIP: LP, MILP, MINLPに対応する高速な非商用ソルバー。アカデミック用途では強力な選択肢です。インストール方法は公式サイトを参照してください。

代表的な商用ソルバー:

Gurobi: 高性能なLP, MILP, QP, MIQPソルバー。アカデミックライセンスは無料。
CPLEX: IBM製の高性能ソルバー。Gurobiと同様の分野で広く使われています。アカデミック用途では無料版あり。

これらのソルバーは、それぞれの公式サイトからダウンロードおよびインストール手順を確認してください。多くの場合、インストール後に環境変数（PATHやライセンスファイルの場所など）の設定が必要です。 Pyomoがソルバーを見つけられるように、ソルバーの実行可能ファイルがあるディレクトリにPATHを通すか、Pyomoのコード内でソルバーの実行可能ファイルのフルパスを指定する必要があります。

注意: ソルバーのインストールは、OSや環境によって手順が異なる場合があります。各ソルバーの公式ドキュメントを参照してください。Condaを使用すると、多くの場合、依存関係を含めて簡単にインストールできます。

基本的な概念とワークフロー

Pyomoを用いた最適化モデリングの基本的な流れは以下のようになります。

モデルオブジェクトの作成: `ConcreteModel` または `AbstractModel` のインスタンスを作成します。
コンポーネントの宣言: モデルに必要な要素（集合、パラメータ、変数、目的関数、制約条件）をモデルオブジェクトの属性として定義します。
モデルのインスタンス化 (AbstractModelの場合): 抽象モデルに具体的なデータを紐付けて、解くことができる具体的な問題インスタンスを生成します。
ソルバーの選択と設定: `SolverFactory` を使用して、利用するソルバーを指定します。
最適化の実行: ソルバーの `solve()` メソッドを呼び出して、モデルを解きます。
結果の確認と分析: ソルバーの終了ステータスを確認し、最適解（目的関数の値や変数の値）を取得して分析します。

ConcreteModel vs AbstractModel

Pyomoでは、主に2つのモデルタイプが提供されています。

ConcreteModel (具象モデル): モデル構造とデータを同時に定義します。小〜中規模の問題や、モデル構造がデータに依存しない場合にシンプルで扱いやすいです。Pythonスクリプト内で直接データを定義します。
AbstractModel (抽象モデル): まずモデルの構造（数式）だけを定義し、後からデータファイル（`.dat`形式など）を読み込んで具体的な問題インスタンスを生成します。大規模な問題や、同じモデル構造で異なるデータセットを使って繰り返し解きたい場合に適しています。モデル構造とデータの分離が可能です。

初心者はまず `ConcreteModel` から始めるのが理解しやすいでしょう。この解説でも主に `ConcreteModel` を用います。

主要なモデリングコンポーネント

Pyomoモデルは、以下の主要なコンポーネント（Pythonクラス）を組み合わせて構築されます。

コンポーネント	クラス名	説明
集合 (Sets)	`Set`	モデルのインデックス（添え字）の集合を定義します。例えば、製品の種類、工場の場所、時間ステップなど。多次元の集合も定義可能です。`RangeSet`は連続した整数の集合を簡単に定義できます。
パラメータ (Parameters)	`Param`	モデルの入力データとなる定数を定義します。例えば、製品の単価、輸送コスト、需要量、リソースの上限など。集合でインデックス付けされたパラメータも定義できます。
変数 (Variables)	`Var`	最適化によって値を決定する決定変数を定義します。定義域（連続、整数、バイナリ）や範囲（上限・下限）を指定できます。集合でインデックス付けされた変数群も定義できます。
目的関数 (Objective)	`Objective`	最大化または最小化したい数式を定義します。通常はモデルに一つだけ定義します。
制約条件 (Constraints)	`Constraint`	決定変数が満たすべき等式または不等式を定義します。集合でインデックス付けされた複数の制約式をまとめて定義することも可能です。
式 (Expressions)	`Expression`	複雑な数式を部品化し、目的関数や制約条件で再利用するために定義します。必須ではありませんが、モデルの可読性や構造化に役立ちます。

これらのコンポーネントを組み合わせることで、数理最適化モデルをPythonコードとして表現します。

モデリングコンポーネント詳解

ここでは、主要なコンポーネントの使い方をもう少し詳しく見ていきます。

Set: 集合の定義

Setコンポーネントは、モデルのインデックスを定義します。リスト、タプル、Pythonのsetオブジェクト、またはジェネレータ関数で初期化できます。

from pyomo.environ import ConcreteModel, Set, RangeSet

model = ConcreteModel()

# リストで初期化
model.Products = Set(initialize=['A', 'B', 'C'])

# タプルで初期化
model.Locations = Set(initialize=('Tokyo', 'Osaka', 'Nagoya'))

# 連続する整数の集合 (1から10まで)
model.TimePeriods = RangeSet(1, 10)

# 多次元集合 (直積集合)
model.Routes = model.Locations * model.Locations

# 集合の集合も可能
model.ProductGroups = Set(initialize=[('A', 'B'), ('C',)])

# 集合演算も可能
model.SetA = Set(initialize=[1, 2, 3])
model.SetB = Set(initialize=[3, 4, 5])
model.UnionSet = model.SetA | model.SetB    # 和集合 {1, 2, 3, 4, 5}
model.IntersectionSet = model.SetA & model.SetB # 積集合 {3}
model.DifferenceSet = model.SetA - model.SetB # 差集合 {1, 2}
model.SymmetricDifferenceSet = model.SetA ^ model.SetB # 対称差 {1, 2, 4, 5}

initializeの他に、filter（要素をフィルタリングする関数）やvalidate（要素を検証する関数）などのオプションも指定できます。

Param: パラメータの定義

Paramコンポーネントは、モデルの定数データを定義します。スカラー値または集合でインデックス付けされた値を持ちます。辞書や関数で初期化できます。

from pyomo.environ import ConcreteModel, Set, Param

model = ConcreteModel()
model.Products = Set(initialize=['A', 'B'])
model.Locations = Set(initialize=['Tokyo', 'Osaka'])

# スカラーパラメータ
model.Budget = Param(initialize=10000)

# インデックス付きパラメータ (辞書で初期化)
demand_data = {('A', 'Tokyo'): 100, ('A', 'Osaka'): 150, ('B', 'Tokyo'): 200, ('B', 'Osaka'): 80}
model.Demand = Param(model.Products, model.Locations, initialize=demand_data, default=0) # default値も指定可能

# インデックス付きパラメータ (関数で初期化)
def cost_rule(model, prod, loc):
    # 何らかの計算ロジック (例)
    costs = {('A', 'Tokyo'): 10, ('A', 'Osaka'): 12, ('B', 'Tokyo'): 9, ('B', 'Osaka'): 11}
    return costs.get((prod, loc), 9999) # 見つからない場合のデフォルト値

model.UnitCost = Param(model.Products, model.Locations, initialize=cost_rule)

# パラメータ値へのアクセス
print(model.Budget.value)
print(model.Demand['A', 'Tokyo'].value)
print(model.UnitCost['B', 'Osaka'].value)

mutable=True とすると、モデル作成後や求解後にパラメータ値を変更できます。

Var: 変数の定義

Varコンポーネントは、決定変数を定義します。定義域(domain)、範囲(bounds)、初期値(initialize)などを指定できます。

from pyomo.environ import ConcreteModel, Set, Var, NonNegativeReals, Binary, Integers

model = ConcreteModel()
model.Products = Set(initialize=['A', 'B'])
model.Locations = Set(initialize=['Tokyo', 'Osaka'])

# スカラー変数 (非負実数)
model.TotalCost = Var(domain=NonNegativeReals)

# インデックス付き変数 (非負実数)
# domain: 変数の種類 (デフォルトはReals: 実数)
#    NonNegativeReals: 非負実数 (0以上)
#    PositiveReals: 正の実数 (0より大きい)
#    NonPositiveReals: 非正実数 (0以下)
#    NegativeReals: 負の実数 (0未満)
#    Integers: 整数
#    NonNegativeIntegers: 非負整数 (0以上)
#    PositiveIntegers: 正の整数 (1以上)
#    Binary: 0または1
model.ProduceAmount = Var(model.Products, domain=NonNegativeReals)

# インデックス付き変数 (バイナリ変数) と範囲指定
# bounds: 変数の範囲をタプル (下限, 上限) で指定
model.UseFactory = Var(model.Locations, domain=Binary)
model.TransportAmount = Var(model.Products, model.Locations, domain=NonNegativeIntegers, bounds=(0, 1000))

# 変数値へのアクセス (求解後)
# print(model.TotalCost.value)
# print(model.ProduceAmount['A'].value)
# for loc in model.Locations:
#    print(f"UseFactory[{loc}] = {model.UseFactory[loc].value}")

Objective: 目的関数の定義

Objectiveコンポーネントは、最適化の目標を定義します。sense引数で最大化(maximize)か最小化(minimize)を指定し、expr引数で目的関数となる数式を指定します。

from pyomo.environ import ConcreteModel, Objective, Var, Param, Set, maximize, minimize, summation

model = ConcreteModel()
# ... (Set, Param, Var の定義は省略) ...
model.Products = Set(initialize=['A', 'B'])
model.Locations = Set(initialize=['Tokyo', 'Osaka'])
model.Price = Param(model.Products, initialize={'A': 100, 'B': 120})
model.Cost = Param(model.Products, initialize={'A': 50, 'B': 60})
model.ProduceAmount = Var(model.Products, domain=NonNegativeReals)

# 目的関数: 総利益の最大化
def total_profit_rule(model):
    # Pyomo の summation を使うと簡潔に書ける
    return summation(model.Price, model.ProduceAmount) - summation(model.Cost, model.ProduceAmount)
    # 通常のPythonループでも書ける
    # profit = 0
    # for p in model.Products:
    #     profit += (model.Price[p] - model.Cost[p]) * model.ProduceAmount[p]
    # return profit

model.ProfitObjective = Objective(rule=total_profit_rule, sense=maximize)

# 目的関数: 総コストの最小化 (例)
# def total_cost_rule(model):
#    return summation(model.Cost, model.ProduceAmount)
# model.CostObjective = Objective(rule=total_cost_rule, sense=minimize)

目的関数の定義には、Pythonの関数（ルール関数）を使うのが一般的です。これにより、複雑な計算も記述できます。

Constraint: 制約条件の定義

Constraintコンポーネントは、モデルの制約式を定義します。等式(==)、不等式(<=, >=)を使って表現します。

from pyomo.environ import ConcreteModel, Constraint, Var, Param, Set, NonNegativeReals

model = ConcreteModel()
model.Products = Set(initialize=['A', 'B'])
model.Resources = Set(initialize=['Material', 'Labor'])
model.ProduceAmount = Var(model.Products, domain=NonNegativeReals)
model.ResourceUsage = Param(model.Products, model.Resources, initialize={
    ('A', 'Material'): 2, ('A', 'Labor'): 3,
    ('B', 'Material'): 4, ('B', 'Labor'): 1
})
model.ResourceLimit = Param(model.Resources, initialize={'Material': 100, 'Labor': 120})

# 単一の制約
# model.SomeConstraint = Constraint(expr=model.ProduceAmount['A'] <= 50)

# インデックス付き制約 (リソース制約)
def resource_limit_rule(model, resource):
    # 各リソースについて、全製品の生産に必要なリソース量の合計が上限を超えない
    usage = summation(model.ResourceUsage, model.ProduceAmount, index=resource)
    # またはPythonループで:
    # usage = sum(model.ResourceUsage[p, resource] * model.ProduceAmount[p] for p in model.Products)
    return usage <= model.ResourceLimit[resource]

model.ResourceConstraint = Constraint(model.Resources, rule=resource_limit_rule)

# 等式制約の例 (需要を満たす制約)
# model.Demand = Param(model.Products, ...)
# def demand_satisfaction_rule(model, product):
#    return model.ProduceAmount[product] == model.Demand[product]
# model.DemandConstraint = Constraint(model.Products, rule=demand_satisfaction_rule)

# 制約式の確認 (例)
# model.ResourceConstraint.pprint()

目的関数と同様に、ルール関数を使って制約式を定義するのが一般的です。インデックス付き制約では、ルール関数は対応するインデックス（上の例ではresource）を引数に取ります。

Expression: 式の定義

Expressionコンポーネントは、再利用可能な数式を定義します。

from pyomo.environ import ConcreteModel, Expression, Var, Param, Set

model = ConcreteModel()
model.Products = Set(initialize=['A', 'B'])
model.Locations = Set(initialize=['Tokyo', 'Osaka'])
model.TransportAmount = Var(model.Products, model.Locations, domain=NonNegativeReals)
model.TransportCost = Param(model.Products, model.Locations, initialize={...}) # コストデータ

# 各製品ごとの総輸送コストを計算する式
def total_transport_cost_per_product_rule(model, product):
    return sum(model.TransportCost[product, loc] * model.TransportAmount[product, loc] for loc in model.Locations)

model.TotalTransportCostPerProduct = Expression(model.Products, rule=total_transport_cost_per_product_rule)

# 目的関数などでこの Expression を利用できる
# def overall_objective_rule(model):
#    ...
#    total_cost = sum(model.TotalTransportCostPerProduct[p] for p in model.Products)
#    ...
#    return total_cost
# model.MyObjective = Objective(rule=overall_objective_rule, sense=minimize)

簡単な例題: 線形計画問題（生産計画）

具体的な例として、簡単な生産計画問題をPyomoでモデル化し、解いてみましょう。

問題設定

ある工場で2種類の製品（製品A、製品B）を生産しています。生産には2種類の資源（資源1、資源2）が必要です。目的は、資源の制約内で総利益を最大化する各製品の生産量を決定することです。

製品情報:
- 製品A: 利益 5 (単位あたり), 資源1を 2単位使用, 資源2を 3単位使用
- 製品B: 利益 4 (単位あたり), 資源1を 1単位使用, 資源2を 1単位使用
資源制約:
- 資源1: 利用可能量 8 単位
- 資源2: 利用可能量 9 単位
決定変数:
- x_A: 製品Aの生産量 (非負実数)
- x_B: 製品Bの生産量 (非負実数)
目的関数: 総利益 = 5 * x_A + 4 * x_B （最大化）
制約条件:
- 資源1制約: 2 * x_A + 1 * x_B ≤ 8
- 資源2制約: 3 * x_A + 1 * x_B ≤ 9
- 非負制約: x_A ≥ 0, x_B ≥ 0

Pyomoコード (ConcreteModel)

import pyomo.environ as pyo

# 1. モデルオブジェクトの作成
model = pyo.ConcreteModel(name="Simple Production Plan")

# 2. コンポーネントの宣言

# 集合 (製品)
model.PRODUCTS = pyo.Set(initialize=['A', 'B'])

# パラメータ (利益、資源使用量、資源上限)
profit_data = {'A': 5, 'B': 4}
model.Profit = pyo.Param(model.PRODUCTS, initialize=profit_data)

resource1_usage_data = {'A': 2, 'B': 1}
model.Res1Usage = pyo.Param(model.PRODUCTS, initialize=resource1_usage_data)

resource2_usage_data = {'A': 3, 'B': 1}
model.Res2Usage = pyo.Param(model.PRODUCTS, initialize=resource2_usage_data)

model.Res1Limit = pyo.Param(initialize=8)
model.Res2Limit = pyo.Param(initialize=9)

# 変数 (生産量) - 非負実数
model.ProduceQty = pyo.Var(model.PRODUCTS, domain=pyo.NonNegativeReals)

# 目的関数 (総利益の最大化)
def total_profit_rule(model):
    return sum(model.Profit[p] * model.ProduceQty[p] for p in model.PRODUCTS)

model.Objective = pyo.Objective(rule=total_profit_rule, sense=pyo.maximize)

# 制約条件
# 資源1制約
def resource1_constraint_rule(model):
    return sum(model.Res1Usage[p] * model.ProduceQty[p] for p in model.PRODUCTS) <= model.Res1Limit

model.ConstraintRes1 = pyo.Constraint(rule=resource1_constraint_rule)

# 資源2制約
def resource2_constraint_rule(model):
    return sum(model.Res2Usage[p] * model.ProduceQty[p] for p in model.PRODUCTS) <= model.Res2Limit

model.ConstraintRes2 = pyo.Constraint(rule=resource2_constraint_rule)

# モデル構造の表示 (任意)
# model.pprint()

# 3. ソルバーの選択と設定
# SolverFactory を使ってソルバーを指定 (ここでは GLPK を使う例)
# 事前に GLPK がインストールされ、PATHが通っている必要がある
try:
    solver = pyo.SolverFactory('glpk')
    solver_available = True
except:
    print("GLPKソルバーが見つかりません。インストールされているか、PATHを確認してください。")
    solver_available = False

# 4. 最適化の実行 (ソルバーが見つかった場合のみ)
if solver_available:
    results = solver.solve(model, tee=True) # tee=True でソルバーのログを表示

    # 5. 結果の確認と分析
    if results.solver.termination_condition == pyo.TerminationCondition.optimal:
        print("\n--- 最適解が見つかりました --- ")
        print(f"最大総利益: {pyo.value(model.Objective):.2f}")
        print("各製品の生産量:")
        for p in model.PRODUCTS:
            print(f"  {p}: {model.ProduceQty[p].value:.2f}")
        # 使用した資源量も確認できる (例)
        print("\n使用資源量:")
        print(f"  資源1: {sum(model.Res1Usage[p].value * model.ProduceQty[p].value for p in model.PRODUCTS):.2f} / {model.Res1Limit.value}")
        print(f"  資源2: {sum(model.Res2Usage[p].value * model.ProduceQty[p].value for p in model.PRODUCTS):.2f} / {model.Res2Limit.value}")
    elif results.solver.termination_condition == pyo.TerminationCondition.infeasible:
        print("\n--- 実行可能解が見つかりませんでした --- ")
    else:
        print(f"\n--- ソルバー終了ステータス: {results.solver.termination_condition} ---")
else:
    print("\nソルバーが利用できないため、最適化を実行できませんでした。")

このコードを実行すると、GLPKソルバーが呼び出され、最適化計算が行われます。計算が成功すれば、最大の総利益とそのときの各製品の生産量が出力されます。

上記の例では、`solver = pyo.SolverFactory(‘glpk’)` の部分でGLPKを指定しています。もしCBCを使いたい場合は `’cbc’` に、IPOPTなら `’ipopt’` に、Gurobiなら `’gurobi’` に変更します（それぞれのソルバーがインストールされ、Pyomoからアクセス可能である必要があります）。

ソルバーとの連携

Pyomo自体はモデリングツールであり、実際に最適化問題を解くのはバックエンドの「ソルバー」の役割です。Pyomoは様々なソルバーとのインターフェースを提供します。

SolverFactory

ソルバーインスタンスを作成する最も一般的な方法は `pyomo.environ.SolverFactory` を使うことです。引数にソルバーの名前（通常は実行ファイル名またはPyomoが認識する名前）を指定します。

import pyomo.environ as pyo

# GLPK ソルバーのインスタンスを作成
solver_glpk = pyo.SolverFactory('glpk')

# CBC ソルバーのインスタンスを作成
solver_cbc = pyo.SolverFactory('cbc')

# IPOPT ソルバーのインスタンスを作成
solver_ipopt = pyo.SolverFactory('ipopt')

# Gurobi ソルバーのインスタンスを作成 (商用、ライセンス必要)
# solver_gurobi = pyo.SolverFactory('gurobi')

# CPLEX ソルバーのインスタンスを作成 (商用、ライセンス必要)
# solver_cplex = pyo.SolverFactory('cplex')

# 実行ファイルのパスを直接指定することも可能 (Windowsの例)
# solver_glpk_path = pyo.SolverFactory('glpk', executable='C:/path/to/glpk/w64/glpsol.exe')

`SolverFactory` は、指定されたソルバーがシステムで見つからない場合、例外を発生させます。`try-except` ブロックで囲むことで、ソルバーが見つからない場合の処理を記述できます。

ソルバーオプションの設定

ソルバーによっては、計算時間の上限、許容誤差、ログの詳細度など、様々なオプションを設定できます。`solve()` メソッドの `options` 引数に辞書形式で渡すか、ソルバーインスタンスの `options` 属性に設定します。

# solve() メソッドでオプションを指定
results = solver_glpk.solve(model, tee=True, options={'tmlim': 60}) # tmlimはGLPKの時間制限オプション (60秒)

# ソルバーインスタンスにオプションを設定
solver_gurobi = pyo.SolverFactory('gurobi')
solver_gurobi.options['TimeLimit'] = 120 # Gurobiの時間制限オプション (120秒)
solver_gurobi.options['MIPGap'] = 0.01   # GurobiのMIPギャップ許容値 (1%)
results_gurobi = solver_gurobi.solve(model, tee=True)

利用可能なオプションはソルバーによって異なります。各ソルバーのマニュアルを参照してください。

ソルバーの利用可能性の確認

特定のソルバーが利用可能かどうかを事前に確認するには、`SolverFactory` を `try-except` で囲む方法の他に、`available()` メソッドを使うこともできます。

solver_name = 'gurobi'
if pyo.SolverFactory(solver_name).available():
    print(f"{solver_name} は利用可能です。")
    solver = pyo.SolverFactory(solver_name)
    # ... 最適化を実行 ...
else:
    print(f"{solver_name} は利用できません。")

NEOS Server

ローカルマシンにソルバーをインストールできない場合でも、Pyomoは NEOS Server を介してリモートで問題を解く機能を提供しています。`SolverManagerFactory` を使って `neos` を指定します。

# NEOSサーバー経由で問題を解く (インターネット接続が必要)
solver_manager = pyo.SolverManagerFactory('neos')
results_neos = solver_manager.solve(model, opt='cbc', tee=True) # NEOS上でCBCソルバーを使用

NEOS Serverでは多くのオープンソースおよび商用ソルバーを利用できますが、計算時間や問題サイズに制限がある場合があります。

求解と結果の分析

モデルを定義し、ソルバーを選択したら、いよいよ問題を解き、その結果を分析します。

`solve()` メソッドの実行

`SolverFactory` で作成したソルバーインスタンスの `solve()` メソッドに、作成したモデルオブジェクトを渡して実行します。

solver = pyo.SolverFactory('glpk')
results = solver.solve(model, tee=True) # tee=True でソルバーのログを標準出力に表示

`solve()` メソッドは `SolverResults` オブジェクトを返します。このオブジェクトには、求解プロセスと結果に関する情報が含まれています。

ソルバーのステータス確認

最適化計算が成功したかどうかを確認することが重要です。`SolverResults` オブジェクトの `solver.status` 属性と `solver.termination_condition` 属性を確認します。

print(f"Solver Status: {results.solver.status}")
print(f"Termination Condition: {results.solver.termination_condition}")

# 最適解が見つかったかどうかのチェック
if results.solver.termination_condition == pyo.TerminationCondition.optimal:
    print("最適解が見つかりました！")
elif results.solver.termination_condition == pyo.TerminationCondition.infeasible:
    print("実行不可能 (Feasible solution does not exist).")
elif results.solver.termination_condition == pyo.TerminationCondition.unbounded:
    print("非有界 (Objective function can be improved indefinitely).")
else:
    print("その他の終了条件:", results.solver.termination_condition)

主な `TerminationCondition` の値:

optimal: 最適解が見つかりました。
infeasible: 制約を満たす解が存在しません（実行不可能）。
unbounded: 目的関数を無限に改善できます（非有界）。モデルの定義に誤りがある可能性があります。
maxTimeLimit: 設定された計算時間制限に達しました。
maxIterations: 設定された反復回数制限に達しました。
other: その他の理由で終了しました。ソルバーログを確認する必要があります。

結果の取得

最適解が見つかった場合、目的関数の値や決定変数の値を取得できます。

目的関数の値: `pyo.value(model.Objective)` または `model.Objective()` で取得できます。
変数の値: `pyo.value(model.VarName[index])` または `model.VarName[index].value` で取得できます。

if results.solver.termination_condition == pyo.TerminationCondition.optimal:
    # 目的関数の最適値
    optimal_objective_value = pyo.value(model.Objective)
    print(f"最適目的値: {optimal_objective_value}")

    # 各変数の最適値
    print("最適変数:")
    # インデックス付き変数の場合
    for p in model.PRODUCTS:
        print(f"  ProduceQty[{p}] = {pyo.value(model.ProduceQty[p])}")
    # スカラー変数の場合
    # print(f"  ScalarVar = {pyo.value(model.ScalarVar)}")

    # display() メソッドで結果をまとめて表示することも可能
    print("\n--- 結果サマリ (display) ---")
    model.display()

    # pprint() でより詳細な情報を表示
    # print("\n--- 結果詳細 (pprint) ---")
    # results.write() # 結果をファイルに書き出す
else:
    print("最適解が見つからなかったため、結果を取得できません。")

model.display() は、求解後のモデルの状態（変数値、目的関数値など）を要約して表示する便利なメソッドです。

感度分析情報（双対変数、被約費用）

線形計画問題(LP)の場合、ソルバーは感度分析情報として双対変数（シャドウプライス）や被約費用を計算することがあります。これらは制約条件や変数の変化が目的関数に与える影響を示します。

これらの情報を取得するには、モデルに `Suffix` コンポーネントを追加し、`solve()` を呼び出す際に `load_solutions=False` とし、結果オブジェクトから明示的に解をロードする必要があります（または、ソルバーとインターフェースによっては自動的にロードされる場合もあります）。

# 感度分析情報を取得する場合 (LPのみ)
model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)  # 双対変数用Suffix
model.rc = pyo.Suffix(direction=pyo.Suffix.IMPORT)    # 被約費用用Suffix

# ... solve() を実行 ...
results = solver.solve(model, tee=True)

if results.solver.termination_condition == pyo.TerminationCondition.optimal:
    print("\n--- 感度分析情報 ---")
    print("双対変数 (Shadow Prices):")
    for c in model.component_objects(pyo.Constraint, active=True):
        print(f"  Constraint '{c.name}':")
        for index in c:
             # 制約 c[index] に対応する双対変数を取得
             print(f"    {index}: {model.dual[c[index]]}")

    # print("\n被約費用 (Reduced Costs):")
    # for v in model.component_objects(pyo.Var, active=True):
    #     print(f"  Variable '{v.name}':")
    #     for index in v:
    #         # 変数 v[index] に対応する被約費用を取得
    #         print(f"    {index}: {model.rc[v[index]]}")
else:
    print("最適解が見つからないため感度分析情報は利用できません。")

注意: 感度分析情報の取得方法は、使用するソルバーやPyomoのバージョンによって若干異なる場合があります。詳細はPyomoのドキュメントや各ソルバーのマニュアルを確認してください。

高度な機能（一部紹介）

Pyomoは基本的な最適化モデリング以外にも、多くの高度な機能を提供しています。

AbstractModelとデータ管理: 前述の通り、`AbstractModel` を使うとモデル構造とデータを分離できます。`.dat` というテキストファイル形式でパラメータ値を定義し、`model.create_instance(‘data.dat’)` のようにしてモデルインスタンスを生成します。これにより、大規模なデータセットの管理や、異なるデータでのシミュレーションが容易になります。
非線形計画 (NLP): Pyomoは非線形な目的関数や制約条件を持つモデルも記述できます。`sin`, `cos`, `exp`, `log` などの数学関数を式の中で使用できます。NLPを解くにはIPOPTなどの非線形ソルバーが必要です。
```
from pyomo.environ import cos, sin, exp

model.x = Var()
model.y = Var()
model.Objective = Objective(expr=exp(model.x) * sin(model.y), sense=maximize)
model.Constraint1 = Constraint(expr=model.x**2 + model.y**2 == 1)
```
混合整数非線形計画 (MINLP): 非線形な要素と整数変数（またはバイナリ変数）が混在する問題です。解くのが非常に難しい問題クラスですが、PyomoはBonmin, Couenne, SCIPなどのMINLPソルバーと連携して対応できます。

一般化離接計画 (GDP – Generalized Disjunctive Programming): 「AまたはBのどちらかの制約を満たす」といった論理的な条件を含む問題をモデル化するためのフレームワークです (`pyomo.gdp`)。例えば、特定の設備を使うか使わないかによって異なる制約が適用される場合などに用います。GDPモデルは通常、内部でMINLPに変換されて解かれます。

from pyomo.gdp import Disjunct, Disjunction

model.Y = BooleanVar() # 離接を制御するブーリアン変数

# YがTrueのときに有効になる制約群
model.DisjunctTrue = Disjunct()
model.DisjunctTrue.Constraint1 = Constraint(expr=...)
model.DisjunctTrue.Constraint2 = Constraint(expr=...)

# YがFalseのときに有効になる制約群
model.DisjunctFalse = Disjunct()
model.DisjunctFalse.Constraint3 = Constraint(expr=...)

# 離接の定義: YがTrueならDisjunctTrueが、FalseならDisjunctFalseが有効
model.MyDisjunction = Disjunction(expr=[model.DisjunctTrue, model.DisjunctFalse])

微分代数方程式 (DAE): `pyomo.dae` パッケージは、時間微分や空間微分を含む最適化問題（動的最適化など）をモデル化するための拡張機能を提供します。連続的な時間領域や空間領域を `ContinuousSet` で定義し、`DerivativeVar` で微分変数を扱います。離散化手法（有限差分法、直交選点法など）を適用して、大規模なNLP問題に変換して解きます。
確率計画法 (Stochastic Programming): パラメータに不確実性が含まれる場合の最適化問題を扱うための機能 (`PySP: Pyomo Stochastic Programming`) も提供されています。複数のシナリオを考慮して、期待値を最適化したり、リスクを管理したりするモデルを構築できます。
二段階計画問題 (Bilevel Programming): 最適化問題の中に別の最適化問題が入れ子になっているような階層構造を持つ問題をモデル化するための機能 (`pyomo.bilevel`) も開発されています。
ソルバーインターフェース (Persistent Solvers): 大規模な問題や、モデルを少しずつ変更しながら繰り返し解く場合に、ソルバーとの通信オーバーヘッドを削減するためのPersistent Solverインターフェースも提供されています。モデル情報をソルバープロセス内に保持し、効率的な更新と再求解を可能にします。

これらの高度な機能により、Pyomoは非常に広範な最適化問題に対応できる強力なプラットフォームとなっています。

応用分野と事例

Pyomoとその背後にある数理最適化技術は、様々な産業や研究分野で活用されています。

製造業: 生産計画、スケジューリング、在庫管理、設備投資計画、品質管理
物流・サプライチェーン: 輸送ルート最適化、倉庫配置、配送計画、サプライチェーンネットワーク設計
エネルギー: 発電計画（ユニットコミットメント）、送電網運用、再生可能エネルギー導入計画、スマートグリッド管理
金融: ポートフォリオ最適化、リスク管理、オプション価格付け、アルゴリズム取引
化学工学: プロセス設計・最適化、反応ネットワーク解析、スケジューリング
航空宇宙: 軌道設計、フライトスケジューリング、クルー割り当て
通信: ネットワーク設計、周波数割り当て、ルーティング
公共サービス: 公共交通機関のスケジュール最適化、施設配置、資源配分
機械学習: ハイパーパラメータ最適化、特徴選択、サポートベクターマシンの学習（二次計画問題として定式化される）

Pyomoは、これらの分野における複雑な意思決定問題を、数学モデルとして定式化し、最適な解（またはそれに近い解）を見つけるための強力なツールとなります。オープンソースであるため、学術研究だけでなく、企業のデジタルトランスフォーメーション(DX)推進においても活用が進んでいます。

まとめ

Pyomoは、Pythonの柔軟性と数理最適化モデリングの強力な表現力を組み合わせた、非常に有用なオープンソースライブラリです。直感的なPython構文で多様な最適化問題をモデル化し、様々なオープンソースおよび商用ソルバーと連携して解を求めることができます。

基本的な線形計画問題から、複雑な混合整数非線形計画問題、微分代数方程式を含む動的最適化まで、幅広い問題に対応できる拡張性を持っています。活発なコミュニティと継続的な開発により、今後も進化が期待されるツールです。

この解説が、Pyomoを使った数理最適化モデリングの世界への第一歩となれば幸いです。ぜひ、ご自身の課題解決にPyomoを活用してみてください！