Emileの備忘録

色々な事をだらだらと

制約付きの微分動的計画法(DDP)

この記事は学ロボアドベントカレンダー2024の最終日の記事です。


さて、2024年も年の瀬が迫って参りました。皆様いかがお過ごしでしょうか？
自分は最近夜の自転車が寒くて辛いです。


adventar.org



この記事では、制約付きの微分動的計画法について概説し、
pythonのjaxライブラリによる実装についても説明します。

また、この記事は一つ前の記事↓とセットになっています。
微分動的計画法(DDP)の基礎 - Emileの備忘録



==============目次==============

===============================



1.はじめに

DDPについて解説する前に、拡張ラグランジュ関数法について解説します。

2.拡張ラグランジュ関数法

拡張ラグランジュ関数法とは、
制約条件をペナルティとしてラグランジュ関数に加えた関数を用いて、
制約付き最適化問題を解く手法です。


元の最適化問題から拡張ラグランジュ関数を構成し、
制約無しの最適化問題を繰り返し解く事で制約付き最適化問題を解きます。



以下の制約付き最適化問題を考えます。

　　　 \displaystyle \min_{x \in \mathbb{R}^n}\ f(x) \\ \ \mathrm{s.t.}\ \ g_j(x) \leq 0 \ \ \ \ (j=1,2,..,m)

※ 今回は不等式制約のみの場合を扱います。(勿論等式制約も扱えます)


この最適化問題に対して、拡張ラグランジュ関数 L_\rhoは以下で定義されます。

　　　 \displaystyle L_\rho (x, \lambda) := f(x) + \sum_{i=1}^m \psi_\rho (g_i (x), \lambda_i)

ただし、 \psi_\rhoはペナルティ関数であり、

