われがわログ

最適化アルゴリズムとかプログラミングについて書きたい

M5Paperで部屋の環境監視用ダッシュボードを作った

M5Stack

最近発売されたM5Paperで、部屋の温湿度・CO2監視用ダッシュボードを作ったのでメモしておく。

成果物は下記の写真の通り。

f:id:estshorter:20210101193314j:plain
M5Paperの画面

ソースコードは以下に置いた。

github.com

M5Paper

M5Paperは2020/11/27に発売されたESP32搭載マイコンボードで、その特徴は大型電子ペーパー（540x960、4.7インチ）を搭載していることだ。発売当初は即品切れになってしまい入手できなかったのだが、年末の再入荷でやっと購入できた。

www.switch-science.com

成果物詳細

今回作ったシステムの構成は以下のとおり。

f:id:estshorter:20210109145633p:plain
システムダイアグラム

M5Paperには環境センサ（ENV II）が接続されており、ここから温湿度情報を取得している。M5Paperにも温湿度センサは内蔵されているが、本体の熱の影響を遮断するため外部センサを用いた。また、ENV IIユニットは気圧も測定できるが、今回は取得していない。CO2データは、以前の記事で紹介したCO2Miniで測定したものだ。Raspberry Pi 4でCO2データを取得し、http通信でM5Paperに渡している。LED UnitはCO2のインジケータとして使用している。

waregawa-log.hatenablog.com

時刻は外部のNTPサーバ（NICT）に接続して同期する。多機能ボタンを押し込むか、一週間経過毎に同期するようにした。また、多機能ボタンを右へ倒すと画面のリフレッシュを、左に2秒以上倒し続けるとシャットダウンするよう設定してある。

画面の文字がギザギザしていることが今後の課題。IPAフォントの40ptをLovyanGFXで3倍に拡大して表示していることが原因だとは思うのだが、滑らかに表示させるにはどうしたらよいかよくわからない、、

感想

数日使ってみた感想としては、 以前購入した、Xiaomi製電子ペーパー時計では30分に一回画面がリフレッシュ（暗転）して視覚的にうるさかったのだが、M5Paperを使うとリフレッシュタイミングを自由に制御できるのがよい。ちゃんと測定していないが、1日くらいの間なら暗転させなくても大丈夫なように思える。

捜索理論勉強メモ

Optimization

完全に自分の勉強メモ。 静止目標に対する捜索活動の最適化についてまとめた。 参考書は下記の2つ。

コロナ社から出ている本の方が入門向けだが、省かれている部分もある。ただ、不完全定距離センサーの横探知確率や、同センサーの区域捜索に関する議論はコロナ社の方が詳しかった。

捜索理論における確率モデル (シリーズ 情報科学における確率モデル 3)

捜索理論における確率モデル (シリーズ 情報科学における確率モデル 3)

全般的にはこっち↓の方が詳しい。両方参照しつつ読み進めると理解の助けになる。

捜索理論 3訂版―捜索オペレーションの数理

捜索理論 3訂版―捜索オペレーションの数理

捜索センサのモデリング

瞬間的な探知能力のモデリング

不完全定距離発見法則

一定距離 R_{0} \in \mathbb{R}^{+}以下では一定確率で探知できるというモデリング。

g(r(t)) =\begin{cases} g_0, & \text{for} \quad r(t) \leq R_0 \\ 0, & \text{otherwise} \end{cases} \\ b(r(t)) = \begin{cases} b_0, & \text{for} \quad r(t) \leq R_0\\ 0, & \text{otherwise} \end{cases}
  • g: \mathbb{R} → \mathbb{R}_0^+は、離散スキャンセンサ（アクティブソナー等）の1回の「べっ見」に対する目標探知確率（べっ見探知確率）
  • b: \mathbb{R} → \mathbb{R}_0^+は、連続スキャンセンサ（目視・パッシブソナー等）の単位時間当たりの目標探知確率
  • r(t) \in \mathbb{R}_0^+はセンサからの距離

なお、 b(r(t))=\inftyのとき、微小時間\Delta tに対してb(r) \Delta t = 1となって探知確率が1となる

逆三乗発見法則

