連続な関数の最適化手法を紹介したい

この記事は TUT Advent Calender 2022 の7日目の記事です。

はじめに

はじめまして、豊橋技科大1系B4のこたっくです。ついこの前まで豊橋技科大のロボコンサークルであるとよはし☆ロボコンズの制御班に所属し、NHK学生ロボコン2022・ABU Robocon 2022に参加していました。

Twitter twitter.com

サイト tutrobo.rm.me.tut.ac.jp

またロボコンが終わってからはOUXTのメンバーとして船の速度制御系の開発を担当しています。OUXTはMaritime Robotx Challengeという完全自動運転な船を作って様々なタスクにチャレンジするロボコンに参加する団体で日本各地から学生・社会人いろんな人が集まって活動しています。

www.ouxt.jp

OUXTではROSを活用してロボットを動かしますがなんとコードは全て公開されています。「なんとなくROSは知ってるけど実際にコード書こうと思ったら分からない」「パッケージの設計どうしたらいいか分からない」みたいに思ったらOUXTのリポジトリを覗いてみると解決するかもしれません。自分も学生ロボコン中にめちゃくちゃ参考にさせていただきました。

github.com

豊橋の大学のアドベントカレンダーなので豊橋のことも軽く触れておきましょう。皆さん中華料理屋「美楽」はご存知でしょうか? とよはし☆ロボコンズでは美楽の唐揚げ定食を食べきれると良いロボットが作れるとされています。詳しくは以下の記事で詳しく解説されています。ロボコンやってる人もやってない人も豊橋に来た際にはぜひ食べてみてください。

shouyuzukenotako.hatenablog.com

本題

最近はロボコンを言い訳に後回しにして遅れに遅れていた研究を巻き返すためにひーひー言ってる毎日です。無事に卒業できるといいですね(白目)
タスクに追われてぎりぎりで生きてると次こそは効率的にやろうとか思うんですけどなかなか上手くいかないもんですね。特に試行錯誤でいい感じにするタスクは時間を使うばかりで進捗効率は悪いことが多いです。

ということでやはりそういうタスクは最適化問題に落とし込んで効率的に行うべきです。

ということで今回は制約付き非線形最適化法である逐次二次最適化(SQP法)について紹介しようと思います。

二次計画法(QP法)

逐次二次最適化に触れる前にいくつかの前提となる手法について紹介します。 まずは二次計画法です。二次計画法とは以下の二次形式で表される問題を最適化する手法を指します。

 \displaystyle
\begin{align}
\underset{x}{\text{min}} \space & \dfrac{1}{2}x^TQx + c^T x\\
\text{s.t.} \space & Ax = b\\
& A_{eq}x \leq b_{eq} \\
& x \in \mathbb{R}^n
\end{align}\\
\\

制約がなければ最急降下法ニュートン法で解くことが出来ます。最急降下法ニュートン法高専や技科大の授業でやった人も多いのではないでしょうか。

制約が在る最適化問題ではKKT条件というものを満たすような解を探します。

まずはスラック変数を導入して不等式制約を等式制約に置き換えます。

 \displaystyle
\begin{align}
\underset{x}{\text{min}} \space & \dfrac{1}{2}x^TQx + c^T x\\
\text{s.t.} \space & Ax = b\\
& A_{eq}x + s = b_{eq} \\
& s > 0
\end{align}\\
\\

この式に問題におけるKKT条件は

 \displaystyle
\begin{eqnarray}
Qx+c + A^T\lambda + A_{eq}\lambda_{eq} &=& 0 \\
\lambda_i s_i &=& \rho\\
Ax-b+s &=& 0\\
A_{eq}x-b_{eq} &=& 0\\
s &>& 0\\
\lambda &>& 0\\
\end{eqnarray}\\
\\

ここで各変数を微少量だけ変化させた場合を考えます

 \displaystyle
\begin{eqnarray}
x &\leftarrow& x + \Delta x\\
s &\leftarrow& s + \Delta s\\
\lambda &\leftarrow& \lambda + \Delta \lambda\\
\lambda_{eq} &\leftarrow& \lambda_{eq} + \Delta \lambda_{eq}\\
\end{eqnarray}\\
\\

上式をKKT条件に代入して微小変化量についてまとめると以下のようになります。

 \displaystyle
\begin{bmatrix}
Q    & 0    & A^T    & A_{eq}^T \\
A    & I         & 0 & 0\\
0    & D_\lambda & D_s    & 0 \\
A_eq & 0    & 0      & 0 \\
\end{bmatrix}