　　　 \psi_\rho (g_i(x), \lambda_i) := \left\{ \begin{array}{lll} \frac{\rho}{2} (g_i(x) + \lambda_i / \rho)^2 - \lambda_i^2 / (2 \rho) & (g_i(x) \leq - \frac{\lambda_i}{\rho}) \\ - \lambda_i^2 / (2\rho) & (\mathrm{otherwise}) \end{array} \right.

で定義されます。

突然出てきた \psi_\rhoについては、不等式制約をスラック変数で等式制約にしてから、
スラック変数についての最適化問題を解く事で求められます。


拡張ラグランジュ関数法の考え方としては、
\nabla L(x,\lambda) = 0を保ちながら、 g(x) \leq 0を満たす解を探すイメージです。
加えて、拡張ラグランジュ関数には制約充足を促す効果があります。(ラグランジュ関数との違いです)


拡張ラグランジュ関数法の手続きは以下のようになります。

　1. 初期化:　 x^0 \leftarrow 0,\ \lambda^0 \leftarrow 0,\ \rho^0 \leftarrow \rho_0
　2. 以下3~6を収束するまで繰り返す( k:反復のindex)
　3. xの更新:　 x^{k+1} \leftarrow \mathrm{argmin}_{x \in \mathbb{R}^n}\ L_{\rho^k} (x, \lambda^k)
　4. \lambdaの更新:　 \lambda^{k+1} \leftarrow \max\{ 0, \lambda_i^k + \rho^k g_i(x^{k+1})\}\ \ \ (i=1,...,m)
　5. \rhoの更新:　 \rho^{k+1} \leftarrow \beta \rho^k\ \ \ \ (\beta > 1)
　6. 収束したら反復を終了

シンプルな手続きで制約付き最適化問題を解くことができます。


3.制約付きの微分動的計画法

DDPは、制約無しの非線形最適制御問題の数値解法であり、
二次の間接法で最適入力を求める方法です。


詳細は前の記事を参照してください。(雑な感じでごめんなさい)
微分動的計画法(DDP)の基礎 - Emileの備忘録



離散時間確定システムの最適制御問題を考えます。

　　　 \displaystyle \min_{u_{0:K}}\ J(x_{0:K}, u_{0:K-1})
　　　 \displaystyle \mathrm{s.t.}\ \ x_{k+1} = f(x_k, u_k)
　　　 \displaystyle \ \ \ \ \ \ \ \ g_k(x_k, u_k) \leq \mathbf{0}

ただし、 x \in \mathbb{R}^p, u \in \mathbb{R}^qとし、 f: \mathbb{R}^p \times \mathbb{R}^q \rightarrow \mathbb{R}^pは十分に滑らかであり、
目的関数は以下で定義される2回微分可能な関数とします。

　　　 J(x_{0:K}, u_{0:K-1}) := \sum_{k=0}^{K-1} l(x_k, u_k) + \phi(x_k)


上記の最適制御問題を、拡張ラグランジュ関数を用いて以下のように書き換えます。

　　　 \displaystyle \min_{u_{0:K}}\ L_\rho (x_{0:K}, u_{0:K-1}, \lambda_{0:K-1}) = J(x_{0:K}, u_{0:K-1}) + \sum_{i,k} L_{\rho_{i,k}}(g_{i,k}(x_k, u_k), \lambda_{i,k})
　　　 \displaystyle \mathrm{s.t.}\ \ x_{k+1} = f(x_k, u_k)


ただし、DDPではコスト関数の二回微分可能性が必要になるので、
ペナルティ関数には、以下の \psi_\rhoを使用します。
(先程のだと有効制約が切り替わる点で二回微分可能性が満たされなくなります)

　　　 \displaystyle \psi_{\rho_{i,k}}(y_{i,k}, \lambda_{i,k}) = \frac{\lambda_{i,k}^2}{\rho_{i,k}} \phi \left( \frac{\rho_{i,k}}{\lambda_{i,k}} y_{i,k} \right)

　　　 \displaystyle \begin{equation} \phi(t):=\left\{\begin{array}{l} \frac{1}{2} t^2+t & (t \geq-\frac{1}{2}) \\ -\frac{1}{4} \log (-2 t)-\frac{3}{8} & (\text { otherwise }) \end{array}\right. \end{equation}

この関数は、上で説明した拡張ラグランジュ関数のペナルティを滑らかに近似した物に対応しています。


拡張ラグランジュ関数を導入した制約無しの最適制御問題に対して、
制約無しのDDPで最適化 → ラグランジュ乗数 \lambda_{0:K-1}及びペナルティ項 \rho_{0:K-1}の更新 を繰り返す事で、
制約を考慮した最適制御則を計算します。

DDPによる最適化の部分は通常のDDPと全く同じ手続きで計算できます。


4.jaxによる実装

こちらのリポジトリに実装をまとめておきました。
対向二輪車での静的障害物回避のexampleも実装してあります。


github.com


目標地点からの定常偏差が残っているような気がしますね...。
まぁ、そういう日もあります。


出典


www.kspub.co.jp


arxiv.org
↑ この論文自体は、内点法を使った制約付きDDPを主軸に提案している論文ですが、
面倒&一般的な実装が難しいので、拡張ラグランジュ法ベースのDDPの実装の方を紹介しています。

微分動的計画法(DDP)の基礎

この記事は学ロボアドベントカレンダー2024の24日目の記事です。


さて、2024年も年の瀬が迫って参りました。皆様いかがお過ごしでしょうか？
自分はあっという間に一年が終わってしまって、少し寂しい気持ちになっています。



adventar.org




この記事では、微分動的計画法の基礎について概説し、
pythonのjaxライブラリによる実装についても説明します。


本当は制約付きのDDPを一つの記事にまとめるはずだったのですが、
時間がなかったので長くなりそうだったので2つに分割しました。



==============目次==============

===============================



1.微分動的計画法(DDP)について

DDPは、制約無しの非線形最適制御問題の数値解法であり、
二次の間接法で最適入力を求める方法です。


離散時間確定システムの最適制御問題を考えます。

　　　 \displaystyle \min_{u_{0:K}}\ J(x_{0:K}, u_{0:K-1})
　　　 \displaystyle \mathrm{s.t.}\ \ x_{k+1} = f(x_k, u_k)

ただし、 x \in \mathbb{R}^p, u \in \mathbb{R}^qとし、 f: \mathbb{R}^p \times \mathbb{R}^q \rightarrow \mathbb{R}^pは十分に滑らかであり、
目的関数は以下で定義される2回微分可能な関数とします。

　　　 J(x_{0:K}, u_{0:K-1}) := \sum_{k=0}^{K-1} l(x_k, u_k) + \phi(x_k)



DDPは、動的計画法に基づいた数値解法であり、
Q関数を二次近似してベルマン方程式中の最適化問題を解析的に解く事で、
最適制御問題を近似的に解く方法です。


ここからはDDPの更新則を導出していきます。
動的計画法の構造から、 k = K, K-1, ...と終端から逆向きに解いていく事に注意して下さい。

まず、上記の最適制御問題に対応する状態価値関数 V_kは、

　　　 \displaystyle V_k(x_k) := \min_{u_{k:K-1}} \sum_{\bar k = k}^{K-1} l(x_{\bar k}, u_{\bar k}) + \phi (x_K)

で与えられ、ベルマン方程式

　　　 \displaystyle V_k(x_k) = \min_{u_k} Q_k(x_k, u_k)
　　　 \displaystyle Q_k(x_k, u_k) = l(x_k, u_k) + V_{k+1} (x_{k+1})

を満たします。ただし、 V_K(x_K) = \phi(x_K)となる事に注意してください。


ここで、与えられたノミナル軌道 \bar{x}_{0:K}, \bar{u}_{0:K-1}に対して、
x_k = \bar x_k + \delta x_k及び、 u_k = \bar{u}_k + \delta u_kを置くと、 Q_kは以下で二次近似できます。

　　　 \displaystyle Q_k\left(x_k, u_k\right) \approx Q_{x, k}^{\mathrm{T}} \delta x_k+Q_{u, k}^{\mathrm{T}} \delta u_k+\frac{1}{2} \delta x_k^{\mathrm{T}} Q_{x x, k} \delta x_k + \delta x_k Q_{x u, k} \delta u_k+\frac{1}{2} \delta u_k^{\mathrm{T}} Q_{u u, k} \delta u_k

ただし、

　　　 \displaystyle \begin{align} Q_{x x, k} & =\ell_{x x}+f_{x}^{\mathrm{T}} V_{x x, k+1} f_{x}+V_{x, k+1} \cdot f_{x x} \\ Q_{x u, k} & =\ell_{x u}+f_{x}^{\mathrm{T}} V_{x x, k+1} f_{u}+V_{x, k+1} \cdot f_{x u} \\ Q_{u u, k} & =\ell_{u u}+f_{u}^{\mathrm{T}} V_{x x, k+1} f_{u}+V_{x, k+1} \cdot f_{u u} \\ Q_{x, k} & =\ell_{x}+f_{x}^{\mathrm{T}} V_{x, k+1} \\ Q_{u, k} & =\ell_{u}+f_{u}^{\mathrm{T}} V_{x, k+1} \end{align}

とします。ここで、 l, f, V等の微分は、 \bar{x}, \bar uで評価しています。
また、 Q_{uu}は Q_{uu} \succ Oとなる様に正則化する必要があります。(単位行列の定数倍を加えて正定行列にします)


ベルマン方程式の \displaystyle V_k(x_k) = \min_{u_k} Q_k(x_k, u_k)について、二次近似を元に右辺を計算すると、

　　　 \displaystyle \delta u^*_k = k_k + K_k \delta x_k
　　　 \displaystyle k_k = -Q^{-1}_{uu,k} Q_{u,k},\ \ K_k = - Q_{uu,k}^{-1} Q_{ux,k}

として近似的な最適制御則が求まります。(ゲインが求まる)
具体的な \delta x_kや \delta u_kは未知な事に注意してください。


更に、求まった制御則を代入して V_kも二次近似する事で、

　　　 \displaystyle \begin{align} V_k & =-\frac{1}{2} Q_{u, k}^{\mathrm{T}} Q_{u u, k}^{-1} Q_{u, k} \\ V_{x, k} & = Q_{x, k}-Q_{u x, k}^{\mathrm{T}} Q_{u u, k}^{-1} Q_{u, k} \\ V_{x x, k} & = Q_{x x, k}-Q_{{u x}, k}^{\mathrm{T}} Q_{u u, k}^{-1} Q_{{u x}, k} \end{align}

として、 V_k, V_x, V_{xx}の更新則も求まります。
ただし、終端時刻 Kでは、 V_K = \phi(x_K),\ V_{x,k} = \phi_x(x_K),\ V_{xx,k} = \phi_{xx}(x_K)となります。


後退パスで、全ての K_K及び k_kを計算したら、次は前方パスで新しい軌道を計算します。

　　　 \displaystyle u_k = \bar u_k + k_k + K_k (x_k - \bar x_k)
　　　 x_{k+1} = f(x_k, u_k),\ \ \ x_0: \mathrm{given}

求まった新しい軌道は、次の反復でノミナル軌道 \bar{x}_{0:K}, \bar{u}_{0:K-1}として使用されます。


以上の手続きを終了条件が満たされるまで繰り返す事で、最適入力を計算する方法をDDPと呼びます。

全体像をまとめておくと、終端から順にゲイン k_k, K_kを求め (backward)、
順方向のrollout (forward)でノミナル軌道を更新する手続きを収束するまで繰り返すという構造になっています。


単純な計算のみから構成されるbackward → forwardを繰り返すだけなので、
非常に高速に制御則を計算できるというメリットがあります。



2.jaxによる実装


こちらのリポジトリに実装をまとめておきました。exampleにLQRとの比較コードが載せてあります。
(理想的には軌道が一致するはずの設定ですが、数値的な要因で少しズレていると考えられます。)


github.com


正直、特筆する様な事は無いのですが、
jitの恩恵を最大限得られるように、forやwhile、条件分岐は出来る限りlaxを使うのがコツだと思っています。


DDPはjaxの練習に丁度よい題材だと思うので、jaxを勉強したい方は是非実装してみて下さい!





雑多な感想

忙しい。jax楽しい。最近寒い。

ブラックボックス(ベイズ)最適化とハイパラの自動調整

この記事は学ロボアドベントカレンダー2023最終日の記事です。
学ロボアドカレ2023、多数のご参加ありがとうございました!!



さて、2023年も年の瀬が迫って参りました。
皆様いかがお過ごしでしょうか？
私はこの年末が学部生活の冬という現実に戦々恐々としています。

adventar.org




この記事では、ハイパーパラメータの自動調整に用いられるblackbox最適化を概説します。
blackbox最適化の中でも特にベイズ最適化を中心に解説し、optunaを用いた実装についても解説します。


環境としては、python 3.11.6を想定しています。




==============目次==============

===============================




1.ブラックボックス最適化とは

式で陽に表現できない関数をblackbox関数と言い、
blackbox関数に対する最適化問題をblackbox最適化と言います。

以降、blackbox関数を f: \mathcal{X} \rightarrow \mathbb{R}と表す事にします。


特にblackbox最適化では以下の様な状況を考える事が多いです。

・ \nabla f(x), \nabla^2 f(x)が観測できない
・ f(x)の評価コストが高い
・ f(x)の観測にノイズが乗る ( y = f(x)+\epsilonが観測される)


勾配やヘッセ行列などが得られないので、一般的な最適化手法を適用する事は出来ません。
blackbox最適化では、この様な関数 fに対して f(x)を最小にする x \in \mathcal{X}を見つける事を考えます。



(例) blackbox最適化の例としては、何らかの問題のハイパーパラメータの自動調整が挙げられます。
この記事の後半では、ハイパラの自動調整に広く用いられるoptunaで実際にblackbox最適化を実装します。





1-1. Blackbox最適化の手法

主なBlackbox最適化の手法として、以下が挙げられます。

・ランダムサーチ/グリッドサーチ
文字通り、予め決めたルール/事前情報で総当りします。
これで問題無いなら一番手っ取り早いです。


・ベイズ最適化
(この記事で主に解説する内容です。)


・進化計算
入力 x\in \mathcal{X}を個体とした進化計算によってblackbox最適化を行います。
GAやCMA-ESといった手法がよく使われます。

CMA-ESについてはこちらの記事の解説が分かり良いです。

horomary.hatenablog.com






2. ベイズ最適化

この記事では、特にblackbox最適化の中でも広く用いられているベイズ最適化を扱います。
ベイズ最適化では、blackbox関数 fに統計モデルを仮定し、 fを分布からのサンプルとして捉えます。


具体的には、以下の手続きで fを最小化する点 xを求めます。

1. 観測データ \mathcal{D}_n := \{(x_i, y_i)\}_{i=1}^nから統計モデルを更新
2. 統計モデルから獲得関数を求め、獲得関数を最大化する点 x'を次の入力として選択
3. 選択した点 x'を入力して、評価値 y'を観測し、 \mathcal{D}_{n+1} \leftarrow \mathcal{D}_n \cup \{(x', y')\}とする。


統計モデルは入力に対する fの分布の近似に使用し、
獲得関数は探索と活用のトレードオフを考慮した入力を定めるのに使用します。


強化学習などでもおなじみ探索と活用のトレードオフはここでも登場します。
勿論UCB系のアプローチも使われます。(今回はminimizeなのでLCB)






3. ベイズ最適化の手法

ベイズ最適化には、統計モデルと獲得関数の組み合わせで様々な種類の手法が存在します。
ここでは、ベイズ最適化でよく使われる手法を統計モデル/獲得関数のそれぞれで列挙します。


3-1. 統計モデル

まず統計モデルから概説します。
統計モデルを用いてblackbox関数 fをモデル化します。

ガウス過程回帰(GPR)

ガウス過程(GP)というのは、入力空間上のランダムな関数を定める確率過程です。

出力の平均値を定める関数 m_n: \mathcal{X} \rightarrow \mathbb{R}とカーネル関数 k_n: \mathcal{X}\times \mathcal{X} \rightarrow \mathbb{R}から構成され、
入力系列 x_1,...,x_n \in \mathcal{X}に対する、関数 g: \mathcal{X} \rightarrow \mathbb{R}の出力系列が
\mathbf{g}:= (g(x_1),.., g(x_n)) \sim \mathcal{N}(\mathbf{m}, \mathbf{K})となる時、関数 gはガウス過程に従うと言います。
ただし、 \mathbf{m}_i = m_n(x_i)、 \mathbf{K}_{ij} = k_n(x_i, x_j)とします。


以降では、 \mathcal{D}_nが得られている時にGPの平均を表す関数を m_n(x)、分散を表す関数を k_n(x, x)として、
関数 gがガウス過程に従う事を g \sim \mathcal{GP}(m_n, k_n)と表す事にします。

更に、観測 yに乗るノイズの分散を \sigma^2として、
\sigma_n^2(x) := k_n(x, x) + \sigma^2を定義します。( y \sim \mathcal{N}(m_n(x), \sigma^2_n(x)))


ガウス過程回帰(GPR)は、ガウス過程と正規分布の再生性・同時分布の性質を用いて、
ガウス過程の m_n, k_nをオンラインで更新する事で関数回帰を行う手法です。

GPRでは、 fの事前分布を \mathcal{GP}(m_0, k_0)とし、観測ノイズ \epsilon \sim \mathcal{N}(0, \sigma^2)とすると、
\mathcal{D}_n = \{(x_i, y_i)\}_{i=1}^nが観測された場合の事後分布は以下の様に表されます。


f \sim \mathcal{GP}(m_n, k_n)
・ m_n(x) = m_0(x) + \mathbf{k}_n(x)^\top (\mathbf{K}_n + \sigma^2 I_n)^{-1} ( \mathbf{y}_n - m_0(x))
・ k_n(x, x') = k_0(x, x') - \mathbf{k}_n(x)^\top (\mathbf{K}_n + \sigma^2 I)^{-1} \mathbf{k}_n(x)

ただし、 \mathbf{k}_n (x) = (k_n(x_1, x),..., k_n(x_n, x))^\topとし (\mathbf{K}_{n})_{ij} = k_n(x_i, x_j)としています。


上記の様にガウス過程としてモデル化する事で、回帰する対象の分散(不確実さ)も表現出来るのがメリットです。
ただし、観測した点数 nについて更新にかかる計算量が O(n^3)と計算量が大きいです。




TPE(Tree-structured Parzen Estimator)

TPEでは、 P(x|y, \mathcal{D}_n)のみをモデル化して P(y|x, \mathcal{D}_n)を求めます。
獲得関数にはEI(詳細は後述)を用います。


具体的には、
P(x|y, \mathcal{D}_n) := \left\{ \begin{array}{lll} l(x | \mathcal{D}_n) & (y < y^*) \\ g(x|\mathcal{D}_n) & (y \geq y^*) \end{array} \right.としてモデル化し、
関数 l(x), g(x)は \mathcal{D}_nを用いたParzen推定(カーネル密度推定法の一種)で求めます。
P(y)は陽にモデル化しない事に注意して下さい。


ただし、 y^*は P(y < y^*) = \gammaとなる様に与えた閾値とします。( \mathcal{D}_nの上位 \gammaとなる yとか)
optunaでは、 \gamma = 0.25がデフォルト値として採用されています。

yが y^*よりも小さいか否かを考えたいので、この様に分割しています。


このモデルを用いたEIを計算すると以下の様になります。
(TPEは獲得関数までセットなので、ここで解説します)

\alpha_{EI}(x) = \int_{-\infty}^{y^*} (y^* -y)P(y | x,\mathcal{D}_n) dy
　　　　 = \int_{-\infty}^{y^*} (y^* - y)\frac{l(x, \mathcal{D}_n) P(y | \mathcal{D}_n)}{\gamma l(x | \mathcal{D}_n) + (1-\gamma)g(x|\mathcal{D}_n)}dy
　　　　 = \left( \gamma + \frac{g(x | \mathcal{D}_n)}{ l(x|\mathcal{D}_n) }(1-\gamma) \right)^{-1} \int_{-\infty}^{y^*} (y^* - y) P(y | \mathcal{D}_n) dy
　　　　 \propto \left( \gamma + \frac{g(x | \mathcal{D}_n)}{ l(x|\mathcal{D}_n) }(1-\gamma) \right)^{-1}


※ \int_{-\infty}^{y^*} P(y|\mathcal{D}_n)dy = \gammaに注意。


以上の結果から、 \alpha_{EI}(x)を最大化するには \frac{g(x|\mathcal{D}_n)}{l(x|\mathcal{D}_n)}を小さくする点 xを考えれば良いと分かります。


特に l(x|\mathcal{D}_n)の定義から、以下の手続きで獲得関数を最大化する点が求まります。
1. x \sim l(x|\mathcal{D}_n)を複数サンプリングする
2. サンプリングした xから、最も \frac{g(x | \mathcal{D}_n)}{l(x | \mathcal{D}_n)}を小さくする点 x \in \mathcal{X}を次の入力として採用

以上がEIとTPEを用いたベイズ最適化法の手続きとなります。


ガウス過程回帰と比べて計算量が小さく、観測した点数が多い場合にも使えるのがメリットです。
Optunaのdefaultの最適化アルゴリズムです。


こちらの記事が分かりやすいです。
gochikika.ntt.com





3-2. 獲得関数

獲得関数についてもまとめておきます。
統計モデルから求めた獲得関数を最大化する点 xを次の入力に用います。

EI(Expected Improvement)

EI(期待改善量)はblackbox最適化で用いられる最も代表的な獲得関数です。
単純かつ比較的高速に動作するので、手っ取り早く試すのに向いている手法です。

EIの獲得関数は、閾値 y^*に対して \alpha_{EI}(x) := \mathbb{E}[\max (0, y^* - y)]として定義されます。


y^*には様々な決め方が存在しますが、 y_1,..., y_nを \mathcal{D}_nに含まれる観測値として、
y^* := \min \{y_1,...,y_n\}と定める方法が用いられたりします。



LCB(Lower Confidence Bound)

多碗bandit問題や強化学習でもおなじみの考え方です。

LCBは信頼区間の下限を用いた獲得関数であり、
GPRなど明示的に分散も求まる統計モデルを用いて以下の様に定義されます。

\alpha_{LCB}(x) := - (m_n(x) - \sqrt{\beta_n} \sigma_n (x))
ただし、 \beta_nはハイパーパラメータで \beta_n = c \log n\ \ (c:const.)などとして定めます。
(最大化問題を考える場合はUCBを考えれば良いです。)

LCBは"不確かな時は楽観的に"というヒューリスティックに基づいていて、
標準偏差を期待値に足し合わせた評価値で探索と活用のトレードオフを取ります。


直感的かつ単純な実装で実現できる事や、リグレット等の理論解析のしやすさがメリットとなりますが、
実用上性能がそこまで良くない事も実験的に知られています。



Thompson Sampling

Thompson Samplingは、目的関数の事後分布から関数を一つサンプリングし、
その関数を最小化する点を次の入力として採用する手法です。

こちらもLCB/UCBと合わせて多碗バンディット問題でよく使われます。
※ blackbox最適化は連続腕bandit問題とも解釈出来るので、banditの手法が幾つか登場します。


ガウス過程を用いる場合、Thompson Samplingの獲得関数は以下の様に定義されます。
\alpha_{TS}(x) = - f^{(n)}(x),\ \ \ \ (f^{(n)} \sim \mathcal{GP}(m_n, k_n))


f^{(n)}は確率的に定まる関数であり、直接 \alpha_{TS}(x)を最大化するのは難しいので、
今回は \mathcal{X}を有限個に分割して近似的に f^n(x)を計算します。


まず、 \mathcal{X} = \sum_{i=1}^\Delta I_iとして互いに素な集合族 \{I_i\}_{i=1}^\Deltaを取り、
各集合 I_iから代表点 p_i \in I_i \subset \mathcal{X}を取ります。
この代表点を用いて、 \mathcal{X} \simeq \{p_i\}_{i=1}^\Deltaとして有限個の点で \mathcal{X}を近似します。

さらに、 \mathcal{GP}(m_n, k_n)を用いて p_iに対応する q_i := f(p_i) \sim \mathcal{GP}(m_n, k_n)をサンプルし、
\{(p_i, q_i) \}_{i=1}^\Deltaによって、関数 f^{(n)}を f^{(n)} \simeq \hat{f}^{(n)}: p_i \mapsto q_iとして有限個の点で近似します。


これらの近似によって、以下の単純な計算で \alpha_{TS}(x)を最大化する点を近似的に求められる様になります。
\mathrm{argmax}_{x \in \mathcal{X}}\ \alpha_{TS}(x) \simeq \mathrm{argmax}_{p_i \in \{x_1,...,x_\Delta\}} \left( - \hat{f}^{(n)}(p_i) \right) = \mathrm{argmin}_{i = \{1,...,\Delta\}}\ q_i

勿論この計算法では予め集合 \mathcal{X}の分割法を定める必要があります。


※ \alpha_{TS}(x)を最大化する xは、 fを最小化する点の事後分布 p_{min}(x^* | \mathcal{D}_n)に従う点としても解釈できます。
この観点から捉えると、banditでのThompson Samplingとの関連性が分かりやすいと思います。



Thompson Samplingは多碗bandit問題では有力な手法なのですが、
ベイズ最適化ではそれほど使われていない様です。(高次元で性能が悪い)



※ banditの手法の比較は以下で試せます。(気になる方は試してみて下さい)
github.com




(補足)獲得関数の最小化について

特にGPRを用いた獲得関数の最小化について補足します。
先程定義した m_n(x)、 \sigma_n(x) = k_n(x, x) + \sigma^2を考えます。


これらの関数の x^{(i)}についての勾配は以下の様に求まります。( m_0(x) =0としています)
\frac{\partial m_n (x)}{\partial x^{(i)}} = \left( \frac{\partial \mathbf{k}_n (x)}{\partial x^{(i)}} \right)^\top (\mathbf{K}_n + \sigma^2 I_n)^{-1} \mathbf{y}_n
\frac{\partial \sigma_n(x)}{\partial x^{(i)}} = \frac{1}{2\sigma_n (x)} \left( \frac{\partial k_0(x, x)}{\partial x^{(i)}} - 2 \left( \frac{\partial \mathbf{k}_n(x)}{\partial x^{(i)}} \right)^\top (\mathbf{K}_n + \sigma^2 I_n)^{-1} \mathbf{k}_n(x) \right)


これらの勾配は \mathbf{K}_n + \sigma^2 I_nを事前計算しておく事で、計算量を O(n^2)まで削減できます。

EIやLCBの獲得関数の勾配は m_n(x),\ \sigma_n(x)の勾配を用いて表現されるので、
上記の行列を事前計算しておく事で、獲得関数の最適化計算の計算量を削減できます。

今回の実装では事前計算を行っておりません






4. 実装(GP-LCB)

GPR+LCBを用いた場合のプログラムを実装してみます。
この手法はGP-LCB(GP-UCB)と呼ばれ、連続腕バンディット問題の方策の1つです。


ここからは、 f(x) = \exp(x/2) + \sin(2\pi x)を最小化する xをblackbox最適化で求める事を考えます。
ただし \mathcal{X} = [-1.0, 1.0]とします。


GPRはscikit-learnのGaussianProcessRegressorを使用し、
kernel関数にはmatern kernelにwhite kernel、Constant kernelを組み合わせた物を用います。

matern kernelの定義は以下を参照して下さい。
sklearn.gaussian_process.kernels.Matern — scikit-learn 1.3.2 documentation


また、LCBの \beta_nには \beta_n := 2.0 * \log nを使用し、
f(x)の観測ノイズ \epsilonは \epsilon \sim \mathcal{N}(0, 0.01)とします。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, WhiteKernel, ConstantKernel
from scipy.optimize import minimize
from functools import partial


def f(x: np.ndarray) -> np.ndarray:
    x = x.reshape(-1)
    return np.exp(x/2.0) * np.sin(2 * np.pi * x)


def obs(value: np.ndarray, noise: bool | float = False) -> np.ndarray:
    return value + np.random.normal(0, noise, size=value.shape) if noise else value


def plot_data(xs: np.ndarray, ys: np.ndarray, x_bound: np.ndarray, gpr: GaussianProcessRegressor, n) -> None:
    x = np.linspace(x_bound[:, 0], x_bound[:, 1], 200).reshape(-1)
    y, sigma = gpr.predict(x.reshape(-1, 1), return_std=True)
    plt.plot(x, f(x), label='f(blackbox)')
    plt.plot(x, y, label='GPR(model)')
    plt.fill_between(x, y - sigma, y + sigma, alpha=0.2, label='sigma')
    plt.fill_between(x, y - sigma * np.sqrt(2.0 * np.log(n)), y + sigma * np.sqrt(2.0 * np.log(n)), alpha=0.2, label='UCB/LCB')
    plt.scatter(xs, ys, marker='x', color='red')
    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.legend()
    plt.show()


def lcb_policy(x_bound: np.ndarray, gpr: GaussianProcessRegressor, n: int) -> np.ndarray:
    def lcb_score(x: np.ndarray, gpr: GaussianProcessRegressor, n: int) -> float:
        mean, sigma = gpr.predict(x.reshape(1, -1), return_std=True)
        score = (mean - sigma * np.sqrt(2.0 * np.log(n)))[0]
        return score

    objective = partial(lcb_score, gpr=gpr, n=n)
    point_num = 30
    x_next_candidates: list = []

    for x in np.random.uniform(x_bound[:, 0], x_bound[:, 1], size=(point_num, x_bound.shape[0])):
        # 多点スタートで最適化 (local minimum対策)
        ans = minimize(objective, x0=x, bounds=x_bound)
        x_next_candidates.append([ans.fun, ans.x])

    x_next = min(x_next_candidates)[1]
    return x_next


def GP_LCB(x_bound: np.ndarray, trial_num: int):
    kernel = Matern(nu=1.5) + WhiteKernel() + ConstantKernel(1.0)  # カーネルの設定
    gpr = GaussianProcessRegressor(kernel=kernel, alpha=0, random_state=0)  # モデルの作成

    # observed data
    xs: np.ndarray = np.empty(0)
    ys: np.ndarray = np.empty(0)

    # 初期値の設定
    x_next: np.ndarray = np.random.uniform(x_bound[:, 0], x_bound[:, 1], size=(x_bound.shape[0],))

    for i in range(trial_num):
        print(f"trial: {i+1}")
        y_next = obs(f(x_next), noise=0.01)  # 観測値の取得
        xs = np.append(xs, x_next).reshape(-1, x_next.shape[0])
        ys = np.append(ys, y_next)
        gpr.fit(xs, ys)  # モデルの学習
        plot_data(xs, ys, x_bound, gpr, i+1)
        x_next = lcb_policy(x_bound, gpr, i + 1)

    return min(ys)


if __name__ == '__main__':
    x_bound_array = np.array([[-1.0, 1.0]])
    GP_LCB(x_bound_array, 20)

※ GPRのkernelの選び方などGP-LCB自体にもハイパラが多く、上手く最小化出来ない時もあります。


最適化計算等めんどうな部分はscikit-learnのモジュールで誤魔化しています。

本来GPRを用いた場合の獲得関数の最適化計算には逆行列の事前計算テクが使えるのですが、
今回は省略させて頂きました。


以下がプログラムの実行結果となります。
青の帯が分散を表し、オレンジの帯の上端がUCBを下端がLCBを表しています。

赤のxが入力した点 xを表しています。

反復に従って、 fを最小化する点が求まっている事が分かります。





5. optunaを用いた最適化

・optuna公式
optuna.org


5-1. install

以下のコマンドでinstall出来ます。

pip install optuna


dashboardで可視化したい人は以下もinstallして下さい。

pip install optuna-dashboard

installで躓く事は無いと思います。



5-2. toy-modelに対する最適化(sample code)

簡単な関数に対して、optunaを用いたblackbox最適化を実装してみます。


先ずは1変数で考えます。
最適化したい変数をtrial objectで記述し、以下のcodeで最適化します。

import optuna

def objective(trial):
    x = trial.suggest_float('x', -10, 10)  # xの範囲を[-10, 10]とする
    return (x - 2) ** 2

study = optuna.create_study()  # problem setting
study.optimize(objective, n_trials=100)  # blackbox最適化の実行
study.best_params  # 結果

結果はstudy objectのパラメータに格納されます。
optuna.study.Study — Optuna 3.5.0 documentation


また、optunaでは以下の様にアルゴリズムが選択出来ます。
Efficient Optimization Algorithms — Optuna 3.5.0 documentation

デフォルトではTPEを使用します。



2変数になっても同じ様に実装出来ます。

import optuna

def func(x, y):
    return x**2+y**2+x*y

def objective(trial):
    x = trial.suggest_float('x', -1, 1)  # xの範囲指定
    y = trial.suggest_float('y', -1, 1)  # yの範囲指定
    return func(x, y)

study = optuna.create_study()
study.optimize(objective, n_trials=100)
study.best_params

簡単にblackbox最適化を試す事が出来ます。



5-3. 結果の可視化

以下のコードで結果を可視化出来ます。

# Empirical Distribution Function Plot
fig = optuna.visualization.plot_edf(study)
fig.write_image('edf.png') 

# Optimization History Plot
fig = optuna.visualization.plot_optimization_history(study)
fig.write_image('optimization_history.png')

# parallel coordinate plot
fig = optuna.visualization.plot_parallel_coordinate(study)
fig.write_image('parallel_coordinate.png')

# slice
fig = optuna.visualization.plot_slice(study)
fig.write_image('slice.png')


こちらの記事に詳しく説明されています。(コードもこちらの記事から引用しました。)
ibetoge.hatenablog.com




5-4. マルチプロセス化

optunaはRDBを用いたマルチプロセス化にも対応しています。
今回はSQLiteを用いて実装します。


以下のコマンドでデータベースを生成する事が出来ます。(pythonからでも生成出来ます)

optuna create-study --study 'optuna_test' --storage 'sqlite:///optuna_testdb.db'

pythonでも生成出来ます。

DATABASE_URI = 'sqlite:///optuna_testdb.db'
study_name = 'optuna_test'
optuna.create_study(study_name=study_name, storage=DATABASE_URI)

データベースが生成できたら、次はmulti-processでの実行を実装します。
今回はmultiprocessingを用いて実装します。


multi-processでoptunaを実行するコードは以下の様になります。
簡単かつ直感的に実装出来る様になっています。

import optuna
import multiprocessing
from multiprocessing import Process


def func(x, y):
    return x ** 2 + y ** 2 + x * y


def objective(trial):
    x = trial.suggest_float('x', -1, 1)  # xの範囲指定
    y = trial.suggest_float('y', -1, 1)  # yの範囲指定
    return func(x, y)


def optimize(study_name, storage, n_trials):
    study = optuna.create_study(study_name=study_name, storage=storage, load_if_exists=True)
    study.optimize(objective, n_trials=n_trials)


n_trials = 10
concurrency = multiprocessing.cpu_count()  # CPUの数を取得
n_trials_per_cpu = n_trials / concurrency

DATABASE_URI = 'sqlite:///optuna_testdb.db'
study_name = 'optuna_test'
optuna.create_study(study_name=study_name, storage=DATABASE_URI)  # dbの生成

# multi-process化
workers = [Process(target=optimize, args=(study_name, DATABASE_URI, n_trials_per_cpu)) for _ in range(concurrency)]

for worker in workers:
    worker.start()

for worker in workers:
    worker.join()

study = optuna.load_study(study_name=study_name, storage=DATABASE_URI)  # 結果の取得

fig = optuna.visualization.plot_edf(study)
fig.write_image('edf.png')

# Optimization History Plot
fig = optuna.visualization.plot_optimization_history(study)
fig.write_image('optimization_history.png')

# parallel coordinate plot
fig = optuna.visualization.plot_parallel_coordinate(study)
fig.write_image('parallel_coordinate.png')

# slice
fig = optuna.visualization.plot_slice(study)
fig.write_image('slice.png')


公式Documentにも記述があります。
Easy Parallelization — Optuna 3.5.0 documentation





6. 終わりに

ベイズ最適化以外の手法やアルゴリズムの使い分けにも触れようと思っていたのですが、
カレンダーの日にちまでに間に合いませんでした。

申し訳ないです...。






参考文献

この記事を書くに当たって参考にした解説/本を列挙しておきます。


・blackbox最適化の幅広いサーベイ
search.ieice.org


・ベイズ最適化
www.kindaikagaku.co.jp


・ガウス過程回帰
www.jstage.jst.go.jp


・bandit問題
www.kspub.co.jp

キャチロボ23まとめ(新規開発物とその知見など)

この記事はポエムです。ｱﾚﾙｷﾞｰのある方は閉じてください。


この記事は学ロボアドベントカレンダー18日目の記事です。
adventar.org

キャチロボ2023で新規開発した物・知見やその他もろもろについてまとめた記事です。


==============目次==============

===============================


1. チームについて

キャチロボ2023では、
「調べてみたけど分かりませんでした！いかがでしたか？」チームとして参加させて頂きました。


メンバー構成は
・ソフト1 (私)
・ハード2 (@m11_Flexxia, @qzrobo_92)
・回路1 (twitter等に居ない)
といった構成で、学ロボ2021・22とデスマにデスマを重ねてきたメンバーです。


当日大事故を起こし、散々な結果ではありましたが楽しめました。
はいちゅーんずさん強かったです! 準優勝おめでとうございます!!




1.1. ハードウェア

詳細は以下のツイートをご覧ください。


結果はともあれ、高い制御性とメンテナンス性を誇る非常に良い機体でした。
(制御の無茶な注文に完璧に応えた機体です。ハードの二人に頭が上がりません....)

開発当初からM3508の強力なトルクを存分に活かせる様に設計されています。
十二分な剛性の確保やバックラッシの低減など細部まで考え抜かれています。
好き放題遊ばせてもらいました。本当に良い機体です!!


対して、まだまだ制御は詰めれたと思います。
ソフトウェアについては2.にまとめました。




2. 新規開発/導入した物まとめ

今回新規開発/導入した物は以下になります。
・ROS2
・micro-ROS
・KONDOサーボ
・ロボマスモーター導入(C620、M3508)
　・部内向けCANライブラリ作成
・データ駆動型制御(FRIT)
・リアルタイム経路生成

ソフトが自分一人というのもあって好き放題動けました。




2-1. 知見

・ROS2
ROS1を知っていれば、2~3日で使える様になると思います。
分からない時は実装読むのが早いかも知れません。

Serial通信には以下のpackageを使いました。
github.com

callbackのデッドロックに注意!
Using Callback Groups — ROS 2 Documentation: Humble documentation





・micro-ROS
hexagon-emile.hatenablog.com

意外と情報が少ないです。環境構築で沼に嵌りました。
FreeRTOSと合わせた時のビルドがなかなか通らず苦戦しました。(特にFPU周り)

F767zi + FreeRTOSで動かしている感覚では、Topicは40Hzが限界という感じでしたが果たして...?
メモリの多いH7を使った方が無難かもしれません。

一度環境構築を突破すれば、簡単にROS<->マイコン間通信が出来る便利なミドルウェアです。





・KONDOサーボ
高トルクで感動!! Amazon謎サーボと大違い!!!
配線を間違えて1つ破壊しました。ごめんなさい....

コマンド送信用packageも作りましたが自前で実装した方が早いと思います。
github.com

Hardware Interfaceで実装し、ROS2 controlで実装するという選択肢もありますが、
ROS2 controlでシステムを実装するのが面倒だったので自前のinterfaceを用意しました。

アクチュエーター制御を全てマイコン側にまとめるのが一番綺麗なシステム構成になると思います。





・ロボマスモーター
C620にトルク指令入れると良い感じに制御してくれます。
ブラシレスのMDを保持/開発していないチームにとっては現状最良の選択肢だと考えています。

トルク制御が掛かるので、位置制御をしたい時は2段階でPIDをかけると手っ取り早く動きます。
(位置制御(PID)で目標速度を出力、目標速度から速度制御(PID)で目標トルクを決定)
次元を考えると位置->速度->トルクで上手く行くのは直感的だと思います。


アンチワインドアップかけるのも良いかもです。
↓分かりやすく、参考になる記事です。
hamachannel.hatenablog.com





・データ駆動型制御(FRIT)
hexagon-emile.hatenablog.com

最近流行っているデータ駆動型制御です。

閉ループ系を安定化できるPIDパラメータを見つけて、
それで制御したデータを1組回収すれば、データを取り直す事無くパラメータの自動調整が出来ます。

分かりやすいデータは準備出来ていないのですが、
割と簡単に実装できて、それなりに良いPIDパラメータが得られます。


FRITと対になってよく使われるVRFTについては、naokichiさんが記事を書いて下さっています。
zenn.dev

FRITは閉ループ系で調整できますが、
系を安定化する初期パラを見つける必要があり、その上非線形最適化となります。(大域最適性が保証されない)


VRFTは二次計画で解ける上に事前調整も要らないですが、
開ループ系なので不安定系には使いづらいという特徴もあります。
(入力に周波数帯の広い信号を入れた方が良い事からも系を選ぶ手法だとは思います)

対象と目的に合わせて適宜使い分ける事で私達をパラチュン地獄から解放してくれるはずです。






・リアルタイム経路生成
hexagon-emile.hatenablog.com

以前に垂直多関節を扱った事があったので、それと比べると遥かに楽でした。
簡単な実装の組み合わせで手っ取り早くそれなりに良い経路が作れます。

ompl自体もモジュールとして扱いやすく、学ロボの経路生成でも使えるかも知れません。





3. チーム内の情報共有について

メンバーの一部がキャンパスが異なる・院試の時期がバラバラと言った背景もあり、
今回のプロジェクトでは特に情報共有のコスト削減に注力していました。

そういった事情もあり、今回は1つのサービスに情報を集約できるNotionというサービスを導入しました。
タスクやスケジュールからドキュメントまで一元管理出来ます。
※ 微妙に使いにくい部分もあります
※ 無料で使えます (一番大事)



議事録等もこれまではGoogle Documentを使う事が多かったのですが、
今回からはNotion内のDocumentとして記録していました。

実際のドキュメント


お互いの進捗状況やハードとソフトの間のパラメータ共有、認識の擦り合わせに使った図等、
slack/discord/Google Documentに散逸してしまいがちですが、
Notionで一つのアプリケーションにまとめられ、効率的に開発が進められたと思います。


「取り敢えずNotion見れば分かるか。取り敢えずNotionにメモっとくか。」と出来るのは非常に大きいです。

ただしスケジュール機能のUIが非常に見づらい、スプリントの機能が使いにくいなど、
微妙に手が届かないむず痒さもあります。

どうせ無料なので暇な時に試してみて下さい。



4. 最後に

ロボマスやデータ駆動型PIDなどはこれからの学ロボでも流行っていくのではと考えています。
少しでも何か参考になれば嬉しいです。

実はまだまだ試したいネタが沢山あったのですが、全ては無理でした。
欲張りすぎるのは良くないですね。

ompl+αでリアルタイム経路生成

この記事は学ロボアドベントカレンダー11日目の記事です。
adventar.org


この記事はompl+αでロボットアームの経路をリアルタイムで生成する手続きをまとめた記事です。


実際の経路

FMT*で生成したglobal pathをスプラインで二重に補間し、台形加減速まで経路側で処理します。
moveitは使いません。

動作環境にはUbuntu 22.04 + ROS2 humbleを想定しています。




==============目次==============

===============================





1. 一般的な経路生成の流れ

実装した手法について説明する前に、一般的な経路生成の流れを説明します。


一般的に、ロボットの経路生成を行う場合、先に大雑把なglobal pathを引いてから
問題設定に合わせて細かいlocal pathを引くという手続きを取る事が多いです。

この様に問題を分割する事で、経路全体を生成するには実行時間が足りない場合でも、
1周期内に必要なパスの生成が間に合うようにlocal pathの計算を調整する事で対応出来るようになります。



特にアームの軌道生成では、global pathの生成サンプリングベースの手法を、
local pathの生成には最適化ベースの手法を用いる事が多いです。


というのも、サンプリングベースの手法は、複雑な系に対してでも実行可能な経路を生成する事が出来ますが、
計算コストが大きく、滑らかな経路が引かれるとも限らない為、後処理が必要となります。(ex. RRT、PRM)

対して、最適化ベースの手法は、高速に滑らかな実行可能解が得られるものの、
目的関数が非凸で局所解に嵌りやすく現実には不適切な経路が引かれる事も多いです。(ex. CHOMP、STOMP)

これらの特徴から、サンプリングベースの手法で予めおおまかな経路を引いた後に、
その経路を元に最適化ベースの手法で滑らかな経路を生成するといった手続きは妥当だと考えられます。
他にも学習ベースで入力からEnd2Endで経路を生成する手法も存在しますが、ここでは扱いません。


moveit!を例に取ると、omplに実装されたサンプリングベースの手法でglobal pathを生成し、内部に
実装されたCHOMP/STOMPでlocal pathを生成するという手続きでロボットの経路/動作を生成しています。
global pathは系の性質に合わせて特化したアルゴリズムも多く、問題に合わせて選べる様になっています。



ここまで長々と説明しましたが、この記事ではlocal pathには最適化ベースの手法を用いずに、
global pathを直接補間する事でlocal pathを補間する方法を紹介します。

この記事のコードはキャチロボ2023年度大会に向けて実装、投入されていたコードです。



2. ompl+αを用いた経路生成

以降では手先座標の始点座標と終点座標を与え、状態制約下でそれらを結ぶ経路を生成する問題を考えます。
複数のwaypointを通過する経路を引く問題も、経由点での連続性を無視すれば、上記の問題に帰着出来ます。

また、アームの座標系が2種類ある事に注意して下さい。(手先座標/パラメータ座標)


今回の実装では、FMT*でstart/goalの2点を繋ぐglobal pathを生成し、
パラメータ空間上で粗めにスプライン補間してから、手先座標上でもう一度補間する事で経路を生成します。

台形加減速も経路側で実装します。


経路生成のコードはこちらにまとまっています。
catch23_robot_controller/src/trajectory/r_theta_optimal_planning.cpp at main · Emile-Aquila/catch23_robot_controller · GitHub

突貫工事で作ったリファクタ無しのコードなのでかなり汚いです!!(ごめんなさい)



プログラム全体の流れとしては、
・本体のNode (Client)からServiceにwaypointsを投げる
・Serviceでglobal / local path含めて全て生成し、Nodeに送信　(経路生成)
・Nodeが受け取ったpathを目標位置として、F7 (micro-ROS)に一定周期で順次送信　(経路追従)
といった流れです。

キャチロボで使ったプログラムはmicro-ROS側もROS2側も両方公開しています。




3. global pathの生成

FMT*を用いてパラメータ空間上でglobal pathを生成します。

3.1 FMT*について

FMT*(Fast Marching Tree*)は、漸近最適なサンプリングベースの経路計画アルゴリズムで、
複雑な制約を伴う高次元空間での動作計画をターゲットとして考案された手法です。

2015年にこちらの論文で提案されました。
arxiv.org


FMT*では、予め生成したサンプル点から、DPを用いて再帰的にグラフの構築と探索を行います。
経路点同士の再接続計算時に、衝突チェックを遅延評価にする事でチェック回数を減らしているのが特徴です。

やたらめったら早くて、それなりに良い経路が得られます。



3.2. omplの導入

ompl公式のInstall手順に従えば良いです。
ompl.kavrakilab.org


omplはaptでinstallできます。

sudo apt install libompl-dev ompl-demos

あとはCMake等で適当に呼び出せばomplを使う事が出来ます。




CMakeであれば、

find_package(Boost REQUIRED)
find_package(Boost COMPONENTS program_options REQUIRED)
find_package(ompl REQUIRED)
include_directories(
        ${EIGEN_INCLUDE_DIR}
        ${OMPL_INCLUDE_DIRS}
        ${Boost_LIBRARIES}
)

でOMPLが使えるはずです。(EigenとBoostが必要)





3.3. FKとIK

FK: パラメータ座標 -> 手先座標
IK: 手先座標 -> パラメータ座標　と定義します。


今回は円筒座標系なので簡単に計算できます。
catch23_robot_controller/src/utils/robot_state.cpp at main · Emile-Aquila/catch23_robot_controller · GitHub

冗長自由度がある場合には、特異姿勢の処理や複数解の存在に気を付けて実装して下さい。




3.4. omplのプログラムを書く

自分で実装する場合には公式のDemoを見るのが早いです。
実装をまとめたリポジトリもあります。

ompl.kavrakilab.org

github.com



対応するStateSpaceを作成 (閉区間を設定)
　-> SpaceInfomationを作成 (状態制約はここに対応)
　-> アルゴリズムを指定し、ProblemDefinitionを作成してstart, goalを設定
　-> planning
の流れが基本になります。



簡単なコードなので詳細は省きますが、以下のコードでglobal pathを生成します。
こちらを参考にしてClassにまとめて実装しました。
ompl.kavrakilab.org


・planning全体

auto state_space(std::make_shared<ob::RealVectorStateSpace>(3));
auto space_info = _space_info_our_area;
if(is_common)space_info = _space_info_common;  // 共通エリアの処理

// start point
ob::ScopedState<> start(state_space);  // (theta, r, phi)
ArmState start_tmp = arm_ik(start_tip);
std::cout<<"start: x,y,theta -> " << start_tip.x <<", "<< start_tip.y <<", "<< start_tip.theta <<std::endl;
std::cout<<"start: r,theta,phi -> " << start_tmp.r <<", "<< start_tmp.theta <<", "<< start_tmp.phi <<std::endl;
start->as<ob::RealVectorStateSpace::StateType>()->values[0] = start_tmp.theta;
start->as<ob::RealVectorStateSpace::StateType>()->values[1] = start_tmp.r;
start->as<ob::RealVectorStateSpace::StateType>()->values[2] = start_tmp.phi;

// goal point
ob::ScopedState<> goal(state_space);  // (1, 1, 1)
ArmState goal_tmp = arm_ik(goal_tip);
std::cout<<"goal: x,y,theta -> " << goal_tip.x <<", "<< goal_tip.y <<", "<< goal_tip.theta <<std::endl;
std::cout<<"goal: r,theta,phi -> " << goal_tmp.r <<", "<< goal_tmp.theta <<", "<< goal_tmp.phi <<std::endl;
goal->as<ob::RealVectorStateSpace::StateType>()->values[0] = goal_tmp.theta;
goal->as<ob::RealVectorStateSpace::StateType>()->values[1] = goal_tmp.r;
goal->as<ob::RealVectorStateSpace::StateType>()->values[2] = goal_tmp.phi;

// problem definition
auto prob_def(std::make_shared<ob::ProblemDefinition>(space_info));  // problem instance
prob_def->setStartAndGoalStates(start, goal);
prob_def->setOptimizationObjective(getBalancedObjective1(space_info));  // 目的関数の設定

// planner
auto planner = allocatePlanner(space_info, PLANNER_FMTSTAR);
planner->setProblemDefinition(prob_def);  // problem instanceを代入
planner->setup();  // plannerのsetup
planner->checkValidity();

ob::PlannerStatus solved = planner->ob::Planner::solve(0.4);ob::PlannerStatus solved = planner->ob::Planner::solve(0.4);  // 0.4sで解く


StateSpace, SpaceInformationについては、Classのコンストラクタに実装しています。
分かりにくくて申し訳ないです。

・classのコンストラクタ

auto state_space(std::make_shared<ob::RealVectorStateSpace>(3));
auto state_space_common(std::make_shared<ob::RealVectorStateSpace>(3));
matrix<double> bounds_pre{
    std::vector<double>{-M_PI_2, M_PI*2.0 - deg_to_rad(10.0)}, // theta
    std::vector<double>{325.0, 975.0},  // r
    std::vector<double>{deg_to_rad(-97.0f), deg_to_rad(110.0f)},
};

ob::RealVectorBounds bounds(3), bounds_common(3);
for(int i=0; i<bounds_pre.size(); i++){
    bounds.setLow(i, bounds_pre[i][0]);
    bounds.setHigh(i, bounds_pre[i][1]);
    bounds_common.setLow(i, bounds_pre[i][0]);
    bounds_common.setHigh(i, bounds_pre[i][1]);
}
state_space->setBounds(bounds);  // bounds for param
state_space_common->setBounds(bounds);  // bounds for param

_space_info_our_area = std::make_shared<ob::SpaceInformation>(state_space);
_space_info_our_area->setStateValidityChecker(std::make_shared<ValidityCheckerRobotArea>(_space_info_our_area));
_space_info_our_area->setup();
_space_info_our_area->printSettings(std::cout);

_space_info_common = std::make_shared<ob::SpaceInformation>(state_space_common);
_space_info_common->setStateValidityChecker(std::make_shared<ValidityCheckerRobotAreaCommon>(_space_info_common));
_space_info_common->setup();
_space_info_common->printSettings(std::cout);


キャチロボのルールに合わせてプログラムをまとめただけで、
基本的にはTutorialやDemoのプログラムと同じです。



4. local pathの生成

スプライン補間を二重にかける事でlocal pathを生成します。
台形加減速も実装します。

勿論、スプライン補間では状態制約を陽に扱う事が出来ないので、
得られた解に対して、制約を守るようにclipした物を最終的な経路として与えます。

正直な話、スプライン補間を選んだのはCHOMP等を実装してる余裕が無かったからです()



4.1. スプライン補完

(3次)スプライン補間は、点同士を3次多項式で補間する補間法です。
高速に計算できる上に、 C^2級の曲線が得られるのが嬉しい所です。

媒介変数を用いて計算する場合には、得られた曲線を区分する時の長さを調整しやすいのもメリットです。


こちらの実装を使わせて頂きました。
github.com


今回はパラメータ座標上で滑らかなpathを引きたいので、先にパラメータ座標上で補間します。
その後、次は手先座標上でもう一度補間する事で、pathの間隔調整(台形加減速)を行います。




4.2. 台形加減速の実装

local path側で台形加減速を実装します。

この様にする事で、固定周期で経路点を目標値として制御器に送るだけで、
滑らかに経路追従させる事ができます。


local pathの点同士の最小間隔と最大間隔、間隔の最大変化幅をパラメータとして与え、
最大間隔を取る点が出来るだけ多くなる様にlocal pathの点を取り直します。

間隔計算は二分探索で、local pathを取り直す処理はスプライン補間で実装します。



まず、最小間隔から最大間隔に変化するステップ数 nを決め打つと、間隔の変化幅が求まり、
入力したステップ数が実行可能か求まります。

実行可能な最小ステップ数を求める事を考えます。


ステップ数が実行可能か判定する手続きを f(n):\mathrm{N} \rightarrow \{0, 1\}として捉えると、
n \in \mathrm{N}について単調増加となるので、二分探索を用いる事でこの問題は容易に解く事が出来ます。

これは競技プログラミングで使われる初歩的なテクニックの一つで、決め打ち二分探索と呼ばれています。


二分探索なので実装も簡単かつ高速です。めぐる式派です。


実装はこちらにまとまっています。
catch23_robot_controller/src/trajectory/spline_utils.cpp at main · Emile-Aquila/catch23_robot_controller · GitHub



5. 実際の動作

実際の動作映像はこちらです。


コントローラーのボタンが押されてから、経路を生成->追従させています。
今回の実装では、FMT*の最大計算時間を0.4 sとしていて、そこがボトルネックになっています。

簡単に実装できる範囲で動いて良かったです...。

STM32 Nucleo-F767ziでmicro-ROSを動かす(CubeMX+CMake)

この記事は学ロボアドベントカレンダーの1日目の記事です。
今年も先駆けを務めさせて頂きます。殿も務める予定です。


adventar.org


この記事はSTM32 Nucleo-F767zi上でmicro-ROSを動かす為の設定をまとめた記事です。
CubeMX+CMakeで動かす事を前提にしています。

動作環境はUbuntu 22.04を想定しています。


==============目次==============

===============================


1. micro-ROSとは

micro-ROSはマイコン等の組み込みシステム上で動作するROS2のclientライブラリです。
real-time性が考慮されており、一般的なROS2よりも低遅延で動作します。

PC上のROS2とマイコンの間で(表面的には)通信レイヤーを意識せずtopic通信等が可能であり、単純に便利です。

micro.ros.org



mROS2と違ってhost PC上でagentを動作させる必要があるのがデメリットですが、
ユーザーが多く、調べれば色々と出てくるのがメリットだとは思います。

※ micro-ROSはμROSと略す事が多いっぽいです。(mROSとするとmROS2と判別しにくい)


今年のキャチロボで実戦投入しました。



2. Agentの準備

PCで動かすmicro-ROS Agentを導入します。

以下のリポジトリをcloneして下さい。
github.com


パス等通したら次を実行して行きます。

・micro-ROS Agentを作成

ros2 run micro_ros_setup create_agent_ws.sh  # Download micro-ROS-Agent packages
ros2 run micro_ros_setup build_agent.sh  # build step
source install/local_setup.bash  # setup env


・micro-ROSのAgentを動かす

ros2 run micro_ros_agent micro_ros_agent serial -b 115200 --dev /dev/ttyACM0  # Run a micro-ROS agent

※ /dev/ttyACM0のデバイスと115200のserial通信でやり取りする場合の設定


Agentは簡単に動かせます。





3. firmware libraryの設定

3-1. CubeMXの設定

CubeMXのpin設定は以下のリポジトリと同様に行います。
github.com


具体的には、以下の様にpinを設定します。

・System Core/RCC
　・HSE: Crystal/Ceramic Resonator

・System Core/SYS
　・Timerbase Source: TIM1

・Connectivity/USART3
　・NVIC Setting: Mode:Asynchronous
　・USART3 global interrupt: Enable

　・DMA SettingsでUSART3_RXを追加
　　・Request SettingのModeをCircular
　　・PriorityをVery High

　・DMA SettingsでUSART3_TXを追加
　　・PriorityをVery High

・Connectivity/USB_OTG_HS
　・External Phy: Device Only

・Middleware and Software Packs/FREERTOS
　・Interface: CMSIS_V2
　・FPU: Enable
　・Taskを追加
　　・stack sizeは必ず3000以上に
　　　

・Middleware and Software Packs/USB_DEVICE
　・Class For FS IP:Communication Device Class(Virtual Port Com)

・Project Manager/Code Generator
　・Generate peripheral initialization as a pair of .c/.h files per peripheralをonにする。(忘れやすい)




3-2. library generationの設定

CubeMX+Makefileで動かす例が公式から上がっているので、こちらを流用します。
github.com



ここで説明する設定の全体像は以下のリポジトリにまとめてあります。
GitHub - Emile-Aquila/micro-ROS_example_for_F767zi: Example Codes for micro-ROS on F767zi with CMake + STM32 CubeMX


まずmicro_ros_stm32cubemx_utilsリポジトリをcloneし、
・extra_sources
・microros_static_library
の２つのフォルダのみ残して、projectフォルダへ移植します。(これ以外のフォルダ/ファイルを削除)


さらに、microros_static_library/library_generation下のlibrary_generation.shを
FPUを使えるように無理矢理書き換えます。

#export RET_CFLAGS=$(make print_cflags)
export RET_CFLAGS="-mcpu=cortex-m7 -mthumb -mfpu=fpv5-sp-d16 -mfloat-abi=hard -specs=nosys.specs -fdata-sections -ffunction-sections"

詳細: micro-ROS_example_for_F767zi/micro_ros_stm32cubemx_utils/microros_static_library/library_generation/library_generation.sh at main · Emile-Aquila/micro-ROS_example_for_F767zi · GitHub


加えて、library_generation下のcolcon.metaのmicroxrcedds_clientに以下を追加します。

"-DUCLIENT_HARD_LIVELINESS_CHECK=ON",
"-DUCLIENT_HARD_LIVELINESS_CHECK_TIMEOUT=1000"

Agentからmicro-ROSの生存確認が出来る様になります。

ﾈｯﾄｻｰﾌｨﾝしてて設定を見つけました。
4.3. Reconnections and liveliness — Vulcanexus 1.0.0 documentation



micro-ROSで使うpackageを追加したい場合は、
library_generation/extra_packages/extra_packages.reposにリポジトリを追記して下さい。



最後に以下を実行してlibraryを生成します。
micro_ros_stm32cubemx_utilsが含まれるディレクトリで実行して下さい。

sudo docker run -it --rm -v $(pwd):/project --env MICROROS_LIBRARY_FOLDER=micro_ros_stm32cubemx_utils/microros_static_library microros/micro_ros_static_library_builder:humble

※ static libraryの設定を書き換える度に上のコマンドを実行します。



3.3. CMakeの設定

CMakeの設定について説明します。

CMakeを選択したのは普段CLionで開発していて都合が良い&Makefileよりも書きやすいからです。
(ClionでCubeMXのProjectを読み込む方法は以下の記事にまとめられています。)
・STM32CubeMX プロジェクト | CLion ドキュメント



設定の要点を抜き出すと以下の様になります。

・CMakeLists.txt

set(MROS_UTILS_DIR "./micro_ros_stm32cubemx_utils")

include_directories(
        ${MROS_UTILS_DIR}/microros_static_library/libmicroros/microros_include
        Core/Inc ... 
)

(略)


file(GLOB_RECURSE SOURCES
        "${MROS_UTILS_DIR}/extra_sources/custom_memory_manager.c"
        "${MROS_UTILS_DIR}/extra_sources/microros_allocators.c"
        "${MROS_UTILS_DIR}/extra_sources/microros_time.c"
        "${MROS_UTILS_DIR}/extra_sources/microros_transports/dma_transport.c"
        "Core/*.*" "Middlewares/*.*" "Drivers/*.*" "USB_DEVICE/*.*"
)

(略)

target_link_libraries(${PROJECT_NAME}.elf ${CMAKE_SOURCE_DIR}/micro_ros_stm32cubemx_utils/microros_static_library/libmicroros/libmicroros.a)
link_directories(${PROJECT_NAME}.elf ${CMAKE_SOURCE_DIR}/micro_ros_stm32cubemx_utils/microros_static_library/libmicroros/libmicroros.a)

(略)

公式のstm32cubemx_utilsに記述されていたMakefileの設定をCMakeに移植しただけです。
細かいCMakeの書き方や中身については省略させて頂きます。


これでCMakeでbuild出来るようになったはずです。
(.ioc -> CMakeはCLionにやらせています。)




4. micro-ROSのCodeを書く

rclcで書きます。頑張って下さい。

micro-ROS公式のドキュメントが参考になります。
micro.ros.org



キャチロボのF767ziのコードを公開しているので、こちらも参考になるかもしれません。

Catch23_F7/Core/Src/freertos.c at main · Emile-Aquila/Catch23_F7 · GitHub

補足すると、上のコードは

rmw_ret_t ping_result = rmw_uros_ping_agent(1000, 5);  // ping Agent

でAgentの接続確認を行っています。(Tutorialだとこの機能が探しにくかったはず)


RTOSのタスクの一つとしてmicro-ROS Clientは動作する様に実装されているので、
他のRTOSのタスクも一応動かす事が出来ます。


何かしら参考になる部分があれば嬉しいです。

エクストリーム帰寮参加記

この記事はポエムです。ｱﾚﾙｷﾞｰのある方は閉じてください。
PCで見る事を想定して書いています。スマホだとレイアウトが崩れてるかも(ごめん)


この記事はエクストリーム帰寮アドベントカレンダー2022の17日目の記事です。
adventar.org




時が経つのは早いもので、気がつけば大学生活も後半です。
大学生になって3年、エクストリーム帰寮も3走目。
2022年も年の瀬が近づいてきました。寂しく賑やかな冬の訪れです。



==============目次==============

===============================
↑ 1~4はこれまで蓄積してきた知見や小話、5がエク帰2023の話。(時間が無い人は5から読んで下さい。)





1. エクストリーム帰寮とは?

熊野寮祭の人気企画の一つで、11月末~12月初頭の金曜の夜中に熊野寮から車で運ばれ、
落とされた所から地図を見ずに歩いて帰ってくるという企画です。

京北や鯖街道、奈良や信楽・伊賀など真っ暗な山の中に落とされる事が多く、本物のスリルを味わえます。



2. なぜ参加するのか?

自分は熊野寮の関係者では無いですし、熊野寮に入るのは年に一度エク帰の時だけです。
ではなぜ熊野料の様な魔境*1に足を踏み入れてまで参加するのかと言うと、
落とされた直後のスリルと達成感が忘れられないからでございます。

ここまで自分の限界に挑める体験は滅多に無いのではないでしょうか*2。


3. これまでの長距離移動歴

本題に入る前に、自己紹介を兼ねてこれまでの長距離移動歴を晒しておきます。
普段はロボティクスや制御工学・強化学習などを勉強したり、何かしら開発してたりします。


・エクストリーム帰寮2020(木津川の山の中、徒歩、45km程度)
・エクストリーム帰寮2021*3(日吉ダム、徒歩、57.9km・893m↑)
・エクストリーム帰寮2022(芦生演習林、徒歩、75.9km・763m↑)
・ビワイチ(クロスバイク、214km)
・ビワイチ(ロードバイク、214km・738m↑)
・神戸(クロス、170km)
・神戸(ロード、175km・640m↑)
・永源寺ダム(ロード、160km・723m↑)
・伊賀(ロード、132km・1362m↑)
・日吉ダム(ロード、123km・1700m↑)
・彦根(ロード/クロス、120km程度) × 3~4
・関ヶ原(クロス、途中で知り合いの車に乗せていただいたので140km程度)
・下宿→大学(2~3km、毎朝通学RTA((交通規則は厳守しています。)))
など。これまでの走行距離のΣは4000km程度。


エクストリーム帰寮2021(貴船でリタイア)
ダムって良いよね(永源寺ダム(コンバインダム(ロックフィル+重力式コンクリート))。記事と関係ない)

4. 準備と装備

最低限の安全の為の装備と、人様に迷惑を掛けない為の装備は準備しましょう。

4.1 これまでの知見

・長時間使えるライトが便利

電灯のある道におろしてもらえるほど現実は甘くない。そして日が出ても明るいとは限らない。
山道は日中でも暗いし、夜間は霧が出ている時もある。

そもそも車が来た時に自分たちの存在を認識してもらえる様にしておかないと危ないです。
誰もいないと思って山道を走っている車ほど怖いものはありません。

ﾌﾞﾚﾌﾞﾚの日吉ダム(2021エク帰)
夜の山道(2020エク帰)
・ポリ袋が便利(回路部品とか入れてるアレ)

補給食やちょっとしたゴミをまとめて管理するのに便利。

ゴミはちゃんと持って帰りましょう。


・補給食は大事

ハンガーノックは本当に怖い。カロリーは正義。カロリーは正義。(大事な事なので2回言いました)
山の中だと数10kmコンビニ・電波無しだったりする。

水がなくても食べやすい羊羹がオススメ。
カロリーメイトは腹持ちはいいけれど、口がﾊﾟｻﾊﾟｻになるのでしんどいです。


・防寒/雨対策

夜の山の中は死ぬほど寒い。2021は雨/霧/雪のオンパレードだったので本当に酷かったのですが、
2022は暖かく天気も良かったので平和でした。
ちなみに寝る時はペットボトルを枕に、ダンボールやゴミ袋を布団にすると暖かいです。(NHKロボコンで学んだ)

雪の芹生峠(エク帰2021)
・レッドブル

翼を授けてくれる。カフェインの翼で数多の修羅場を駆け抜けてきました。
デスマの時に使うもので便利なものはエクストリーム帰寮でも使えると思います。


・土地勘があると便利

運が良いと酷道・腐道を回避したり、逆に突撃したりできます。

物好きな方はぜひ持越(日吉)や芹生(上黒田)などに挑戦してみてください。
477は最近道が綺麗になってるので物足りないかもしれない。

ちなみに芹生峠には教習車も走っていました。頭おかしい。

この先酷道(サイクリングにて)
・不審者扱いされる

ゆめゆめ忘るるなかれ。真夜中に歩き回っている我々は、傍から見れば不審者そのものなのだ。
そもそも車で知らない場所に運ばれて歩いて帰ってくる行為そのものが不審な行為なのかもしれない。



4.2 装備

今回のエク帰寮では以下の装備での出走を予定しています。
・タオル
・リュック(18L、560g)
・羊羹(セブンイレブン、kCal)
・モバイルバッテリー(Ankerはいいぞ)
・ライト(RN800、GVolt70、テールライト)
・ｼﾞｯﾌﾟﾛｯｸ × 2
・キャラメル(森永)




5. 出走後の経過

歩行ルート全体(芦生演習林 --> 熊野寮) : 75.9km・763m↑・総移動時間21時間10分
※ 重要なスポットにはgoogle mapのリンクを貼ってるので、位置を確認してみて下さい。

エクストリーム帰寮2022

5.1 輸送中の話

花背峠(酷道477号) --> 佐々里峠(腐道38号) --> 芦生演習林のルートで運ばれました。

輸送ルート(青)

とんでもない場所に運ばれている中、自分はと言えばひたすら車酔いしてました...。
花背付近では既にｸﾞﾛｯｷｰ状態になっていて、なかなか気持ち悪かったです。


40kmの所(広河原)で別のチームが落とされ、そこで一度トイレ休憩を挟みました。
新幹線以外のボットン便所は久しぶりだったので、楽しかったです。
↓ファンが回る音が室内にも伝わっているらしく、用を足してて落ち着かない。

広河原のトイレ(張り紙)

5.2 投下

自分は美山に運ばれてるんだと思っていたのですが、まさかの芦生演習林でした。
運んで下さった方にお礼を言って真っ暗な道へ出発です。

芦生演習林


芦生演習林にも色々と面白い話がある様です。こうした人間の営みって良いですよね...。
www.geoid.gr.jp

https://takeshi-ise.tumblr.com/post/171158859784/%E8%8A%A6%E7%94%9F%E7%A0%94%E7%A9%B6%E6%9E%97%E3%81%AE%E6%AD%B4%E5%8F%B2%E3%81%A8%E5%88%86%E6%B0%B4%E5%B6%BA%E4%BA%AC%E9%83%BD%E6%BB%8B%E8%B3%80%E7%A6%8F%E4%BA%95
takeshi-ise.tumblr.com



ちなみに降ろされた時点で、以下の事が分かっていました。*4
・佐々里峠を超えた美山の奥に落とされている
・帰るルートは2通り*5
　　　府道38号(南へ) --佐々里峠--> 国道477 --花脊峠--> 鞍馬寺 --> 熊野寮 (鞍馬までずっと山道)
　　　府道38号(西へ) --> 美山町 --> 国道162(無限162編) --> 高雄 --> 熊野寮 (遠回り・人が存在)　
・クマが居るらしい
・佐々里峠付近はほとんど人通りが無い
(京都北部の道では162と367が車通りが多い事で知られる。477も意外と車が居る。)


今回は「出来るだけ勾配が少なく、補給しやすいルートを選ぶ」という方針だったので後者を歩く事にしました。
自己位置の同定やルート選択については、地の利があった事も相まって割と簡単でした。
まさか15kmも遠回りするとは...



ちなみにロボットなどで経路を経過する時は、
「大域的な経路 --> 局所的な経路」の順番に計画する事が多いです。


これはエクストリーム帰寮でも同じ事が言えて、現在の位置と京大の位置関係(方角など)から、
　1. 経由地点(特定のスポットや市や街でも良い)を一通り決める
　2. 各経由地点を結ぶ経路を青看板や地域の案内地図から選択
という手続きを繰り返す事で熊野に帰る経路が計画できる事が多いです。

5.3 道中

75kmは長すぎて自分の文章力では一纏めに記述出来無いので、3パートに分けて記述します。

5.3.1 明るくなるまで(芦生演習林〜 知井の里 情報発信館)

まずは夜明けまでの15kmです。芦生演習林から38号をひたすら西へ進みます。

芦生演習林〜知井の里 情報発信館まで


夜の山道は真っ暗で、めちゃくちゃ怖いです。
鹿やクマ・イノシシが来ないことを祈りながら歩きます。

土砂崩れで道が崩落していた


土砂崩れで道が崩落して迂回路が作られていました。最近崩落したっぽいですね。
(詳細 : 京都府道路情報管理・提供システム、ID382)


人間が作った道を自然が破壊しては人間が対応する、
自然の雄大さとそれを変形する人間の力強さにはロマンを感じます。


情報過多


迂回路が煌々と照らされていてやたら明るいので、
補給を取りつつ看板から位置の認識が合っている事を確認します。


エリアマップ発見

さらに歩くこと数十分、エリアマップを発見しました。ようやく街に着きそうです。
(実はこの地図、距離が圧縮されていて街は当分先です。)






ここから数時間歩き続けます....。
途中で犬に吠えられたり、鹿に遭遇したりしました。


田歌周辺で犬に吠えられた時は死を覚悟しましたが、どうやらお家の番犬だったようで命拾いしました。
勿論写真を取る余裕はありませんでした。まぁそういう日もありますよね。

このエリアでは番犬を飼っているお家が多い様で、2軒のお家のﾜﾝｺに吠えられました。
ﾜﾝｺさんや...起こしてすまんな...。


交番


更に歩くこと数時間、5:00頃にようやく村(?)に到着し、初めて車&通行人に遭遇しました。
ちょうど交番があったので、その明かりを頼りに休憩です。


ここからは街っぽくなり、時たま人間を見かける事もありました。
更に1時間ほど歩くと「知井の里 情報発信館」と書かれた道の駅(?)っぽい所に到着です。

5.3.2 明るくなってから

次は日が出ている間の30km弱です。ここまではまだ記憶があります。

知井の里〜ウッディー京北


知井の里 情報発信館で休憩し、補給を取った所で出発です。
ここまでお世話になったRN800(自転車ライト)をリュックに戻して、夜明けを迎えた美山を進んでいきます。

綺麗な町並み

美山は綺麗な茅葺屋根の家が連なる事で有名です。
鹿対策のネットで街の山際が囲われる構造になっていてｳｫｰﾙﾏﾘｱを彷彿とさせます。

看板とヴィンテージクロモリ

gazelleというメーカーのクロモリロードが置いてありました。
オランダの老舗自転車メーカーらしいです。自分も初めて聞きました。

お饅頭とコロッケ

道の駅で売店のおばちゃんに鹿肉コロッケを頂きました。滅茶苦茶美味しかったです。

京都まで52km
お洒落な建物

お洒落なお店もありました。


5.3.3 陽が落ち始めてから(ウッディー京北〜熊野寮)
ウッディー京北〜熊野寮

正直な話、ここらへんから記憶がありません...。

自転車で走ったことのあるルートだったので、
余計に自転車の距離感と徒歩での距離感の差に苦しめられていました。



6. まとめ

よく40kmまでは健全、50kmから苦行、60kmを超えると悟りの領域とはよく言ったものです。
70kmを超えた段階では痛みのあまり、信号の度に手すりにもたれたり、地べたに座ったりしていました。

100km超えの皆様は尊敬しか無いです。



・後日談
先日チャリでもう一度美山に行ってきました。自転車だと一瞬ですね。





注釈

*1:批判の意図はありません

*2:他にはNHKロボコンや二郎系ラーメンなどが知られる

*3:実は神明峠に抜ければすぐに帰れる距離だったのだが、間違えて京北に抜け、チキって162をスキップしたので地獄になった。

*4:以前友達と京北の峠をチャリで走破しようぜw的なノリでルートを組んだ事があった。

*5:実は3通りあって、鯖街道へ抜けるルートが存在する