単位時間当たりの目標探知率b: \mathbb{R}\rightarrow \mathbb{R}_{0}^{+} が、 センサからの距離r(t) \in \mathbb{\text{R}}^{+} の三乗に反比例するという連続スキャンセンサ向けのモデリングであり、広く用いられている。

b\left( {r\left( t \right)} \right) = \left( \frac{k}{r\left( t \right)} \right)^{3}
  • k \in \mathbb{R}^+は定数

目視による海上の捜索では立体角をもとに

b(r(t)) = \frac{\alpha A' h}{( {r^{2}( t ) + h^{2}} )^{\frac{3}{2}}} \approx \frac{\alpha A'h}{r^{3}\left( t \right)} = \left( \frac{k}{r( t )} \right)^{3}

とモデリングされる。最後の式変形では k^ 3 = \alpha A'h とおいた。

  • \alpha \in \mathbb{R}^+は定数
  • A'は物体の面積（XY平面への射影面積）
  • hは飛行高度 ( h \ll r(t)を仮定)

目標の相対移動に対する探知能力のモデリング

以下では、センサを直交座標系の原点に固定し、目標が相対運動している場合を考える。 センサと目標との相対距離が時間の関数 r: \mathbb{R}_ 0 ^ + \to \mathbb{R} _ 0 ^ +で変化する場合には、探知確率 p(t) \in [0,1]は

p(t) = 1 - \exp\left( {- F\left( r \right)} \right)

である。ここで、探知ポテンシャルF(r) \in \mathbb{R} _ 0 ^ + は次のように与えられる。

F(r) = \begin{cases} \sum_{i=1}^n \{ -\log (1-g(r(t_i))) \} & \text{離散スキャンセンサの場合}\\ \int_0^t b(r(\tau)) d\tau & \text{連続スキャンセンサの場合} \end{cases}
  • t_{i}( {i = 1,2, \cdots,n})はべっ見を行う時刻

無限直線状を等速運動する目標で、センサとの最接近距離が xの目標の探知確率は横距離探知確率PL(x) \in [0,1]と呼ばれている。ただし、べっ見探知確率や横距離探知確率はセンサの探知能力の指標としては複雑であり、 実務では、次式で定義される有効探索幅 W \in \mathbb{\text{R}} ^ {+}が用いられている。

W = \int _ {-\infty} ^ {\infty} PL ( x) dx

この式は、横距離探知確率 PL\left( x \right)を持つセンサを、幅 W内で目標に行き会えば確率1で目標を探知し、その外では決して探知しないとする完全定距離センサに近似していることを表している。

不完全定距離センサ

\begin{aligned} PL(x) &= \begin{cases} 1-\exp \left(-\dfrac{b_0 \sqrt{R _ 0^2-x^2}}{u}\right) & \text{for} \quad |x| \leq R_0\\ 0 & \text{otherwise} \end{cases} \\ W &= 2R _ 0 \left[ 1 - \int _ 0 ^1 \exp (-2Z \sqrt{1- y ^ 2 })dy \right] \end{aligned}

ここで

Z = \begin{cases} -\dfrac{R_0 \log (1-g_0)}{ut_0} & \text{離散スキャンセンサの場合} \\ \dfrac{b_0R_0}{u} & \text{連続スキャンセンサの場合} \end{cases}

であり、uは相対速度の絶対値である。

逆三乗センサ

\begin{aligned} PL(x) &= 1- \exp \left(- \dfrac{2k^3}{ux^2}\right)\\ W&= 2 \sqrt{\dfrac{2\pi k^3}{u}} \end{aligned}

なお、 有効捜索幅 Wの経験的な値（商船の見張りから見た有効捜索幅など）は、国際航空海上捜索救難マニュアルに掲載されている。

平行捜索の探知確率

面積Aの地域を、掃引幅Sの平行経路で捜索するときの探知確率 P\left( S \right)を求める。 Sは捜索時間 T \in \mathbb{R}で面積Aを捜索し終わるように設定されるので、 S = \frac{A}{uT}の関係がある。

不完全定距離センサ

まず、 p _ 0 = PL (0)とおき、 2ip_0S \leq W \lt (2i+1) p _ 0 あるいは、 (2i+1)p_0S \leq W \lt (2i+1) p_0S を満たす i \in \mathbb{\text{N}}_{0}（一つしか存在しない）を求める。 この iを用いて、探知確率 P\left( S \right)は次のように書ける。

