Unityで太陽系の簡易シミュレーション

過去2回に渡ってUnityで万有引力による惑星の運動を実装してきましたが、それらを応用して太陽系を製作しました。UnityプロジェクトをGithubで公開します。シミュレーションというには簡易過ぎると思いますが、ご自由に活用ください。時間が取れれば製作過程についても少しずつ書いていきたいと思います。

■Githubリポジトリ
ronron-gh / unity-space

■過去記事
・Unityで万有引力による運動の確認 - mizu-mi blog
・Unityで万有引力による運動の確認(2) ～月も回してみた - mizu-mi blog

概要

各惑星の自転と公転をUnityの物理エンジンでシミュレーションしています。シンプルに作りたかったため、惑星の大きさや軌道は現実に忠実ではありません。

◎こだわった点
・各惑星の運動のスクリプトは共通化した。
・月も実装した（衛星用のスクリプトを別途作成）。
・各惑星にはテクスチャを貼り付けてリアルな見た目にした（NASAが公開するフラット画像を使用）。また、土星の輪はBlenderで作成した。
・光は太陽の位置から全方位に照射。
・スマホのジャイロでカメラを動かせる（カメラ用スクリプト内の＃ENABLE_GYROを有効化してスマホ向けにビルドすることで使用可能。※Androidのみで確認）。

開発に用いたUnityバージョン

2018.2.15f1

使い方

GithubからcloneしたフォルダをUnityで開いてください。Assets->Scenes->SampleScene をHierarchyにドラッグ＆ドロップすると、デバッグ実行できる状態になります。

また、demoフォルダに、Windows向けにビルドしたexeファイル（Windows10で動作確認済み）と、android用のインストーラ(.apk)が入っています。

テクスチャ画像の配布元

以下のサイトから集めましたが、配布元はNASAのようです。NASAの画像は通常著作権フリーのようですが、商用利用は自己責任でお願いします。

・地球
　Blue Marble Next Generation

・水星、金星、火星、木星、天王星、海王星
　ＣＧ講座POV-Ray地球

・土星、土星の輪
　JHT's Planetary Pixel Emporium

・月
　CGで月を作る | CGBeginner

ライセンスについて

著作権を主張するほどのものではありませんが、とりあえずMITライセンスとしています。ご自由にご利用ください。

ニューラルネットワークの表現力を確かめる

こちらで、「ニューラルネットワークは普遍性定理により任意の関数を表現できる」ということに触れました。今回は、実際に複雑な非線形の関数をニューラルネットワークに学習させ、表現力を確認したいと思います。

確認のために作成したプログラムはGithubに置きます。フレームワークはchainerを利用しています（chainerはメジャーアップデートが終了してしまいましたね。勉強目的で使っているので気にしないことにします^^;）。
ronron-gh / neural-regress-sinc

■普遍性定理について
以下のリンクで直感的にイメージできる形で説明されています（「ニューラルネットワークと深層学習」というオンライン書籍の翻訳のようです）。
ニューラルネットワークが任意の関数を表現できることの視覚的証明

