数式処理ソフトウェアで立式した運動方程式をSimulinkで実行する

この記事は制御工学アドベントカレンダー15日目の記事です．

制御工学で避けては通れないのが運動方程式の立式です．
プラントモデリングではもちろんのこと，現代制御に代表されるモデルベース制御では制御ロジックを構築するために運動方程式（状態方程式）が必要になります．

運動方程式の立式というと，ほぼラグランジュ法（オイラー＝ラグランジュ方程式）一択のような気がしています．

いまいちこのあたりの用語を掴みきれていないのですが，汎函数と呼ばれる“函数を変数とした函数”を変分（汎函数での微分）することで運動方程式を求めるといった流れです．

たとえば，$f(t) = at$は函数と呼ばれ，$\frac{df(t)}{dt}$という操作が微分と呼ばれることはご存知のとおりです．
それと似たような関係で，$I[f(t)] = Af(t)$を汎函数と呼び，$\frac{dI[f(t)]}{df(t)}$を変分と呼ぶようです．

物体の運動の場合，最も必要エネルギの少ない方向に物体が運動するという原理があるらしく（最小作用の原理？），以下の形式で与えられます．

$$I[ y ] = \int^{b}_{a}f( x, y, y' ) dx$$

特に運動の場合は，右辺$f(x,y,y')$をラグランジアン$L$というエネルギ量で与えるようです．
合っている自信はないですが，以下の式のようなイメージです．（時間微分なのでドットで書いています）

$$I[ x ] = \int^{b}_{a}L( t, x, \dot{x} ) dt$$

これはある時刻$a$から$b$までのエネルギの合計となります．
上の式が最小になるようにすれば良いということで，$I[x]$の変分（微分）がゼロとなるような$x$（停留値）を求めれば良いことになります．

この解（停留値）は決まっていて，これをオイラー＝ラグランジュ方程式と呼びます．
$$\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{x}} \right) - \left( \frac{\partial L}{\partial x} \right) = 0$$

解の形式が決まっているのでシーケンシャルに解けて便利なのですが，非常に計算量が多いです．

前置きが長くなりましたがこれを数式処理ソフトウェアで楽にしようというのが今回のコンセプトです．

---

この記事では数式処理ソフトウェアとして以下を使用しています．

• wxMaxima（Windows版Maxima，無償）
• Symbolic Math Toolbox（MATLAB，有償）

恐らく運動方程式を求めたあとはMATLAB/SimulinkとかPyControlとかでシミュレーションすることになると思うのですが，結論としてはwxMaxima優位といった印象です．

---

今回，2重振り子を例として，Simulinkでのシミュレーションまでやってみようかと思います．

ラグランジュ法では物体の運動エネルギ，位置エネルギの関係を与えてあげることが必要になります．
逆に言うとエネルギまで出せてしまえばあとはソフトウェア任せにできます．

事前にやってみたところ，wxMaximaのインストールから結果確認まで45分程度で終わりました．

最適化函数を使って船を横移動させてみる

この記事は制御工学アドベントカレンダー9日目の記事です．
さいきん，仕事の方で船に関わることが多いので船の動きも面白いよ，と布教してみる．
船，とひとくちに言ってもタンカーなどの超大型船から客船，個人所有のクルーザー，漁船までいろいろ種類があります．
また，推進方法によっても船内にエンジンを積んだ船内機艇，エンジンごと取り外し可能な船外機艇，ジェット推進の艇などの種類があります．

ここではクルーザーなどに多く使われる下図のような船外機艇(Outboard motor)に焦点を当ててみます．

Mercury社のウェブページ
図の船は船外機が2つ取り付けられています（2機掛けと呼びます）が，最大で5つ取り付ける方もいるようです．
250馬力を5機掛け．パワーの面でも価格の面でも恐ろしいですね．

---

その船外機艇のなかでも面白いなあと思ったものがMercury社のSkyhookというJoystick操船システムです．
1分34秒あたりを見るとわかりやすいですが，着岸状態から船がほぼ真横に移動しています．

Mercury Joystick Piloting for Outboards
通常，船外機は何機掛けであろうと同じ方向を向かせて舵を切るのですが，この操船システムではそれぞれの船外機を独立させて反対方向に動かすことで横移動を実現しています．

このSkyhookを再現するため，最適化函数を使って遊んでみようというのがこの記事の趣旨です．
環境はMATLAB/Simulink R2015aです．

C/GMRESによる非線形モデル予測制御1

今回は珍しく2部構成で，本記事ではシステムの立式から評価函数の設定まで，次の記事では具体的な解き方について述べます．

第1部となる本記事ではC/GMRESの核となる部分には言及しないため，そちらについては2部以降を見ていただければと思います．

---

非線形モデル予測制御(Nonlinear Model Predictive Control; NMPC)とは，その名の通り”状態方程式が非線形なシステム”に対するモデル予測制御(Model Predictive Control; MPC)です．

最近流行りの制御則で，車の自動運転やロケットの再着陸システムにも使用されているそうです．

MyEnigmaさん
モデル予測制御(Model Predictive Control:MPC)の応用例

そもそも制約有りのMPCは非線形な制御系ですが，”非線形MPC”といった場合には“制御対象が非線形”な系を指すように思います．

線形のモデル予測制御については別途，記事に起こすつもりです．

通常，線形MPCでは，システムの表現に以下のような線形な状態方程式を使うのに対して，
$$\dot{x}(t) = Ax(t)+Bu(t)$$
非線形MPCでは状態方程式の形式を問わず使用でき，以下のような形式で表現されます．
$$\dot{x}(t)=f(x(t),u(t)) = \left[\begin{array}{c} x_{2}(t) \\ ax_{1}(t)+bx_{2}(t)u_{1}(t) \end{array}\right]$$

これはセミアクティブダンパの状態方程式ですが，式中に$x_{2}(t)u_{1}(t)$という“状態変数と操作量の積”が含まれています．少し考えるとわかりますが，線形の状態方程式に書き直すことはできません．

実際，このシステムを制御する際に単純に線形MPCを適用することはできません．
このような，非線形なシステムの制御を目的としたMPCが非線形MPCです．