P(S) = \begin{cases} 1-(1-p_0)^{2i}\left\{ 1+p_0 \left(2i- \dfrac{W}{p_0S}\right) \right\} & 2ip_0S \leq W < (2i+1) p_0\\ 1-(1-p_0)^{2i+1}\left\{ 1+p_0 \left(2i+1- \dfrac{W}{p_0S} \right) \right\} & (2i+1)p_0S \leq W < (2i+1) p_0S \end{cases}

逆三乗センサ

P(s) \approx \sqrt{1 - \exp \left\{ - \left( \frac{W}{S} \right)^{2} \right\}}

最適化問題としての定式化

「捜索資源（捜索努力）」の配分問題を解き、最適な捜索計画を求める。 ただし、捜索者の運動の連続性や捜索経路は考慮しづらいため、ここでは考慮しない。

目標探知確率を最大にする最適努力配分

仮定

  1. 目標空間は Q個の離散地域（セル）から成り、セルのインデックス集合を \mathcal{Q}とおく。
  2. 捜索計画全体（捜索資源の配分）は \Phi = [\phi_1, \phi_2, \cdots \phi _ Q ] ^T \in \mathbb{R} _ 0^{+Q}なる連続量で表される。
  3. 1つの静止目標が Q個の地域のいずれかに存在し、各地域の目標存在確率 p_{i} \in \left\lbrack {0,1} \right\rbrack\left( {i \in \mathcal{\text{Q}}} \right)は捜索者が推定可能である。
  4. 捜索者が使用できる総捜索コストは C \in \mathbb{R}（例えば、捜索に要する金額や、捜索者の被るリスク）であり、また、 Cを目標地域へ分配する場合、任意の大きさに分割可能である。
  5. 捜索資源を1単位だけ地点 iに投入したときのコストは c _ i \in \mathbb{R} _ 0 ^ + である。
  6. 地域 iに目標が存在し、その地域に資源 \phi _ i \in \mathbb{R} _ 0 ^ + が投入された場合、目標の発見確率は関数 f _ i ( \phi _ i )で表され、次の性質を持つ。 f _ i(0) = 0, \quad f(\infty)\leq 1, \quad 0 \lt f'(\psi _ i) \lt \infty, \quad f''(\psi_i) \lt 0。 例えば、逆三乗センサを用いて平行捜索する場合の f _ i ( \phi _ i )は f _ i ( \phi _ i ) = \sqrt{1 - \exp \left\{ - \left( \frac{W \phi _ i}{S} \right)^{2} \right\}} である（これだと f(0)で微分不可なので、理論的には好まれない？）。また、ランダム捜索する場合（センサの種別は問わず）は f _ i \left( {\phi_{i}} \right) = 1 - \exp\left( - \frac{W \phi _ i}{S} \right) である。
  7. 捜索者は総捜索コスト Cの制約下で、目標探知確率を最大とする捜索を計画する。

以上の前提の下、最適捜索のための資源配分問題は次の通り定式化される。

\min_{\Phi} \quad \sum_{i \in \mathcal{Q}} p _ i f _ i (\phi _ i) \quad \text{subject to} \quad \sum_{i \in \mathcal{Q}} c _ i \phi _ i \leq C, \quad \phi _ i \geq 0 \quad \text{for all} \quad i \in \mathcal{Q}

Near-identity diffeomorphismsを用いた差動二輪ロボットの制御

Python Robotics

Turtlebot3をシミュレータ上で動作させて以来、「差動二輪ロボットの駆動則（速度指令値から進行方向速度・回転角速度への変換則）って何があるんだろう」と思っていたのだが、よさげなやつを最近発見したのでメモしておく。

東工大＆Georgia Techによる論文[1]を読んでいたところ、Georgia Tech側のロボットは差動二輪ロボットのため、論文[2]の技術で運動方程式を

\dot{x} = u

なるsingle integrators（なんて訳したらよいのかわからん）としてモデル化しているとの記述があった。

