前章：モデル予測制御の考え方
MPCの本題に入る前に制約付きの非線形最適化問題を解く手法として主双対内点法を紹介いたします。
宣伝
拙著「PythonとCasADiで学ぶモデル予測制御 (KS理工学専門書) 」が発売されました！
是非こちらも読んでいただければ幸いです。
表記法について
ベクトルxの各成分x^{(j)}を対角に並べて拡張した行列をX \triangleq diag(x)、
全成分が1のベクトルをe \triangleq [1, 1, \dots]^Tとします。
制約付き最適化問題の標準形
一般に、x \in \mathbb{R}^n上での制約付きの最適化問題は以下の形で表せます。
x^* = \argmin_x f(x) \\
s.t.
\left\{
\begin{align*}
g_i(x) & = 0 ~~~~ i = 1,2,\dots \\
h_j(x) & \leq 0 ~~~~ j = 1,2,\dots
\end{align*}
\right. \tag{1}
ここで、補助変数として新たにx^+, x^-,s \geq 0を導入して（sはスラック変数と呼ばれます）
x = x^+ - x^- \\
h_j(x) + s_j = 0
と置きます。ここで改めてx \triangleq [s, x^-, x^+]と再定義して制約条件もまとめてg=[g^{(1)}, g^{(2)}, \dots]^Tと書き直すことで上記の問題は以下のように簡単な形式に帰着させられます。（nも新たなxの次元として再定義します）
x^* = \argmin_x f(x) \\
s.t.
\left\{
\begin{align*}
g(x) & = 0 \\
x & \geq 0
\end{align*}
\right. \tag{2}
バリア関数の導入
バリア関数法とは、x>0という不等式制約を無理やりコスト関数の方に押し付けてやる方法で、正の実数値パラメーター\muを用いて次のような問題に変形します。
x^* = \argmin_x \varphi_\mu(x) \\
\varphi_\mu(x) \triangleq f(x) - \mu \sum_{j=1}^n \log{(x^{(j)})} \\
s.t. ~~~~~ g(x) = 0 \tag{3}
最適化問題の解が存在する場合、上記のバリア問題(3)の解は\mu \rightarrow 0でオリジナルの問題(2)の解と一致します。
ここでラグランジュ関数L=\varphi_\mu- \sum_i \lambda^{(i)} g^{(i)}のx^{(j)}微分は
\partial_j L = \partial_j f - \mu \frac{1}{x^{(j)}} - \sum_i \lambda^{(i)} \partial_j g^{(i)}
となるため、z^{(j)} \triangleq \mu / x^{(j)}とすると、ラグランジュの未定乗数法より解の条件は以下のようになります。
\begin{align*}
\nabla f(x) + \lambda \nabla g(x) - z & = 0 \\
g(x) & = 0\\
XZe - \mu e & = 0 \tag{4}
\end{align*}
この時、\mu \rightarrow 0とともにx,z \geq 0が成立し、(2)のKKT条件と一致します。
この(4)の解を直線探索付きニュートン法を用いて求めます。
直線探索付きニュートン法によるバリア問題の解
直線探索付きニュートン法とは、例えばある関数F(x):\mathbb{R}^n \rightarrow \mathbb{R}を最小化するとき、すなわち\nabla F(x) = 0 を求めるときに次の方程式を解いて得られた方向ベクトル d_k \in \mathbb{R}
\nabla^2 F(x_k) d_k = -\nabla F(x_k) \tag{5}
を用いて、
x_{k+1} = x_k + \alpha_k d_k ~~~ (0 < \alpha_k \leq 1)
という形で反復的にx_kを更新していく方法です。（\alpha_k \equiv 1と固定したものを単にニュートン法と呼ぶこともあります。）
\alpha_kは各kステップにおいてさらに探索され最適なものが選択されます。この時、方向ベクトルd_k上で\alpha_kが探索されるため、「直線探索付き」という名前がついています。
この\alpha_kを求める方法としてはフィルター法が良く使われています。
バリア問題の話に戻ると
\nabla F =
\left(
\begin{gathered}
\nabla f(x_k) + \lambda_k \nabla g(x_k) - z_k\\
g(x_k) \\
X_k Z_k e - \mu e
\end{gathered}
\right) \tag{6}
を(5)に代入して
\begin{bmatrix}
W_k & A_k & -I\\
A_k^T & 0 & 0 \\
Z_k & 0 & X_k \\
\end{bmatrix}
\begin{pmatrix}
d^x_k\\
d^\lambda_k\\
d^x_k\\
\end{pmatrix}
= -
\left(
\begin{gathered}
\nabla f(x_k) + A_k \lambda_k - z_k\\
g(x_k) \\
X_k Z_k e - \mu e
\end{gathered}
\right) \tag{7}
が得られます。ただしA_k \triangleq \nabla g(x_k), ~~~ W_k \triangleq \nabla_{xx}^2 \mathcal{L}, ~~~ \mathcal{L} \triangleq f(x) + g(x)^T \lambda - zとして与えました。
最適化ステップ
ここまでの準備で(1)の最適化を行うことが可能になりました。
- 初期値x_0, \lambda_0, z_0とバリアパラメータ\mu>0を与える。
- (7)を解いて方向d_kを決める。
- フィルター法を用いて\alpha_kを与える
-
x_{k+1} = x_k + \alpha_k d^x_k, ~ \lambda_{k+1} = \lambda_k + \alpha_k d^\lambda_k,~ z_{k+1} = z_k + \alpha_k d^z_k で更新する。
-
\muを小さくする
- 2.に戻る
このステップを事前に決められた閾値を下回るまで繰り返すことで最適化を実行することが可能になります。
次の章では実際にMPCで簡単な制御を行ってみましょう。
次章：離散時間の非線形MPC
Discussion