%load_ext lab_black
準備
以下では，基本的には数理最適化ソルバーとしては Gurobi（パッケージはgurobipy）を用いるが， これは商用パッケージ（アカデミック・非商用は無料）である． Gurobiを使用できない場合には，mypulp をかわりに読み込んで使う（以下では，コメントアウトしてある）． mypulpはPuLP (https://coin-or.github.io/pulp/index.html) のラッパーであり，Gurobiの基本的なAPIをサポートする． PuLPの標準ソルバーはCBC (https://github.com/coin-or/Cbc) であり，2次（錘）最適化問題は解くことができない．
%matplotlib inline
import numpy as np
from gurobipy import Model, quicksum, GRB, multidict
# from mypulp import Model, quicksum, GRB, multidict
import intvalpy
最適化問題とは
集合 $\mathbf{F}$, $\mathbf{F}$ から実数への写像 $f: \mathbf{F} \rightarrow \mathbf{R}$ （$\mathbf{R}$ は実数全体の集合） が与えられたとき $$ \begin{array}{l l } minimize & f(x) \\ s.t. & x \in \mathbf{F} \end{array} $$
を与える $x \in \mathbf{F}$ を求める問題を最適化問題(optimization problem)とよぶ．ここで minimize と s.t. (subject to) は，それぞれ「最小化」と「制約条件」を表す記号である． $\mathbf{F}$ が有限集合のときには，"minimize"は要素の中の最小値を表す $\min$ を意味し， 無限集合の場合には下界の最大値を表す $\inf$ を意味する．$\inf$ は $\min$ を一般化した関数である． たとえば $\inf \{ (0,1] \}$ は $0$ である． また，$\inf \{ \emptyset \}$ は $\infty$ と定義する． $\mathbf{F}$ を実行可能解の集合，その要素 $x \in \mathbf{F}$ を実行可能解(feasible solution)または単に解(solution)とよぶ． 上式における関数 $f(x)$ を 目的関数(objective function)とよぶ． 目的関数の最小値を最適値(optimal value)とよぶ．最適値 $z^*$ は以下のように定義される． $$ z^* = \inf \{ f(x) | x \in \mathbf{F} \} $$ $z^* = \infty$ のとき，問題は実行不可能(infeasible)もしくは実行不能とよばれ， $z=-\infty$ のとき非有界(unbounded)とよばれる．
最適値 $z^*$ を与える解 $x \in \mathbf{F}$ を大域的最適解(globally optimal solution)とよび，その集合を $\mathbf{F}^*$ と記す． 大域的最適解の集合 $\mathbf{F}^*$ は以下のように定義される. $$ \mathbf{F}^* =\{ x \in \mathbf{F}| f(x) = z^* \} $$
ある実行可能解 $x \in \mathbf{F}$ の「近所」にある解の集合を $x$ の近傍(neighborhood)とよぶ． 近傍の定義は対象とする問題によって異なるが，一般には， 実行可能解の集合 $\mathbf{F}$ を与えたとき，近傍 $N$ は以下の写像と定義される． $$ N: \mathbf{F} \rightarrow 2^\mathbf{F} $$ すなわち，実行可能解の集合から，そのべき集合（部分集合の集合）への写像が近傍である．
実行可能解 $x \in \mathbf{F}$ で $$ f(x) = \inf \{ f(y)| y \in N(x) \} $$ を満たすものを（近傍 $N$ に対する）局所的最適解(locally optimal solution)とよぶ．
目的関数が実数値関数 $f: \mathbf{R}^n \rightarrow \mathbf{R}$ のときには，$\epsilon >0$ に対する解 $x$ の $\epsilon$ 近傍 $$ N(x,\epsilon) = \{ y \in \mathbf{R}^n | \|x-y\| \leq \epsilon \} $$ を用いる．ここで $\| x-y \|$ はEuclidノルム $\sqrt{\sum_{i=1}^n (x_i-y_i)^2}$ を表す．
問題によっては大域的最適解でなく，局所的最適解を求めることを目的とすることもある．
最適化問題の分類
以下では，最適化問題を幾つかの基準にしたがい分類する．実際には，これらの分類基準は正確なものではなく， 解法を考えていくときに便利だからということに過ぎないことを付記しておく．
連続か離散か
問題に内在する変数が連続値をとるか，離散値をとるかによって最適化問題を分類する． 連続値をとる最適化問題は，連続最適化(continuous optimization)，離散値をとる最適化問題は離散最適化(discrete optimization)とよばれる． 多くの場合，連続値は適当な精度で打ち切られた有理数として計算されるので，理論的には離散最適化によって連続最適化を解くことができる． 実際には連続最適化では，変数が連続値であることを利用した最適化手法を用いて効率的に解くことができるので， 連続最適化と離散最適化を分けて考えるのである．
離散最適化問題は，変数が数値か否かによってさらに2つに分類される．数値に限定した離散最適化問題を整数最適化(integer optimization)とよび， 数値でなく組合せ構造をもつ有限集合から解を選択する問題を組合せ最適化(combinatorial optimization)とよぶ， 整数を用いて組合せ構造を表すことができ，逆に整数値をとる解の組も組合せ構造として表すことができるので， この2つの問題は実質的に同値である．実際には整数最適化では，整数値を実数に緩和した連続最適化問題を基礎として解法を組み立て， 組合せ最適化では組合せ構造を利用して解法を組み立てるので，分けて考えるのである． なお，整数に限定された変数と連続最適化で取り扱う実数をとる変数の両者が含まれる問題を特に， 混合整数最適化(mixed integer optimization)とよぶ．多くの数理最適化ソルバーはこの混合整数最適化に特化したものである．
変数が数値で表されている問題を特に数理最適化(mathematical optimization) とよぶ．数理最適化には連続最適化と（混合）整数最適化が含まれる．
線形か非線形か
連続最適化は，大きく線形と非線形に分類される． 目的関数や制約条件がすべてアフィン関数（線形関数を平行移動して得られる関数）として記述されている問題を 線形最適化(linear optimization)， そうでない問題を非線形最適化(nonlinear optimization)とよぶ．
非線形最適化には様々なクラスがある． 目的関数だけが2次関数の2次最適化(quadratic optimization)， 制約にも2次関数が含まれる2次制約最適化(quadratic constrained optimization)の他にも，より一般的な 2次錐最適化(second-order cone optimization)，半正定値最適化(semi-definite optimization)， 多項式最適化(polynomial optimization)などがある．
変数が整数に限定された線形最適化問題を整数線形最適化(integer linear optimization)とよぶ． 任意の非線形関数は，整数変数を用いて区分的線形関数として近似できるので，理論的には（効率を度外視すれば） 非線形最適化は整数線形最適化で解くことができる．一方，$x(1-x)=0$ という形の多項式制約によって変数 $x$ が $0$ か $1$ かに限定でき， このような2値の整数変数を用いれば任意の整数変数を表現できるので，理論的には整数最適化は多項式最適化で解くことができる． 繰り返しになるが，原理的には同値である問題に対して分類が必要なのは，個々の問題に特化した解法が必要であるからである．
凸か非凸か
非線形最適化は，目的関数ならびに制約領域が凸(convex)であるか否かによって問題が分類できる．
関数 $f: \mathbf{R}^n \rightarrow \mathbf{R}$ は，すべての $x,y \in \mathbf{R}^n$ と $\lambda \in [0,1]$ に対して $$ f( \lambda x +(1-\lambda) y) \leq \lambda f(x) + (1-\lambda) f(y) $$ が成立するとき，凸関数(convex function)とよばれ， $$ f( \lambda x +(1-\lambda) y) \geq \lambda f(x) + (1-\lambda) f(y) $$ が成立するとき，凹関数(concave function)とよばれる．
また，集合 $S$ は，すべての $x,y \in S$ と $\lambda \in [0,1]$ に対して $$ \lambda x +(1-\lambda) y \in S $$ が成立するとき，凸集合(convex sex)とよばれる．
凸集合で表される制約領域内で凸関数を最小化（もしくは凹関数を最大化）する場合には， 大域的最適解と局所的最適解が一致するため，探索が容易になる．そのような問題を凸最適化(convex optimization)とよび， そうでない問題を非凸最適化(nonconvex optimization)とよぶ．
整数最適化や組合せ最適化も（特殊な場合を除いては）非凸最適化の範疇に含まれるが， 一般に大域的最適化という用語は，非線形な非凸最適化問題を対象として用いられる． 非凸最適化に対して大域的最適解を求める際には何らかの列挙法か近似解法に頼ることが多い．
ちなみに線形関数は凸でもあり凹でもあるので，線形最適化は凸最適化に分類される． 他にも半正定値最適化，2次錐最適化，凸2次制約最適化も凸最適化であるので，その性質を利用した効率的な解法が構築できる．
大域的最適解か局所的最適解か
ほとんどの非線形最適化の目的は，局所的最適解を求めることである．凸最適化問題に対しては局所的最適解が大域的最適解になるので，それで十分である． 非凸な問題の場合には，本来ならばすべての実行可能解の中で最も目的関数が良い大域的最適解を求めたいのであるが， それが困難であるので，大域的最適解の必要条件を満たす解（局所的最適解）を求める訳である．
非線形最適化問題に対する主流の解法は，初期解を改善していく反復法であるので， 初期解を変えて何度も局所的最適解を求めることによって大域的最適解に近い解を求めることができる． 一方，厳密な意味での大域的最適解を求める問題を大域的最適化(global optimization)とよぶ． 大域的最適化は本質的に難しい問題となるが，問題の構造を利用した厳密解法の研究が進められている．
不確実性の有無
問題に内在するパラメータに不確実性が含まれている問題を確率最適化(stochastic optimization)とよぶ． また，不確実性をパラメータの範囲に変換して，範囲内での最悪のパラメータに対する最適解を求める問題を ロバスト最適化(robust optimization)とよぶ． 将来発生する事象に対して何の情報ももたない問題をオンライン最適化(online optimization)とよぶ． 一方，すべてのパラメータが確定値の問題を，特に確定最適化(deterministic optimization)とよんで区別する場合もある．
ネットワーク構造をもつか否か
ネットワーク構造をもつ最適化問題に対しては，その構造を生かした解法がより高速に設計できる場合がある． そのような最適化問題をネットワーク最適化(network optimization)とよぶ．
ネットワーク上の幾つかの最適化問題は，線形最適化に帰着できる． しかし，多くの場合，ネットワークの構造を利用した，より高速なアルゴリズムを用いることができる．
線形最適化問題
線形最適化問題の一般形は以下のように書ける．
$$ \begin{array}{ll} minimize & c^T x \\ s.t. & A x \leq b \\ & x \in \mathbf{R}^n \end{array} $$ここで，$x$ は変数を表す $n$次元実数ベクトル， $A$ は制約の左辺係数を表す $m \times n$ 行列，$c$ は目的関数の係数を表す $n$次元ベクトル， $b$ は制約式の右辺定数を表す $m$次元ベクトルである．
簡単な例題と実行可能領域（実行可能解の範囲）を示す．
$$ \begin{array}{l c c c c} maximize & 15 x_1 & + 18 x_2 & & \\ s.t. & 2x_1 & + x_2 & \leq & 60 \\ & x_1 & + 2 x_2 &\leq & 60 \\ & & x_1,x_2& \geq & 0 \end{array} $$ここで，"maximize" は最大化（負号をつけて最小化）を表す．
model = Model("lo1")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
model.update() # Gurobiの怠惰な更新(lazy update)という仕様（忘れずに！）
model.addConstr(2 * x1 + x2 <= 60)
model.addConstr(x1 + 2 * x2 <= 60)
model.setObjective(15 * x1 + 18 * x2, GRB.MAXIMIZE)
model.optimize()
if model.Status == GRB.Status.OPTIMAL:
print("Opt. Value=", model.ObjVal)
for v in model.getVars():
print(v.VarName, v.X)
A = np.array([[-2, -1],
[-1, -2],
[1, 0],
[0, 1]])
b = np.array([-60, -60, 0, 0])
intvalpy.lineqs(A, b, title='Feasible Region', color='gray', alpha=0.5, s=10, size=(10,10), save=False, show=True);
錐最適化問題
錐最適化の枠組みを用いると，2次錐最適化や半正定値最適化などの様々な最適化問題を統一的に記述できる．
凸集合 $C \subseteq \mathbf{R}^n$ は，任意の $x \in C$ と任意の $\lambda >0$ に対して $\lambda x \in C$ となるとき，凸錐(convex cone)とよばれる． 凸錐 $C$ を用いると，錐線形最適化問題(cone linear optimization problem)は以下のように定義される．
$$ \begin{array}{ll} minimize & c^T x \\ s.t. & A x \leq b \\ & x \in C \end{array} $$ここで，$x$ は変数を表す $n$次元実数ベクトル， $A$ は制約の左辺係数を表す $m \times n$ 行列，$c$ は目的関数の係数を表す $n$次元ベクトル， $b$ は制約式の右辺定数を表す $m$次元ベクトルである．
$C= \mathbf{R}_+^n$ のとき線形最適化問題になり， $C$ が以下に定義される $n$次元空間の2次錐(second-order cone) $K_2^n$ のとき2次錐最適化問題になる．
$$ K_2^n = \{ x \in \mathbf{R}^{n} | \| (x_2,x_3,\ldots,x_n) \|_2 \leq x_1 \} $$ここで $$ \| (x_2,x_3,\ldots,x_n) \|_2 = \sqrt { x_2^2 + x_3^2 + \cdots + x_n^2 } $$ である．
しばしば，以下の回転つき2次錐(rotated second-order cone) $K_r^n$ で考えた方が，自然な定式化ができる．
$$ K_r^n = \left \{ x\in \mathbf{R}^{n} \mid x_3^2 + \cdots + x_n^2 \leq 2x_1 x_2 , \: x_1, x_2 \geq 0 \right \} $$以下の直交変換 $T_n$ によって $x \in K_2^n \Leftrightarrow T_n x \in K^n_r$ となるので，回転つき2次錐も2次錐に他ならない．
$$ T_n = \left[ \begin{array}{ccc} 1/\sqrt{2} & 1/\sqrt{2} & 0 \\ 1/\sqrt{2} & -1/\sqrt{2} & 0 \\ 0 & 0 & I_{n-2} \end{array} \right] $$また，$C$ が以下に定義される $3$次元空間の指数錐(exponential cone) $K_{exp}^3$ のとき指数錐最適化問題になる．
$$ K_{exp}^3 = \left\{(x,y,z) \in \mathbf{R}^{3} | y e^{x/y} \leq z, y >0 \right\} $$次に，半正定値最適化問題(semidefinite optimization problem)を定義する．
$p \times p$ の実対称行列の集合を以下のように定義する． $$ S^p = \{ X \in \mathbf{R}^{p \times p} | X=X^T \} $$
$p \times p$ の正方行列 $X$ は任意のベクトル $\lambda \in \mathbf{R}^p$ に対して $\lambda^T X \lambda \geq 0$ のとき半正定値(positive semidefinite)とよばれ，$X \succeq 0$ と書かれる．
以下に定義される対称半正定値行列の集合 $S_+^p$ は凸錐になる． $$ S_+^p= \{X \in S^p | X \succeq 0 \} $$
対称性から $p(p+1)/2$個の実数値を定めれば実対称行列が定まる． $n= p(p+1)/2$ としたとき，1つの行列は $n$次元実数ベクトル $x$ とみなすことができる． 錐 $C$ を半正定値行列の集合 $S_+^p$ としたとき，錐最適化問題は半正定値最適化問題とよばれる．
2次錐最適化問題はGurobiで解くことができるが，より一般的な錐に対してはSCS (https://www.cvxgrp.org/scs/ オープンソース)やMOSEK(https://www.mosek.com/ 商用)などのソルバーを用いる必要がある．
以下に，2次錐最適化問題の簡単な例題を示す．最初の制約式は２次錐制約 $\| (x,y) \|_2 \leq z$ を表している．
$$ \begin{array}{l c c c c c} maximize & 2x & +2y & + z & & \\ s.t. & x^2 & + y^2 & & \leq & z^2 \\ & 2x & + 3y & + 4z & \leq & 10 \\ & & & x,y,z & \geq & 0 \end{array} $$model = Model()
x = model.addVar(name="x")
y = model.addVar(name="y")
z = model.addVar(name="z")
model.update()
model.addConstr(x**2+y**2 <= z**2)
model.addConstr(2*x+3*y+4*z <= 10)
model.setObjective(2*x + 2*y + z, GRB.MAXIMIZE)
model.optimize()
if model.Status == GRB.Status.OPTIMAL:
print("Opt. Val.=", model.ObjVal)
print("(x,y,z)=", (x.X, y.X, z.X))
整数最適化問題
整数最適化問題は，以下のように書ける．
$$ \begin{array}{lll} minimize & f(x) & \\ s.t. & g_i(x) \leq 0 & i=1,2,\ldots,m \\ & x \in \mathbf{Z}^n \end{array} $$ここで $\mathbf{Z}$ は整数全体の集合，$x$ は変数を表す$n$次元ベクトル， $f: \mathbf{R}^n \rightarrow \mathbf{R}$ は目的関数，$g_i: \mathbf{R}^n \rightarrow \mathbf{R} (i=1,2,\ldots,m)$ は制約の左辺を表す関数である．
一般の整数最適化問題は解くことが困難であるが，目的関数と制約関数が線形である整数線形最適化問題ならびにその拡張に対しては，汎用のソルバーが準備されている． 一般の整数線形最適化問題は以下のように書ける．
$$ \begin{array}{ll} minimize & c^T x \\ s.t. & A x \leq b \\ & x \in \mathbf{Z}^n \end{array} $$ここで，$x$ は変数を表す $n$次元ベクトル，$A$ は制約の左辺係数を表す $m \times n$ 行列， $c$ は目的関数の係数を表す $n$次元ベクトル， $b$ は制約式の右辺定数を表す $m$次元ベクトルである．
変数が一般の整数でなく，$0$ もしくは $1$ の2値をとることが許される場合， 整数最適化問題は特に $0$-$1$整数最適化問題 ($0$-$1$ integer linear optimization problem)とよばれる， 整数変数に $0$以上 $1$以下という制約を付加すれば$0$-$1$整数最適化問題になり， さらに任意の整数は2進数で表現することができる（つまり$0$-$1$整数最適化問題を解ければ一般の整数最適化問題を解くことができる）ので， この2つは（少なくとも計算量の理論においては）同値である．
変数の一部が整数でなく実数であることを許した問題を混合整数最適化問題とよぶ． 多くの数理最適化ソルバーはこの問題を対象としている．また，モダンな数理最適化ソルバーは凸2次の目的関数や凸2次制約， さらには2次錐制約を取り扱うことができる．
簡単な例題を示す．
$$ \begin{array}{l c c c c c} minimize & & y & + z & & \\ s.t. & x & + y & + z & = & 32 \\ & 2x & + 4y & + 8z & = & 80 \\ & & & x,y,z & \in & \mathbf{Z}_+ \end{array} $$model = Model()
x = model.addVar(vtype="I", name="x")
y = model.addVar(vtype="I", name="y")
z = model.addVar(vtype="I", name="z")
model.update()
model.addConstr(x + y + z == 32)
model.addConstr(2 * x + 4 * y + 8 * z == 80)
model.setObjective(y + z, GRB.MINIMIZE)
model.optimize()
if model.Status == GRB.Status.OPTIMAL:
print("Opt. Val.=", model.ObjVal)
print("(x,y,z)=", (x.X, y.X, z.X))
ロバスト最適化
確率最適化では，パラメータ（問題に含まれる数値データ）の不確実性を特定の確率分布が既知として扱うが， 実際には，パラメータの確率分布を推定することは難しく，単にパラメータの動きうる範囲だけが与えられている場合が多い． ここでは，データがある範囲内で変化しても，実行可能解になることを保証したロバスト最適化（robust optimization）について考える．
以下の線形最適化問題を考える．
$$ \begin{array}{l l l} maximize & \sum_{j \in N} c_j x_j & \\ s.t. & \sum_{j \in N} a_{ij} x_j \leq b_i & \forall i \in M \end{array} $$ここで，変数 $x_j$ は（負の値もとれる）実数変数であることに注意されたい．
いま，制約の係数 $a_{ij}$ は不確実性もつ確率変数 $\tilde{a}_{ij}$ であり，区間 $[a_{ij}-\hat{a}_{ij}, a_{ij}+\hat{a}_{ij}]$ 内で変化するものとする． ロバスト最適化の目的は，係数が変化しても実行可能性を保証するような解の中で，目的関数を最大にする解を求めることである． 係数 $\tilde{a}_{ij}$ が独立に変化するものと仮定すると，ロバスト最適化の意味での最適解は，以下の線形最適化問題によって得ることができる．
$$ \begin{array}{l l l} maximize & \sum_{j \in N} c_j x_j & \\ s.t. & \sum_{j \in N} a_{ij} x_j + \sum_{j \in N} \hat{a}_{ij} u_j \leq b_i & \forall i \in M \\ & - u_j \leq x_j \leq u_j & \forall j \in N \\ & u_j \geq 0 & \forall j \in N \end{array} $$ここで，$u_j$ は非負の実数変数であり，最適解 $x^*$ においては $u_j= |x^*_j|$ が成立する．
上のロバスト最適化のフレームワークは，パラメータの変化が最悪の場合を想定したものであるが， 現実的には最悪のシナリオを考慮して最適化を立てるのではなく，ばらつきをパラメータによって制御できるものが望ましい． ここでは，制約 $i$ に対して，高々 $\Gamma_i$ 個の変数が変化するという制限を付けたロバスト最適化を考え， 上と同様に線形最適化に帰着されることを示す． 以下では簡単のため，制御パラメータ $\Gamma_i$ は正数であると仮定する．
上の仮定の下でのロバスト最適化は，以下のように書ける．
$$ \begin{array}{l l l} maximize & \sum_{j \in N} c_j x_j & \\ s.t. & \sum_{j \in N} a_{ij} x_j + \max_{S \subseteq N, |S| \leq \Gamma_i} \left\{ \sum_{j \in S} \hat{a}_{ij} |x_j| \right\} \leq b_i & \forall i \in M \end{array} $$上の定式化において，制約 $i$ における $$ \max_{S \in N, |S| \leq \Gamma_i} \left\{ \sum_{j \in S} \hat{a}_{ij} |x_j| \right\} $$ の部分を（$|x_j|$ を定数と見なして）線形最適化として記述すると， $$ \begin{array}{l l l} maximize & \sum_{j \in N} \hat{a}_{ij} |x_j| z_{ij} & \\ s.t. & \sum_{j \in N} z_{ij} \leq \Gamma_i & \\ & 0 \leq z_{ij} \leq 1 & \forall j \in N \end{array} $$ となる． 線形最適化問題の実行可能領域は有界であり，実行可能解（$z=0$）をもつので，強双対定理より，双対問題と同じ最適目的関数値をもつことが言える． 最初の式 $\sum_{j \in N} z_{ij} \leq \Gamma_i$ に対する双対変数を $\theta_i$，2 番目の式 $z_{ij} \leq 1$ に対する双対変数を $y_{ij}$ とすると， 上の問題の双対問題は $$ \begin{array}{l l l} minimize & \Gamma_i \theta_i + \sum_{j \in N} y_{ij} & \\ s.t. & \theta_i + y_{ij} \geq \hat{a}_{ij} |x_j| & \forall j \in N \\ & y_{ij} \geq 0 & \forall j \in N \\ & \theta_i \geq 0 \end{array} $$ となる．
ロバスト最適化の定式化における $\max_{S \subseteq N, |S_i| \leq \Gamma_i} \{ \sum_{j \in S} \hat{a}_{ij} |x_j| \}$ の部分を 上で導いた双対問題に置き換え，さらに $|x_j|$ の絶対値を補助変数 $u_j$ を用いて外すことによって，ロバスト最適化と同値な以下の線形最適化問題を得る．
$$ \begin{array}{l l l} maximize & \sum_{j \in N} c_j x_j & \\ s.t. & \sum_{j \in N} a_{ij} x_j + \Gamma_i \theta_i + \sum_{j \in N} y_{ij} \leq b_i & \forall i \in M \\ & \theta_i + y_{ij} \geq \hat{a}_{ij} u_j & \forall i \in M, j \in N \\ & -u_j \leq x_j \leq u_j & \forall j \in N \\ & y_{ij} \geq 0 & \forall i \in M, j \in N \\ & \theta_i \geq 0 & \forall i \in M \\ & u_j \geq 0 & \forall j \in N \end{array} $$今度は，制約 $i$ に対して係数のベクトル $\tilde{a}_i =(\tilde{a}_{i1},\tilde{a}_{i2},\ldots, \tilde{a}_{i|N|} )$ が， ベクトル $a_i$ を中心とした楕円 $\{ a_i + P_i u \mid \| u \|_2 \leq 1 \}$ 内にあるという制約を付加した場合を考える．
ここで $$ \max_{\| u \|_2 \leq 1} \{ (a_i + P_i u)^T x \} \leq b_i $$ は $$ a_i^T x + \| P_i^T x \|_2 \leq b_i $$ と同値であるので，楕円の不確実性を考えたロバスト最適化は，以下の2次錘最適化問題に帰着される．
$$ \begin{array}{l l l} maximize & \sum_{j \in N} c_j x_j & \\ s.t. & a_i^T x + \| P_i^T x \|_2 \leq b_i & \forall i \in M \end{array} $$栄養問題
線形最適化問題の例として栄養問題(diet problem)を考え，実行不可能性とその対処法について述べる．
あなたは，某ハンバーガーショップの調査を命じられた健康オタクの諜報員だ． あなたは任務のため，毎日ハンバーガーショップだけで食事をしなければならないが， 健康を守るため，なるべく政府の決めた栄養素の推奨値を遵守しようと考えている． 考慮する栄養素は，カロリー(Cal)，炭水化物(Carbo)，タンパク質(Protein)，ビタミンA(VitA)，ビタミンC(VitC)，カルシウム(Cal)，鉄分(Iron)であり， 1日に必要な量の上下限は，以下の表の通りとする． 現在，ハンバーガーショップで販売されている商品は， CQPounder, Big M, FFilet, Chicken, Fries, Milk, VegJuice の7種類だけであり， それぞれの価格と栄養素の含有量は，以下のようになっている． さらに，調査費は限られているので，なるべく安い商品を購入するように命じられている． さて，どの商品を購入して食べれば，健康を維持できるだろうか？ ただし，ここでは簡単のため，商品は半端な数で購入できるものと仮定する．
| 栄養素 $N$ | Cal | Carbo | Protein | VitA | VitC | Calc | Iron | 価格 |
|---|---|---|---|---|---|---|---|---|
| 商品名 $F$ | $n_{ij}$ | $c_j$ | ||||||
| CQPounder | 556 | 39 | 30 | 147 | 10 | 221 | 2.4 | 360 |
| Big M | 556 | 46 | 26 | 97 | 9 | 142 | 2.4 | 320 |
| FFilet | 356 | 42 | 14 | 28 | 1 | 76 | 0.7 | 270 |
| Chicken | 431 | 45 | 20 | 9 | 2 | 37 | 0.9 | 290 |
| Fries | 249 | 30 | 3 | 0 | 5 | 7 | 0.6 | 190 |
| Milk | 138 | 10 | 7 | 80 | 2 | 227 | 0 | 170 |
| VegJuice | 69 | 17 | 1 | 750 | 2 | 18 | 0 | 100 |
| 上限 $a_i$ | 3000 | 375 | 60 | 750 | 100 | 900 | 7.5 | |
| 下限 $b_i$ | 2000 | 300 | 50 | 500 | 85 | 660 | 6.0 |
一般には最適化問題は常に最適解をもつとは限らない． 特に，現実的な問題を考える場合には，（制約条件がきつすぎて）解が存在しない場合も多々ある．
実行可能解が存在しない場合を実行不可能(infeasible)もしくは実行不能と呼ぶ． たとえば，以下の線形最適化問題は，すべての制約を満たす領域（実行可能領域）が 空なので，実行不可能である．
$$ \begin{array}{l c c c} maximize & x_1 + x_2 & & \\ s.t. & x_1 - x_2 & \leq & -1 \\ & -x_1 + x_2 &\leq & -1 \\ & x_1,x_2 & \geq & 0 \end{array} $$また，目的関数値が無限に良くなってしまう場合を非有界(unbounded)と呼ぶ． たとえば，以下の線形最適化問題は，目的関数値がいくらでも大きい解が存在するので， 非有界である．
$$ \begin{array}{l c c c} maximize & x_1 + x_2 & & \\ s.t. & x_1 - x_2 & \geq & -1 \\ & -x_1 + x_2 & \geq & -1 \\ & x_1,x_2 & \geq & 0 \end{array} $$非有界の場合の実行可能領域を以下に示す（intvalpyパッケージ https://pypi.org/project/intvalpy/ を利用している）． 実行不可能な場合は，実行可能領域は外側になるので，空になる．
A = np.array([[1, -1],
[-1, 1],
[1, 0],
[0, 1]])
b = np.array([-1, -1, 0, 0])
intvalpy.lineqs(A, b, title='Feasible Region', color='gray', alpha=0.5, s=10, size=(10,10), save=False, show=True);
実際に実行不可能もしくは非有界の問題をGurobiで解いてみよう．
model = Model()
x ={}
for i in range(1,3):
x[i] = model.addVar(vtype="C", name= f"x({i})")
model.update()
#infeasible
#model.addConstr(x[1]-x[2]<=-1)
#model.addConstr(-x[1]+x[2]<=-1)
#unbounded
model.addConstr(-x[1]+x[2]>=-1)
model.addConstr(x[1]-x[2]>=-1)
model.setObjective(x[1]+x[2], GRB.MAXIMIZE)
model.optimize()
print("Status=", model.status)
いずれの場合も，statusが $4$ の Infeasible or unbounded model という答えが返ってくる． これは，Gurobiの前処理によって問題に解がないことが判明したが，実行不可能か非有界かは分からないということを意味する．
実際に，実行不可能か非有界なのかを判別するのは，前処理の可否を表す DualReductions パラメータを $0$ に設定して， 再求解する必要がある．
また，非有界の場合には，目的関数が無限に良くなるベクトルを得ることができる（これはBendersの分解法などの特殊解法を実装するときに役に立つ）． それには，パラメータ InfUnbdInfo を $1$ に設定し，最適化後にモデルの UnbdRay 属性を参照すれば良い．
model.Params.DualReductions=0
model.Params.InfUnbdInfo=1
model.optimize()
print("Status=", model.status)
print("無限に目的関数が大きくなるベクトル=", model.UnbdRay)
実行不可能性の対処法（制約の逸脱を許すモデル化）
実行不可能な問題に対する現実的な対処法を，栄養問題を用いて解説する．
古典的な栄養問題は線形最適化問題であるので，半端な数の商品を購入することも許されている．
商品の集合を $F$ （Foodの略），栄養素の集合を $N$（Nutrientの略）とする． 栄養素 $i$ の1日の摂取量の下限を $a_i$， 上限を $b_i$ とし， 商品 $j$ の価格を $c_j$，含んでいる栄養素 $i$ の量を $n_{ij}$ とする． 商品 $j$ を購入する個数を非負の実数変数 $x_j$ で表すと，栄養問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{j \in F} c_j x_j & \\ s.t. & a_i \leq \displaystyle\sum_{j \in F} n_{ij} x_j \leq b_i & i \in N \\ & x_j \geq 0 & j \in F \end{array} $$上の問題例は，実行不可能である． ここでは，制約の逸脱を表す以下の2つの変数を追加することによって，対処する．
- $d_i$: 不足変数
- $s_i$: 超過変数
制約は，以下のように変更する．
$$a_i - d_i \leq \displaystyle\sum_{j \in F} n_{ij} x_j \leq b_i +s_i \ \ i \in N$$
また，目的関数にも制約の逸脱ペナルティを追加する．ここで， $M$ は十分に大きな数とする．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{j \in F} c_j x_j + \displaystyle\sum_{i \in N} M (d_i+s_i) & \end{array} $$Gurobiでモデル化してみよう．
F, c, n = multidict(
{
"CQPounder": [
360,
{
"Cal": 556,
"Carbo": 39,
"Protein": 30,
"VitA": 147,
"VitC": 10,
"Calc": 221,
"Iron": 2.4,
},
],
"Big M": [
320,
{
"Cal": 556,
"Carbo": 46,
"Protein": 26,
"VitA": 97,
"VitC": 9,
"Calc": 142,
"Iron": 2.4,
},
],
"FFilet": [
270,
{
"Cal": 356,
"Carbo": 42,
"Protein": 14,
"VitA": 28,
"VitC": 1,
"Calc": 76,
"Iron": 0.7,
},
],
"Chicken": [
290,
{
"Cal": 431,
"Carbo": 45,
"Protein": 20,
"VitA": 9,
"VitC": 2,
"Calc": 37,
"Iron": 0.9,
},
],
"Fries": [
190,
{
"Cal": 249,
"Carbo": 30,
"Protein": 3,
"VitA": 0,
"VitC": 5,
"Calc": 7,
"Iron": 0.6,
},
],
"Milk": [
170,
{
"Cal": 138,
"Carbo": 10,
"Protein": 7,
"VitA": 80,
"VitC": 2,
"Calc": 227,
"Iron": 0,
},
],
"VegJuice": [
100,
{
"Cal": 69,
"Carbo": 17,
"Protein": 1,
"VitA": 750,
"VitC": 2,
"Calc": 18,
"Iron": 0,
},
],
}
)
N, a, b = multidict(
{
"Cal": [2000, 3000],
"Carbo": [300, 375],
"Protein": [50, 60],
"VitA": [500, 750],
"VitC": [85, 100],
"Calc": [660, 900],
"Iron": [6.0, 7.5],
}
)
model = Model("modern diet")
x, s, d = {}, {}, {}
for j in F:
x[j] = model.addVar(vtype="C", name=f"x({j})")
for i in N:
s[i] = model.addVar(vtype="C", name=f"surplus({i})")
d[i] = model.addVar(vtype="C", name=f"deficit({i})")
model.update()
for i in N:
model.addConstr(quicksum(n[j][i] * x[j] for j in F) >= a[i] - d[i], f"NutrLB({i})")
model.addConstr(quicksum(n[j][i] * x[j] for j in F) <= b[i] + s[i], f"NutrUB({i})")
model.setObjective(
quicksum(c[j] * x[j] for j in F) + quicksum(9999 * d[i] + 9999 * s[i] for i in N),
GRB.MINIMIZE,
)
model.optimize()
status = model.Status
if status == GRB.Status.OPTIMAL:
print("Optimal value:", model.ObjVal)
for j in F:
if x[j].X > 0:
print(j, x[j].X)
for i in N:
if d[i].X > 0:
print(f"deficit of {i} ={d[i].X}")
if s[i].X > 0:
print("surplus of {i} ={s[i].X}")
逸脱最小化
Gurobiには，モデルを緩和（整数条件を取り払い線形最適化問題にする）し，さらに制約の逸脱を自動的に最小化してくれるメソッド feasiRelaxS が準備されている． これを使って，実行不可能性に対処することができる．
引数は以下の通り．
- relaxobjtype: 逸脱量の評価方法を指定する．$0$ のとき線形，$1$ のとき2乗，$2$ のとき逸脱した制約の本数
- minrelax: Falseのとき制約の逸脱量を最小化し，Trueのとき逸脱最小化の条件下で目的関数を最小化する
- vrelax: 変数の上下限を緩和するとき True
- crelax: 制約の逸脱を許すとき True
栄養問題の制約の逸脱量の2乗和を最小化し，さらに逸脱を許した解の中で目的関数を最小化するものを求めてみる．
model = Model("modern diet (overcome infeasibility using feaseRelaxS)")
x = {}
for j in F:
x[j] = model.addVar(vtype="C", name=f"x({j})")
model.update()
for i in N:
model.addConstr(quicksum(n[j][i] * x[j] for j in F) >= a[i], f"NutrLB({i})")
model.addConstr(quicksum(n[j][i] * x[j] for j in F) <= b[i], f"NutrUB({i})")
model.setObjective(
quicksum(c[j] * x[j] for j in F),
GRB.MINIMIZE,
)
model.feasRelaxS(relaxobjtype= 1, minrelax=True, vrelax=False, crelax=True)
model.optimize()
violation = 0.
for i in N:
val = sum(n[j][i] * x[j].X for j in F)
violation+= ( max(a[i]-val,0.))**2 + (max(val-b[i],0.))**2
print(a[i]," <= ",val," <= ",b[i])
print("violation=", violation)
print("Obj. func. val=", sum(c[j] * x[j].X for j in F))
model = Model("modern diet (overcome infeasibility using IIS)")
x = {}
for j in F:
x[j] = model.addVar(vtype="C", name=f"x({j})")
model.update()
for i in N:
model.addConstr(quicksum(n[j][i] * x[j] for j in F) >= a[i], f"NutrLB({i})")
model.addConstr(quicksum(n[j][i] * x[j] for j in F) <= b[i], f"NutrUB({i})")
model.setObjective(
quicksum(c[j] * x[j] for j in F),
GRB.MINIMIZE,
)
model.Params.OutputFlag=0
model.computeIIS()
model.write("model.ilp")
while True:
model.computeIIS()
for con in model.getConstrs():
if con.IISConstr:
print("Remove constraint:", con.ConstrName)
model.remove(con)
break
model.optimize()
status = model.Status
if status == GRB.OPTIMAL:
break
混合問題
ここでは，古典的な問題である混合問題(product mix problem}を考え，問題のパラメータが変動するロバスト最適化が， 2次錐最適化問題に帰着できることを示す．
4種類の原料を調達・混合して1種類の製品を製造している工場を考える． 原料には，3種類の成分が含まれており，成分1については $20 \%$ 以上に， 成分2については $30 \%$ 以上に，成分3については $40 \%$ 以上になるように混合したい． 各原料の成分含有比率は，以下の表（各原料に含まれる成分 $k$ の比率） のようになっている． 原料の単価が1トンあたり $5,6,8,20$万円であるとしたとき，どのように原料を混ぜ合わせれば， 製品1トンを最小費用で製造できるだろうか？
| 成分 | 1 | 2 | 3 |
|---|---|---|---|
| 原料 $1$ | 25 | 15 | 30 |
| 原料 $2$ | 30 | 30 | 10 |
| 原料 $3$ | 15 | 65 | 0 |
| 原料 $4$ | 10 | 5 | 80 |
原料 $i$ の価格を $p_i$，原料 $i$ に含まれる成分 $k$ の比率を $a_{ik}$， 製品に含まれるべき成分 $k$ の比率の下限を $LB_k$ とする． 原料 $i$ の混合比率を表す実数変数を $x_i$ としたとき， 今回の混合問題は，以下のように記述できる．
$$ \begin{array}{rcl} minimize & \sum_{i=1}^4 p_i x_i & \\ s.t. & \sum_{i=1}^4 x_{i} =1 & \\ & \sum_{i=1}^4 a_{ik} x_i \geq LB_k & k=1,2,3 \\ & x_i \geq 0 & i = 1,2,3,4 \end{array} $$ところが実際の問題では，このように比率が正確に決まっていることは珍しく， いくらかの誤差を伴っているのが普通である． 例えばカタログ上では 原料 $1$ に成分 $3$ が $30\%$含まれることになっているが， 実際の原料ではそれが $29\%$ かもしれないし，$30.5\%$ とか多いかもしれない． このように誤差を考慮し，最悪の誤差でも制約を満たすという条件の下で 最適化することをロバスト最適化(robust optimization)という．
ここでは成分の比率 $a_{ik}$ が誤差 $e_{ik}$ をもっている状況を考え， この誤差に関して「実際の値が$a_{ik}$より常に上であったり，常に下であったりということはなく，$a_{ik}$が中心である」と仮定する． つまり誤差 $e_{ik}$ は正の値も負の値も同じように取りうると仮定する．
また，誤差はそれほどは大きくないと考えなければこのような問題は 解けるはずもないので，誤差の量に関して何らかのモデル化が必要である． ここでは次が成り立っているモデルを考える． $$ \sum_{i=1}^4 e_{ik}^2 \leq \epsilon^2 \qquad k = 1,2,3 $$ ただし， $\epsilon>0$ は誤差の量を押さえるパラメータである． この $\epsilon$ が大きければ大きな誤差を見込むことになる．
さて，上記の誤差を考慮の上，製品に対しては成分の比率を保証しなければならない． もともとの制約が $ \sum_{i=1}^4 a_{ik} x_i \geq LB_k $ であったとすれば，誤差が含まれた制約は以下のように書ける． $$ \sum_{i=1}^4 (a_{ik} + e_{ik}) x_i \geq LB_k $$ ここで誤差 $e_{ik}$ は $\sum_{i=1}^4 e_{ik}^2 \leq \epsilon^2 $ を満たす どんな値でも取りうることを考慮し，そのいずれに対しても 成分は $LB_k$ 以上なければならないので， ロバスト性に関する制約は次のようになる． $$ \min\left\{\sum_{i=1}^4 (a_{ik} + e_{ik}) x_i \ : \ \sum_{i=1}^4 e_{ik}^2 \leq \epsilon^2 \right\} \geq LB_k $$ あるいは誤差 $e_{ik}$ に関連する項だけ左辺にまとめれば， 以下のようにも書ける． $$ \min\left\{\sum_{i=1}^4 e_{ik} x_i \ : \ \sum_{i=1}^4 e_{ik}^2 \leq \epsilon^2 \right\} \geq LB_k - \sum_{i=1}^4 a_{ik} x_i $$
ここで左辺の $\min$ を取っているところは $4$次元の半径 $\epsilon$ の球上で 線形関数 $\sum_{i=1}^4 e_{ik} x_i$ を最小化している．ただし，変数は $e_{ik}$ である．
球上での線形関数の最小化問題の最適解は以下となることが知られている． $$ e_{ik} = - \epsilon\frac{x_i}{\| x \|} \ \ \ \ i=1,2,3,4 $$ ただし，$\| x \| = \sqrt{\sum_{i=1}^4 x_i^2}$ である． よって， $$ \min\left\{\sum_{i=1}^4 e_{ik} x_i \ : \ \sum_{i=1}^4 e_{ik}^2 \leq \epsilon^2 \right\} = -\epsilon \sum_{i=1}^4 \frac{x_i^2}{\| x \|} = -\epsilon \| x \| $$ となり，結局， 以下の2次錐制約で表現できる． $$ \epsilon \| x_i \| \leq -LB_k + \sum_{i=1}^4 a_{ik} x_i $$
まとめると，混合問題のロバスト最適化問題は以下の2次錐最適化問題として定式化できる．
$$ \begin{array}{rcl} minimize & \sum_{i=1}^4 p_i x_i &\\ s.t. & \sum_{i=1}^4 x_{i} =1 & \\ & \sqrt{\epsilon^2 \sum_{i=1}^4 x_i^2} \leq -LB_k + \sum_{i=1}^4 a_{ik} x_i & k=1,2,3 \\ & x_i \geq 0 & i = 1,2,3,4 \end{array} $$このモデルを使って，誤差の上限 $\epsilon$ を $0$ から $0.05$ まで変えたときの目的関数値の変化を計算してみよう．
def prodmix(I, K, a, p, epsilon, LB):
"""prodmix: robust production planning using soco
Parameters:
I - set of materials
K - set of components
a[i][k] - coef. matrix
p[i] - price of material i
LB[k] - amount needed for k
Returns a model, ready to be solved.
"""
model = Model("robust product mix")
x, rhs = {}, {}
for i in I:
x[i] = model.addVar(vtype="C", name=f"x({i})")
for k in K:
rhs[k] = model.addVar(vtype="C", name=f"rhs({i})")
model.update()
model.addConstr(quicksum(x[i] for i in I) == 1)
for k in K:
model.addConstr(rhs[k] == -LB[k] + quicksum(a[i, k] * x[i] for i in I))
model.addConstr(
quicksum(epsilon * epsilon * x[i] * x[i] for i in I) <= rhs[k] * rhs[k]
)
model.setObjective(quicksum(p[i] * x[i] for i in I), GRB.MINIMIZE)
model.update()
model.__data = x, rhs
return model
def make_data():
a = {
(1, 1): 0.25,
(1, 2): 0.15,
(1, 3): 0.2,
(2, 1): 0.3,
(2, 2): 0.3,
(2, 3): 0.1,
(3, 1): 0.15,
(3, 2): 0.65,
(3, 3): 0.05,
(4, 1): 0.1,
(4, 2): 0.05,
(4, 3): 0.8,
}
epsilon = 0.01
I, p = multidict({1: 5, 2: 6, 3: 8, 4: 20})
K, LB = multidict({1: 0.2, 2: 0.3, 3: 0.2})
return I, K, a, p, epsilon, LB
I, K, a, p, epsilon, LB = make_data()
obj_list = []
for i in range(5):
epsilon = i * 0.01
model.Params.OutputFlag = 0
model = prodmix(I, K, a, p, epsilon, LB)
model.optimize()
print("obj:", model.ObjVal)
x, rhs = model.__data
for i in I:
print(i, x[i].X)
obj_list.append(model.ObjVal)
pd.Series(obj_list).plot()
Gurobiの主なパラメータ
パラメータの変更方法
- setParam(“パラメータ名”,値 )
- model.Params.パラメータ名 =値
一般
- TimeLimit :制限時間（秒）
- OutputFlag : 出力制御
- 0 : ログファイルを出力しない
- 1 : ログを出力する（既定値）
- Threads 使用するスレッド数
- LogFile : ログファイル名の設定 （出力した結果は grblogtoolsパッケージのparse関数で読み込み，解析することができる）
線形最適化
Method 線形最適化手法の選択（連続モデルか混合整数最適化の根ノード）
- -1 ：自動選択 (既定値）
- 0 ：主単体法
- 1 ：双対単体法（分枝ノードの中での既定値）
- 2 ：内点法（2次，2次錐での既定値）
- 3 ：非決定性並流（線形最適化での既定値）
- 4 ：決定性並流
並流最適化：1つのスレッドで双対単体法， 1つのスレッドで主単体法， 残りのスレッドで並列内点法
- Crossover (線形最適化のための）内点法から単体法への移行方法
- -1 ：自動設定
- 0 ： クロスオーバーを行わない
- 1 ：双対変数の固定 => 主変数の固定 => 主単体法
- 2 ：双対 => 主 =>双対
- 3 ：主 => 双対 => 主
- 4 ：主 => 双対 => 双対
- BarHomogenuous 同次自己双対内点法の選択
- -1 ：自動選択（分枝限定法のノードで使われたときのみ使用）
- 0 ：使用しない
- 1 ：使用する（実行不可能性や非有界の判定をしたいとき）
- InfUnbInfo 実行不可能性や非有界の情報を得たいとき
- 0 ：使用しない(既定値）
- 1 ：使用する（端線情報 UnbdRayや実行不可能性の証拠 FarkasDual, FarkasProodを得ることができる）
混合整数最適化
- MIPGap : 相対誤差（0以上の実数）
- 1e-4 : 既定値
- Presolve 前処理
- -1 : 自動選択（既定値）
- 0 : オフ
- 1 : 控えめ
- 2 : 目一杯
- GomoryPasses Gomoryの切除平面を追加する回数
- -1 : 自動選択（既定値）
- 0 : オフ
- 正数 : 回数
- Cut 切除平面の追加の制御
- -1 : 自動選択（既定値）
- 0 : 切除平面を追加しない
- 1 : 多少追加
- 2 : 多めに追加
- 3 : 目一杯追加
- MIPFocus 混合整数最適化の戦略
- 0 : バランス型（既定値）
- 1 : 実行可能解探索優先
- 2 : 最適性検証優先
- 3 : 限界値改良優先
- Heuristics ヒューリスティクスに使う時間の割合（0以上1以下）
- 0.05 : 既定値
- ZeroObjNodes 目的関数を $0$ とし，実行可能解の探索を行うヒューリスティクスを適用する分枝ノードの個数；
- -1 : 他のヒューリスティクスで実行可能解が得られないとき，根ノードで実行（既定値）
VarBranch 分枝変数選択戦略
- -1 : 自動選択（既定値）
- 0 : 擬似被約費用分枝
- 1 : 擬似潜在価格分枝
- 2 : 最大実効不能分枝
- 3 : 強分枝
#### チューニング
モデルインスタンス
modelのtuneメソッドで，パラメータのチューニングができ， model.write("ファイル名.prm")でチューニング結果のファイルへの出力ができる．
- TuneTimeLimit チューニング用の計算時間上限（秒）
- -1 : 自動選択（既定値）
- TuneCriteria チューニングの尺度（基本は最適解を得るまでの計算時間．最適解を得られないときの評価尺度を決めるためのパラメータ）
- -1 : 自動選択（既定値）
- 0 : 第2評価尺度を用いない
- 1 : 誤差
- 2 : 最良の実行可能解
- 3 : 限界値