私の頭でもすんなり理解できたので詳細は論文[2]の方を読んでほしいが、 論文[2]のキモはx=[x_1\;x_2]^ T（2次元位置）を並進速度v_1と角速度\omegaで制御するとき、差動二輪ロボットでは\dot{x}_ 1と\dot{x} _ 2を両方同時に制御できないので、 それだったらz=x+ \lambda (R e_1)という新しい変数を考えてやって、これを制御すればよいというもの（違ってたら指摘ください）。 ここで、e_1はロボット固定座標系（進行方向が+X、紙面上向きが+Zの右手系）の一軸方向の単位ベクトルの、ロボット固定座標系での表現、Rはロボット固定座標系から慣性座標系への座標変換行列である。 実際にzを慣性系で時間微分してみると

\dot z = R_\lambda v + (R_\lambda e_1) \dot{\lambda} \quad \text{where} \quad v=[v_1\; \omega]^ T, R_\lambda = [Re_1 | \lambda Re_2 ]

となって並進速度v_1 \in \mathbb{R}と角速度\omega \in \mathbb{R}で\dot{z}を自由に制御できることがわかる。 \lambda(t)=0のときは\omegaが現れないことに注意してほしい。 そして、\lambda (t) \to 0となるようにうまく（論文では \dot{\lambda} = -c _ \lambda (\lambda - \varepsilon), \lambda(0) > \varepsilon >0, c _ \lambda> 0と） 設定したのちにz（これ自体はsingle integratorにできる）を制御してやれば、zを通じてxも所望の値の近傍に到達する、というわけである。なお、zがnear-identity diffeomorphismsの一つである。

制御結果は以下の通り。ダイナミクス制御（制御量が力・トルクレベル）版とキネマティック制御（制御量が速度・角速度レベル）版の両方をPythonで実装した。

dynamics:

https://raw.githubusercontent.com/estshorter/nir/outputs/dynamics.gif

kinematics:

https://raw.githubusercontent.com/estshorter/nir/outputs/kinematics.gif

リポジトリ：

github.com

なお、論文[3]によれば、Georgia Techは\lambdaが固定のものを使っているようだ。

参考文献

  • [1] R. Funada et al., "Coordination of Robot Teams Over Long Distances: From Georgia Tech to Tokyo Tech and Back-An 11,000-km Multirobot Experiment," in IEEE Control Systems Magazine, vol. 40, no. 4, pp. 53-79, Aug. 2020, doi: 10.1109/MCS.2020.2990515.
  • [2] R. Olfati-Saber, "Near-identity diffeomorphisms and exponential ε-tracking and εstabilization of first-order nonholonomic SE(2) vehicles," Proceedings of the 2002 American Control Conference (IEEE Cat. No.CH37301), Anchorage, AK, USA, 2002, pp. 4690-4695 vol.6, doi: 10.1109/ACC.2002.1025398.
  • [3] P. Glotfelter and M. Egerstedt, "A Parametric MPC Approach to Balancing the Cost of Abstraction for Differential-Drive Mobile Robots," 2018 IEEE International Conference on Robotics and Automation (ICRA), Brisbane, QLD, 2018, pp. 732-737, doi: 10.1109/ICRA.2018.8461234.

温湿度計付き置時計 MHO-C303 へのラズパイからの時刻同期

Programming C++

210411追記

MHO-C303は廃版になったもよう、、

CO2計測ついでに湿度もロギングしたいなぁと思って色々調べた結果、 Xiaomi製、温湿度計付き置時計 MHO-C303を買った。*1

旧機種のlywsd02とは違い、電池が単4二つになっているので長持ちするのではと思う。ただ、要らない機能（タイマー？使ったないのでよくわからん）のためにタッチ式ボタンが上側についており、持ち上げるたびに誤作動するのが不満。 また、湿度がEMPEXのアナログ式のものに比べて10%程度高く、どっちが真値なのかよくわからない。体感的にはEMPEXのに近いが、体のセンサはあてにならんことだし。

