最適制御向け最適化ライブラリOpEnに入門する

※記事内に商品プロモーションを含むことがあります。

はじめに

Rust製の最適制御向け最適化ライブラリOpEnに入門するためチュートリアルの非線形計画問題を解いたので、備忘録を兼ねてまとめた。
Opengen basics · OpEn

OpEnは非凸最適化問題を解くためのライブラリである。
PANOC (Proximal-Averaged Newton-type method for Optimal Control) という最適化手法が実装されている。
最適化問題のソースコードをRustで出力・コンパイルして解くため、計算の実行速度は速い。
また、OpEnはPythonやMATLABのインターフェイスを備えているため、Rustの知識がなくても利用できる。

環境

本記事ではOpEnのインターフェイスはPython+CasADiとした。
環境は以下の通り。また、OSはWindows 10である。

Rust 1.55.0
Python 3.8.8
opengen 0.6.9
OpEn 0.7.1
CasADi 3.5.5

2022/01/30追記 opengen v0.6.8以前ではPythonからOpEnを実行するとエラーが発生します。
最適制御向け最適化ライブラリOpEnのRust build of TCP interface failedエラーについて – Helve Tech Blog

OpEnのインストール方法は以下の記事を参照。
Rust製最適化ライブラリOpEnのインストール – Helve Tech Blog

また、CasADiは以下でインストールできる。

1

pip install casadi

最適化問題の定義

ここからはチュートリアルの問題に沿って進める。
まず、CasADiとopengenをインポートする。

1
2


import casadi.casadi as cs
import opengen as og

次に、変数と目的関数を定義する。目的関数はRosenbrock関数である。

$$\phi(u) = \sum_{i=0}^{4} \left\{ p_1 (u_{i+1}-u_i^2)^2 + (p_0 - u_i)^2 \right\}$$

1
2
3


u = cs.SX.sym("u", 5)               # decision variable
p = cs.SX.sym("p", 2)               # parameter
phi = og.functions.rosenbrock(u, p) # cost function

uは変数、pはパラメータである。パラメータは最適化計算の実行時に指定する。
なお、OpEnは最適化制御を解くために開発されたため、制御工学で入力を表す記号uを未知変数としていると思われる。

phiは目的関数である。ここでは用意された関数og.functions.rosenbrockを使っているが、自分で定義することも可能である。
例：

1

phi = sum([p[1]*(u[i+1]-u[i]**2)**2 + (p[0]-u[i])**2 for i in range(4)])

さらに、制約条件を定義する。

$$ u_1^2 + u_2^2 < 1.5^2, \\ -1 < u_3 < 0, -2 < u_4 < 10, -3 < u_5 < -1, $$

1
2
3
4
5
6


ball = og.constraints.Ball2(None, 1.5)  # ball centered at origin
rect = og.constraints.Rectangle(xmin=[-1,-2,-3], xmax=[0, 10, -1])
segment_ids = [1, 4]
bounds = og.constraints.CartesianProduct(segment_ids, [ball, rect])

problem = og.builder.Problem(u, p, phi).with_constraints(bounds)

Ball2(None, 1.5)は半径1.5の超球内に変数を収める制約である。最初の引数でNoneの代わりに数値を入れると、球の中心を指定できる。
Rectangle(xmin, xmax)は最小値、最大値を指定する制約である。
また、segment_idsで制約式を付加する変数のインデックスを指定している。

ビルドとソルバの設定

生成されるオプティマイザのメタ情報を設定する（省略可能）。

1
2
3
4
5


meta = og.config.OptimizerMeta()                \
    .with_version("0.0.0")                      \
    .with_authors(["P. Sopasakis", "E. Fresk"]) \
    .with_licence("CC4.0-By")                   \
    .with_optimizer_name("the_optimizer")

生成されたソースコードをRustでビルドときの設定を指定する。
with_build_modeは"debug"または"release"を指定できる。"release"を指定するとビルドは遅くなるが、最適化の実行は速くなる。
また、.with_tcp_interface_configはソルバをPythonまたはMATLABから呼び出すときに必要。

1
2
3
4


build_config = og.config.BuildConfiguration()  \
    .with_build_directory("python_build")      \
    .with_build_mode("debug")                  \
    .with_tcp_interface_config()

ソルバの設定を行う。

1
2
3
4


solver_config = og.config.SolverConfiguration() \
    .with_lbfgs_memory(15)                      \
    .with_tolerance(1e-5)                       \
    .with_max_inner_iterations(155)

最適化問題、メタ情報、ビルドとソルバの設定をbuilderにまとめる。

1
2
3
4
5


builder = og.builder.OpEnOptimizerBuilder(problem,
    metadata=meta,
    build_configuration=build_config,
    solver_configuration=solver_config)
builder.build()

最適化計算の実行

PythonとオプティマイザはTCPで通信するため、TCPサーバを起動する。

1
2
3
4
5


mng = og.tcp.OptimizerTcpManager('python_build/the_optimizer')
mng.start()

pong = mng.ping() # check if the server is alive
print(pong)

OptimizerTcpManagerには、metaとbuild_configで指定したオプティマイザのパスを指定する。
TCPサーバとの通信が確立できれば、print(pong)で{'Pong': 1}と表示される。

次に、最適化計算を実行する。callにはパラメータの値を指定する。

1

response = mng.call([1.0, 50.0]) # call the solver over TCP

最適化の結果を確認する。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11


if response.is_ok():
    # Solver returned a solution
    solution_data = response.get()
    u_star = solution_data.solution
    exit_status = solution_data.exit_status
    solver_time = solution_data.solve_time_ms
else:
    # Invocation failed - an error report is returned
    solver_error = response.get()
    error_code = solver_error.code
    error_msg = solver_error.message

response.get()で最適化の結果が得られる。主なプロパティは以下の通り。

プロパティ	説明
`solution`	最適解
`exit_status`	状態（`Converged`など）
`solve_time_ms`	実行時間（ミリ秒）
`cost`	評価関数の値
`num_outer_iterations`	反復回数

最後にTCPサーバを停止する。

1

mng.kill()

最適化結果の確認

最後に、最適化結果を確認する。

1
2
3
4
5
6


>>> print(u_star)
[0.4940235276165833, 0.23381718884060976, 0.0, 0.006622525726116431, -1.0]
>>> print(exit_status)
Converged
>>> print(solver_time)
3.0322

u_starは最適解である。また、計算時間は3.0322ミリ秒であった。