非線形MPCの中にも種々ありますが，ここではC/GMRESによる非線形MPCについて取り上げます．

ただ，このシステムは減衰が操作量なのでダンパとして動作させること自体はPID等でも可能です．というよりずっとONにしておけばよいのです．しかし，できるだけ少ない操作量で最大の減衰を得よう，などと考え始めるとちゃんと制御する必要があります．

カルマンフィルタの基礎式を代数とベイズ定理から見る

カルマンフィルタの歴史ももう長いことと思いますが，各所で各人なりにまとめられているあたり，なかなか理解の難しいものなんだなあと感じています．
例に漏れず自分もカルマンフィルタの原理がいまいちピンとこなかったので自分なりにまとめてみました．

カルマンフィルタ(Kalman Filter)はフィルタと名付けられてはいますが，使われ方としてはフィルタよりも現代制御のオブザーバ（Observer; 状態推定器）に近いような気がします．
観測した情報から，状態ベクトルが真値に最も近くなるように推定するものです．

概念は以下のところでわかりやすく説明されています．

モノを作りたくなるブログ

シンプルなモデルとイラストでカルマンフィルタを直感的に理解してみる

ただ，数式として少し省略されているところもあったので，改めてまとめるに至りました．
基本的に，代数とベイズの定理のみから丁寧に導出していこうという感じです．

カルマンゲインの式の導出，予測値の算出，分散の更新は天下り的な部分も多いので，今回はこれらの導出に焦点を当ててみました．
統計的にみて，予測値と実プラントの誤差が最小となるように演繹的に求めていくとカルマンフィルタの式に行き着くあたり，統計学の応用先の代表例になっているだけのことはあるなあ，という感想です．

ここでは，状態空間と観測空間が同一の場合のみ扱います．

4ストロークエンジンの物理モデル

4ストロークのガソリンエンジンの物理モデリングをしてみました．
シミュレーションはMATLAB/Simulink，可視化ツールはProcessingです．

参考元はこのあたり．
双方ともに高知工科大学の磯村先生の文献です．

4 サイクルガソリンエンジンの動特性について

並列計算機によるエンジンのシミュレータ開発

---

以前にも倒立振子とかパーティクルフィルタでのグリッドマップマッチングとか2リンク平面マニピュレータとかのモデリングをやってきたのですが，そのどれとも色が違うなあ，ということで．

例によって，ラグランジュ法で運動方程式を立てます．
このあたりは参考文献の上の方が詳しいかな．

シンプルな外乱オブザーバ

外乱オブザーバ(Disturbance observer; DOB)はフィードバックシステムに外乱が加わるような状況で，比較的簡単に外乱を除去し安定化できる強力なツールですが，ウェブにはあまり記述が無いようだったので書いてみました．

• 外乱オブザーバとは？
• なぜ外乱オブザーバが必要なのか
• 基本原理
• 実際の外乱オブザーバ（ノミナルモデル，ローパスフィルタ）

---

そもそも外乱オブザーバとは？

“外乱”でググってみると，以下のような説明が得られます．

ある通信系に、所定の信号系以外から加わり妨害となる信号。 --- 大辞林

大体内容としては間違いないですが，今回扱うのは”通信系”に限らず，制御システムに対してはたらく，妨害となる信号くらいに留めておきます．

外乱オブザーバとは，このような外乱を除去するために，制御系のコントローラに組み込む補償器のひとつです．

---

なぜ外乱オブザーバが必要なのか

次のような1次遅れ系のフィードバックループを考えてみます．

ここで，$C(s)$はコントローラ，$G(s)$は制御対象です．また，目標値を$r$，出力を$y$で表します．
外乱というと，図の1, 2, 3番の箇所に加わることが想定されますが，本稿では1について記述します．

また，以下の理由で2, 3番を省いています．

2番の部分に入力される外乱は，フィードバックループを構成する主な目的のひとつのため，比較的簡単に除去（補償）できる．
3番については，システム自体は目標値通り動作している（外乱によって我々使用者が目標値どおりに動作しているか確認できないだけの状態）ため．

1番は，コントローラは制御対象を最適に動作させるための操作量を正常に生成できているが，それが外乱によって変えられ制御対象に入力されてしまっている状態です．

つまり，PIDなどのコントローラを適切に設計しても，設計通りの操作量が制御対象に入力されないということです．

今回はこの外乱を除去する外乱オブザーバについて述べていきます．

ここでは，"最適な操作量"は設計時に想定している操作量を指します．
最適レギュレータなどの最適とは根本的に異なる内容ですのであしからず．

ニューラルネットで雑な手書き文字認識

以前にもニューラルネットについて言及しましたが，かなり原初的な内容ばかりだったので今回は簡単な応用について．

C言語でニューラルネット（関数近似）

C言語でニューラルネット（2次元平面での分離）

といっても，修士の授業で作った内容そのままなので，いざ実装してみると動かないかも．

ニューラルネットの応用先として，真っ先に上がるのがパターン認識の分野でしょう．

とくに文字認識の分野ではすでにかなりの数，実用化されていると聞きます．

ここでは，最も単純な文字の認識について述べていきたいと思います．

具体的には以下の2点．

• ニューラルネット（誤差逆伝播法）の基本的な部分
• 特徴量（ニューラルネットへの入力）

パーティクルフィルタによるグリッドマップマッチング

今回，著作権的にグレーな感じがするので問題があればご一報ください．

専門ではないですが，ロボットの自己位置推定に興味を持ったのでMATLAB/simulinkでシミュレーションしてみました．
ところどころ用語が間違っていると思いますがご容赦ください．

手を付けたきっかけは，“グリッドマップのマッチングに基づく未知障害物にロバストな自己位置推定”という論文の発表を聞いたことです．

面白そうだったので今回は，これを実装したいと思います．

また，今回の記事の大部分はMyEnigmaさんの記事を参考にしています．

パーティクルフィルタやグリッドマップそれぞれについては，こちらがとてもに詳しいのですが，これらを組み合わせた実装について記載がなかったため，本記事で述べていきます．

