読者です 読者をやめる 読者になる 読者になる

MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

Python製凸最適化モデリングツールCVXPYの使い方

A First Course in Optimization Theory

A First Course in Optimization Theory

目次

はじめに

先日、最適化技術の一つである

線形計画法や二次計画法の紹介をしましたが、

myenigma.hatenablog.com

myenigma.hatenablog.com

最適化の式が、そこで紹介した標準形と同じでないと

最適化が実現できないという問題があります。

 

実はかなりの数の最適化の問題は、

上手く数式変換することで

上記の標準形に変換することができます。

このような変換を自動的に実施するソフトを

最適化モデリングツールと呼ぶのですが、

今回はこの最適化モデリングツールの一つであり、

Python製のCVXPYというツールの紹介をしたいと思います。

 

CVXPYとは?

CVXPYはPython製の最適化モデリングツールです。

github.com

 

このモデリングツールを使うことにより、

二次計画法などの標準形ではない形式の最適化問題を

標準形に変換することなく、計算することができます。

つまり元の数式の記述に近いまま、

最適化計算が可能だということです。

 

その他のモデリングツールに関しては、

下記の資料を参照下さい。

myenigma.hatenablog.com

 

インストール方法

下記に各プラットフォームの

インストール方法が記述されています。

 

Macの場合、pipだとエラーが出てしまったので、

上記のようにAnacondaでインストールする必要がありました。

 

CVXPYによる最小二乗法のサンプルコード

下記のコードは最小二乗法をCVXPYで解く、

サンプルコードです。

コスト関数と制約条件を設定することで、

標準形ではない形式の最適化問題を解くことができます。

#! /usr/bin/python 
# -*- coding: utf-8 -*- 

from cvxpy import *
import numpy
import matplotlib.pyplot as plt

m = 30
n = 20
numpy.random.seed(1)
A = numpy.random.randn(m, n)
b = numpy.random.randn(m)

x = Variable(n)
# コスト最適化
objective = Minimize(sum_squares(A*x - b))
# 制約
constraints = [0 <= x, x <= 1]
prob = Problem(objective, constraints)

result = prob.solve()
# 最適値
print "optimal parameter:\n",x.value

# ラグランジュパラメータ
print "Lagrange parameter\n",constraints[0].dual_value

#最適化の状態
print ("status:"+prob.status)

上記のコードのように、

最適化の結果を取得することができます。

 

シンプルな最適化問題を解いてみる(その2)

下記の最適化問題を解いてみます。

f:id:meison_amsl:20161110144114p:plain

 

下記のコードで解くことができます。

#! /usr/bin/python 
# -*- coding: utf-8 -*- 

from cvxpy import *
import numpy
import matplotlib.pyplot as plt


z1 = Variable()
z2 = Variable()

objective = Minimize(abs(z1+5)+abs(z2-3))

constraints = [2.5 <= z1, z1 <= 5]
constraints += [-1.0 <= z2, z2 <= 1]
prob = Problem(objective, constraints)

result = prob.solve()

# 最適値
print "z1",z1.value
print "z2",z2.value
print "cost:",prob.value

f:id:meison_amsl:20161110144206p:plain

 

シンプルなモデル予測制御をCVXPYで解いてみる

下記の記事を参照下さい。

myenigma.hatenablog.com

 

CVXPYを使う時の基礎

下記のマニュアルの翻訳です。

名前空間

CVXPYでは、maxやmin,sumなど、

既存のビルドイン関数と

同じ名前の関数が多く設定されています。

そんな時に

from cvxpy import *

とすると、関数がオーバライドされてしまうので注意が必要です。

最適問題の変更

最適化問題(prob)を設定して、

一回計算した後でも、コスト関数を再設定したり、

制約条件を追加しても大丈夫です。

# Replace the objective.
prob.objective = Maximize(x + y)
print "optimal value", prob.solve()

# Replace the constraint (x + y == 1).
prob.constraints[0] = (x + y <= 3)
print "optimal value", prob.solve()

最適化問題が解けない場合の状態管理

最適化問題が解けなかったかったり、

コスト関数が無限方向に発散している場合、

prob.statusが"infeasible"や"unbounded"

になります。

この場合prob.valueはinf or -infになり、

最適化変数は変更されません。

 

また、推定精度が低い場合は

それぞれ下記のようなステータスが表示されます。

  • “optimal_inaccurate”

  • “unbounded_inaccurate”

  • “infeasible_inaccurate

また、Solverが完全に最適化に失敗すると、

SolverErrorという例外を投げます。

この場合はソルバーを変更することを検討した方が良いようです。

 

ちなみにcvxpyでは、それぞれのステータスを簡単に評価するための

