実装：cvxpyで凸問題を解く

🎓 レベル：基礎　|　重要度：B（標準）

📎 前提：凸集合と凸関数・二次計画（QP）　|　関連：実装ツールの使い分け

要点（BLUF）

cvxpy は DCP（規律ある凸計画） で問題を書く。凸性が崩れる式はエラーで弾いてくれる。
変数・目的・制約を宣言するだけで、LP・QP・SOCP・SDP を 同じ書き方で解ける。
双対変数（影の価格、凸最適化問題と双対理論）も constraint.dual_value で取り出せる。

概念 ── モデリング言語の発想

scipy の linprog/minimize は「行列で問題を渡す」低レベル API。cvxpy は 数式に近い形で書け、内部で適切なソルバー（LP→シンプレックス/内点法、SOCP/SDP→内点法）へ自動変換する。人間は「何を最適化したいか」だけ書けばよい。

DCP（Disciplined Convex Programming）

cvxpy は、式が 凸性を保つ規則（凸集合と凸関数の演算規則）で組み立てられているかを検査する。たとえば：

cp.sum_squares(x)、cp.norm(x, 2)、cp.quad_form(x, Q)（ $Q\succeq0$ ）は凸 → 最小化 OK。
凸関数同士の非負和・最大値・アフィン合成は凸。
cp.log(x)、cp.sqrt(x) は凹 → 最大化なら OK、最小化はエラー。

DCP に従わない式（例：凸×凸、非凸な積）は その場でエラーになる。「解けたが実は非凸で局所解だった」という事故を未然に防ぐ ── これが DCP の安全性。

graph LR
  A["数式で記述 (Variable/Minimize/制約)"] --> B["DCP検査: 凸性が保たれるか"]
  B -->|"凸ならOK"| C["標準錐形式へ自動変換"]
  C --> D["ソルバー(内点法等)で求解"]
  B -->|"非凸ならエラー"| E["DCPError で警告"]

同じ書き方で LP・QP・SOCP

LP（生産計画、線形計画の標準形と幾何）・QP・SOCP を cvxpy の統一構文で解く。

import numpy as np
import cvxpy as cp

# --- LP: 生産計画 max 40xA+30xB ---
x = cp.Variable(2, nonneg=True)
lp = cp.Problem(cp.Maximize(40*x[0] + 30*x[1]),
                [2*x[0] + x[1] <= 120, x[0] + x[1] <= 100])
lp.solve()
print(f"[LP]   xA={x.value[0]:.1f}, xB={x.value[1]:.1f}, 利益={lp.value:.1f}")

# --- QP: min (x-1)^2+(y-2)^2 s.t. x+y<=2 ---
y = cp.Variable(2)
qp = cp.Problem(cp.Minimize(cp.sum_squares(y - np.array([1.0, 2.0]))),
                [y[0] + y[1] <= 2, y >= 0])
qp.solve()
print(f"[QP]   x={y.value[0]:.3f}, y={y.value[1]:.3f}, f={qp.value:.4f}")

# --- SOCP: min ||z||_2 s.t. sum(z)=3 ---
z = cp.Variable(3)
socp = cp.Problem(cp.Minimize(cp.norm(z, 2)), [cp.sum(z) == 3])
socp.solve()
print(f"[SOCP] z={np.round(z.value,3)}, ||z||={socp.value:.4f}")

実行結果：

[LP]   xA=20.0, xB=80.0, 利益=3200.0
[QP]   x=0.500, y=1.500, f=0.5000
[SOCP] z=[1. 1. 1.], ||z||=1.7321

3つの異なる錐の問題（錐計画（SOCP・SDP）入門）が、Variable・Minimize/Maximize・制約リストという 同じ3点セットで書けた。これまでの章で scipy/PuLP を使って別々に解いた問題（LP=3200、QP=0.5、最小ノルム=√3）が、すべて再現している。

双対変数と実務の作法

# 制約の双対変数(影の価格)を取り出す
print(lp.constraints[0].dual_value)   # 機械時間の影の価格

求解後に status を確認：optimal 以外（infeasible・unbounded）なら結果を信じない。
prob.value が目的の最適値、var.value が最適解。
ソルバー選択：prob.solve(solver=cp.ECOS / cp.SCS / cp.CLARABEL) で切替。問題規模・精度で使い分け。
パラメータ化（cp.Parameter）で、係数だけ変えた再解を高速化できる。

数式の直観的意味

cvxpy（DCP）が体現するのは「凸最適化は 合成可能（composable）」という事実。凸性を保つ演算（凸集合と凸関数）だけを部品として与えれば、それらをどう組み合わせても凸性が保たれる ── だからユーザーは個々の式の凸性を証明せずとも、規則に従うだけで凸問題を構築できる。内部では問題を標準錐形式（錐計画（SOCP・SDP）入門）に変換し、内点法（内点法入門）で解く。「モデリング（人間）」と「求解（ソルバー）」を分離するこの設計が、最適化を専門家以外にも開いた。低レベル API（scipy/PuLP）との使い分けは実装ツールの使い分け。

⚠️ よくある誤解・落とし穴

DCP エラー＝非凸の警告：x*y や cp.sqrt(cp.sum_squares(x)) の素朴な記述は弾かれる。凸な等価表現（cp.norm）に直す。
status を必ず確認：infeasible/unbounded で value が None/inf になる。
整数変数も書けるが凸でなくなる：cp.Variable(boolean=True) は MILP になり、別系統のソルバー（分枝限定法）が走る。
大規模 SDP・SOCP はソルバーと前処理で性能が大きく変わる。精度（eps）と速度のトレードオフを調整する。