自己位置推定は，ロボットが現在，どこに居るのかを知るための手法です．
見晴らしの良い屋外であればGPSで事足りることが多いかと思いますが，そうでない場合にはロボット自身が現在の位置を知ることは難しくなります．

自己位置推定はよく知るものだとカーナビがありますが，これは普段はGPSで位置を測定していて，トンネル内などのGPSの電波が入らないところでは，タイヤの回転数（回転角）や加速度センサの2回積分（$\simeq$位置）を使って，自己位置を推定しています（多分）．

カーナビでは，トンネル内であまり自己位置の精度が要求されない（トンネル内での分岐は少ない）ことと，道路自体がカーナビ内部にあるので推定された位置が多少ズレていても道路に引き寄せられること，最終的な運転は人間が行うことから，あまり自己位置推定は難しくない/重要でないのでは，と以前は思っていました．

ロボットの場合だと，精度の高い自己位置の推定が要求されるため，いろいろな手法が開発されてきたとのことです．あまり詳しくないのでこの辺りは省略します．

今回はその中でもパーティクルフィルタによるグリッドマップマッチングについてシミュレーションをしてみました．

多入力多出力系とは

自分も初めて聞いた時に，

多入力多出力系ってなんぞや？入出力ごとに制御系組めばいいんじゃないの？

とおもったので，覚え書き程度にまとめておきます．

多入力多出力系（システム）とは，制御工学において，入力および出力が1つだけではないことを指します．

一体何言ってるかわかんないよね．

ということで制御工学でときおり出てくる，多入力多出系について．

多入力多出力系（Multi Input Multi Output; MIMO; マイモ）系とは，上で述べているように入出力が1つではない制御対象の事を言います．

多入力多出力系は，よく制御が困難であると言われます．

この記事では，多入力多出力系の制御のどこがどう困難であるのかについて，簡単に述べたいと思います．

倒立振子シミュレーション(MATLAB/Simulink)

せっかく個人用（学生用）のMATLABを買ったので，それを使ってみましょう的な記事．
制御系の研究ではお馴染みの倒立振子のシミュレーションをやってみます．

動画はシミュレーションの様子．
実は，ToolBoxにSimscapeというものがありまして，その中のSimMechanicsではMATLAB上で3Dソリッドモデルを作るだけで物理シミュレーションができたりします．

ただし，今回は普通のSimulinkのみで行います（笑）

やっぱり解析的に色々するには，どちらにしろ運動方程式も立てないといけないので．
大体の場合において，制御をしようと思うと，運動方程式のパラメータから
制御系の設計パラメータを求める必要が出てきます．
そんな理由から，運動方程式を立てて，制御パラメータを決めていく方法で今回は制御系の設計を行います．

STM32F4 DISCOVERYことはじめ

なんとも良い時代になりました．
32bit 高速CPUの評価ボードが2,000円の時代です．

その名もSTM32F4 DISCOVERY．

STMicroelectronics社の出している組み込みデバイス用の32bit CPU，STM32F407VGT6の評価用ボードです．
主な特徴は以下．

CPUの特徴

• 168MHz動作の高性能CPU，STM32F407VGT6（ARM Cortex-M4）を搭載．
• flashが1MB，RAMが192kBのデータ領域．
• 3つの12bit ADC，2つの12bit DACを搭載．
• 浮動小数点演算装置（FPU）搭載．

ボードの特徴

• 3軸MEMSモーションセンサー（LIS302DL）搭載．
• 無指向性MEMSマイクロホン搭載．
• DAコンバータ内蔵D級アンプ搭載．
• その他，LED，プッシュスイッチ，ST-LINK用インターフェースなどなど．

ST-LINK経由で書き込みができ，統合開発環境もメーカーのページから無料でDLできるので，導入が簡単です．
特に実数の処理に非常に有利なFPUが搭載されているため，現代制御などでのフィジカルコンピューティングを目的としている時に気兼ねなく実数演算を使えます．
値域を気にして固定小数点を作る時代は終わったのだ…

タイマも簡単で，数も豊富なため，PWMが気軽に使えるところもなかなかのメリットだと感じています．

他には，ひとつのピンに対して複数のAlternate Finction（AF）が割り当てられているので，
ピン配置が思いのほか自由なところが気に入っています．

例えば，PA1とPA2にはタイマ5とADC123が割り当てられているので，
PA1をADCに，PA2をPWM（タイマ）に使用することが簡単にできるようになっています．

micro USB経由での書き込みができるようになっている（ST-Link）ので，
ボードさえ買ってしまえば即座に使用可能です．
（ピンヘッダもついてるし）

また，SysTick_Handlerという割り込みハンドラが標準で定義されていて，
周期割り込みの設定にほとんど手間取らず，わかりやすいところもよいです．

秋月で2,000円ちょうどで販売しているので，入手性も申し分ないため一度買ってみるのもアリかとおもいます．
1ランク下のSTM32F3 DISCOVERY（1,250円）もあります．

主な違いはコアの性能の違いと，周辺機器の違いです．
個人的にはF3 DISCOVERYのほうが好きなんですが（個人で購入したのはこっちだけ．安いから），
研究室の諸事情によりF4のほうを使っています．

残念ながら，F3とF4，完全互換というわけにはいかず，多少コードを書き換える必要があったように思います．
ピン割り当てはモチロンのこと，Peripherala の設定が多少違ったような…

正直FPU搭載に目が眩んで買ってしまいましたが，用途はかなり幅広いと思われます．
（雑なまとめ）

STM32F4 DISCOVERYでのPID制御

STM32F4 DISCOVERYで簡易なPID制御をやってみました．

PID制御は，制御入力が以下の式で表される制御のことです．

$$u = K_{p}(e(t) + \frac{1}{T_{I}}\int e(t) \ \mathrm{dt} + T_{D}\mathrm{\frac{d}{dt}}e(t)) \\ \because K_{p}: \mathrm{比例ゲイン, const.}, \\ T_{I}: \mathrm{積分時間, const.},\\ T_{D}: \mathrm{微分時間, const.},\\ e(t): \mathrm{誤差}$$
上記はIdeal形式の表記で，Parallelでは以下のようになりますが，本質的には同じです．