エイリアスが設定されています。

  • OPTIMAL
  • INFEASIBLE
  • UNBOUNDED
  • OPTIMAL_INACCURATE
  • INFEASIBLE_INACCURATE
  • UNBOUNDED_INACCURATE

上記のエイリアスを使って、

下記のように簡単に最適化が成功したかどうかを判断することができます。

if prob.status == OPTIMAL:

ベクトルと行列

最適化用の変数はスカラーやベクトル、行列に対応しています。

a = Variable()

# vector .
x = Variable(5)

# Matrix  4 x 7.
A = Variable(4, 7)

またcvxpyではNumpyのndarraysやmatrices,

SciPyのsparse matricesにも対応しており、

最適化変数以外のパラメータはこれらの型が使えます。

制約条件

等式制約や不等式制約は、==,<=,>=で指定することができます。

例えば0<=x, x<=1とすると、xのすべての要素は0と1の間に設定されます。

もし、行列の不等式制約において,半定値(semi-definite)凸制約を

追加したい場合は、こちらの資料を参照下さい。

また、0<=x<=1, x==y==2のように、

複数の等号、不等号を同時に設定することはできません。

パラメータ

CVXPYでは、定数をパラメータとして設定することができます。

これにより最適化問題を再構築すること無く、

最適化のパラメータを変更することができます。

 

パラメータは変数のように、スカラー、ベクトル、行列で設定できます。

またパラメータには符号を設定することができます。

この符号はDeiscipline Convex Programmingという

最適化問題に使うことができます。

# スカラ
m = Parameter(sign="positive")

# ベクトル
c = Parameter(5)

# 行列
G = Parameter(4, 7, sign="negative")

# numpyの行列をコピー
G.value = -numpy.ones((4, 7))

パラメータは下記のように初期化することができます。

# Create parameter, then assign value.
rho = Parameter(sign="positive")
rho.value = 2

# Initialize parameter with a value.
rho = Parameter(sign="positive", value=2)

 

このパラメータを上手く使うことにより、

パラメータを変更しただけの同一の最適化問題を

並列処理することができます。

並列処理の方法としては、

最適化を解く関数を作り、

multiprocessingモジュールのPoolを使うことで、

簡単に実行することができます。

from multiprocessing import Pool

# Assign a value to gamma and find the optimal x.
def get_x(gamma_value):
    gamma.value = gamma_value
    result = prob.solve()
    return x.value

# Parallel computation (set to 1 process here).
pool = Pool(processes = 1)
x_values = pool.map(get_x, gamma_vals)

 

CVXPYにおける関数

こちらの説明を元にまとめています。

演算子

+, -, *, /の二項演算子は、関数として扱われます。

+, -はアフィン関数として、

*と/もCVXPYもアフィン関数として扱われます。

なぜなら、CVXPYでは a*bはどちらか一つがスカラーの時しか計算されず、

a/bではbがスカラーの時のみ計算されるからです。

インデクシングとスライシング

すべてのスカラーでない項は、a[i,j]という形でインデックス指定することができます。

このインデックス指定もアフィン関数です。

ある項が列ベクトルの時、expr[i,0]でインデックス指定することもできますが、

expr[i]でもできます。行ベクトルも同じように簡略化したインデックス指定ができます。

 

また、Pythonの標準のスライシングも対応しています。(ex a[i:j:k, r])

加えてCVXPYでは、インデックスリストやブール配列によるデータ取得にも対応しています。

詳しい方法は下記を参照下さい。

Indexing — NumPy v1.11 Manual

転置

すべての項における転置はa.Tで取得できます。

転置もアフィン関数です。

べき乗

すべての項におけるべき乗は、a**p または power(a,p)で計算できます。

べき乗はアフィン関数ではありません。

スカラー関数

スカラー関数は入力に対して、スカラー値を返す関数になります。

cvxpyには様々なスカラー関数が準備されています。

詳しいスカラー関数の説明は下記を参照下さい。

Functions — CVXPY 0.4.8 documentation

上記のスカラー関数の中で、

sum_entriesやnorm,max_entries, min_entriesは各軸に対して実行することができます。

(ex: func(a,axis=0))

要素関数

cvxpyでは各要素毎に計算を実施する関数も多数用意されています。

詳しい要素関数の説明は下記を参照下さい

Functions — CVXPY 0.4.8 documentation

ベクター/行列関数

ベクターや行列を返す関数も多数用意されています。

詳しい各関数の説明は下記を参照下さい。

Functions — CVXPY 0.4.8 documentation

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

A First Course in Optimization Theory

A First Course in Optimization Theory

Convex Optimization

Convex Optimization