工学的アプローチですと大抵は自然現象を素性の良い関数のだけでモデル化を構成して取り扱う事情もあり、実数や関数の距離について厳密に追いかけなくても困らないことが多いです。

しかし、最適制御では関数の性質を扱う前に、$ L^p $ 空間、はたまた $ \Vert \cdot \Vert_{\infty} $ さらに、ハーディ空間やら $ H^p $ノルムやら出てきますし、 MPCでも定義域の表現に集合の記法を使うこともあります。 こういった言葉を目にしながら制御に取り組むと数学的表記でいちいち思考が途切れてしまいます。 今でこそ、特に困らず流して読めますが、やはり体系立てて身につけてこなかったので穴があり、 しばらく前に話題になったような宇宙際タイヒミュラー理論の雰囲気を今後の人生で楽しむには全くの力不足ですし、 ここまで来るのに遠回りしてきたように思います。

また、数年前、TLに数学から民間企業に移られた方の作った「現代数学解説(入門)」が流れてきたのを拝見して視界が開ける思いをしたのですが、それ以降そのままになっていました。

（もう一度読みたいので心当たりがあったら教えてほしいです）

こういった伏線を回収しようと(抽象)数学に入門することにしました。

集合・位相・圏 - 数学の言葉への最短コース

• 作者: 原 啓介
• 出版社/メーカー: 講談社
• 発売日: 2020/1/26

数学屋さんからは異論もあるでしょうが、私が以前見かけた「現代数学解説」のスライドと比べると数学的定義を追いかけている分、より厳密な議論で構成されています。 絵本よりも厳密で、プロの本よりも易しく学びたい方にはピッタリです。

当書で「圏」の入り口まで連れて行ってもらえると、計算機科学のための圏論の基礎の基礎といった資料をまず苦しまずに眺めることができるようになります。

同様に制御にかかわる高度な(?)数学的表現をWikipediaなどで追いかけたとき(例としてハーディ空間 - Wikipedia)、議論が何を話しているのか要略して解釈できるよう「言葉」の解像度をあげるためにきっかけの一冊にもなると思います。

最近発売された柔らかめの本らしくKindleのリフローに対応した電子書籍の扱いもあります。残念ながらKindle版だと脚注と本文の距離があり、補注を追いかけようとするとかなり辛みがありますので、 原先生の議論を一回で深くまで理解されたい方は紙の本を、 何度も読み返しても良いという方はKindleを選ばれると良いかなと思います。

PID制御のチューニング手法にはZiegler Nicholsの限界感度法はじめ様々な手法があるが、これらの手法がどのように異なるか簡単なケースについて調べた。

前提

むだ時間あり1次遅れ系(First-order plus deadtime; FOPDT)のPI制御を対象にチューニング手法のパラメータ比較。特に

\[ \rm{P}(s) = \frac{1}{1+\tau s}e^{-Ls} \]

の $ \tau = 1, L=1 $ のケースについてチューニング手法の結果をプロットする。システムをブロック線図で表すと下図。

ベンチマーク環境

現代制御的な手法はむだ時間なしの制御対象を

\[ \rm{\tilde{P}}(s) = \frac{1}{1+\tau s} \]

と表し、モデル化誤差が含まれたプラント$ \tilde{P} $に対してチューニングした上でロバスト性に期待する前提もおいてみる。

モデル化「誤差」が含まれたシステム

パラメータ観測用おもちゃ

複素関数解析を組み合わせて安定条件を解くむだ時間を持つ1次遅れ系のPI制御パラメータの解析的な分析で求めた安定領域(下図 緑の領域)を示したおもちゃがあるので遊んでいってください。 シミュレーションのタイムステップと積分誤差で境界より少し内側で振動し、境界上で発散するけど許して。

制御応答(左)と制御パラメータ(右)

黒円をマウスでグリグリしてみてね

結果とウンチク

古典制御については過去エントリーの二番煎じ。図に書き加えた $ r, q $ は最適制御の重み付けパラメータ。1状態1入力システムの最適制御に導出経緯を掲載した。

パラメータの比較

$ H_\infty $ ノルム

モデル化誤差を含んだプラントのシステムの相補感度関数 $ \frac{\tilde{P}C}{1+\tilde{P}C} $ すなわち

\[ G(s) = \frac{K_p s + K_i}{s^2 + (K_p + 1)s + K_i} \]

をロバスト制御の考えに基づいて $ \vert G \vert_\infty < 1 $ となる $ K_p, K_i $ の領域を求めた。システムの$ L_2 $ゲイン

\[ \vert G(\omega j) \vert = \sqrt{ \frac{K_p^2 \omega^2 + K_i^2}{(K_i - \omega^2)^2 + (K_p + 1)^2 \omega^2} } \]