\begin{bmatrix}
\Delta x \\
\Delta s \\
\Delta \lambda \\
\Delta \lambda_{eq} \\
\end{bmatrix}

= 
-
\begin{bmatrix}
Qx + c + A^T\lambda + A_{eq}^T \lambda_{eq}\\
Ax - b + s\\
-\rho I + \lambda \bigodot s \\
A_{eq} x - b_{eq}\\
\end{bmatrix}\\
\\

ここで

 \displaystyle
\begin{array}
\\
\lambda \bigodot s &= \text{(アダマール積)}\\
D_\lambda &= \text{diag}(\lambda)\\
D_s &= \text{diag}(s)\\
\\
\end{array}

です。

この連立方程式を解いて各変数の更新量を求めて更新する作業を収束するまで繰り返すことで解を求めることが出来ます。

なにやら複雑な式に見えますが制約の数が0の時を考えると左辺は係数行列の左上の Qだけ残り右辺は一番上の成分の Qx + cが残る形になり、ニュートン法と一致することが分ります。ここでは制約を考慮しながら解くニュートン法ぐらいに捉えてもらって大丈夫です。

BFGS法

次にBFGS法について紹介します。

BFGS法とは準ニュートン法等で使用される目的関数のヘッシアンを近似する行列を求める手法です。 連続関数の最適化問題が収束するためには目的関数のヘッシアンが半正定値である必要がありますが、任意の非線形な目的関数のヘッシアンが半正定値である保証はないのでヘッシアンを半正定値の行列で近似してやります。

以下の更新式を使ってイテレーションごとに近似ヘッシアンを更新します

 \displaystyle
B_{k+1} = B_{k} - \dfrac{B_k s_s \left(  B_k s_k \right)^\top}{ (s_k)^\top B_k s_k } + \dfrac{y_k (y_k)^\top}{(s_k)^\top y_k}\\
\\

近似ヘッシアンの初期値としては単位行列などが使われます。

また実際にはSQP法では {(s_k)^\top y_k} > 0を満たすとは限らないためパウエルの修正BFGS公式が使われます。

逐次二次最適化(SQP法)

ようやくですが本題の逐次二次最適化について紹介します。 逐次二次最適化とは以下のような形式の問題を最適化する手法の一つです。 目的関数や制約関数が非線形の問題に対して使用することが出来ます。

 \displaystyle
\begin{align}
\underset{x}{\text{min}} \space & f(x)\\
\text{s.t.} \space & g_i(x) = 0, \space i = 0, \dots , m_e\\
&  g_i(x) \leq 0, \space i = m_e+1, \dots , l\\
& x \in \mathbb{R}^n\\
\\
\end{align}

「逐次」で「二次最適化」な名前から分かるようにこの最適化手法は任意の最適化問題を最適化変数 x周りで二次近似して問題を解き最適化変数を更新する手順を解が収束するまで繰り返し行います。

SQP法の手順はおおまかに以下になります。

  1. 最適化変数 x、近似ヘッシアン Bの初期値を決定する
  2. 現在の最適化変数 x周りで目的関数と制約関数のヤコビアン、制約関数を評価する。
  3. サブ問題であるQP問題を解き探索方向を決定する。
  4. メリット関数の直線探索でステップ幅を計算する
  5.  xを更新する
  6. メリット関数の制約重みを更新する
  7. BFGS法で Bを更新する
  8. 収束していなければ2.に戻る

ここで一番大事なのは3. のサブ問題を解いて探索方向を決定することです。 二次形式に近似して二次計画法で解くことを考えるとサブ問題であるQPは以下のように書くことが出来ます。

 \displaystyle
\begin{align}
\underset{x}{\text{min}} \space & \dfrac{1}{2}d^TBd + \nabla f(x)^T x\\
\text{s.t.} \space & \nabla g_i(x)^Td + g_i(x) = 0, \space i = 0, \dots , m_e\\
& \nabla g_i(x)^Td + g_i(x) \leq 0, \space i = m_e+1, \dots , l\\
\\
\end{align}

この問題を解くと探索方向 dが求まるのでステップ幅を直線探索で求めて xを更新します。
この作業を収束するまで繰り返します。

ここではサブ問題をニュートン法ベースの二次計画法で解く方法を紹介しましたが、修正コレスキー分解と最小二乗法によって解くものをSLSQP(Sequential Least SQuares Programming)と言います。 SLSQPは論文とFORTRANのルーチンが公開された後、PythonのScipyやCのnloptにも移植されています。

気になる方は以下のタイトルで論文を調べてみてください

