JavaScriptを有効にしてください

CasADiとIPOPTで非線形計画問題を解く

 ·  ☕ 4 min read  ·  ✍️ Helve

はじめに

CasADiは自動微分と非線形最適化のためのライブラリである。C++で実装されており、C++, Python, Matlab, Octaveのインターフェースを備えている。
本記事ではPythonを使い、CasADiで非線形最適化問題を解く方法をまとめた。

CasADiにはIPOPTという非線形ソルバが含まれているため、これを用いる。IPOPTは主双対内点法を利用したソルバであり、大規模な非線形問題を高速に解くことができる。ただし、問題は連続である必要があり、大域的最適解を求めるには問題が凸である必要がある。

環境は以下の通り。

ソフトウェア バージョン
Python 3.7.4
CasADi 3.5.5

ソフトウェアのインストール

pipでCasADiをインストールする。CasADiには最適化ソルバも含まれており、同時にインストールされる。

pip install casadi

対象とする非線形計画問題

以下の制約付き非線形最小化問題を考える。

$$ \begin{array}{ll} \text{minimize} & f(x_1, x_2) = x_1^2 + x_2^2 \\ \text{subject to} & x_1 x_2 \ge 1, \\ & x_1 \ge 0, x_2 \ge 0 \end{array} $$

図示すると以下のようになり、$(x_1, x_2)=(1, 1)$で最小値$2$をとる。青色で図示した領域は実行可能領域を示す。

pyomo-nlp

CasADiのソースコード

上記の問題をCasADiを使って記述し、最適化を実行したコードは以下のようになる。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import casadi

opti = casadi.Opti()

# 変数を定義
x1 = opti.variable()
x2 = opti.variable()

# 目的関数を定義
OBJ = x1**2 + x2**2
opti.minimize(OBJ)

# 制約条件を定義
opti.subject_to( x1*x2 >= 1 )
opti.subject_to( x1 >=0 )
opti.subject_to( x2 >=0 )

opti.solver('ipopt') # 最適化ソルバを設定
sol = opti.solve() # 最適化計算を実行

print(f'評価関数:{sol.value(OBJ)}')
print(f'x1: {sol.value(x1)}')
print(f'x2: {sol.value(x2)}')

実行結果は以下のようになり、最適値を得られている。

This is Ipopt version 3.12.3, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
(中略)
EXIT: Optimal Solution Found.
      solver  :   t_proc      (avg)   t_wall      (avg)    n_eval
       nlp_f  |        0 (       0)        0 (       0)        19
       nlp_g  |        0 (       0)        0 (       0)        19
  nlp_grad_f  |   1.00ms ( 76.92us) 997.00us ( 76.69us)        13
  nlp_hess_l  |        0 (       0)        0 (       0)        11
   nlp_jac_g  |        0 (       0)        0 (       0)        13
       total  |  22.00ms ( 22.00ms)  21.94ms ( 21.94ms)         1
評価関数:1.9999999825046224
x1: 0.9999999956261556
x2: 0.9999999956261556

ソースコードの解説

以下、ソースコードの詳細を述べる。

最適化問題のオブジェクトの作成

まず、Optiクラスを用いて最適化問題のオブジェクトを作成する。

1
opti = casadi.Opti()

変数の定義

次に、variableクラスを用いて変数を定義する。

1
2
x1 = opti.variable()
x2 = opti.variable()

variable()に引数を与えない場合、変数はスカラになる。一方、
x = opti.variable(2)
のように整数を引数にとることで、ベクトルにできる。
ベクトルにした場合、個々の要素にはx[0], x[1]のようにアクセスする。

目的関数の定義

次に、目的関数を定義し、opti.minimizeメソッドに引数として与える。

1
2
OBJ = x1**2 + x2**2
opti.minimize(OBJ)

もしくは、以下のように目的関数を直接与えることも可能である。

1
opti.minimize(x1**2 + x2**2)

しかし、直接与えた場合、最適化後の目的関数の値を取得するには、以下のようにする必要があるため、目的関数を別途定義したほうが便利である。
sol.stats()['iterations']['obj'][-1]

制約条件の定義

subject_toメソッドで制約条件を設定する。

1
2
3
opti.subject_to( x1*x2 >= 1 )
opti.subject_to( x1 >=0 )
opti.subject_to( x2 >=0 )

等式制約は==で、不等式制約は<=, >=, <, >でそれぞれ設定する。
また、変数に上限と下限の両方がある場合、以下のように定義できる。

1
opti.subject_to( opti.bounded(0,x,1) )

ソルバの指定と最適化

最後に、ソルバを指定して最適化を行う。

1
2
opti.solver('ipopt') # 最適化ソルバを設定
sol = opti.solve() # 最適化計算を実行

solverメソッドでソルバを指定する。ここではIPOPTを用いる。次に、solveメソッドを実行すると、最適化が行われる。solには最適化の結果が格納される。

最適化ソルバのオプション設定

最適化ソルバ(ここではIPOPT)にオプションを設定するためには、solverメソッドに辞書形式で渡す。
ここで、solverメソッドの2番目の引数はCasADiプラグインのオプション(本記事では説明を省略)、3番目の引数が最適化ソルバのオプションになる。

1
2
3
4
5
p_opts = {}
s_opts = {'print_level' : 4,
          'tol' : 1E-3}
opti.solver('ipopt', p_opts, s_opts)
sol = opti.solve()

この例では、以下のように設定した。

  • 最適化計算のログ出力の詳細度合い (print_level) を4(0~12の範囲。数字が大きいほど詳細になる)
  • 最適解の許容誤差 (tol) を1E-3

IPOPTで設定可能なオプションについては以下のページを参照。
Ipopt: Ipopt Options

一覧の形式で確認したい場合には以下のページを参照。
https://www.coin-or.org/Bonmin/option_pages/options_list_ipopt.html

最適化結果の評価

最適化後の変数や目的関数の値を取得するには、solvalueメソッドに変数名を渡す。

1
2
3
4
5
6
>>> sol.value(OBJ)
1.9999999825046224
>>> sol.value(x1)
0.9999999956261556
>>> sol.value(x2)
0.9999999956261556

最適化問題の概要を表示するには、sol.disp(), sol.stats()などを実行する。

1
2
3
4
5
6
7
8
9
>>> sol.disp()
Opti {
  instance #7
  #variables: 2 (nx = 2)
  #parameters: 0 (np = 0)
  #constraints: 3 (ng = 3)
  CasADi solver allocated.
  CasADi solver was called: Solve_Succeeded
}
1
2
3
4
5
6
7
8
9
>>> sol.stats()
{'iter_count': 11,
 'iterations': {'alpha_du': [0.0,
   1.0,
   0.00011355335462534469,
(中略)
 't_wall_nlp_jac_g': 0.0,
 't_wall_total': 0.011968,
 'unified_return_status': 'SOLVER_RET_UNKNOWN'}

参考

CasADiの公式リファレンス
CasADi - Docs

シェアする
ブログランキングに参加中です。クリックして頂けると励みになります。

Helve
WRITTEN BY
Helve
関西在住、電機メーカ勤務のエンジニア。Twitterで新着記事を配信中です