$$u = K_{P} \ e(t) +\frac{1}{T_{I}}\int e(t) \ \mathrm{dt}+T_{D}\mathrm{\frac{d}{dt}}e(t)$$

今回はIdeal形式の方で，制御則を組んでいます．
目的としては，ポテンショメータの値をAD変換し，PWMでモータを駆動して，位置サーボをするというものです．

---

まずは，GPIOの設定です．

STM32のシリーズでは，ピンとして用意されているGPIO（PA0, PA1, …）のそれぞれに色々な機能がついています．
このため，何か機能を使いたいときには，使いたい機能を持っているGPIOを設定してから，個別の機能を設定する必要があります．

GPIO（General Purpose IO）には，デジタル（1か0）の入力，出力の他に，それぞれのピンに割り振られた特殊機能を使うことができます．
特殊機能の詳細はベンダーホームページのデザイン・リソースから，UM1472: Discovery kit for STM32F407/417 linesをみるとわかります．

この表を見ると，PA1にはADC_123_IN1が，PA2にはTIM2_CH3があることが分かります．

このため，PA1をADCとして，PA2をPWM（タイマ機能の一部であるため）の出力として利用します．

GPIO No.	Alternate function	details
PA1	ADC1	AD変換器1．ポテンショメータの値を読む用
PA2	TIM2 CH3	タイマ2のチャンネル3．モータへのPWM出力（正回転）

　
このときのGPIOの設定はこちら．

GPIO_InitTypeDef GPIO_InitStructure;

/* GPIO for ADC */
RCC_AHB1PeriphClockCmd(RCC_AHB1Periph_GPIOA, ENABLE);
GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AN;
GPIO_InitStructure.GPIO_Speed = GPIO_Speed_50MHz;
GPIO_InitStructure.GPIO_PuPd = GPIO_PuPd_NOPULL;
GPIO_InitStructure.GPIO_OType = GPIO_OType_PP;
GPIO_InitStructure.GPIO_Pin = GPIO_Pin_1;
GPIO_Init(GPIOA, &GPIO_InitStructure);

/* GPIO for TIM2 */
GPIO_PinAFConfig(GPIOA, GPIO_PinSource2, GPIO_AF_TIM2);
GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AF;
GPIO_InitStructure.GPIO_Speed = GPIO_Speed_100MHz;
GPIO_InitStructure.GPIO_PuPd = GPIO_PuPd_NOPULL;
GPIO_InitStructure.GPIO_OType = GPIO_OType_PP;
GPIO_InitStructure.GPIO_Pin = GPIO_Pin_2;
GPIO_Init(GPIOA, &GPIO_InitStructure);

まずは，GPIOA全体の有効化として，以下を記述．
RCC_AHB1PeriphClockCmd(RCC_AHB1Periph_GPIOA, ENABLE);

次に，GPIOAの1番ピンはADCとして使うので，
GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AN;
としてAN（Analog）モードに．

2番ピンはPWM（タイマ）として使うので，
GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AF;
としてAF（Alternate function）モードに設定します．

---

次に，ADCの設定．

　
個別にADCの設定をしていきます．

ADC_CommonInitTypeDef ADC_CommonInitStructure;
ADC_InitTypeDef ADC_InitStructure;

RCC_APB2PeriphClockCmd(RCC_APB2Periph_ADC1, ENABLE);

ADC_StructInit(&ADC_InitStructure);
ADC_DeInit();

ADC_CommonInitStructure.ADC_DMAAccessMode = ADC_DMAAccessMode_Disabled;
ADC_CommonInitStructure.ADC_Mode = ADC_Mode_Independent;
ADC_CommonInitStructure.ADC_Prescaler = ADC_Prescaler_Div4;
ADC_CommonInitStructure.ADC_TwoSamplingDelay = ADC_TwoSamplingDelay_5Cycles;
ADC_CommonInit(&ADC_CommonInitStructure);

ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;
ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right;
ADC_InitStructure.ADC_ExternalTrigConvEdge = ADC_ExternalTrigConvEdge_None;
ADC_InitStructure.ADC_NbrOfConversion = 1;
ADC_InitStructure.ADC_Resolution = ADC_Resolution_12b;
ADC_InitStructure.ADC_ScanConvMode = DISABLE;
ADC_Init(ADC1, &ADC_InitStructure);

ADC_RegularChannelConfig(ADC1, ADC_Channel_1, 1, ADC_SampleTime_15Cycles);
ADC_Cmd(ADC1, ENABLE);

ADC_SoftwareStartConv(ADC1);
while(!ADC_GetFlagStatus(ADC1, ADC_FLAG_EOC));

　
まずは，今回使うADC1の有効化．
RCC_APB2PeriphClockCmd(RCC_APB2Periph_ADC1, ENABLE);

　
次にADC全般の設定．

ADC_CommonInitStructure.ADC_DMAAccessMode = ADC_DMAAccessMode_Disabled;
ADC_CommonInitStructure.ADC_Mode = ADC_Mode_Independent;
ADC_CommonInitStructure.ADC_Prescaler = ADC_Prescaler_Div4;
ADC_CommonInitStructure.ADC_TwoSamplingDelay = ADC_TwoSamplingDelay_5Cycles;
ADC_CommonInit(&ADC_CommonInitStructure);

今回は，DMAを使わず，個別にADCするので，DMAAccessModeはDisableにしています．

次に，もう少し詳細なADCの設定．

ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;
ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right;
ADC_InitStructure.ADC_ExternalTrigConvEdge = ADC_ExternalTrigConvEdge_None;
ADC_InitStructure.ADC_NbrOfConversion = 1;
ADC_InitStructure.ADC_Resolution = ADC_Resolution_12b;
ADC_InitStructure.ADC_ScanConvMode = DISABLE;
ADC_Init(ADC1, &ADC_InitStructure);

前述のように個別でADCするので，ContinuousConvModeはDisableに．
また，今回はADCはこの1chしか使わないので，NbrOfConversionは1に設定しています．
DataAlignは右揃えに，Resolutionは12bit精度にしています．

　
設定が終わったら，有効化して，

