notebooks/OPT_robust_basic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ロバスト最適化の基本\n",
    "\n",
    "参考：\n",
    "\n",
    "* [Robust Optimization](https://www2.isye.gatech.edu/~nemirovs/FullBookDec11.pdf) \n",
    "* [Robust optimization - A comprehensive survey](https://www.honda-ri.de/pubs/pdf/1808.pdf) \n",
    "\n",
    "\n",
    "何かのシステムを最適化することを考えましょう．\n",
    "例えば工場の製造プロセスのコスト最小化や，運送業者の最適な運搬経路の計画などが相当します．\n",
    "ここでのエンジニアの任務は目的関数$f(x)$を最小化するような最適化変数$x$の設計です．\n",
    "\n",
    "このプロセスは数学的には次のような最適化問題によって表現されます：\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\text { optimize: } & f(\\mathbf{x}) \\\\\n",
    "\\text { subject to: } & g_i(\\mathbf{x}) \\leqslant 0, \\quad i=1, \\ldots, I \\\\\n",
    "& h_j(\\mathbf{x})=0, \\quad j=1, \\ldots, J\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "この問題を最適化すること自体，しばしばNP困難だったりと難しいことがあります．\n",
    "しかし，そもそもこの最適化ができたとして，それは嬉しいでしょうか？\n",
    "いくつか考えるべきことがあります：\n",
    "\n",
    "1. この関数$f$や$g$は実世界の問題やモデルを表現しますが，完全に正しい$f$や$g$をを求めることに意味はあるでしょうか？たいてい実世界は何か不確実な要素を含みます．そのため，$f, g$に精度を求めすぎることに意味はなさそうです．\n",
    "2. もし完全な最適解$x^*$が求まったとして，それは実際にデプロイできるでしょうか？例えば浮動小数点の精度などの影響で，完全な$x^*$はデプロイできそうにありません．\n",
    "3. 上の最適化問題は「定常」な問題を考えていますが，実世界の問題はたいてい変化します（温度や部品の摩耗など）．\n",
    "\n",
    "上の問題の最適解は，大概の場合でモデルや解のデプロイ時に生じる小さな変化に敏感です．\n",
    "そこで，ロバスト最適化では，何らかの不確実性にさらされたとしてもパフォーマンスがあまり変化しないような解を求めることを目標にします．\n",
    "\n",
    "ここで重要になるのは，「どのようにロバスト性をモデル化するか？」と「どうやって関連する最適化問題を解くか？」の２つです．\n",
    "この２つを紹介する前に，日本の田口玄一博士が考えたタグチメソッドを紹介します．一般に，品質を保証する製品の製造は次の流れで実現されるべきです：\n",
    "1. システム設計：製品に関する特性を調査し，モデル化するフェーズ\n",
    "2. 変数設計：製品の質が一定以上に達するように，特定の変数を最適化するフェーズ\n",
    "3. Tolerance設計：ステップ２で設計された変数を調整するフェーズ\n",
    "\n",
    "タグチメソッドは自動車のバックミラーの最適なロバスト設計などで使われてます：[Robust optimization of an automobile rearview mirror for vibration reduction](https://link.springer.com/article/10.1007/s001580100107)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ロバスト線型計画問題\n",
    "\n",
    "通常の線形計画問題は次の最適化問題をときます：\n",
    "\n",
    "$$\n",
    "\\min_x \\{c^T x + d: Ax \\leq b\\}\n",
    "$$\n",
    "\n",
    "ここで，$x, c \\in \\mathbb{R}^n$，$d\\in \\mathbb{R}$，$A\\in \\mathbb{R}^{m\\times n}$，そして$b\\in \\mathbb{R}^m$とします．\n",
    "（通常の線形計画問題では定数ベクトル$d$は無視しても問題ないですが，ロバスト最適化を考えるときはこれも考慮すると嬉しいことがあります．また，線形な等式は不等式に変形できることにも注意しましょう（[Convex Optimization](https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf)のp.132）．）\n",
    "\n",
    "簡略化のため，使用するデータ$(c, d, A, b)$はまとめて$(m+1)\\times (n+1)$の行列$D$としておきましょう．\n",
    "\n",
    "$$\n",
    "D=\\left[\\begin{array}{c|c}\n",
    "c^T & d \\\\\n",
    "\\hline A & b\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "とします．\n",
    "\n",
    "まず，このデータが揺れてしまう状況を定義します．\n",
    "\n",
    "---\n",
    "\n",
    "**不確実線形計画問題**\n",
    "\n",
    "次のような線形計画問題の集合を，不確実線形計画問題（$\\text{LO}_\\mathcal{U}$）と呼びます：\n",
    "\n",
    "$$\n",
    "\\left\\{\\min _x\\left\\{c^T x+d: A x \\leq b\\right\\}\\right\\}_{(c, d, A, b) \\in \\mathcal{U}}\n",
    "$$\n",
    "\n",
    "ここで，$\\mathcal{U} \\subset \\mathbb{R}^{(m+1) \\times(n+1)}$を不確実集合と呼びます．\n",
    "この不確実集合は次の形式でかけるとします：\n",
    "\n",
    "$$\n",
    "\\mathcal{U}=\\{\\left[\\begin{array}{c|c}\n",
    "c^T & d \\\\\n",
    "\\hline A & b\n",
    "\\end{array}\\right]=\\underbrace{\\left[\\begin{array}{c|c}\n",
    "c_0^T & d_0 \\\\\n",
    "\\hline A_0 & b_0\n",
    "\\end{array}\\right]}_{\\text {nominal data $D_0$}}+\\sum_{\\ell=1}^L \\zeta_{\\ell} \\underbrace{\\left[\\begin{array}{c|c}\n",
    "c_{\\ell}^T & d_{\\ell} \\\\\n",
    "\\hline A_{\\ell} & b_{\\ell}\n",
    "\\end{array}\\right]}_{\\text {basic shifts $D_\\ell$ }}: \\zeta \\in \\mathcal{Z} \\subset \\mathbb{R}^L\\}\n",
    "$$\n",
    "\n",
    "つまり，何らかの基本的な摂動の集合$D_1, D_2, \\dots, D_\\ell$があり，それがノミナルデータ$D_0$に加わるとします．\n",
    "\n",
    "---\n",
    "\n",
    "この不確実線形計画問題$\\text{LO}_\\mathcal{U}$はまだ線形計画問題の集合に過ぎません．\n",
    "この集合に対して，ロバスト最適化問題を定義していきます（なぜそんなことを考えるのか？は教科書p.9参照）．そのためにいくつか定義を導入します．\n",
    "\n",
    "---\n",
    "\n",
    "**ロバスト実行可能解**\n",
    "\n",
    "$$\n",
    "A x \\leq b \\quad \\forall(c, d, A, b) \\in \\mathcal{U}\n",
    "$$\n",
    "\n",
    "を満たす解$x$をロバスト実行可能解と呼びます．\n",
    "\n",
    "---\n",
    "\n",
    "---\n",
    "\n",
    "**ロバスト実現値**\n",
    "\n",
    "不確実集合の中で最悪ケースの実現値をロバスト実現値と呼ぶ．\n",
    "\n",
    "$$\n",
    "\\widehat{c}(x)=\\sup _{(c, d, A, b) \\in \\mathcal{U}}\\left[c^T x+d\\right]\n",
    "$$\n",
    "\n",
    "---\n",
    "\n",
    "---\n",
    "\n",
    "**Robust Counterpart**\n",
    "\n",
    "次の問題を，$\\text{LO}_\\mathcal{U}$のRobust Counterpartと呼ぶ．\n",
    "\n",
    "$$\n",
    "\\min _x\\left\\{\\widehat{c}(x)=\\sup _{(c, d, A, b) \\in \\mathcal{U}}\\left[c^T x+d\\right]: A x \\leq b \\forall(c, d, A, b) \\in \\mathcal{U}\\right\\}\n",
    "$$\n",
    "\n",
    "つまりこれは，ロバスト実行可能解の中から，最もロバスト実現値が小さくなる解を探しています．\n",
    "\n",
    "**等価な形式**\n",
    "\n",
    "このRobust Counterpartは次の問題と等価であることに注意しましょう．\n",
    "\n",
    "$$\n",
    "\\left.\\min _{x, t}\\left\\{t: \\begin{array}{rl}\n",
    "c^T x-t & \\leq-d \\\\\n",
    "A x & \\leq b\n",
    "\\end{array}\\right\\} \\forall(c, d, A, b) \\in \\mathcal{U}\\right\\}\n",
    "$$\n",
    "\n",
    "このようにスラック変数$t$を導入すると，$d$を導入した意味がわかります．\n",
    "もし$d$が決定的ならば，$t\\to t-d$とすればよいですが，$d$を非決定的な場合はこれはできません．\n",
    "また，このように目的関数に関するデータが非決定的でも，$t$を目的関数にした決定的な問題に還元することができます．\n",
    "\n",
    "**決定的な目的の場合**\n",
    "\n",
    "$c$と$d$が決定的な場合は，Robust Counterpartは次の形になります．\n",
    "\n",
    "$$\n",
    "\\min _x\\left\\{\\widehat{c}(x)=\\sup _{(A, b) \\in \\mathcal{U}}\\left[c^T x+d\\right]: A x \\leq b \\forall(A, b) \\in \\mathcal{U}\\right\\}\n",
    "$$\n",
    "\n",
    "すると，目的と制約が独立しているのがわかります．つまり，制約の行をそれぞれ独立させた\n",
    "\n",
    "$$\n",
    "(A x)_i \\leq b_i \\Leftrightarrow a_i^T x \\leq b_i\n",
    "$$\n",
    "\n",
    "を考え，そのRobust Counterpartを\n",
    "\n",
    "$$\n",
    "a_i^T x \\leq b_i \\forall\\left[a_i ; b_i\\right] \\in \\mathcal{U}_i\n",
    "$$\n",
    "\n",
    "とします．ここで，$\\mathcal{U}_i=\\left\\{\\left[a_i ; b_i\\right]:[A, b] \\in \\mathcal{U}\\right\\}$は元の問題の$i$行だけ抜きたしたものです．\n",
    "\n",
    "すると，元の問題は$\\widehat{\\mathcal{U}}=\\mathcal{U}_1 \\times \\ldots \\times \\mathcal{U}_m$　の形を考えても等価になります．\n",
    "\n",
    "* $x$がロバスト実行可能界ならば，不確実集合の凸包についても実行可能\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 不確実線形計画問題への取り組み方\n",
    "\n",
    "不確実線形計画問題のRobust Counterpartはどのように解けばいいでしょうか．\n",
    "これを考えるために，次の定義を導入します：\n",
    "\n",
    "---\n",
    "\n",
    "集合$X^{+} \\subset \\mathbb{R}_x^n \\times \\mathbb{R}_u^k$は次を満たすとき，$X \\subset \\mathbb{R}_x^n$を「表す」という．\n",
    "\n",
    "$$\n",
    "X=\\left\\{x: \\exists u:(x, u) \\in X^{+}\\right\\}\n",
    "$$\n",
    "\n",
    "---\n",
    "\n",
    "これを使うと嬉しい例として，例えば次のような最適化問題を考えてみます．\n",
    "\n",
    "$$\n",
    "\\min _x\\left\\{f(x) \\text { s.t. } x \\text { satisfies } \\mathcal{S}_i, i=1, \\ldots, m\\right\\}\n",
    "$$\n",
    "\n",
    "ここで，$\\mathcal{S}_i$は変数$x$についての制約を持つシステムとします．\n",
    "一方で，$\\mathcal{S}_i^+$を変数$x, v^i$についての制約を持つシステムとして，$\\mathcal{S}_i$を表すとします．このとき，\n",
    "\n",
    "$$\n",
    "\\min _{x, v^1, \\ldots, v^m}\\left\\{f(x) \\text { s.t. }\\left(x, v^i\\right) \\text { satisfies } \\mathcal{S}_i^{+}, i=1, \\ldots, m\\right\\}\n",
    "$$\n",
    "\n",
    "は上の問題と等価です．\n",
    "\n",
    "これを踏まえて，次のRobust Counterpartを考えてみましょう（目的関数が決定的なら，このような形だけ考えても十分です）\n",
    "\n",
    "$$\n",
    "a^T x \\leq b \\quad \\forall\\left([a ; b]=\\left[a^0 ; b^0\\right]+\\sum_{\\ell=1}^L \\zeta_{\\ell}\\left[a^{\\ell} ; b^{\\ell}\\right]: \\zeta \\in \\mathcal{Z}\\right)\n",
    "$$\n",
    "\n",
    "---\n",
    "\n",
    "**簡単な例**\n",
    "\n",
    "$$\n",
    "\\mathcal{Z}=\\operatorname{Box}_1 \\equiv\\left\\{\\zeta \\in \\mathbb{R}^L:\\|\\zeta\\|_{\\infty} \\leq 1\\right\\}\n",
    "$$\n",
    "\n",
    "を考えます．このとき，\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccr} \n",
    "& {\\left[a^0\\right]^T x+\\sum_{\\ell=1}^L \\zeta_{\\ell}\\left[a^{\\ell}\\right]^T x \\leq b^0+\\sum_{\\ell=1}^L \\zeta_{\\ell} b^{\\ell}} & \\forall\\left(\\zeta:\\|\\zeta\\|_{\\infty} \\leq 1\\right) \\\\\n",
    "\\Leftrightarrow & \\sum_{\\ell=1}^L \\zeta_{\\ell}\\left[\\left[a^{\\ell}\\right]^T x-b^{\\ell}\\right] \\leq b^0-\\left[a^0\\right]^T x \\quad \\forall\\left(\\zeta:\\left|\\zeta_{\\ell}\\right| \\leq 1, \\ell=1, \\ldots, L\\right) \\\\\n",
    "\\Leftrightarrow \\quad & \\max _{-1 \\leq \\zeta_{\\ell} \\leq 1}\\left[\\sum_{\\ell=1}^L \\zeta_{\\ell}\\left[\\left[a^{\\ell}\\right]^T x-b^{\\ell}\\right]\\right] \\leq b^0-\\left[a^0\\right]^T x &\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "が成立します．これの最大を与える$\\zeta_\\ell$は明らかに$\\sum_{\\ell=1}^L\\left|\\left[a^l\\right]^T x-b^{\\ell}\\right|$です．\n",
    "\n",
    "よって，制約としては\n",
    "\n",
    "$$\n",
    "\\left[a^0\\right]^T x+\\sum_{\\ell=1}^L\\left|\\left[a^{\\ell}\\right]^T x-b^{\\ell}\\right| \\leq b^0\n",
    "$$\n",
    "\n",
    "を考えれば良いことがわかります．\n",
    "さらに変形すれば，次の不等式が元の不等式を表現するのがわかります．\n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{array}{l}\n",
    "-u_{\\ell} \\leq\\left[a^{\\ell}\\right]^T x-b^{\\ell} \\leq u_{\\ell}, \\ell=1, \\ldots, L \\\\\n",
    "{\\left[a^0\\right]^T x+\\sum_{\\ell=1}^L u_{\\ell} \\leq b^0}\n",
    "\\end{array}\\right.\n",
    "$$\n",
    "\n",
    "これより，元は無限の制約だったものが，$L$個の有限の制約になっていることがわかります．\n",
    "これを行列形式に変形しましょう．\n",
    "今回は$x$の変わりに$x'=[x, [u_1, u_2, \\cdots, u_L]]$を目的変数とします．\n",
    "\n",
    "1. $\\left[a^{\\ell}\\right]^T x-b^{\\ell} \\leq u_{\\ell}$の実現\n",
    "\n",
    "変形して$\\left[a^{\\ell}\\right]^T x - u_{\\ell} \\leq b^{\\ell}$です．\n",
    "また，$i$番目だけ$1$のベクトル$e_i=[0, 0, \\cdots, 1, \\cdots, 0]$を考えます．\n",
    "すると，$a'_\\ell=[a^\\ell\\mid -e_\\ell]$とすれば，$[a'_\\ell]^T x'\\leq b^\\ell$で$\\ell$番目の制約を満たせます．\n",
    "\n",
    "2. $-u_{\\ell} \\leq \\left[a^{\\ell}\\right]^T x-b^{\\ell}$の実現\n",
    "\n",
    "変形して$-\\left[a^{\\ell}\\right]^T x-u_{\\ell} \\leq -b^{\\ell}$なので，$a''_\\ell=[-a^\\ell \\mid -e_\\ell]$として，\n",
    "$[a''_\\ell]^T x'\\leq -b^\\ell$で$\\ell$番目の制約を満たせます．\n",
    "\n",
    "3. ${\\left[a^0\\right]^T x+\\sum_{\\ell=1}^L u_{\\ell} \\leq b^0}$の実現\n",
    "\n",
    "全部$1$のベクトル$u=[1, 1, \\cdots, 1]$を使って，\n",
    "$a'''_\\ell=[a^0\\mid u]$とすれば，$[a'''_\\ell]^T x'\\leq b^0$で制約を満たせます．\n",
    "\n",
    "試しに簡単な不確実線形計画問題を解いてみましょう．p.7の問題を実装します．\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ノミナルモデルの解:  -8819.657744624841 [  0.         438.78894252  17.5515577    0.        ]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import linprog\n",
    "\n",
    "\n",
    "c0 = np.array([100, 199.9, -5500, -6100])\n",
    "d0 = 0.0\n",
    "b0 = np.array([0.0, 1000, 2000, 800, 100000, 0, 0, 0, 0])\n",
    "A0 = np.array([\n",
    "[-0.01, -0.02, 0.500, 0.600],\n",
    "[1, 1, 0, 0],\n",
    "[0, 0, 90, 100 ],\n",
    "[0, 0, 40, 50],\n",
    "[100, 199.9, 700, 800],\n",
    "[-1, 0, 0, 0],\n",
    "[0, -1, 0, 0],\n",
    "[0, 0, -1, 0],\n",
    "[0, 0, 0, -1],\n",
    "])\n",
    "\n",
    "\n",
    "# ノミナルモデルを先に解いてみます\n",
    "\n",
    "res = linprog(c0, A_ub=A0, b_ub=b0)\n",
    "print(\"ノミナルモデルの解: \", res.fun, res.x)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ロバストモデルの解:  -8272.619999790008 [8.77659368e+02 4.39005286e-02 1.74144617e+01 4.39005286e-02\n",
      " 4.38829684e-02 1.75602114e-05]\n"
     ]
    }
   ],
   "source": [
    "# Robust Counterpartをときます\n",
    "# 摂動集合は上の簡単な例の場合と同じです．\n",
    "\n",
    "L = 2\n",
    "A1 = np.zeros_like(A0)\n",
    "b1 = np.zeros_like(b0)\n",
    "A1[0, 0] = 5e-5\n",
    "\n",
    "A2 = np.zeros_like(A0)\n",
    "b2 = np.zeros_like(b0)\n",
    "A2[0, 1] = 4e-4\n",
    "\n",
    "As = [A1, A2]\n",
    "bs = [b1, b2]\n",
    "new_A = []\n",
    "new_b = []\n",
    "\n",
    "for l, AB in enumerate(zip(As, bs)):\n",
    "    A, B = AB\n",
    "    e = np.zeros(L)\n",
    "    e[l] = 1.0\n",
    "    for a, b in zip(A, B):\n",
    "        # 1の実現\n",
    "        al1 = np.hstack((a, -e))\n",
    "        bl1 = b\n",
    "        new_A.append(al1)\n",
    "        new_b.append(bl1)\n",
    "\n",
    "        # 2の実現\n",
    "        al2 = np.hstack((-a, -e))\n",
    "        bl2 = -b\n",
    "        new_A.append(al2)\n",
    "        new_b.append(bl2)\n",
    "\n",
    "\n",
    "# 3の実現\n",
    "for a, b in zip(A0, b0):\n",
    "    a3 = np.hstack((a, np.ones(L)))\n",
    "    b3 = b\n",
    "    new_A.append(a3)\n",
    "    new_b.append(b3)\n",
    "\n",
    "new_A = np.array(new_A)\n",
    "new_b = np.array(new_b)\n",
    "new_c = np.hstack((c0, np.zeros(L)))\n",
    "\n",
    "# 最終的に2L+1倍の制約ができます\n",
    "assert new_A.shape == ((2*L + 1) * A0.shape[0], A0.shape[1] + L)\n",
    "\n",
    "res = linprog(new_c, A_ub=new_A, b_ub=new_b)\n",
    "print(\"ロバストモデルの解: \", res.fun, res.x)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "教科書のp.10と似た答えは出てますが，ちょっと違いますね..."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}