szmlb.net

tips for robotics

PICOSメモ (半正定値計画問題)

半正定値計画問題(Semi definite programming: SDP)を解きたくてライブラリを探していたところ, この記事にたどり着いた.
記事では, Pythonの最適化モデリングライブラリPICOSを使っていて, 使いやすそうだったので使ってみたくなった.
ユーザーフレンドリーな形で設計されているためより簡単に問題を記述できるらしい.
サンプルは本家のサイトにたくさんあるが, 複雑な例題が多かったのでシンプルな問題を解いてみた.
以下, いくつかサンプルコードをメモしておく.

基本的なコードの書き方

基本的な流れは以下の通り.

1. インポート

import picos as pic
import cvxopt as cvx
import scipy as sp
import scipy.linalg as sl
import numpy as np

2. 係数行列の定義

cnp = np.array([[2,  1]])
c = pic.new_param('c', cnp)

3. picosインスタンスの生成

sdp = pic.Problem()

4. 説明変数の定義

x = sdp.add_variable('x', (cnp.shape[1], 1))

5. 制約の記述

sdp.add_constraint(-x[0] + x[1] <= 1.0)
sdp.add_constraint(x[0] + x[1] >= 2.0)
sdp.add_constraint(x[1] >= 0.0)
sdp.add_constraint(x[0] - 2.0 * x[1] <= 4.0)

6. 目的関数の記述

sdp.set_objective('min', c*x)

7. 解く

sdp.solve()

線形計画問題 (Linear programming: LP)を解く

SDPは, LPの対称行列空間への拡張と考えることができるので, 最適化問題のクラスとしては, LPを内包した形になっている. よって, LPも問題なく解ける (ライブラリの処理としてはLPのソルバを呼んでいるだけ?)

import picos as pic
import cvxopt as cvx
import scipy as sp
import scipy.linalg as sl
import numpy as np

def testSolveLP():
    """
    線形計画問題を解く
    CVXOPTの線形計画問題のサンプル:

    min 2.0 * x1 +  1.0 * x2
    s.t. -x1 + x2 <= 1.0
         x1 + x2 >= 2.0
         x2 >= 0.0
         x1 - 2.0 * x2 <= 4.0
    """
    cnp = np.array([[2,  1]])
    c = pic.new_param('c', cnp)

    sdp = pic.Problem()

    x = sdp.add_variable('x', (cnp.shape[1], 1))
    sdp.add_constraint(-x[0] + x[1] <= 1.0)
    sdp.add_constraint(x[0] + x[1] >= 2.0)
    sdp.add_constraint(x[1] >= 0.0)
    sdp.add_constraint(x[0] - 2.0 * x[1] <= 4.0)

    sdp.set_objective('min', c*x)

    print(sdp)
    sdp.solve()
    print(x.value)

二次計画問題 (Quadratic programming: QP)を解く

SDPはQPも内包している. 例えば以下のQPの標準形をSDPに変換することを考える.

f:id:szmlb:20180101225006j:plain
QPの目的関数と制約条件

{\bf Q}が正定行列であれば, その平方根行列を求めることができ, \frac{1}{2}{\bf Q} = {\bf M}^T{\bf M}のように分解できる. さらにシュールの補題を駆使して問題を変形すると, 以下のようになる.

f:id:szmlb:20180101230629j:plain
QPの標準形と等価なSDPの問題設定

問題の式変形は以下の書籍を参考にした.

システム制御のための最適化理論 (システム制御工学シリーズ)

システム制御のための最適化理論 (システム制御工学シリーズ)

上記の問題を, PICOSでは以下のように簡潔に記述できる.