ADC_RegularChannelConfig(ADC1, ADC_Channel_1, 1, ADC_SampleTime_15Cycles);
ADC_Cmd(ADC1, ENABLE);

　
一度だけADC1をテスト動作しています．テスト動作が終わる（ADC_FLAG_EOCが立つ）まで，待っています．

ADC_SoftwareStartConv(ADC1);
while(!ADC_GetFlagStatus(ADC1, ADC_FLAG_EOC));

これでADCの設定は終わり．

---

次はタイマの設定．

　
タイマ個別の設定は以下．

#define PWM_FREQ 20e3
#define PWM_DIV 200

TIM_TimeBaseInitTypeDef TIM_TimeBaseInitStructure;

RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE);

/* timer 2 for PWM */
TIM_TimeBaseInitStructure.TIM_ClockDivision = 0;
TIM_TimeBaseInitStructure.TIM_CounterMode = TIM_CounterMode_Up;
TIM_TimeBaseInitStructure.TIM_Period = (uint32_t)(PWM_DIV - 1); //20kHz, 200div, 0.5percent per div, period = 84MHz / (prescaler + 1) / 20kHz - 1
TIM_TimeBaseInitStructure.TIM_Prescaler = (uint16_t)((float)84e6 / (PWM_FREQ * PWM_DIV) - 1); //desired : 200steps in 20kHz, prescaler = 84 MHz (Core clock) / (20kHz*200step) - 1
TIM_TimeBaseInitStructure.TIM_RepetitionCounter = 0;
TIM_TimeBaseInit(TIM2, &TIM_TimeBaseInitStructure);

/* configuration for PWM */
TIM_OCInitStructure.TIM_OCMode = TIM_OCMode_PWM1;
TIM_OCInitStructure.TIM_OutputState = TIM_OutputState_Enable;
TIM_OCInitStructure.TIM_OCPolarity = TIM_OCPolarity_High;
TIM_OCInitStructure.TIM_Pulse = PWM_DIV * 25 / 100; //25% duty
TIM_OC3PreloadConfig(TIM2, TIM_OCPreload_Enable);
TIM_Cmd(TIM2, ENABLE);
TIM_CtrlPWMOutputs(TIM2, ENABLE);
TIM_OC3Init(TIM2, &TIM_OCInitStructure);

上部defineのPWM_FREQはPWMの周波数，PWM_DIVはPWMの分解能を表しています．
PWM_DIVを分解能と呼ぶのは少し正規の表現から外れていますが，便宜上こう呼びます．
PWM_DIVは，PWMのデューティ比を何分割までできるかを表しており，
今回は200分割，つまり0.5%単位でのデューティ比の設定が可能となっています．

PWMの場合，モータドライバの電圧によって同じデューティ比でも実際にモータにかかる電圧が変わってくるので，
こういった設定のほうがいいんじゃないかということで．

たとえば，モータドライバの最大出力電圧が10Vで，
TIM_OCInitStructure.TIM_Pulse = PWM_DIV * 25 / 100;
とすると，出力電圧は10Vの25%である，2.5V相当となります．

次にタイマ2の有効化．
RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE);

PWMを200分割で使用するための設定．
TIM_TimeBaseInitStructure.TIM_Period = (uint32_t)(PWM_DIV - 1);

プリスケーラの設定．
TIM_TimeBaseInitStructure.TIM_Prescaler = (uint16_t)((float)84e6 / (PWM_FREQ * PWM_DIV) - 1);

タイマ2の設定の反映．
TIM_TimeBaseInit(TIM2, &TIM_TimeBaseInitStructure);

次に，タイマをPWMで使用することや，タイマのピン出力をONにすることなどを設定し，反映しています．

TIM_OCInitStructure.TIM_OCMode = TIM_OCMode_PWM1;
TIM_OCInitStructure.TIM_OutputState = TIM_OutputState_Enable;
TIM_OCInitStructure.TIM_OCPolarity = TIM_OCPolarity_High;
TIM_OCInitStructure.TIM_Pulse = PWM_DIV * 25 / 100; //25% duty
TIM_OC3PreloadConfig(TIM2, TIM_OCPreload_Enable);
TIM_Cmd(TIM2, ENABLE);
TIM_CtrlPWMOutputs(TIM2, ENABLE);
TIM_OC3Init(TIM2, &TIM_OCInitStructure);

途中，TIM_Pulseに値を代入しているのは，ただのテストです．

これで設定は終わり．

---

次は時間割り込みの処理について記述します．

PIDでは積分器および微分器を使うので，PIDの計算が周期的に実行されてほしいです．
このため，PIDが20ms周期で実行されるように，時間割り込みの処理について記述しています．

STM32F407では，周期割り込みが関数として記述されているみたいなので，それを使います．
SysTick_Handlerという関数です．

別スクリプトの設定で，SysTick_Handlerには1ms周期で入るようになっています．

#define STROKE 25.0
#define PID_SMPL_PERIOD 20

void SysTick_Handler(void){
    static float c = 0, u = 0, pos = 0;
    static int i = 0;

    i++;

    /* ADC */
    if(!(i % 20)){
        ADC_SoftwareStartConv(ADC1);
        while(!ADC_GetFlagStatus(ADC1, ADC_FLAG_EOC));
        pos = STROKE * (float)ADC_GetConversionValue(ADC1) / (float)((1<<12) - 1);
    }

    /* PID */
    if(!(i % PID_SMPL_PERIOD)){
        u = PID(c, pos);

        /* saturation */
        if(u > PWM_MAX_VOLTAGE){
            u = PWM_MAX_VOLTAGE;
        }else{
            u = 0;
        }


        TIM_OCInitStructure.TIM_Pulse = (uint32_t)(PWM_DIV * (u / PWM_MAX_VOLTAGE));
        TIM_OC3Init(TIM2, &TIM_OCInitStructure);
    }

    if(i == 10000) i = 0;
}

この関数内では，AD変換は20msごとに実行され，
ポテンショメータの電圧から，実際のストロークに変換されています．