の最大値が$ \omega \in \mathbb{R}^+ $ によらず1未満であるのでsolve[{(a^2*x+b^2)/((b-x)^2+(a+1)^2*x)<1, a>0, b>0,x>0}] - Wolfram|Alphaから

\[ 0 < K_i \leq K_p + \frac{1}{2} \]

このままだとプラントの前提知識がうすすぎるのでむだ時間をパデ近似で埋めてノルムを求める方法もある。

先日、線分交差判定の方法をまとめました。こういった当たり判定を多用する物理エンジンでは、Sweep and Pruneアルゴリズムなどを用いて当たり判定の枝刈りをしたり、そもそもネイティブPythonではなく、C++で実装して高速化するなどしています。

そういった2D向け物理エンジンでPython向けのインタフェースがあるものにBox2DやChipmunkがありますが、OpenAI GymのBipedalWalkerではBox2Dが採用されています。このBox2Dがネイティブ実装と比較してどのくらいの性能が出るのか比較しました。

Sweep and Pruneアルゴリズムの概念図

実行環境

• python 3.8
• Box2D==2.3.10
• numba==0.54.1

ベンチ環境

大量の隔壁に向けてLidarを照射し、距離を測ります。ネイティブ実装では純粋に"壁の数×レーザー数"回だけの当たり判定処理を行います。アルゴリズムについては先日のエントリーをご参照ください。

ベンチマーク環境

結果

ノートをこちらで公開。

縦10x横1000枚の隔壁にレーザーを照射した結果がこちらで、Box2Dが圧倒的でした。Pythonインタフェースのドキュメントが充実していないので使うのがしんどいですが、今度から手組み辞めようかなと思うくらい速い…。

	方法	結果[msec]	指数
1	Box2D	0.037	1.0
2	Numba実装	202.878	5,483.0
3	Native Python実装	549.742	14,857.9

そして、試した範囲内では理論通りO(n)になっているっぽい(※ Box2Dのグラフは細かく読んでない)

実行結果

Box2D実装メモ

Box2D Pythonインタフェースのドキュメントは体系的なものがないため毎度困るので実装メモを残す。

シミュレーション世界の初期化

import Box2D
from Box2D.b2 import (edgeShape, fixtureDef)

world = Box2D.b2World()
# v2.1.2でedgeShapeが導入された[1]
walls_2d = [
  fixtureDef(shape = edgeShape(
    vertices= [(w[0][0], w[0][1]), (w[1][0], w[1][1])])
  ) for w in walls
]
world.CreateStaticBody(fixtures = walls_2d)
# 上記とこれに大きな違いはない
# for f in walls_2d:
#   world.CreateStaticBody(fixtures = f)
class LidarCallback(Box2D.b2.rayCastCallback):
    # BipedalWalkerでは基底クラスに含まれている扱いだが明示的に宣言必要
    p1 = [0., 0.]
    p2 = [0., 0.]
    fraction = 1.0
    def ReportFixture(self, fixture, point, normal, fraction):
        self.p2 = point
        self.fraction = fraction
        # ReportFixtureの戻り値で振る舞いが変わる
        # [2] b2RayCastCallback Class Referenceを参照
        return fraction
lidar = [LidarCallback() for _ in range(ray_n)]

1. Fixtures - Box2D tutorials - iforce2d
2. Box2D: b2RayCastCallback Class Reference

RayCast

Lidarのレーザー起点と有効圏で初期化してRayCastする。当たり判定に伴う処理はReportFixtureで定義済み。

import numpy as np
# ロボットの定義(再掲)
o = np.array([50., 150.]) # ロボット位置
ray_l = 200 # Lidar有効距離
ray_n = 9   # Lidar線数
ray_angles = [math.pi/2 * i / (ray_n - 1.) - math.pi/4 for i in range(ray_n)]
ray_ends = o + ray_l * np.array([[math.cos(theta), math.sin(theta)] for theta in ray_angles])

for l, p in zip(lidar, ray_ends):
    l.p1 = o
    l.p2 = p
    world.RayCast(l, l.p1, l.p2)

シミュレーション内のオブジェクト取得

C++/Javaのインタフェースでは定義されているworld.getBodies()はPythonでは提供されておらず、world.bodiesでアクセスする。

import matplotlib.pyplot as plt
for b in world.bodies:
    for f in b.fixtures:
        w = f.shape.vertices
        # 隔壁が2点で定義されている前提の描画処理
        plt.plot([w[0][0], w[1][0]], [w[0][1], w[1][1]], color = 'black')
for l in lidar_b2d:
    plt.plot([l.p1[0], l.p2[0]], [l.p1[1], l.p2[1]], color = 'blue')
plt.show()

参考

• gym/bipedal_walker.py at master · openai/gym
• Box2D: Overview
• erincatto/box2d: Box2D is a 2D physics engine for games
• slembcke/Chipmunk2D: A fast and lightweight 2D game physics library.