このMHO-C303であるが、Bluetooth LEによる通信機能が付いており、スマートホン向け専用アプリ(Mi Home）で接続＆時刻同期できる。 そして、本機種は例によってハックされており、ラズパイからでもデータ取得＆時刻同期が可能である。 具体的には以下のPythonライブラリを使えばよい。

github.com

旧機種向けだが、MHO-C303に対しても問題なく動いた。 なお、温湿度データ、バッテリレベルはadvertiseもされているとのこと。

tasmota.github.io

さて、C++に移植するか、ということで成果物がこちら。

github.com

データ取得機能はなく（思ったより湿度が不正確だったので作るモチベを失った）、時刻をラズパイと同期するだけの簡単なプログラムである。 これでめんどくさい、スマホを使った同期から解放される。 なお、BluetoothのCライブラリとしてはgattlibを使っている。 ライブラリを使う際にハマったのが、スキャンした直後でないと機器に接続できないということ。 MACアドレスさえわかれば接続できるでしょ、と思っていたのだが、どうやら違ったようだ（バグじゃないよな？）。 というわけで上記プログラムでは、一度接続してみてダメだったらスキャン＆再接続する、という微妙なロジックになっている。

最近色々C++でプログラムを書いて思ったのは、C++はCの資産をシームレスに使えるのが便利だ、ということ。 Python、Golang、Rustだとひと手間必要だし。 加えて、C++17ならかなりモダンな感じで書けるのがよい。 パッケージ管理がないとか、ビルドファイル書くのがめんどいとか、標準ライブラリのカバー範囲が狭いとか、テンプレートに関するエラーがわかりづらいとか、辛い点は色々あるのだけれど。。

追記

MHO-C303と同じセンサを搭載したENV IIで測定した湿度は、EMPEXのものとほぼ同等の値だった。 やはりXiaomiの値がズレている？何か変なロジックでも入っているのだろうか。

*1:実はMHO-C401も買った。完全に無駄だったが。

マイクラっぽい地図ビューワを作った

Programming

https://raw.githubusercontent.com/estshorter/VoxelMapViewer/images/screnshot.png

↑こんなんを作った。

ここ↓からグリグリ動かせる

https://estshorter.github.io/VoxelMapViewer/

業務でマイクラ風3次元地図ビューワが欲しくなったのだが、Matplotlibのvoxels()は重すぎて使用に耐えないし、どうしたものかと思っていたところ、以下のサイトを見つけたのが発端。

threejsfundamentals.org

Three.jsを使った上記サイトの実装では、ボクセル描画に関して見えないところで手を抜いて高速化している。 これをそのまま使ってもよかったのだが、どうせなら少し改良してみるかと思って、worldの構成要素を立方体（ボクセル）ではなく直方体にした。 なので、最初の画像は、高さの違う細長いブロックをXY平面に敷き詰めたworldとなっており、上記サイトのものより汎用性は落ちるが処理自体は軽くなっているはず（計測はしていない）。 ソースは以下のリポジトリに置いた。

github.com

ちなみに、画像上の赤い線はA*で求めた最短経路である。A*はTypeScriptではなくC++で実装しており、A*の出力結果をTypeScript側で読み込み描画している。

github.com

www.redblobgames.com

なお、ノイズを用いた地図の生成には以下のサイトを参考にした。

www.redblobgames.com

CO2Mini用C++ライブラリを作った

Programming C++

今年はリモートワークの機会が多かったのだが、どうもパフォーマンスが上がらず困っていた。そこで、CO2Mini（OEM元はこれ）を購入して室内環境の改善を図った。

www.monotaro.com

CO2MiniをPCに接続するとHIDデバイスとして認識され、暗号化を解けば温度・CO2濃度が取得できる。先人によって様々な言語に対してライブラリがつくられているので（Python実装、JS実装、C実装など）、手軽にCO2濃度を計測できるのだ*1。さて、私のマイブーム言語はC++であり、なぜかC++にはCO2Miniのライブラリがなかったので自分で作ってみた。ライブラリといっても250行程度の簡単なものである。 Windows、Linuxで動作を確認している。

github.com

READMEにも書いてある通り、以下のような感じで使うことができる。

#include <chrono>
#include <iostream>
#include <stdexcept>

#include "co2mon.hpp"

int main(int argc, char* argv[]) {
    using namespace std::literals::chrono_literals;

    co2meter::Co2meter dev;
    try {
        dev.Open();
    }
    catch (std::exception& e) {
        std::cerr << e.what() << std::endl;
        return 1;
    }
    // non-blocking monitoring
    dev.StartMonitoring(4s); // pass monitoring cycle
    while(true){
        std::this_thread::sleep_for(30s);
        auto temp = dev.GetTemp();
        auto co2 = dev.GetCo2();
        if (temp) std::cout << "TMP: " << temp.value().value << std::endl;
        if (co2)  std::cout << "CO2: " <<  co2.value().value << std::endl;
    }
    dev.StopMonitoring();
    co2meter::Co2meter::Exit();
}

StartMonitoring()に与えるデータ取得周期を長くしすぎると、何故か応答性が悪くなるので注意。内部でバッファリングされているのかもしれない。

私は、以下のリポジトリのソフトウェアをRaspberry Pi 4 (Ubuntu 20.10)上で運用している。

github.com

CO2Mini本体のディスプレイは温度とCO2を交互に表示しており少し見づらいため、ブラウザから両者を一度に見られるインジケータをVue.js + Bulmaで作った。 換気をするときにはこの表示を見ている。 https://raw.githubusercontent.com/estshorter/co2logger/images/display.png

グラフはAmbientから見られるようにした。 cpp-httplibでhttp通信をしている。 https://raw.githubusercontent.com/estshorter/co2logger/images/ambient.png

で、在宅環境はどうなったかというと、CO2濃度に気を付け始めて以来すこぶる調子が良い。ただ、DELL製32インチモニタ（U3219Q）に買い替えた影響が大きいかもしれない。金は正義。。。

*1:最近のファームウェアは暗号化しないようだが、2020/11に購入した私のものは暗号化されていた。

Bruno Wen-li は誰か / Who is Bruno Wen-li?

Misc

まとめ

Bruno Wen-li = 市川淳 (Jun Ichikawa)　（敬称略）

本文

Bruno Wen-liという人名を聞いてピンとくる人は何人いるだろうか。サントラに興味がある人以外はわからないと思うが、アニメ版ヨスガノソラのBGM作曲者の一人である。ヨスガノソラというとその内容（エロ表現、近親相姦）のみが注目されがちであるが、実はBGMが素晴らしい。下記のように、様々なサイトで高く評価されている。個人的には、ノスタルジックかつ哀愁漂うメロディーにとても惹かれる。

thanksyunno.blog.fc2.com

3kmusic.jugem.jp

haitubo.hatenablog.com

Amazon、moraなどの配信サイトから購入できるので、気に入ったら購入してほしい。Spotifyでの配信はされていない。

ヨスガノソラ　オリジナルサウンドトラック／オープニング主題歌＆番組挿入歌アルバム

ヨスガノソラ　オリジナルサウンドトラック／オープニング主題歌＆番組挿入歌アルバム

  • 発売日: 2016/09/21
  • メディア: MP3 ダウンロード

mora | ヨスガノソラ オリジナルサウンドトラック/オープニング主題歌&番組挿入歌アルバム

さて、ここからが本題。Brunoさんはヨスガノソラ以外一つも作品を手がけておらず、誰かの別名義ではと囁かれていた。私もその正体がずっと気になっていたのだが、答えを遂に得たのでここに報告する。

経緯としては「未確認で進行形」のアニメを見ていたときふと、この作曲者（市川淳さん）って他にどんな作品に関わっているのだろうと思い、同氏のWikipediaを見たことが始まりだった。すると作品リストの項目に「ヨスガノソラ（Bruno Wen-li名義）」としれっと出典なく書いてあるではないか。そこで「市川淳 "Bruno Wen-li"」とググって他にソースがないか探してみるも、出てきたのは2014年の5chの書き込み（これ書き込んだ人は何者なんだ…）と、VGMdbの本人ページ（出典なし）のみ。これはもう本人に直接聞いてみるしかないなと思い、Twitterで以下のように問い合わせてみた。

回答：

というわけで、Bruno Wen-li = 市川淳 との答えを得た。ついに表名義を知れた興奮から、回答を頂いた夜はなかなか寝付けなかった。

ちなみに、市川さん作曲の「未確認で進行形」のサントラもオススメである。

未確認で伴奏形 O.S.T.

未確認で伴奏形 O.S.T.