AD変換をスタートし，
ADC_SoftwareStartConv(ADC1);

終了するまで待っています．
while(!ADC_GetFlagStatus(ADC1, ADC_FLAG_EOC));

そのAD変換値ADC_GetConversionValue(ADC1)から，実際の位置posに変換しています．
pos = STROKE * (float)ADC_GetConversionValue(ADC1) / (float)((1<<12) - 1);
ここで，AD変換が12bitなため，その最大値((1<<12) -1)で割っています．

---

PIDの計算は，PID_SMPL_PERIODごとに行われます．
今回はADCの周期と同じ20msごとです．
u = PID(c, pos)

PID()は，目標値cと現在の位置posを引数として与えて，返り値として計算された操作量uを得る関数です．

ここでは，計算された操作量がドライバで出力できるかどうかの判断を行っています．

/* saturation */
if(u > PWM_MAX_VOLTAGE){
    u = PWM_MAX_VOLTAGE;
}else{
    u = 0;
}

あとは，操作量（電圧）からPWMのデューティ比を計算し，
TIM_OCInitStructure.TIM_Pulse = (uint32_t)(PWM_DIV * (u / PWM_MAX_VOLTAGE));

実際の出力に反映させています．
TIM_OC3Init(TIM2, &TIM_OCInitStructure);

---

最後に，PID()関数の中身について．

PID()関数は非常にシンプルで，以下のようになっています．

float PID(float c, float pos){
    static float e = 0, ei = 0, ed = 0, e_pre = 0;

    e = c - pos;
    ei += e * (PID_SMPL_PERIOD * 1e-3);
    ed = (e - e_pre) / (PID_SMPL_PERIOD * 1e-3);
    e_pre = e;

    return Kp * (e + 1 / Ti * ei + Td * ed);
}

ここでcは目標の位置，posは現在の位置です．

まず，現在の偏差$e(t)$を計算し，
e = c - pos;

偏差の積分値$\int{e(t)} \ \mathrm{dt}$，
ei += e * (PID_SMPL_PERIOD * 1e-3);

偏差の微分値$\mathrm{\frac{d}{dt}}e(t)$，を計算します．
ed = (e - e_pre) / (PID_SMPL_PERIOD * 1e-3);

微分計算では，前回の誤差と現在の誤差を使って計算するので，前回の誤差e_preを更新します．
e_pre = e

これらの値を元に，$K_{p}(e(t) + \frac{1}{T_{I}}\int e(t) \ \mathrm{dt} + T_{D}\mathrm{\frac{d}{dt}}e(t))$を計算し，returnします．
return Kp * (e + 1 / Ti * ei + Td * ed);

---

これらの記述で，STM32F4DISCOVERYでPID制御ができます．
以下はこれらを基にしたプログラムリスト(main.c, timer.c, timer.h)です．

このままビルドするとSysTick_Handler()が重複しているとかでエラーが出るので，
stm32f4xx_it.c内のSysTick_Handler()を消してあげてください．

下の方のCall back関数は消すとリンクエラーが出るので，
これを消したい場合にはリンク先の箇所を消せば動くと思います．

今回は正方向の出力しかできないので，両方向にモータを回転させたい場合は，
PWMのチャンネルを1つ追加する必要があると思います．

C言語でニューラルネット（2次元平面での分離）

前回に引き続き，ニューラルネットワークです．

前回は，ニューラルネットによる関数近似を行いましたが，今回はちょっとパターン認識寄り．
というのも，2次元平面上のある点が，AとBどちらに分類されるのかという問題なためです．

たとえば，以下の様な状況の場合．

これを2種類に分類したいとします．
理想的な場合，このように分類できると想像できると思います．

先ほどの4点を与えたとき，大体の人はこのように線を引いて分けると思います．
このような1本の線では分離ができない問題のことを，線形分離不可能な問題と言ったりします．

線形分離可能な例はこんな感じ．

1本の直線で，分類できています．

線形分離可能な問題は，単純パーセプトロンなどでも解けますが，線形分離不可能の場合には，隠れ層を持つようなニューラルネットワークの構成である必要があります．

このような隠れ層を持つニューラルネットワークの学習方法のひとつが，誤差逆伝播法（Back Propagation）です．

というわけで，2次元での分類を行っていきます．

まず，2次元の入力（x軸，y軸）に対して出力（●，×）があるので，入力次元は2，出力次元は1です．

data.dat		rho.dat
0	0	1
0	1	0
1	0	0
1	1	0

まずは，線形分離可能なパターンから試してみます．
ここで，data.datの左側がx軸，右側がy軸で，rho.datの0が●，1が×を表しています．
よって，このデータは先ほどの直線で分離できるデータを数字で表したものとなっています．

これをニューラルネットワークに入力として与えると，以下のようになります．

このように，確りと分離できていることがわかります．
プログラム中では，positive.datに0.5以上（厳密には超過）の出力群を，negative.datに0.5以下の出力群をファイル出力しています．

次に，線形分離不可能の場合．

data.dat		rho.dat
0	0	1
0	1	0
1	0	0
1	1	1

これは最初の図に対応します．
このときの出力は，以下のようになります．

左下と右上で，×になっています．
ここからわかるように，非線形の分離問題（EX-OR問題）に対しても分離を行えていることがわかります．

プログラムリストは，前回の物から少し変えてあります．（主にTest関数内）
以下リスト．

C言語でニューラルネット（関数近似）

C言語で関数近似を行うためのニューラルネットワーク（誤差逆伝播法，B ack P ropagation N eural N etwork）のプログラムです．

プログラム中で，Input_sizeは学習用データの入力の長さ，Input_dimは学習用データの次元，test_sizeは検証用データの入力の長さです．

たとえば，以下のような点の集合を関数で近似したい場合を考えます．

このとき，学習用データは次のように与えられます．

data.dat	rho.dat
0	1
1	3
2	2
3	4

というわけで，上のように実行ファイルの場所に，data.datとrho.datを作成します．
見て分かる通り，data.datがx軸，rho.datがy軸を表しています（一次元の場合）．

