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とかだと解けたりするんでしょうか？イマイチわかりきっていないところがあります。

2018-10-07

Latex 表組みのサイズを微調整

tex

表のサイズの微調整
GUIに依る表の作成

表のサイズの微調整

従来，Latexの表が大きすぎた時 \smallなどのコマンドで表の文字を小さくしてサイズを収めるといったことをしていました。しかしこの調整は非常に表を醜くします。

正しくは\resizebox{}環境を用いるべきです。 \columnwidthで横幅ぴったりに合わせます。2カラムぶち抜くときは\textwidthを用いると良いです。

% \usepackage{graphics}
\resizebox{\columnwidth}{!}{%
%table 環境を中に
}

https://tex.stackexchange.com/questions/10863/is-there-a-way-to-slightly-shrink-a-table-including-font-size-to-fit-within-th

GUIに依る表の作成

Latexの表の作成は最もだるい事の一つですが，Tablegeneratorを用いることで簡略化出来ます。 Excelとかで表を作ってしまってGUIで整形，その後にTexに変換して貼り付けるのが最近のやり方です。

www.tablesgenerator.com

お試しあれ。

2018-09-30

Google ColaboratryでGPU（TPU）付の機械学習環境が使えるらしい

雑記 NVIDIA GoogleDrive Keras

趣味人としてはそろそろ独自のコードを書いて公開したいところ。（なお本業…） [: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マシンを持っていなくても（環境構築しなくても）機械学習の学習が出来るのは初学者にはかなり良いと思います。特に環境構築や機材は初学者にとって大きな壁なので回線さえあればきちんとした環境で実行できるというのは教育効果が非常に高そうです。

とりあえず，コードを試して満足したい方はここから初めて良いと思いますが，一方で今後ちゃんとやっていくならどこかで環境構築はしておくべきだと思います。結局自分でデータ集めてモデル組まないと意味ないので僕もこんなサンプルコードばかり実行してないで自分のやりたいプロジェクトを進めていこうとおもいます。

2018-09-21

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')

2018-09-21

「Linear＝線形」ではない

勉強勉強-数学

俺たちは雰囲気で数学をやっている… ”Linear” は線形ないし一次（Affine）の意味でよく出ますがごっちゃになっちゃうので気をつけましょう。

f:id:ossyaritoori:20180921090015p:plain

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]ではキチンと「アファインな関数」と言っていますね。 $https://wikimedia.org/api/rest_v1/media/math/render/svg/9146e5a3c0820b3321a9624832c9724496582f15$

[1]でもこう言っています。線形とは…？？？ f:id:ossyaritoori:20180921095033p:plain

行列一次不等式って言ったほうが良いのでは？よくわからなくなってきました。

[1]制御におけるRiccati方程式：https://www.jstage.jst.go.jp/article/sicejl1962/35/12/35_12_971/_pdf 計測と制御第35巻12号 p971-973(1996)

2018-09-20

Windowsを使って組み込みっぽい案件をやった話

Windows Python 雑記

組み込みっぽい案件とは
- やったこと
Pythonスクリプトをbatから実行
- 管理者権限が必要な操作をbatでやる例
Boot時に特定のプログラムを実行する
- 管理者権限が不要の時
- 管理者権限が必要の時
その他の作業
感想

組み込みっぽい案件とは

ここでの「組み込み」とは「単一の機能を実現する」の意味合いです。

単一の機能を実行する
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ファイルを作成し，起動時に動くようにした
権限が必要な操作とかでちょっと困った

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アドレスを固定化します。

Windows 10 IPアドレス固定方法

感想

もっと賢くできないですか…ね？？？？

2018-09-12

MATLABでGPSデータから距離を計算

matlab

GPSの緯度経度の行列から移動距離を計算するアルゴリズムがmatlab公式に用意されてないので自分で書きました。探せばあるのかもしれません。

今回の手法ではヒュベニの公式というのを使います。詳細な説明等はここのサイトがよいです。

画像も借りてきました。

関数は以下の通り。整備したらGithubにでも投げます。

function D = calcGPS_dist(LON,LAT)

len = length(LON);

D = 0;
for i = 1:len-1
    D = D+hubeny([LON(i),LAT(i)],[LON(i+1),LAT(i+1)]);
end

end

function d = hubeny(L1,L2)
% input: pair of [LON,LAT]
% output : distance [m]
% e2 = 0.00669437999019758;
a = 6378137;
b = 6356752.314140;
e2 =   (a^2 - b^2) / a^2;

mu_y = (L1(2) + L2(2))/2;
W = sqrt(1-e2*(sind(mu_y))^2);
M = a*(1-e2)/(W^3);
N = a/W;

dx = (L1(1)-L2(1))/180*pi;
dy = (L1(2)-L2(2))/180*pi;

d = sqrt( (dy*M)^2+(dx*N*cosd(mu_y))^2 );

end

ハマリポイントは角度をradに変換するのを忘れていたところです。

Software Installation

SeDuMi

SDP3

YALMIP