普遍性定理により、中間層が1層あればあらゆる関数を近似できるとされています。層がディープだから何でも近似できるのだと思っていたので、ちょっと衝撃です(^^;

ニューラルネットワークに学習させる関数は次のグラフに示すsinc関数とします。

中間層1層の場合

まずは、普遍性定理を確認するために、中間層を1層としました。中間層のニューロン数は500としました。ソースコードはリポジトリのnn_regress_sinc.pyを参照ください。

学習結果は次のグラフの通りです。やや歪ですが、sinc関数に近い結果が得られました（ニューロン数やエポック数などのハイパーパラメータを調整すれば、もっと良い結果が得られるかもしれませんが、ここで力尽きました^^;）。

中間層2層の場合

次に、中間層を2層に増やして学習させてみました。ニューロン数は1層目2層目共に100としました。ソースコードはリポジトリのnn_regress_sinc_mid2.pyを参照ください。

学習結果は次のグラフの通りです。中間層1層のときよりもニューロン数が少ないにもかかわらず、sinc関数をより忠実に表現できています。

まとめ

中間層が1層しかないニューラルネットワークでも、普遍性定理によりsinc関数のような複雑な関数を表現できることを、実際にプログラムを組んで確認しました。また、中間層を2層に増やすと、1層のときよりも少ないニューロン数でより良い学習結果が得られることを確認しました。

ニューラルネットワークと最小二乗法の違いについて考えてみた

ニューラルネットワークと最小二乗法は、”損失関数（二乗和誤差）が最小となる重みを探す”という点で似ているように感じ、違いをはっきりと理解したいと思ったため、両者の違いについて自分なりに考察してみました（機械学習初心者なので間違いがあるかもしれません）。

考察のために作成したプログラムはGithubに置きます。フレームワークはchainerを利用しています（chainerはメジャーアップデートが終了してしまいましたね。勉強目的で使っているので気にしないことにします^^;）。
ronron-gh / neural-regress

ニューラルネットワークで回帰問題を解く場合

次のグラフに示すサンプルデータの近似関数を導出する回帰問題を考えます。

まずは、最小二乗法は使わずに、ニューラルネットワークで解いてみます。非線形の問題なので次の図のように中間層を1層設けました。

余談ですが、中間層が1層だけでも、ニューラルネットワークはあらゆる非線形関数を表現できます。これを普遍性定理というようです。以下のWebページで直感的にイメージできる形で説明されています（「ニューラルネットワークと深層学習」というオンライン書籍の翻訳のようです）。
ニューラルネットワークと深層学習
それでも、層を深くしたディープラーニングが注目されているのにはいろいろと訳があるようです（自分も勉強中のため、これ以上はもっと詳しくなってから書きます^^;）

結果は次のグラフのようになりました。中間層のニューロン数が増えると表現力が増していきます。

中間層のニューロン数：5

中間層のニューロン数：10

中間層のニューロン数：50

最小二乗法で回帰問題を解く場合

それでは、最小二乗法で解く場合を考えます。最小二乗法ではまずどんな関数で近似するかを決める必要があります。ここでは次のような三次関数で近似することにします。

　f(x) = a₃x³ + a₂x² + a₁x + a₀

最小二乗法ではサンプルデータとf(x)の二乗和誤差が最小となる係数a₀～a₃を求めることになります。表計算ソフトで最小二乗法を実行した結果を次のグラフに示します。

最小二乗法のアルゴリズムについては説明しませんが、自分が参考にしているWebページを紹介させていただきます。（Webで二次以上の関数の最小二乗法について説明しているのは、自分が知っている中ではここだけです）。
 ITエンジニアの数学 最小二乗法

このように、最小二乗法ではまず三次関数で近似できると仮定していますが、ニューラルネットワークではそのような仮定はしないため、両者は回帰問題へのアプローチの仕方が違います。

ニューラルネットワークで最小二乗法と同じアプローチを取ってみる

もう少し踏み込んで、今度はニューラルネットワークで最小二乗法とまったく同じアプローチを取るとどうなるかを考えてみました。それが次のニューラルネットワークです。

このように、最小二乗法は中間層なしで表現できました。結果は次のグラフの通り、先の表計算ソフトと同じような結果が得られました。このようにニューラルネットワークで表現すると、最小二乗法は実は（三次関数近似であっても）中間層なしで表現できるため、線形の問題を解いているに過ぎないということがわかります。

ニューラルネットワークで最小二乗法を解いた結果。重みは次の通り算出された(表計算ソフトの最小二乗法にほぼ一致)。 w3 = 0.298, w2 = 0.527, w1 = -4.010, b = 4.83

念のため補足しますが、ニューラルネットワークは複雑な問題に対応するために、損失関数（二乗和誤差）が最小となる重みの算出に勾配降下法を用いているため、計算時間は最小二乗法と比べ物にならないくらい大きくなります。一方、最小二乗法は解析的に二乗和誤差関数を偏微分して、勾配=0となる重みを求めるため高速です（最小二乗法の詳細は先に紹介したリンクを参照ください）。

まとめ

最小二乗法とニューラルネットワークの違いについて、自分なりに考察しました。最小二乗法は、ニューラルネットワークで解くことを考えると、線形の回帰問題を解くことに帰着します。一方ニューラルネットワークは普遍性定理により、どんなに複雑な非線形関数でも表現できるというポテンシャルを秘めています。
ゼロから独学で学び始めると、このような基本的な部分の理解にたどり着くまでに結構な時間がかかったため、記事にしてみようと思いました（どこか間違っていたらすみません）。

Unityで万有引力による運動の確認(2) ～月も回してみた

前回、地球が太陽から受ける万有引力によって楕円軌道で運動することをUnityで再現しましたが、今回は、さらに月を加えてみたいと思います。そんなの、月と地球の間に万有引力を発生させればいいじゃないかと思うかもしれませんが（自分も最初はそうしました）、実は太陽による万有引力も重要な役割を持っています。今回もパラメータやC#スクリプトを公開します。

月の公転の条件

冒頭で述べた通り、太陽の周りを公転している地球の周りを月が公転するためには、月に対して地球による万有引力と太陽による万有引力の2つを加える必要があります（前回同様に単純化したいので、地球と太陽には万有引力は与えません）。

月も、地球と一緒に太陽の周りを公転していると考えると、太陽による万有引力を合成するのは自然な考えですね。

月の公転を実装

太陽、地球、月それぞれのパラメータと、C#スクリプトを示します。

パラメータ

■万有引力定数
・太陽　Gs = 6.67E-06 Nm^2/kg^2  (実際の1E+05倍)
・地球　Ge = 6.67 Nm^2/kg^2          (実際の1E+11倍)

　※スケールがでたらめなので、万有引力定数は太陽に引かれる場合と
　　地球に引かれる場合とで分ける。

■質量
・太陽　Ms = 2E+07 kg
・地球　Me = 50 kg
・月　　Mm = 1 kg

■距離
・太陽と地球の距離 Rse = 20 m
・地球と月の距離 　Rem   = 5 m

■初速
・地球 Ve = sqrt(Gs×Ms / Rse) = 2.582634 m/s (太陽に対して垂直方向に与える)
・月　 Vm = sqrt(Ge×Me / Rem) = 8.167007 m/s (地球に対して垂直に与える)

　※等速円運動となるvを求める。(式の導出は前回記事参照)

Unityコンポーネント設定

■太陽

コンポーネント			設定	備考
Transform	Position	X	0 m
		Y	0 m
		Z	0 m
Rigidbody	Mass		2E+07 kg
	Use Gravity		チェックなし

■地球

コンポーネント			設定	備考
Transform	Position	X	20 m
		Y	0 m
		Z	0 m
Rigidbody	Mass		50 kg
	Use Gravity		チェックなし
C# Script	初速（Z方向）		2.582634 m/s (=Ve)	C#スクリプトのStart（）で設定する。
	力		太陽から受ける万有引力を計算して設定	C#スクリプトのFixed Update（）で設定する。

■月

コンポーネント			設定	備考
Transform	Position	X	25 m
		Y	0 m
		Z	0 m
Rigidbody	Mass		1 kg
	Use Gravity		チェックなし
C# Script	初速（Z方向）		10.74964 m/s (=Ve+Vm)	C#スクリプトのStart（）で設定する。 地球から見た相対速度がVmとなるようにする。
	力		太陽と地球それぞれから受ける万有引力を計算し、合成した結果を設定	C#スクリプトのFixed Update（）で設定する。

C#スクリプト

月のオブジェクトにアタッチするC#スクリプトを示します。
地球については前回記事を参照ください。

public GameObject のsunとearthはUnity画面のInspectorビューに現れるので、太陽と地球のオブジェクトをドラッグ＆ドロップしてください。万有引力や初速もpublicにするとInspectorビューから変更できるので便利かもしれません。

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class MoonBehaviourScript : MonoBehaviour {
    //オブジェクト
    public GameObject sun;
    public GameObject earth;
    //パラメータ
    private float G_s = 6.67E-6f;    //万有引力定数(太陽)
    private float G_e = 6.67f;       //万有引力定数(地球)
    private float initVelocityZ = 2.582634f + 8.167007f;    //初速
    private float M_m;      //月の質量
    private float M_s;    //太陽の質量
    private float M_e;    //地球の質量
    private float r_s;    //太陽と月の距離
    private float r_e;    //地球と月の距離
    private float f_s;    //太陽から受ける万有引力
    private float f_e;    //地球から受ける万有引力

    void Start()
    {
        //質量の情報を取得
        M_m = GetComponent<Rigidbody>().mass;
        M_s = sun.GetComponent<Rigidbody>().mass;
        M_e = earth.GetComponent<Rigidbody>().mass;

        //初速を与える
        Vector3 initVelocity = new Vector3(0f, 0f, initVelocityZ); ;
        GetComponent<Rigidbody>().velocity = initVelocity;
    }

    void FixedUpdate()
    {
        //月から太陽に向かうベクトルを計算
        Vector3 direction_s = sun.transform.position - GetComponent<Transform>().position;
        //月から地球に向かうベクトルを計算
        Vector3 direction_e = earth.transform.position - GetComponent<Transform>().position;

        //距離rを計算
        r_s = direction_s.magnitude;
        r_e = direction_e.magnitude;

        //単位ベクトルに変換（後で万有引力の方向として使う）
        direction_s.Normalize();
        direction_e.Normalize();

        //万有引力を計算
        f_s = G_s * M_m * M_s / (r_s * r_s);
        f_e = G_e * M_m * M_e / (r_e * r_e);

        //加速度を与える
        GetComponent<Rigidbody>().AddForce((f_s * direction_s) + (f_e * direction_e), ForceMode.Force);
    }
}

実行結果

うまくいくとこのようになります。ちなみに地球の公転が速すぎると月がどこかへ飛んで行ってしまいました。いろいろパラメータをいじってみると面白そうですね。

実行結果～太陽による万有引力を省略した場合

太陽の必要性に気付くまで、しばらくこの状態でハマってた^^;

スマホから取り込んだ動画をFFmpegで圧縮するバッチファイルを作ってみた

最近家族からスマホ（Android）内の写真・動画をPCに取り込んで欲しいと頼まれたのですが、全部で40GB程あり、自分のPCのHDD容量をかなり圧迫するサイズでした。この際RAID対応NASを導入しようかとも思いましたが、少々値が張るので、とりあえずサイズの大部分を占める動画（mp4）を圧縮してみることにしました。動画のサイズを見てみると、1分程度の動画で100MBを超えており、FFmpegで圧縮すると5分の1程度のサイズにすることができました。

バッチファイルの勉強を兼ねて、指定したフォルダ内のmp4ファイルをまとめて圧縮するバッチファイルを作成したので、せっかくなので公開したいと思います。

前置き① スマホからPCへの取り込みについて

詳細は割愛しますが、スマホをUSBケーブルでPCに接続し、スマホ側でUSB転送モードをMTPに設定すると、PCのエクスプローラからスマホ内（内臓ROM及びSDカード）のフォルダを参照できます。例ですが、SDカードのDCIM > 100ANDRO といった場所にスマホで撮影した写真・動画が格納されています。自分はフォルダごとPCのHDDにコピーしました（デジカメの付属ソフトのように日付毎にフォルダ分けしてくれたりするツールは、探しても見つからないですね）。

前置き② FFmpegについて

FFmpegのダウンロードや使い方はインターネットで検索するとたくさん出てくるので割愛します。自分はこちらのサイトを参考にさせていただきました。

creatorsblog.nijibox.jp

前置き③ Google フォトの「バックアップと同期」について

Androidスマホを使っている方なら、Googleフォトの「バックアップと同期」を使えば、わざわざ自前でバックアップを取る必要がないことにお気付きかと思います。「バックアップと同期」は、容量無制限でオンラインストレージに写真・動画をバックアップでき、スマホ本体からはバックアップ済みの写真・動画を削除できます。また、バックアップした写真・動画は別のデバイスからも見ることができます。自分も「バックアップと同期」を使っていますが、念のために自分のPCや外付けHDDにバックアップを取っています。

作成したバッチファイル

引数で指定したフォルダパスにある動画（mp4）をすべて圧縮するバッチファイルを作成しました。対象フォルダをバッチファイル（.bat）にドラッグ＆ドロップすることでも、フォルダパスを引数として渡せます。また、次回スマホから追加で取り込んだときに、同じファイルを圧縮しないようにするために、圧縮完了したmp4ファイルのファイル名をテキストファイル（compressed_file_list.txt）に保存し、次回はcompressed_file_list.txtにあるファイルは圧縮しない仕組みにしました。以下に、バッチファイル内容と解説を載せます。

※元ファイルは削除する（ゴミ箱にも残らない）ので、試す場合は念のためバックアップをお願いします。

バッチファイル内容

解説

@ECHO OFF setlocal enabledelayedexpansion IF "%1" == "" GOTO ERROR ECHO 入力されたフォルダパス：%1 CD %1 SET is_comp=0 FOR %%a IN (*.mp4) DO ( 　　FOR /F %%b IN (compressed_file_list.txt) DO ( 　　　　IF %%b==%%a ( 　　　　　　SET is_comp=1 　　　　) 　　) 　　IF !is_comp!==0 ( 　　　　ECHO %%aを圧縮 　　　　 RENAME %%a %%a.tmp 　　　　ffmpeg -i %%a.tmp -c:v libx264 -c:a aac %%a 　　　　DEL %%a.tmp 　　　　ECHO %%a >> compressed_file_list.txt 　　) ELSE ( 　　　　ECHO %%aは圧縮済み 　　) 　　SET is_comp=0 ) GOTO END :ERROR ECHO 入力エラー SET /P inp = "何かキーを押して終了" :END SET /P inp = "何かキーを押して終了"

環境変数（is_comp）の値を処理の途中で変更可能とするために遅延環境変数を使用する。 カレントディレクトリを入力されたフォルダパスに変更する。 フォルダ内のmp4ファイル名についてループ処理を行う。 compressed_file_list.txtにファイル名が存在する場合はフラグ(is_comp)を立てて圧縮処理をスキップする。 ファイル名を”ファイル名”.tmpにリネームし、ffmpegに入力する。ffmpegの出力は元ファイルの名前にすることで、圧縮前後でファイル名が変わらないようにする。 圧縮が終わったら、.tmpは削除し、元ファイルの名前をcompressed_file_list.txtに記録する。

補足

圧縮率が期待外れな場合の対処法

スマホの機種によっては、圧縮後のサイズがほとんど変わらない場合があるようです。その場合は、-crfオプションで動画の品質を指定すると圧縮後のサイズを調整できます。

    ffmpeg -i input.mp4 -crf 30 -c:v libx264 -c:a aac output.mp4

-crfオプションの値の範囲とデフォルト値はコーデックにより違うようです。小さい値ほど品質は良くなります。

mp4以外も圧縮できる

入力ファイルはmp4以外でも問題ありません。以下はmovを入力してmp4を得る場合。

    ffmpeg -i input.mov -c:v libx264 -c:a aac output.mp4

また、出力ファイルの形式も変更可能であり、指定した出力ファイル名の拡張子を見て自動で切り替えてくれます。

組込みソフト開発で学んだ知識まとめ

プロフィールの通り、自分は組込みソフトエンジニアであり、就職してからいつの間にか10年近く経ちました。10年目にしてはあまり自信はないのですが、思い起こせばいろいろと考えたり学んだりしてたので、ここでまとめてみたいと思います。学生の頃にマイコンのプログラミングは経験していましたが、就職してから3～4年目くらいまでは結構ついていくのに必死でした。その頃を思い出しながら書きたいと思います。

RepublicaによるPixabayからの画像

はじめに～組込みソフト開発ってどんなことするの？

組込み系と一口に言っても、システムの規模によって使う技術も様々ですが、自分は主にマイコンF/Wの開発に携わってきました（システムの規模は小さく、OSもなし）。どんなことをしているのかについては、最近たまたま同じような境遇の方のブログを見つけ、とても共感したので、自分で書き直すことはせずリンクを貼らせていただきます(^^)

nmmmk.hatenablog.com

組込みソフト開発で学んだ知識

MISRA-C

マイコンF/W開発のチームに入ってまず渡されたのが、MISRA-Cの本でした。

　『組込み開発者におくるMISRA‐C:2004―C言語利用の高信頼化ガイド』
　　MISRA‐C研究会 日本規格協会

MISRA-Cは自動車業界で生まれた、C言語で安全なプログラムを書くためのコーディング規約です。本を読んでみた感想として、ほとんどのルールがコンパイラの基本的な仕様に関するものであり、自動車業界に限らず知っておいたほうがよいと思える内容でした。

実際、MISRA-Cを学ぶまで知らなかったコンパイラ仕様もあるので、少し紹介します。

・汎整数拡張と平衡化
異なる型の間での演算を行う場合、暗黙的に同じ型に変換されますが、この暗黙の変換は汎整数拡張と平衡化により実行されます。まず汎整数拡張により、式中の整数型（char、short）はint型に昇格します。その後、平衡化により、より精度が高いオペランド（longやfloat、double）に合わせられます。自分は汎整数拡張のことを知らなかったので、今更C言語でこんな発見があるのかと、目から鱗でした。

・プロトタイプ宣言がない場合のコンパイル結果
プロトタイプ宣言がなくてもコンパイルは通るため、プロトタイプ宣言が抜けていないかのチェックは疎かになりがちですが、実はプログラムの動作結果に影響します。プロトタイプ宣言がない関数は、コンパイラが勝手に戻り値と引数をint型と解釈します。そのため、もし定義側の戻り値や引数が浮動小数点数だったとしても、動作時は小数点以下が丸められてしまいます。（コンパイルエラーにしてほしいですね。warningは出ますが、warningがいくつも放置されているプロジェクトだと見落としてしまうので。）

マイコンの仕組み

マイコンには学生の頃から触れていましたが、就職してから新たに学んだことも多いので、挙げていきたいと思います。なお、すべて以下のサイトで気軽に学ぶことができます（会社HPの技術コラムのようですが、かなり丁寧に書かれています）。

[技術コラム集]　組込みの門

・ヒープとスタック
学生のときは、サイズが大きな配列はローカル変数ではなくmallocで動的確保しましょうと習いましたが、メモリ領域におけるヒープとスタックの違いについて理解するとその理由がわかります。

・メモリのセクション
PCアプリのプログラミングではプログラムコードやグローバル変数、定数などがメモリのどこに配置されているか、あまり気にしないと思いますが、マイコンはメモリサイズが小さいため、すべて考えて配置する必要があります。　

・キャッシュ
マイコンでキャッシュを使う場合は、キャッシュの有効／無効をメモリ領域ごとに検討する必要があります。また、組込みではDMAなどを行った場合にキャッシュに最新の値が反映されず、不具合の原因になることがあります（コヒーレンシの問題）。

・スタートアップルーチン
main関数が呼ばれる前のマイコンの初期化処理であり、アセンブラで書かれています。大抵はマイコンのベンダーがサンプルを用意しているため、自分でゼロから作ることはありませんが、簡単なカスタマイズは自分でする必要があります。

デバッガ（ICE）

学生のときは、マイコンプログラムのデバッグはLEDやprintf(出力先はシリアル通信のターミナル)で行っていましたが、仕事ではICEによりIDE上でデバッグすることが多いです。ICEによって、Visual Studioのようにプログラムを1ステップずつ実行したり、メモリダンプしたりできます。また、CPUレジスタやペリフェラルのレジスタを参照・設定したり、キャッシュの内容を表示したり、スタックトレースできたりと、デバッグが大変はかどります。

また、IDEがマイコンのCPUやメモリをエミュレートする機能を持っている場合があり、開発初期に実機が入手できないといった状況でもデバッグを行えます(ペリフェラルはさすがにエミュレートできないため、主にアルゴリズムのデバッグや単体テストで活躍しました)。

科学技術計算

マイコンプログラムは基本的にC言語で開発するため、科学技術計算の標準的なライブラリはありません。そのため、客先から要求されたときは、まずは自力で実装方法を調べることになります。

過去経験したものに、最小二乗法があります。なんとなく原理は知っていましたが、自分で実装したことはなく、実装できたとしても処理速度や計算精度といった性能面で要求を満たせるのかが心配でした（もちろん、インターネットに転がっているプログラムでも同じ心配があります）。そんなとき、たまたま持っていた本が、『C言語による最新アルゴリズム事典（奥村 晴彦 著）』でした。久しぶりに本を開いてみると、最小二乗法について計算効率まで考慮したアルゴリズムがソースコード付きで載っており、本を買った時の期待以上に助けられました。もちろん、本に載っているからといっても使うときは自己責任ですが、自作やインターネットに比べたら安心感は段違いです。

高速フーリエ変換（FFT）を要求されたこともありますが、これはマイコンのベンダーがライブラリを提供していたため、自分で組まずに済みました（信号処理向けのマイコンだったので、当然といえば当然ですが）。ただ、コンパイルを成功させるまでに時間がかかったり、その後使いこなすまでにも時間がかかったりはします。

難しいアルゴリズムの実装はプレッシャーもありますが、動いた時の喜びも大きいので自分は好きです。

ロジカルライティング　(追記 2021/12)

組込みソフトとは直接関係はありませんが、仕事を進める上で役に立った知識です。メールやプレゼン資料で（もちろん、口頭でのやり取りでも）　相手に伝えたいことを正確に伝えるにはどのような文章構成にすればよいかを、体系的に示したものがロジカルライティングです（ロジカルシンキングが有名ですが、それを文書作成に展開したものです）。

いろいろな書籍や講座がありますが、自分が読んだのはこちらの書籍です。

『ロジカル・ライティング―論理的にわかりやすく書くスキル』

　照屋 華子(著)

全体を説明してから詳細を説明するとか、メールなら何を相手に依頼したいのかを初めに書くとか、ある程度経験を積んでいれば無意識にやっていることもありますが、なぜそのようにしているのかを体系的に理解できるので一読の価値はあると思います。

組込みソフト開発におけるキャリアアップ　(追記 2021/12)

この記事を最初に投稿してから早くも2年以上が経ちましたが、実は最近転職しました。

前職は受託開発の会社で、ハードウェアはすべて客先が製作し、ソフトウェアの要求仕様も基本的にはすべて客先から提示されます。学生時代のソフトウェア系の友人も皆、受託開発の会社に就職していたため、ソフトウェアは受託開発が当たり前だと思い込んでいました。

でも10年以上働いていると、次第に物足りなさを感じるようになりました。元々ハードウェアにも興味があり、もっと上流の、ハードウェアを含めた構想設計の段階から関わってみたいという思いが強くなりました。

また、全ての受託開発がそうとは限りませんが、前職では10年目くらいになるとプロジェクト管理を任されるようになり、チームメンバーの進捗管理や品質管理、コスト管理等がメインになるのが当たり前でした。

そんなこんなで、自分の年齢でも組込みソフト技術者として自社製品の開発に携われる会社はないのだろうかと、転職サイトで検索してみると、そのような会社は思っていたよりもたくさんあることがわかりました。（そこからは、転職サイトのエージェントからのアプローチを皮切りに、どんどん転職活動が進んでいきました。）

とはいえ、受託開発で経験を積み上げてきたからこそ転職が成功したのも事実です。また、新卒で今の会社に採用されたとしても、メカや電装など幅広い会社であるため、最初から組込みソフトをやらせてもらえたかはわかりません（そもそも、好きな会社を選べるほどの学歴や成績でもなかった^^;）。学生時代にあまり深く考えずに就職活動をしていたことを悔やんだこともありましたが、今となっては、転職で徐々にキャリアアップしていくのが自分には合っていたのではないかと感じています。

Unityでロボットアームの逆運動学

前回記事でロボットアームの順運動学をUnityで実装しましたので、今回は次のステップとして逆運動学を実装してみました。逆運動学は三角関数で解析的に解く方法と、ヤコビアンを利用して数値的に解く方法がありますが、今回は数値的に解いています。本記事ではアルゴリズムの概要を解説し、C#スクリプトのソースコードも載せますが、基礎から説明するととても長くなってしまうので、ある程度書籍などで学んだ知識を前提とします。私は梶田秀司編著「ヒューマノイドロボット」を参考にしています。

前回記事をご覧いただいた方へ
前回の順運動学は、Matrix4x4で同次変換行列を作り、回転と並行移動をまとめて計算しましたが、今回は同次変換行列を用いていません。前回のソースコードを利用できないことについてご了承願います。

ロボットアームの構造

ロボットアームの構造については前回記事「Unityでロボットアームの順運動学」に同じとします。逆運動学の計算で利用するパラメータを示します。

・pw_j ：ワールド座標系での位置
・p_j：ローカル座標系での位置（親リンク相対）
・Rw_j：ワールド座標系での姿勢
・R_j：ローカル座標系での姿勢（親リンク相対）
・angle_j：関節の回転角度
・a_j：関節軸ベクトル（親リンク相対）

逆運動学（数値解法）のアルゴリズム

逆運動学ではアーム先端の位置と姿勢を指定できますが、今回は位置のみにして簡略化しました（初心者なので^^;）。逆運動学を数値計算で解くアルゴリズムは次のようになります。

①逆運動学で求める、各関節の角度を並べたベクトルをangleとし、初期値を設定する。
②アーム先端の目標位置pw_{4_target}を指定する。
③順運動学で、angleから各関節の位置と姿勢を計算する。
④目標位置からの誤差Δpw₄を計算する。(Δpw₄ が十分に小さければ終了。)
　Δpw₄ = pw_{4_target} – pw₄

⑤Δpw₄を小さくする角度修正量Δangleを計算する。
⑥Δangleを関節角度に反映する。
　angle = angle + Δangle

⑦③から繰り返す。

 角度修正量Δangleをどのように求めるかが問題となりますが、ここでヤコビアンを使います。各関節を微小に変化（Δangle）させたときのアーム先端の位置の微小変化量（Δpw₄）をヤコビアンJを用いてΔpw₄ = J Δangleと表現できれば、Jの逆行列を使ってΔangleを求めることができます。

　Δangle = J^-1 Δpw₄

ヤコビアンの計算方法

ある関節が微小に回転したとき、アーム先端の位置がどのように変化するかを考えます。例えば、第2関節の微小回転角度をΔangle₂とすると、アーム先端の位置の微小変化量Δpw₄は次のように計算できます（×は外積）。

　Δpw₄ = Rw₂a₂×(pw₄ - pw₂)Δangle₂

全ての関節についてこれを計算し、足し合わせれば、Δpw₄を求めることができます。それを行列で表すと、

 

先に示した式Δpw₄ = J Δangleと比較すると、

ということになります。（Jは3x3行列となる。）

C#スクリプト

上記のアルゴリズムを実装したソースコードを載せます。一つ面倒なこととして、Unityには4x4行列しか用意されていないため、行列演算にはMath.Net Numericsというライブラリを使っています。このライブラリは.NETFrameworkのバージョン4.0以降ですが、Unityは3.5なので使うにはひと手間必要です。こちらの記事が参考になりました。

Unityで数値計算ライブラリMath.NET Numericsを使う方法 - arXiv探訪

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using MathNet.Numerics.LinearAlgebra.Single;    //Unityの行列は4x4しかないため、行列演算にはこのライブラリを用いる

public class RobotControllerScript : MonoBehaviour {

    public Transform axis1_tf, axis2_tf, axis3_tf;      //各軸のTransform(逆運動学の結果を出力してロボットを動かす)
    public Transform hand_transform;                    //アーム先端のTransform(逆運動学の結果確認用)
    public Vector3 pw4_target;                          //アーム先端の目標位置（Inspectorビューから入力）

    //---リンクのパラメータ---//
    private Vector3 pw1, pw2, pw3, pw4;    //ワールド座標系での位置
    private Quaternion Rw1, Rw2, Rw3;      //ワールド座標系での姿勢
    private Vector3 p1, p2, p3, p4;        //ローカル座標系での位置（親リンク相対）
    private Quaternion R1, R2, R3;         //ローカル座標系での姿勢（親リンク相対）
    private float angle1, angle2, angle3;  //関節角度
    private Vector3 a1, a2, a3;            //関節軸ベクトル（親リンク相対）

    private float lambda = 0.1f;    //逆運動学の収束を調整する係数（0～1で設定）

    //UI
    public UnityEngine.UI.Text transform_label;

    void Start () {
        //固定値のパラメータを設定
        p1 = new Vector3(0f, 1f, 0f);   //第1軸のローカル座標（親リンク相対）
        p2 = new Vector3(0f, 1f, 0f);   //第2軸のローカル座標（親リンク相対）
        p3 = new Vector3(0f, 3f, 0f);   //第3軸のローカル座標（親リンク相対）
        p4 = new Vector3(0f, 3f, 0f);   //アーム先端のローカル座標（親リンク相対）
        a1 = Vector3.up;                //第1軸の関節軸ベクトル（親リンク相対）
        a2 = Vector3.right;             //第2軸の関節軸ベクトル（親リンク相対）
        a3 = Vector3.right;             //第3軸の関節軸ベクトル（親リンク相対）

        //関節角度の初期値を設定
        //※0度だとヤコビアンの逆行列が存在しないため少し角度を付けておく
        angle1 = 10f;                   //第1軸の回転角度
        angle2 = 10f;                   //第2軸の回転角度
        angle3 = 10f;                   //第3軸の回転角度
    }

    void FixedUpdate () {

        for (int i = 0; i < 50; i++)
        {
            //// 順運動学 ////
            pw1 = p1;
            Rw1 = Quaternion.AngleAxis(angle1, a1);

            pw2 = Rw1 * p2 + pw1;
            R2 = Quaternion.AngleAxis(angle2, a2);
            Rw2 = Rw1 * R2;

            pw3 = Rw2 * p3 + pw2;
            R3 = Quaternion.AngleAxis(angle3, a3);
            Rw3 = Rw2 * R3;

            pw4 = Rw3 * p4 + pw3;


            //// ヤコビアンを作成 ////
            Vector3 j1 = Vector3.Cross(Rw1 * a1, pw4 - pw1);
            Vector3 j2 = Vector3.Cross(Rw2 * a2, pw4 - pw2);
            Vector3 j3 = Vector3.Cross(Rw3 * a3, pw4 - pw3);

            var J = DenseMatrix.OfArray(new float[,]     // ※Math.Netを使用
            {
                { j1.x, j2.x, j3.x },
                { j1.y, j2.y, j3.y },
                { j1.z, j2.z, j3.z }
            });

            //// 逆運動学 ////

            //アーム先端の目標位置からの誤差を計算
            Vector3 err = pw4_target - pw4;
            float err_norm = err.sqrMagnitude;

            //err_norm < 1E-5なら計算終了
            if (err_norm < 1E-5) break;

            // ※Jの逆行列と掛け算するためにerrをMath.Netの行列に入れなおす
            var err2 = DenseMatrix.OfArray(new float[,]
            {
                { err.x },
                { err.y },
                { err.z }
            });

            var d_angle = lambda * J.Inverse() * err2;

            //逆運動学で得た角度修正量を各軸に反映する(d_angleはラジアンなので度に変換)
            angle1 += d_angle[0, 0] * 180f / 3.14159f;
            angle2 += d_angle[1, 0] * 180f / 3.14159f;
            angle3 += d_angle[2, 0] * 180f / 3.14159f;
        }

        //逆運動学で計算した関節角度をロボットに反映
        Quaternion axis1_rot = Quaternion.AngleAxis(angle1, Vector3.up);
        Quaternion axis2_rot = Quaternion.AngleAxis(angle2, Vector3.right);
        Quaternion axis3_rot = Quaternion.AngleAxis(angle3, Vector3.right);

        axis1_tf.localRotation = axis1_rot;
        axis2_tf.localRotation = axis2_rot;
        axis3_tf.localRotation = axis3_rot;

        //UIにアーム先端の位置を表示
        transform_label.text = "Transform X:" + hand_transform.position.x.ToString("F2") +
                                " Y:" + hand_transform.position.y.ToString("F2") +
                                " Z:" + hand_transform.position.z.ToString("F2");


    }
}

実行結果

アーム先端の目標位置（pw4_target）をx = 3、y = 3、z = 3 として逆運動学を実行した結果を示します。初期状態は各関節の角度を10度としています。0度だとヤコビアンの逆行列が存在しないためです。画面にはアーム先端のTransform.positionを表示しており、指定した通りになっています。

また、関節角度がどのように収束しているのか気になったため、グラフ化してみました。ヤコビアンのおかげでうまく収束していることがわかります（第2関節（angle2）が一旦大きくマイナスに振れているのが気になりますが）。

追記（2020.01.24）

当ブログをtwitterで参照してくださり、6軸の位置・姿勢制御に拡張までされている方を発見したので紹介させていただきます。twitter活動はしていないのでこのような形で恐縮ですが、参照してくださりありがとうございました。とても励みになりました(^^)

こちらのソースを試させてもらって順・逆運動学テスト。
同じ教科書のヒューマノイドロボット（オーム社 梶田秀司 著）使ってる！
これはタイムリー！　4月の記事で、やはりこの教科書は絶版でもまだメジャーなんだなぁ。

Unityでロボットアームの逆運動学-mizu-mi bloghttps://t.co/G9tC5a7D6C pic.twitter.com/6AfqUo1xNU

— でべ (@devemin) 2019年12月10日

っしゃー！
逆運動学 (IK) のgithub MATLAB コードを unity に移植して実行できた！
６軸アームをIKで位置・姿勢の制御。嬉しい～
赤いオブジェ（位置・姿勢）をターゲットにして、アーム先端を動かしています。

参考：ヒューマノイドロボット（オーム社 梶田秀司 著）githubhttps://t.co/SmHO1mXrwt pic.twitter.com/r1K8zrGH4F

— でべ (@devemin) 2019年12月16日