この2種類のファイルから，Input_sizeとInput_dimを設定します．
この場合，Input_sizeは，データが横軸に4つなので，Input_size=4とします．
また，入力は1次元なのでInput_dim=1です．

また，test_sizeは，学習が終わったあとのデータ点の数で，x軸の0~3を何分割するかを決める滑らかさの値です．
理想的な場合学習が起こった場合では，test_size=4とすると，上の図の点と同じ位置に点が打たれるはずです．
ここでは，どのように近似されたかを詳細に見るために，test_size=100とします．

これらのデータから実行した結果が以下の図．（学習回数1,000,000回，test_size=100）

このように，学習データとして与えた点を全て通るような曲線が描かれています．
これが，ニューラルネットワークによる関数近似です．

先ほど与えた点の値と，現在のニューラルネットワークの出力の誤差がどんどん小さくなっていくことで，与えた点を通る曲線が描かれています．
この曲線を拡大してみると，100個の点からなる直線の集合です（=test_sizeの値）．

test_sizeを小さくするとこんな感じに．（学習回数1,000,000回，test_size=10）

また，学習回数が少ない場合など．（学習回数10,000回，test_size=100）

学習回数が少ないと，与えた点との近似率が落ちています．（点を通っていない）

学習回数が更に少ない場合．（学習回数1,000回，test_size=100）

ここまでくると，もはや点を通ってすらいないです．

このように，ニューラルネットワークでの関数近似は，与えられた点と現在のニューラルネットワークによる出力の誤差が小さくなるように繰り返し学習することで，与えられた点を通る出力を得られるというものです．

というわけで以下リスト．

プログラム内の変数の詳細についてまとめておくと以下のようになっています．

label	details
Input_size	入力の長さ（大きさ）
Input_dim	入力の次元
Hidden_dim	中間層の数
alpha	シグモイド関数の傾き（後述）
eta	学習係数（他サイト参照）
gain	ゲイン（後述）
test_size	学習後の検証用データの量（大きさ）
th	未使用

ここで，ニューラルネットワークに使用しているシグモイド関数は，0~1の範囲の実数しか取ることができないので，このままでは0~1の出力の範囲の点の近似しかできません．

このため，シグモイド関数にgainを掛けることで，擬似的にシグモイド関数の出力範囲を広げています．
たとえば，gain=100の場合では，シグモイド関数の出力範囲は-50~50の実数となります．
一応，負の値も取れるように，負方向にgainの半分だけずらしています．

これらの処理により，ニューラルネットワークで-50~50までの範囲の点を近似することができるようになっています．

MATLAB学生版の値下げ，個人ライセンスの開始

（朗報）なんとMATLABの個人向けライセンスが2014年4月に始まっていたそうです．

全然知らなかった…

これで趣味でMATLAB使いたい人も簡単に導入できるようになりそうです．

主な料金は以下．（MATLAB Student R2015a，2015年4月11日現在）

	MATLAB 本体	MATLAB本体 +Simulink	Toolbox （MATLAB本体と同時購入の場合）
同梱 ソフト	Symbolic Math Toolbox
	Statistics and Machine Learning Toolbox
	Optimization Toolbox
	Control System Toolbox
	Signal Processing Toolbox
	DSP System Toolbox
	Image Processing Toolbox
	Data Acquisition Toolbox
	Instrument Control Toolbox
	Simulink
	Simulink Control Design
価格	8,950円		各895円

なぜか，MATLAB本体のみ（Simulink無し）を選んでもSimulinkが付いてきてます．
どうせ使うだろうしいいんだろうけど．

また，個人向けライセンス料金は以下．

	MATLAB本体	Toolbox（MATLAB本体と同時購入の場合）
同梱ソフト	無し
価格	14,950円	各4,500円

個人向けライセンスでは，一切Toolboxは付いてきません．
こう見ると学生版ライセンスはやっぱり安いですね．

それでも，これまで個人は基本料金（約20万円！）を支払うしかなかったことを考えると大分安くなりました．
実は学生版ライセンスも，以前と比べて半額くらいになっているそうです．

ただし，学生ライセンスで注意することとして，

1. 有効期限内の学生証のコピーが必要（大学の生協経由でも買えるかも？）
2. 期限は無期限．ただし，購入から4年以上後に再インストールする場合は再度学生証の提出が必要．

これくらいでしょうか．まだ購入手続きが完了していないので買ってからいろいろ書き足そうと思います．

自分の場合は，Control System ToolboxとかCurve Fitting Toolboxとか入れても13,000円程度だったのですぐにでも買いたいです．

残念ながらクレジットの上限にかかっているので来週まで待たなければ…

StackEditとBloggerの連携

StackEditとBloggerの連携について．

何度かStackEditでMarkdown記述→Bloggerに投稿の流れを試しましたが，いろいろ問題が出たのでそのまとめでも．
内容については以下．（この記事はStackEdit→Bloggerで書いてます．）

1. StackEditの始め方
2. StackEditの記述をBloggerに投稿
3. Bloggerとのテンプレートのすり合わせ

まずは1. StackEditの始め方から．

StackEditは基本的にブラウザで使用します．（というかこれ以外知らない）
まずはStackEditのサイトにアクセス．