Kraft, D. A software package for sequential quadratic programming. 1988. Tech. Rep. DFVLR-FB 88-28, DLR German Aerospace Center -- Institute for Flight Mechanics, Koln, Germany.

例題

実際にSQP法を使ってみましょう。

以下の個人開発しているライブラリの中にSQP法を実装してあります。 github.com

cpp_roboticsは自分が実装したいものを実装する趣味用のC++ライブラリです。

問題設定

以下の問題を最適化してみましょう。

 \displaystyle
\begin{align}
\underset{x}{\text{min}} \space & 100 (x_2-x_1^2)^2 + (1 - x_1)^2    \\
\text{s.t.} \space & x_1^2 + x_2^2 \leq 1 \\
& x_1 \geq 0 \\
\\
\end{align}

目的関数はグラフ化してみるとこんな感じの関数です。

図は下記より引用

qiita.com

制約なしの時の最適解は (x_1, x_2) = (\pm 1,  1) です。

制約の1つ目は解が原点から半径1の円の中に収まっていること、2つ目が2つある最適解のうち正の方に収束することを制約にしています。

コード

#define CR_USE_MATPLOTLIB
#include <iostream>
#include <cpp_robotics/optimize/optimize.hpp>
#include "cpp_robotics/third_party/matplotlib-cpp/matplotlibcpp.h"

int main()
{
    using namespace cpp_robotics;
    SQP solver;
    SQP::Problem prob;

    //////////////////// 問題設定 ////////////////////
    // 目的関数
    prob.func = [](Eigen::VectorXd x) -> double
    {
        return 100*( std::pow( (x(1) - std::pow(x(0),2)), 2) ) + std::pow(1 - x(1), 2);
    };
    
    // 制約1
    prob.con.push_back({
        Constraint::Ineq,
        [](Eigen::VectorXd x)
        {
            // |x| <= 1
            return (x(0)*x(0) + x(1)*x(1)) - 1;
        },
    });

    // 制約2
    prob.con.push_back({
        Constraint::Ineq,
        [](Eigen::VectorXd x)
        {
            // x1 >= 0
            return -x(0);
        },
    });

    // 解の初期値
    Eigen::VectorXd x0(2);
    x0 << 0, 0;

    //////////////////// 解く ////////////////////
    // 反復によって生成される点列の可視化用
    std::vector<double> x1_, x2_;

    // 解く
    auto result = solver.solve(prob, x0, [&](auto x){ x1_.push_back(x(0)); x2_.push_back(x(1)); });

    //////////////////// 結果表示 ////////////////////
    std::cout << "iter: " << result.iter_cnt << std::endl;
    std::cout << "x1 = " << result.x(0) << std::endl;
    std::cout << "x2 = " << result.x(1) << std::endl;
    std::cout << "x norm: " << result.x.norm() << std::endl;
    
    // 点列のグラフ化
    namespace plt = matplotlibcpp;
    plt::plot(x1_, x2_, "o--");
    // 補助線
    plt::plot({1,-1}, {1,1}, "xk"); // 最適解
    std::vector<double> edge_x(100), edge_y(100);
    for(size_t i = 0; i < 100; i++) // 制約によるxの取る範囲の境界
    {
        double theta = i/99.0 * 2 * 3.141592;
        edge_x[i] = std::max(0.0, std::cos(theta));
        edge_y[i] = std::sin(theta);
    }
    plt::plot(edge_x, edge_y, "--");

    plt::xlim(-1.1, 1.1);
    plt::ylim(-1.1, 1.1);
    plt::set_aspect(1.0);
    plt::show();
}

実行結果

実行するとこんな感じで解とグラフが表示されます。グラフは反復による解の収束を見ることが出来ます。

ターミナルの表示

iter: 19
x1 = 0.785066
x2 = 0.618036
x norm: 0.999149

青が反復によって生成された点列、オレンジが制約による実行可能領域の境界、黒のバツが制約なしの時の最適解です。 原点からスタートして最適解に近づきつつも制約を満たすところで止まってるのが分ります。実際にターミナルの出力を確認すると最適化した変数のノルムがほぼ1になっています。

最後に

今回は最適化手法の紹介と自前実装のソルバーを使って実際に最適化問題を解いてみました。 整数の最適化問題とか整数混合問題とかも勉強してみたいですね。

記事にミスがあった場合はコメントしてもらえるとありがたいです。

最後に宣伝ですが来年は学生ロボコン2023で後輩達がかっこいいロボット作って優勝してくれると思うので応援よろしくおねがいします!

Twitter twitter.com

サイト tutrobo.rm.me.tut.ac.jp