EEIC後期実験演習例 Python Control版
前にjupyterで書いてみたは良いが重くてGithubで公開をためらった代物です。 Markdown変換してみました。
はじめに：Environment settings
Using anaconda you just need to activate following command.
pip install python-control
For further information see here. folllowing command is additional.
- Textbook for this experiment is here.
# !python -m pip install --upgrade pip # !pip install slycot # I failed to install these command
Initialization
Import matlab like package from control.matlab
from control.matlab import * # MATLAB-like functions
basic feature
At first you can use tf function like matlab.
- 下記は演習2の(1)の問題を例に取り上げます。
# Q2-1 num1 = [1,2] den1 = [3,1] sys = tf(num1,den1) print(sys) # this is the system output
s + 2
-------
3 s + 1
Load matplotlib for drawing figure
import matplotlib.pyplot as plt %matplotlib nbagg # import scipy from scipy import arange
# step responce # Generate step responses plt.figure(1) (y1a, T1a)=step(sys,T = arange(0, 10, 0.01)) plt.plot(T1a,y1a)
<IPython.core.display.Javascript object>
[<matplotlib.lines.Line2D at 0x15ae2ea4438>]
Error メモ
sysの関数を変更すると？ときどき'tuple' object is not callableというエラーが起きます。
イマイチ自信がないですがどこか変更してはいけない部分を変更しようとしているようです。
一度カーネルをRestartすると治るのでとりあえずはこれで勘弁してください。
# bode plot
bode1 = bode(sys)
<IPython.core.display.Javascript object>
演習4：2次系の応答
数式も打ち込みたいけど時間がないので割愛。 減衰比（ダンピング）を変えたらどうなる？という話です。
共振周波数を変えるVerもありましたが，ダンピングの例だけ示します。
Gs = [] K = 1 w_n = 10 Ys = [] Ts =[] legend = [] # changing zeta plt.figure(3) for z in range(0,8): zeta = z * 0.2 Gs.append(tf([K*w_n*w_n],[1,2*zeta*w_n,w_n*w_n])) (y1, T1)=step(Gs[-1],T = arange(0, 1, 0.001)) Ys.append(y1) Ts.append(T1) # plot plt.plot(Ts[-1],Ys[-1]) legend.append('\zeta is '+str(zeta)) plt.legend(legend) plt.show()
<IPython.core.display.Javascript object>
- どうせならbode線図も見てみましょう
plt.figure(4) for G in Gs: bode(G)
<IPython.core.display.Javascript object>
演習5：ゼロ点の性質
最終値の定理から行き着く値は1/15で0.06くらいになることはわかりきっている。 ではゼロ点はどこに寄与するのか？
よしんば過渡応答であったとしてどのように？
という設問であるが，私が思うに分子がs+1になるような極端な例がないとこの設問の意図は組みにくいであろう。
早急に実験レポートに反映されるべきである。
# sys5 = [tf([1],[1,8,15]), tf([1,4],[4,32,60]) , tf([1,50],[50,400,750]), tf([1,1],[1,8,15])] plt.figure(5) for sys in sys5: (y1, T1)=step(sys,T = arange(0, 1, 0.001)) plt.plot(T1,y1) plt.legend(['G1','G2','G3','Additional']) plt.show()
<IPython.core.display.Javascript object>
レポートの解答で無難に言うならば，「遅いゼロ点はオーバーシュートを引き起こす。早い零点はまぁ無視しても良い。」ということになるが，これはどうしてか。 せっかくなのでBode線図も見てみよう
plt.figure(6) for sys in sys5: print(sys) bode(sys) plt.legend(['G1','G2','G3','Additional'])
<IPython.core.display.Javascript object>
1
--------------
s^2 + 8 s + 15
s + 4
-----------------
4 s^2 + 32 s + 60
s + 50
--------------------
50 s^2 + 400 s + 750
s + 1
--------------
s^2 + 8 s + 15
<matplotlib.legend.Legend at 0x15ae3a986a0>
遅いゼロ点ではゲインが不用意に盛り上がっていることが確認できる。 不要な周波数成分を増幅することにより変な振動が起きてしまう → オーバーシュートが出てくる といった感じである。
直感的に説明すると，「分子が1の場合，安定極であれば分母は単調増加して高周波ほどゲインが下がるのに対し，分子にsの項が加わると分子も増加する要素が加わるためsの寄与が相対的に高い遅い零点ではある周波数周辺でゲインが盛り上がる → オーバーシュート」という感じであろうか。
演習８
結構難しい気がする。問題の意図が汲めない… 一応低域の位相を遅らせることと引き換えに高域の位相を少し早めているって言うことなんだろうが。誰か解説手伝って。
余談：zpk記述
因数分解されたシステムは記述が面倒くさい。 zpk2tfというscipyの関数を用いることでゼロ点，極，DCゲインの３つからtfを生成可能だが，生成結果はtupleであり２変数を受け取るtfに直接食わせられない。
そんなときは*を用いた展開表現を用いることで，tupleを複数argumentとして認識させることが出来る。
Gp =tf(*zpk2tf([],[0,-2,-10],40)) Gc =tf(*zpk2tf([-0.2,-1],[-0.02,-5],5)) plt.figure(7) bode(Gp) bode(Gp*Gc) bode(feedback(Gp*Gc)) plt.legend(['Plant','Plant+C', 'Feedback']) plt.show()
<IPython.core.display.Javascript object>
根軌跡
rlocusやってrlockfindがない！ってOctaveで毎年なるんですがpythonにもなかった…
tau = 0.01 plant = tf([1],[1,3,2,0]) gc = tf([-1.5],[-1.0/tau]) x1,y1 = rlocus(plant*gc)
<IPython.core.display.Javascript object>
Saving
データとか保存用になにかする（予定）
MATLAB＆YALMIPでLMIを解けるようになるまでの手順
Software Installation
ソルバはYALMIP以外にCVXgen等もありますが，どっちがいいんですかね？議論ありますが有識者の意見を聞きたいです。
Yalmipまわりのインストールは以下のスライドが比較的参考になります。2018年9月の時点ですでに情報は古いですがURL等は参考になるでしょう。 http://www.maizuru-ct.ac.jp/control/kawata/iscie/doc/install_yalmip_new.pdf
SeDuMi
単にダウンロード MATLAB直下に解凍。
SDP3
SDPT3 Ver3-4をダウンロード。 MATLAB直下に解凍。
その後，解凍したディレクトリに移動して以下のコマンドを実行します。
Installmex(1)
YALMIP
実はYALMIP単体でもそれなりにソルバはありますが，SDP（半正定値計画問題）を解くには上記のソルバが必要です。
addpath(genpath('C:\Users\<ユーザ名>\Documents\MATLAB\SeDuMi_1_3'))
addpath(genpath('C:\Users\<ユーザ名>\Documents\MATLAB\SDpt3-4.0'))
addpath(genpath('C:\Users\<ユーザ名>\Documents\MATLAB\YALMIP-master'))
yalmiptest
実行結果は以下のとおりです。とりあえずSDPは解けているようなのでOKです。
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | Test| Solution| Solver message| +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | Core functionalities| N/A| Successfully solved (YALMIP)| | LP| Correct| Successfully solved (LINPROG)| | LP| Correct| Successfully solved (LINPROG)| | QP| Correct| Successfully solved (QUADPROG)| | QP| Correct| Successfully solved (QUADPROG)| | SOCP| Correct| Successfully solved (SeDuMi-1.3)| | SOCP| Correct| Successfully solved (SeDuMi-1.3)| | SOCP| Correct| Successfully solved (SeDuMi-1.3)| | SDP| Correct| Successfully solved (SeDuMi-1.3)| | SDP| Correct| Successfully solved (SeDuMi-1.3)| | SDP| Correct| Successfully solved (SeDuMi-1.3)| | SDP| Correct| Successfully solved (SeDuMi-1.3)| | MAXDET| Correct| Successfully solved (SeDuMi-1.3)| | MAXDET| Correct| Successfully solved (SeDuMi-1.3)| | Infeasible LP| N/A| Infeasible problem (LINPROG)| | Infeasible QP| N/A| Infeasible problem (QUADPROG)| | Infeasible SDP| N/A| Infeasible problem (SeDuMi-1.3)| | Moment relaxation| Correct| Successfully solved (SeDuMi-1.3)| | Sum-of-squares| Correct| Successfully solved (SeDuMi-1.3)| | Bilinear SDP| N/A| No suitable solver| +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
LMIの試験プログラム
状態FBでのレギュレーター設計問題を考えます。 LMIでは極配置では出来ない入出力制約等も考えられます。（解があるとは言っていない）
今回はただ単に極が-αより小さいという条件を設定します。
式の解説
ちょっと書くのが面倒なのでしばらくは（省略）。 2次形式のリアプノフ関数の安定化のための条件を解きます。
以下の教科書4章p80~を参照のこと。
LMIによるシステム制御 - ロバスト制御系設計のための体系的アプローチ
- 作者: 蛯原義雄
- 出版社/メーカー: 森北出版
- 発売日: 2012/03/08
- メディア: 単行本（ソフトカバー）
- クリック: 1回
- この商品を含むブログを見る
プログラム
-パラメータ設定（状態方程式）
k = 0.2; f = 0.01; A = [ 0 0 1 0; 0 0 0 1; -k k -f f; k -k f -f]; B = [ 0; 0; 1; 0];
- LMI 立式
αは0.1とします。あんまり大きくすると解けなくなるっぽい。極をある一定以上早くなるよう設定した上でγでフィードバックゲインKをなるべく小さくなるように縛ります。
X = sdpvar(4,4,'symmetric'); Y = sdpvar(1,4,'full'); gamma = sdpvar(1); LMI = [X>eye(4), A*X+B*Y+(A*X+B*Y).'+ 2*alpha*X<0,[gamma,Y;Y.', eye(4)]>0]; sol = solvesdp(LMI,gamma) X = double(X); Y = double(Y); gamma = double(gamma); K = Y/X; % LMI2 result pole = eig(A+B*K); display('Solved LMI pole') display(pole);
配置結果
ゲインを上げずにという制約があるので極の制約のギリギリを攻めるのは予想が付きます。 はいこんな感じ。
pole = -0.1000 + 0.6672i -0.1000 - 0.6672i -0.1000 + 0.0979i -0.1000 - 0.0979i
X>eyeとしているのは安全のためなのでこれを本来のX>0としてしまえばもう少しゆるい条件で解けそうと思っているのですがどうでしょう？。
やってみた所感
とりあえず動作を確認したまでで，LMIの利点を活かした設計などはまだ触れられていません。 ここまでのプログラムで一応ロバスト安定化制御器などは導けますが。
例えばalphaを3とかにするだけでInfeasibleになります。
あと，Successfully Solvedと出ても解けてなかったりして意外といろいろ気を配ります。
そもそもソルバに全部投げているのでCVXgenとかだと解けたりするんでしょうか？ イマイチわかりきっていないところがあります。
Latex 表組みのサイズを微調整
表のサイズの微調整
従来，Latexの表が大きすぎた時 \smallなどのコマンドで表の文字を小さくしてサイズを収めるといったことをしていました。 しかしこの調整は非常に表を醜くします。
正しくは\resizebox{}環境を用いるべきです。
\columnwidthで横幅ぴったりに合わせます。2カラムぶち抜くときは\textwidthを用いると良いです。
% \usepackage{graphics} \resizebox{\columnwidth}{!}{% %table 環境を中に }
GUIに依る表の作成
Latexの表の作成は最もだるい事の一つですが，Tablegeneratorを用いることで簡略化出来ます。 Excelとかで表を作ってしまってGUIで整形，その後にTexに変換して貼り付けるのが最近のやり方です。
お試しあれ。
Google ColaboratryでGPU（TPU）付の機械学習環境が使えるらしい
趣味人としてはそろそろ独自のコードを書いて公開したいところ。（なお本業…） [:contents]
Google Colaboratry事始め
どうも今年の初めにGoogleがCloudで機械学習環境を提供しはじめたようで，それに関する記事が林立しています。 流行りには乗り遅れましたが，気になったのでフォローしてみます。
引用
以下の記事を基にしてテストしました。図なども借りています。
[1] Google Colaboratoryの無料GPU環境を使ってみた – MISO
[2] 【秒速で無料GPUを使う】深層学習実践Tips on Colaboratory
仕様
文献[2]によると重要そうな仕様は以下のとおり。 Ubuntu環境です。
- n1-highmem-2 instance
- Ubuntu 17.10
- 2vCPU @ 2.2GHz
- 13GB RAM
- 40GB Free Space
- GPU NVIDIA Tesla K80
- アイドル状態が90分続くと停止
- 連続使用は最大12時間
- Notebookサイズは最大20MB
Teslaというやつは浮動小数点に強い大規模計算機用のより強いGeforceみたいなやつです。単純にRAMだけでも1080Tiの倍以上の24Gあります。 容量が40Gという点以外は個人が趣味で買えるスペックは大幅に超えていますね。
Jupyterで開発をするという仕様上，コマンドライン操作は!をつけて実行します。
!pip freeze
等でパッケージを確認できますが，TensorflowやKeras，matplotlibなんかは当然のように入っています。環境構築も不要です。
Ubuntuベースにしているのでaptgetも使えるのですが今日は試しませんでした。
準備することとか
以下のURLにGoogleアカウントを持った状態でアクセスすれば環境が手に入ります。 https://colab.research.google.com/notebooks/welcome.ipynb
「ファイル」からPython3のノートを作成して実行していきましょう。 最初に上げたURLの記事からたどればよいでしょう。
CPU → GPU環境への移行
「ランタイム」の欄から「ランタイムのタイプ」を開き，図のように選択します。
今はTPUも選択できるそうな。
作成したファイルの保存
ファイルの保存方法は
- ダウンロード
- Githubへ保存
- GoogleDriveへ保存
がありました。Githubに保存するときはポップアップが出るので事前にポップアップブロックを無効にしておいてください。
KerasでのMNISTのテストと所感
KerasでのMNISTサンプルを試してみました。 ノートブックを上げておきます。 https://colab.research.google.com/gist/YoshiRi/b7a984902222a14d3f4fba31625060ca/kerastest.ipynb
6万Stepを12Epocやって自前のGeforce1050iを載せたPCで計算した時は13s/Epocくらいでしたが，この環境では11s/Epocなのでちょっと早くなっていますね。コードの書き方で変わるんでしょうか。
[http://ossyaritoori.hatenablog.com/entry/2018/03/27/Keras%26_Tensorflow%28GPU%E6%9C%89%29%E3%81%AE%E7%92%B0%E5%A2%83%E6%A7%8B%E7%AF%89_on_Windows_with_Anaconda:embed:cite]
結果
このあたりはサンプルプログラムなのでどうせほぼ同じですが一応。
x_train shape: (60000, 28, 28, 1) 60000 train samples 10000 test samples Test loss: 0.02876817006157935 Test accuracy: 0.9913
所感
GPUマシンを持っていなくても（環境構築しなくても）機械学習の学習が出来るのは初学者にはかなり良いと思います。 特に環境構築や機材は初学者にとって大きな壁なので回線さえあればきちんとした環境で実行できるというのは教育効果が非常に高そうです。
とりあえず，コードを試して満足したい方はここから初めて良いと思いますが，一方で今後ちゃんとやっていくならどこかで環境構築はしておくべきだと思います。 結局自分でデータ集めてモデル組まないと意味ないので僕もこんなサンプルコードばかり実行してないで自分のやりたいプロジェクトを進めていこうとおもいます。
MATLAB背景色変更
Texstudioの背景色変更に加え，matlabも背景色変更します。 ossyaritoori.hatenablog.com
ここで配ってます。 github.com
windows環境では以下のコマンドをmatlabのホームで実行しておしまい。私はdarksteel推しです。
!git clone https://github.com/scottclowe/matlab-schemer.git
cd matlab-schemer\
schemer_import('schemes\darksteel.prf')
「Linear＝線形」 ではない
俺たちは雰囲気で数学をやっている… ”Linear” は 線形ないし一次（Affine）の意味でよく出ますがごっちゃになっちゃうので気をつけましょう。
Linearの訳
Linearの訳ですがシチュエーションによりだいたい2パターンあります。
- 線形：Linear transformationとか
- 一次：Linear equationとか
線形性とは
線形性とは以下の2要素を満たす場合に言います。
- 加法性：f(x+y)=f(x)+f(y)
- 斉一次性 ：f(cx)=cf(x)
スカラー関数では f(x)=ax ということ。
Affineとは
1次方程式という意味でのLinear equationでは f(x)=ax+b のような形を容認します。
これは線形写像に原点の平行移動を含むAffine写像になります。
ただ，Linear solutionとかLinear equationとか言う場合はAffine性よりは「一次方程式」という意味合いが強いです。
「Linear = first-order」 ということです。
具体例
いろんな語句を"線形"と"一次"に分類してみます。なんか思いついたら追記します。
Non linear → 場合による
先程の考察からf(x)=ax+bは線形ではないためnon-linearな写像と言う事ができる一方で1次方程式という点でlinearな方程式ということは出来ます。
ただ，linearな写像とこれを区別するため，affine-linearといった言い回しが用いられる例があるようです。
https://www.researchgate.net/post/y3x_4_is_this_equation_linear_or_nonlinear
Linear inequality → 一次
Linear equationから来てます。線形不等式と訳すなかれこれは一次不等式と訳すべきでしょう。
Linear matrix inequality → 一次？
LMIとして有名な線形行列不等式ですが以下の形をみるに方程式自体はAffineに見えます。
追記：原辰次先生の論文[1]ではキチンと「アファインな関数」と言っていますね。
[1]でもこう言っています。線形とは…？？？
行列一次不等式って言ったほうが良いのでは？ よくわからなくなってきました。
[1]制御におけるRiccati方程式：https://www.jstage.jst.go.jp/article/sicejl1962/35/12/35_12_971/_pdf 計測と制御 第35巻12号 p971-973(1996)
Windowsを使って組み込みっぽい案件をやった話
組み込みっぽい案件とは
ここでの「組み込み」とは「単一の機能を実現する」の意味合いです。
- 単一の機能を実行する
- PCに電源を入れたときに起動し機能し続ける
「何故，DOSなのか」などはout of controlなので聞かないでください。
Intel NUC Celeron搭載 小型PC BOXNUC6CAYSAJR【Win10/ 2GB DDR3L/ 32GB eMMC プレインストール】
- 出版社/メーカー: インテル
- 発売日: 2017/03/08
- メディア: Personal Computers
- この商品を含むブログを見る
Windows Embedded とかIoTとかあるらしいですがよくわからんし時間的な問題でいろいろと諦めました。各所との連携とか難しいですね。 Windows 10 IoT Enterprise/Windows Embedded OSのサポートと供給期間|Microsoft|製品情報|東京エレクトロンデバイスEmbedded Solution
やったこと
そもそもいろいろガバですがこんな感じに進めました。
Pythonスクリプトをbatから実行
Anacondaで環境作ってましたが，単一動作しかしないので結局Anacondaフォルダ以下のpython.exeとScriptsフォルダにパスを通してしまいました。 batファイルは非常に簡単で以下のものを作成します。
@echo off python XXX.py
管理者権限が必要な操作をbatでやる例
netshでIPアドレスを固定化したかったのですがこのコマンドは管理者権限でやる必要があります。 WindowsのnetshコマンドでTCP/IPのパラメータを設定する：Tech TIPS - ＠IT
こんな感じのbatファイルを管理者権限で実行する必要があるわけです。
@echo off netsh interface ip set address "イーサネット" static 172.16.1.101 255.255.255.0 172.16.1.1
以下のブログにあるように，powershellを使ってやりました。 管理者権限でbatを実行したい時にやッた事
@echo off powershell start-process 実行対象.bat -verb runas
このbatを実行するときには管理者権限が必要とされます。つまり，クリックでOKを選択する必要があることに注意。 電源を入れただけでは実行出来ないのです。
Boot時に特定のプログラムを実行する
「組み込み」なので電源を入れさえすればプログラムが勝手に動き出すようにしてほしいです。
Boot時に特定のプログラムを起動する方法はいくつかありますが，管理者権限の必要不要に応じて２つの方法を紹介します。
管理者権限が不要の時
以下のフォルダにショートカットファイルを置けば起動時に実行してくれるようです。
C:\Users[ユーザ名]\AppData\Roaming\Microsoft\Windows\Start Menu\Programs\Startup
Debian系なら.bashrcの他に/etc/rc.localに置くことに該当します。 linux起動時に自動的に実行するコマンド（プログラム）の設定
管理者権限が必要の時
上記の管理者権限をKickするbatファイルを実行するときに必要です。 タスクスケジューラを使うようですが，今回は結局やりませんでした。 管理者権限が必要な常駐アプリはタスクスケジューラでスタートアップさせる
一応メモとして。
その他の作業
UltraVNCを用いたファイル操作
あと，物理的な接続なしでwifiで今後の作業を進める必要があったので，Winscpのようにファイル転送できるUltraVNCというソフトをインストールしました。 UltraVNCでコンピュータをリモート制御する（サーバ編）：Tech TIPS - ＠IT
Windowsの自動更新をOFF
Windows10が勝手に自動更新をかけてきます。組み込みでいきなり例の青画面が続くのは心臓に悪いので確実に息の根を止めます。
IPアドレスを固定化
IPアドレスを固定化します。
感想
もっと賢くできないですか…ね？？？？