def testSolveQP():
    """
    二次計画問題をSDP形式でわざわざ解く
    CVXOPTの二次計画問題のサンプル:

    min 2.0 * x1**2 + x2**2 + x1 * x2 + x1 + x2
    s.t. x1 >= 0.0
         x2 >= 0.0
         x1 + x2 = 1.0
    """

    Anp = np.array([[4, 1],
                    [1, 2]])
    bnp = np.array([[1], [1]])

    Mnp = sl.sqrtm(Anp/2.0) # 行列Anpの平方根

    b = pic.new_param('b', bnp)
    M = pic.new_param('M', Mnp)
    I = pic.new_param('I', sp.eye(bnp.shape[0]) )

    sdp = pic.Problem()
    f0 = sdp.add_variable('f0', (1, 1))
    x = sdp.add_variable('x', (bnp.shape[0], 1))
    sdp.add_constraint((( f0 - b.T*x & (M*x).T ) // (M*x & I))>>0 ) # シュールの補題より
    sdp.add_constraint(x[0] >= 0.0)
    sdp.add_constraint(x[1] >= 0.0)
    sdp.add_constraint(x[0] + x[1] == 1.0)

    sdp.set_objective('min', f0)

    print(sdp)
    sdp.solve()
    print(x.value)

最小固有値問題を解く

任意の行列の固有値のなかで最小のものを求める問題を, PICOSで解く.
最小固有値問題はSDPじゃなくても解けるが, 固有値が全て正になるような制約条件下で特定の関数を最小化(最大化)したいようなケースが, 工学の分野では見受けられる.

最小固有値問題を解くためには, 以下のSDPを解けば良い.

f:id:szmlb:20180101232358j:plain
最小固有値問題の目的関数と制約条件

PICOSで記述すると以下のようになる.

def testSolveMinEigVal():
    """
    最小固有値問題をSDP形式でわざわざ解く:
    A=[4, 1; 1, 2]の最小固有値1.59を求める
    """

    Anp = np.array([[4, 1],
                    [1, 2]])

    A = pic.new_param('A', Anp)
    I = pic.new_param('I', sp.eye(Anp.shape[0]) )

    sdp = pic.Problem()
    mu = sdp.add_variable('mu', (1, 1))
    sdp.add_constraint(A - mu * I >> 0 )
    sdp.set_objective('min', -mu)

    print(sdp)
    sdp.solve()
    print(mu.value)

LMI制約付き問題を解く

最後に一つ制御っぽいサンプルを. LMI制約付き問題を解く.
制御の分野では, ロバスト安定化問題や多目的制御問題をLinear Matrix Inequality: LMIと呼ばれる(半)正定性を制約とした問題設定に落とし込む設計法がよく見受けられる. SDPソルバに投げれば確実に解けるため (BMIと呼ばれるちょっと簡単に解けない問題のクラスもあるが, ここでは言及しない). MatlabにはYALMIPと言う有名なLMIソルバが存在するが, YALMIPのサンプルコードを眺めているとPICOSでも解けそうだったのでポーティングしてみた.

具体的に解いた問題は「LMIによるシステム制御」と言う本の例題4.3.

LMIによるシステム制御 - ロバスト制御系設計のための体系的アプローチ

LMIによるシステム制御 - ロバスト制御系設計のための体系的アプローチ

具体的な問題の内容について興味がある方は, 書籍を読むことをお勧めします.
(ものすごく)ざっくりと説明すると, 制御対象のパラメータ変動がある場合の安定性について判定するために, LMIを解くと言う問題設定になっている. 問題の数式表現は以下のようになる.

f:id:szmlb:20180102000124j:plain
ロバスト安定性判別のためのLMI

上式が解けて, 解が存在すれば安定.
PICOSを使うと以下のように解ける.

def testSolveLyapnovEq():
  """
  "LMIによるシステム制御" p.85の例4.3(ロバスト安定性解析)を解いてみる

  P - I >> 0
  P * A1 + A1.T * P << 0
  P * A2 + A2.T * P << 0
  P * A3 + A3.T * P << 0
  """

  A1a = np.array([[-1.0, 2.9], [1.0, -3.0]])
  A2a = np.array([[-0.8, 1.5], [1.3, -2.7]])
  A3a = np.array([[-1.4, 0.9], [0.7, -2.0]])

  A1 = pic.new_param('A1', A1a )
  A2 = pic.new_param('A2', A2a )
  A3 = pic.new_param('A3', A3a )
  I = pic.new_param('I', sp.eye(2) )

  sdp = pic.Problem()
  Xp = sdp.add_variable('X', (2, 2), vtype='symmetric')
  sdp.add_constraint(Xp-I>>0)
  sdp.add_constraint(Xp*A1+A1.T*Xp << 0)
  sdp.add_constraint(Xp*A2+A2.T*Xp << 0)
  sdp.add_constraint(Xp*A3+A3.T*Xp << 0)

  sdp.set_objective("min", pic.trace(Xp))

  print(sdp)
  sdp.solve()
  print(Xp.value)