こんな感じの画面になるはず．
最初はHello!か何かの記事がデフォルトで表示されていると思うので，右上のフォルダ（ディレクトリ）アイコンからNew Documentを選択．
このとき，Please consider sponsoring StackEdit for $5/year (or sign in if you're already a sponsor).とか出ますが，無料で使えるので気にせずOKを．

すると，さっきのようにごちゃごちゃしていない，きれいな状態のページが出るはず．
ここにMarkdown形式で記述していきます．
記述方法に関しては他のサイトを参照してください．

右上の部分に記事のタイトルが，画面の左にMarkdownで記述するエディタが，画面の右にMarkdownのHTML上での実際の記事が映し出されます．

次に，2.StackEditの記述をBloggerに投稿します．

投稿は簡単で，画面左上のシャープマーク#を押し，出てくるリストの中から，Publish (Export to the web)を選択します．

こんな画面が出てくるので，Blog URLに自分のBlogger URLを記入．
Update existing post IDは空でいいです．
FormatはHTMLを選びましょう．

さっきと同じ，スポンサーになってください的なメッセージが出るけど，そのままOKを押せば大丈夫です．

すると，Bloggerの認証を求められるので，流れのまま認証してください．

これで簡単にBloggerでMarkdownで書いた記事が更新……できない！

厳密には更新できているんですが，数式やらが表示できません．
試しに，
$$\frac{a}{b}$$$
と書いて更新してみましょう．
StackEdit上では数式 $$\frac{a}{b}$$ が表示されているんですが，Bloggerではまったく表示されていません．
また，見出しのオサレな背景↓

見出し

も表示できていません．

これは，記事のスタイルがBloggerのCSSスタイルになってしまっているから．
なのでStackEditのMarkdown用CSSをlinkしましょう．

Blogger編集画面のテンプレートから，HTMLの編集を選択．
上のほうに<head>~</head>で囲まれた部分があるので，その中に，

<link rel="stylesheet" href="https://stackedit.io/res-min/themes/base.css" />
<script type="text/javascript" src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML"></script>

を追加しましょう．
これで変更を保存すると，BloggerでもMarkdownのようなオサレな表示ができるようになります．

ただし，条件付きで．

このテンプレートをそのまま使うと，Bloggerの方で設定した背景色やら見出し文字の大きさやらが崩れてしまいます．

これを防ぐには，3.Bloggerとのテンプレートのすり合わせをする必要が．

実はBloggerで設定したCSSと，linkしたCSSが競合してしまっています．
なので，Markdownのスタイルシートであるhttps://stackedit.io/res-min/themes/base.cssを拾ってきて，自前で書き換える必要があります．

これが一番面倒．

とりあえず今回は，linkを外して，拾ってきたCSSファイルの見出しに関する部分と背景色に関する部分だけ削除することにしました．

削除したものが以下．いくらなんでも長いので開けるときは注意してください．
<style>~</style>を<link>の行の代わりに置いてください．
これでMarkdownのオサレな記述がBloggerでも問題なく表示できるはず．

以上，本文でした．

以下，スタイル部分の中身です．

Stack EditでのMarkdown記法テスト

Stack EditorでのMarkdown記法テストの続き．

Template形式でのアップロードでは，Bloggerのテンプレート自体に影響が出てしまうようで，使わないほうが無難かも．

Stack EditでのMarkdown記法テスト

--- tags: ["other"] --- Markdown ========= `Stack Editor`というMarkdown形式で書いた文章をHTMLに自動で変換して，Bloggerに投稿してくれるサービスのテスト． リストも簡単だったり，Syntax Highlightもあるとかなんとか． * リストのテスト * bulletで表示されるのかな？ > 小見出しも簡単に作れるのがよさそう． ちょっと前に，`Markdown→LaTeX`をやってみたりしたけど，TeXの場合は直接書くほうが楽だったんで，すぐにお蔵入りに． >数式も簡単に書けるみたい． $$ e^{i \theta}=\cos{\theta}+i\sin{\theta}$$ 便利な世の中になりました． Stack Editorでエクスポートする際の選択肢に， * Markdown形式 * HTML形式 * Template形式 の三つがあるので，それぞれどうなるか試してみる． > これはTemplate形式．

Stack EditでのMarkdown記法テスト

Markdown

Stack EditorというMarkdown形式で書いた文章をHTMLに自動で変換して，Bloggerに投稿してくれるサービスのテスト．

リストも簡単だったり，Syntax Highlightもあるとかなんとか．

• リストのテスト
• bulletで表示されるのかな？

小見出しも簡単に作れるのがよさそう．

ちょっと前に，Markdown→LaTeXをやってみたりしたけど，TeXの場合は直接書くほうが楽だったんで，すぐにお蔵入りに．

数式も簡単に書けるみたい．

$$e^{i \theta}=\cos{\theta}+i\sin{\theta}$$

便利な世の中になりました．

Stack Editorでエクスポートする際の選択肢に，

• Markdown形式
• HTML形式
• Template形式

の三つがあるので，それぞれどうなるか試してみる．

これはTemplate形式．

近況

久しぶりの更新です．

久しぶりという次元ではないかもしれないです．（前回の更新が1年以上前の2013年11月24日）

怒涛の2015年3月が終わったので，

一年間の近況（もはや近況ではないのかもしれない）について少し．

2013年

12月

当時B4（学部4年）なので，卒論の執筆．

1月の中旬提出なのでこの時期は大変でした．

2014年

1月

卒論提出と卒論の発表会．

実はそこまで追い込まれてはなかったですけど，意外とブログとか忘れていたり．

2-3月

国際会議の論文執筆．津

前回のオーストリアの国際会議は，著者は先生だったので今回は自分で書いてみようということに．

これまで英語で文章を書いたことない学生が必死で書きました．

4-5月

授業とか実験とかいろいろ．

M1（修士1年）になり，授業も始まって新年度．

6月

前に出した国際会議に参加．

プレゼンの準備とかいろいろ．

場所はスウェーデンのカールスタードでした．

7-8月

8月の中ごろにはちょっとした短期留学でアメリカのカリフォルニアに．

夏場でありながらなかなか涼しかった思い出が．

このあたりのことはまた記事にするかも．

9-10月

大学院生らしく実験とか．

11月

次の国際会議の論文執筆．

実は出すたびにテーマが変わっていたりします．

12月

後輩の卒業論文手伝い．

1月が提出&プレゼンなのでなかなか忙しかったり．

2015年

1月

後輩の卒論提出&発表．

自分もショートペーパーと論文抱えていたのでずっと書きっぱなし．

2月

研究やら実験やら

3月

1週目の土日：妹の引越し先探しに東京まで．

2週目の土日：機会学会の関西の学生委員なので，卒業研究発表会というイベントを開催．

3週目の土日：11月に書いた国際会議論文の発表でシンガポールへ．

4週目：修士の中間報告会＋就職活動の説明会＋エントリーシート提出．

とまあ，1年を通してこんな感じでした．

4月になって，やっと休みができたのでちょっと更新．

自分の備忘録にでもなれば．