読者です 読者をやめる 読者になる 読者になる

MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

C++, PythonユーザのためのSwift3基本文法まとめ

swift

Swift実践入門 ── 直感的な文法と安全性を兼ね備えた言語 (WEB+DB PRESS plus)

Swift実践入門 ── 直感的な文法と安全性を兼ね備えた言語 (WEB+DB PRESS plus)

目次

  • 目次
  • はじめに
  • hello world
  • コメント
  • 定数
  • 変数
  • Swiftの基本的な型
  • 演算
  • if文
  • switch文
  • while文
  • for文
  • nilとオプショナル型
  • Array(配列)
  • tuple
  • Set
  • Dictionary
  • function
  • Class and method
  • Computed Property
  • Observer property
  • Inheritance
  • Static function and Static property
  • cast
  • Protocol
  • extension
  • struct
  • 列挙型
  • Generics
  • Subscript
  • 参考資料

 

はじめに

最近、新しい言語の勉強として

Swiftを勉強しているのですが、

developer.apple.com

自分がC++とPythonユーザなので、

C++, Pythonユーザの観点から

Swiftの基本的な文法をまとめておきます。

 

続きを読む

GitHubのサマリを自分のブログやホームページに表示させる方法

javascript

f:id:meison_amsl:20170321090708p:plain

目次

はじめに

最近、すっかり

公開して問題無いコードは

GitHubに公開するようになったのですが、

GitHubでの活動のサマリー

(アカウント名や、スターの多いリポジトリなど)

を自分のブログやホームページに埋め込みたくなり、

色々調べた所、下記のjavascriptライブラリを見つけました。

github.com

 

いくつか修正を加えて、使ってみた所

冒頭の写真ようなGitHubのサマリー(GitHub Widget)を

自分のホームページに簡単に表示させることができたので、

その方法を紹介したいと思います。

  

改良版github-widgetの使い方

前述のjavascriptライブラリですが、

いくつか問題があったので、

下記のリンクのforkし改良したものを現在使用しています。

github.com

  

修正点としては、

  • 表示されるTop repositories (スターの多いリポジトリ)を3から5つに

  • Widgetの末尾に総取得スター数を表示

するようにしました。

 

使い方としては、

上記のリポジトリをWebサーバのどこかにcloneして、

下記のようにsrc/widget.jsをhtmlの中でロードし、

divタグのdata-usernameを自分のGitHubのアカウント名にするだけです。

<div class="github-widget" data-username="AtsushiSakai"></div>
<script src="static/github-widget/src/widget.js"></script>

すると冒頭のスクリーンショットのような

GitHubのアカウント名や、

リポジトリ数、フォロー数、フォロワー数、

スターの多いリポジトリリスト、

フォローボタン、総スター取得数

などを表示したウィジェットが表示されるはずです。

 

ちなみに自分は下記のような個人プロジェクトのホームページに

上記の方法で埋め込んでいます。

atsushisakai.github.io

 

最後に

できれば、

GitHubの緑の草も表示させるようにしたいです。

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

わかばちゃんと学ぶ Gitの使い方入門〈GitHub、Bitbucket、SourceTree〉

わかばちゃんと学ぶ Gitの使い方入門〈GitHub、Bitbucket、SourceTree〉

忘れがちな&間違えがちなPythonデフォルト機能メモ

Python

入門 Python 3

入門 Python 3

目次

  • 目次
  • はじめに
  • printの表示フォーマット
  • クラス変数
  • パッケージの下のモジュールのインポート
  • 集合型(Set)の使い方
  • ジェネレータ
  • リスト内包表記でfilter
  • 参考資料

はじめに

いつも忘れて、ググったり、

長い間上手く使えていなかった

Pythonのデフォルト機能をメモとしてまとめておきます。

 

これらの機能は主に下記の資料を元に勉強しました。

入門 Python 3

入門 Python 3

myenigma.hatenablog.com

 

続きを読む

JITコンパイラライブラリNumbaを使ってPythonコードを劇的に高速化する方法

Python

ふつうのコンパイラをつくろう 言語処理系をつくりながら学ぶコンパイルと実行環境の仕組み

ふつうのコンパイラをつくろう 言語処理系をつくりながら学ぶコンパイルと実行環境の仕組み

目次

  • 目次
  • はじめに
  • Python JITコンパイラライブラリ Numba
  • Numbaのインストール
  • Numbaの速度比較
  • 引数の指定
  • GPUの利用
  • 最後に
  • 参考資料

はじめに

Pythonは言語的にも、

エコシステム的にも大好きなのですが、

いかんせん、複雑で大規模なシステムを作ると、

速度が遅いのが気になる時があります。

 

下記の記事のように、

Pythonの一部をC/C++コードに置き換えてもいいのですが、

やっぱりコードを置き換えるのは面倒ですね。。

(またコードを変換している時に、

バグを埋め込んでしまう可能性もあります)

myenigma.hatenablog.com

 

そこで、色々調べた所、

JIT (Just In Time) コンパイラという技術を利用した

Numbaというライブラリを使うことで、

ある特定のPythonコードを修正無しで、

かなり高速化できたので、

その方法を紹介したいと思います。

 

続きを読む

Raspberry PiにTime Machineサーバを立てる方法

RaspberryPi

APPLE AirMac Time Capsule - 2TB ME177J/A

APPLE AirMac Time Capsule - 2TB ME177J/A

目次

  • 目次
  • はじめに
  • 0. セットアップの方法
  • 1. ハードディスクの設定
  • 2. Raspberry Piでsshログインできるようにする。
  • 3. Raspberry Piで外付けハードディスクをマウントする
  • 4. Netatalkのインストール
  • 5. Netatalkの設定
  • 6. Time Machineの設定
  • 改善したい所
  • 参考資料

はじめに

先日、Raspberry Piを購入して、

セットアップした記事を公開しましたが、

myenigma.hatenablog.com

 

そのRaspberry Piの用途として、

Macのバックアップ機能であるTime Machineのサーバを

Raspberry Piと外付けUSBハードディスクで作ってみたので、

その方法をメモとして残しておきたいと思います。

 

続きを読む

Raspberry PiをMacから操作できるようにセットアップする方法

RaspberryPi

目次

  • 目次
  • はじめに
  • OSのインストール
  • 日本語キーボードレイアウトにする
  • vimを入れる
  • IPアドレスを固定する
  • デフォルトユーザpiのパスワードを変更する
  • sshサーバの設定
  • VNCサーバのインストールと起動
    • VNCサーバのセットアップ
    • VNCクライアントのセットアップ
  • セキュリティの初期設定をする
  • 最後に
  • 参考資料

 

はじめに

ネットを見てたら、

Raspberry Pi 3のセットが、

かなり安くで売っていたので、

つい買ってしまいました。

 

ブームからはかなりズレていますが、

OSをインストールし、

ディスプレイやキーボード, マウスに繋げなくても、

Macから操作できるようにセットアップしたので、

その方法をメモしておきたいと思います。

 

続きを読む

飛行機移動を快適にするためにしていること

Trip

パラダイス山元の飛行機のある暮らし―――年間最多搭乗1022回「ヒコーキの中の人」が贈る空の過ごし方

パラダイス山元の飛行機のある暮らし―――年間最多搭乗1022回「ヒコーキの中の人」が贈る空の過ごし方

目次

 

はじめに

先日、知り合いから

飛行機移動する時ってどんな工夫してる?

と聞かれました。

その時にはあまりすぐに答えられなかったのですが、

その後、落ち着いて考えてみると

約10年前初めて飛行機に乗ってから

色々試行錯誤の結果、

様々な工夫をしていることに気がついたので、

今後同じことを聞かれた時用に、

メモとして残しておきたいと思います。

 

iPadに映画やドラマをダウンロードしておく

自分は機内では、

iPadで映画を見ていることが多いです。

 

使っているiPadはiPad Proの9.7インチです。

9.7インチですと大抵の場合、

充電無しで使い続けることができます。

 

映画を見る時は、手で持つか、

テーブルにアップル公式のカバーを使って固定して見ています。

 

映画や動画のダウンロードは、

以前はDVDをリッピングして、

それをiTunesでiPadに入れて見ていましたが、

最近はAmazon Prime VideoとNetflixのiPadアプリが、

動画のダウンロード&オフライン視聴に対応したので、

事前に映画を大量にダウンロードして、

飛行機の中で見続けています。

Amazonプライム・ビデオ
Amazonプライム・ビデオ
開発元:AMZN Mobile LLC
無料
posted with アプリーチ
Netflix
Netflix
開発元:Netflix, Inc.
無料
posted with アプリーチ

(HuluとYoutubeは早くオフラインダウンロード機能つけて欲しいです。。。)

 

国際線の機内では、前の席のモニターで映画が見れることも多いですが、

解像度が悪かったり、見たい映画が無かったり、

再生時間を選べなかったりなど、色々不満が多いので、

最近はもっぱらiPadで見るようになりました。

 

また、以前は機内でプログラミングなども試してみましたが、

ネットが無いとどうしても効率が上がらなかったので、

止めてしまいました。

 

前から見たかった映画を2,3本見れば、

長い飛行機移動も結構あっという間に過ぎてしまう気がします。

 

ちなみに大量の動画を入れておくので、

iPad Proは120GBのモデルを使っています。

 

iPadに電子書籍をダウンロードしておく

上記の動画に疲れた時は、

iPadで本を読むことも多いです。

大抵はKindleか、

Kindle電子書籍リーダー:人気小説や無料漫画、雑誌も多数
Kindle電子書籍リーダー:人気小説や無料漫画、雑誌も多数
開発元:AMZN Mobile LLC
無料
posted with アプリーチ

事前にダウンロードしたpdfの論文などを

Good Readerで見ています。

 

また、機内は狭いですし、

旅の疲れで集中できないことが多いので、

Kindleの気になる漫画を一気に読んでいることも多いです。

前回の長距離移動では、NARUTOを全巻読破しました。

 

また、新しい場所への旅行の場合は、

機内で地球の歩き方などのガイドブックを

Kindleで読んでいることも多いです。

 

iPadにPodcastをダウンロードしておく

動画や読書に疲れた時は、

iPadでPodcastを事前にダウンロードしておいて、

それを聞いています。

(大抵はそのまま寝てしまいますが笑)

 

Podcastを聞く時は、

iPhoneでも使っている

Overcastというアプリを使っています。

Overcast: Podcast Player

Overcast: Podcast Player

  • Overcast Radio, LLC
  • ニュース
  • 無料

 

BoseのQC20で、飛行機の騒音からの疲労を軽減する

自分は、飛行機の席に座った瞬間から、

飛行機を降りるまで、

下記の記事でも紹介したBoseのQuiet Comport20を使って、

常にノイズキャンセリングをしています。

myenigma.hatenablog.com

 

飛行機は常にエンジンや飛行時の振動音がしていますが、

それを聞き続けているとかなり体と頭が疲労するようです。

自分はこのBoseのQuiet Comport20を使い始めてから、

確かに、飛行機移動の後の疲労が少なくなった気がしています。

 

上記で説明した、映画やPodcastなども

すべてこのQC20でノイズキャンセリングしながら、

音声を聞いています。

 

QC20の詳しい説明は上記の記事を参照下さい。

 

機内食は食べない

自分は、どうしても飛行機の機内食が美味しいと感じられません。。

気圧の問題もあると思いますが、

そもそも食欲もあまりわきませんし、

逆に機内食の匂いで気持ち悪くなってしまいます。

 

そこで、自分はほとんど機内食を食べないようにしています。

食べる時も果物や甘いお菓子のみにしています。

 

そこで、出発前にかならず何かしら食べておきますし、

目的地についた時には、

現地の美味しいものをすぐに食べるようにしています。

 

乾燥と匂い対策用のマスクを準備する

飛行機の中は、

飛行機の機体のメンテナンスの問題で、

かなり湿度を低く保っているようです。

そのせいか、飛行機に長時間のると、

すぐに喉が痛くなってしまいます。

 

また前述の通り、

自分は機内食の匂いがあまり得意ではありません。

そこで、下記のような

加湿機能+匂い機能の2つを有したマスクを必ず付けています。

これを使うことで、

喉の痛みや、風邪を機内でうつされる可能性を軽減できますし、

いい匂いが長時間持続するので、

機内の快適さがかなり向上しました。

 

搭乗前に飲み物を二本買っておく

機内では、ドリンクはCAの人に言えばもらえますが、

ジュースなどはコップに入れてもらう形になるため、

少しずつ飲みたい時などに、

狭い機内では邪魔になります。

 

そこでかならず、

飛行機に乗る前に飲み物を2本、

自分で買うようにしました。

また、できるだけ体力を温存するために、

日本発の時は、ポカリスエットを買っています。

外国発の時は大抵、水とオレンジジュースが多いです。

 

好きな時に好きな量、

飲み物を飲めるのはかなり快適だと思います。

 

長時間の場合はスウェットに着替えておく

飛行機が長時間になる時は、

ジーパンなどを履いていると、

疲れが蓄積する気がしていたので、

最近は、搭乗の寸前に下だけスウェットに着替えています。

これにより飛行機の疲労が軽減されている気がします。

 

また一応、他の人の目もあるので、

あまりダボダボのスウェットではなく、

下記のような、

多少スキニータイプのスウェットを履いています。

最近は家から、スウェットで空港に向かうようになりました笑

 

寝る時はアイマスクをする

となりの人が、

読書灯などを付けて、

仕事をしている可能性もあるので、

必ずアイマスクを準備しています。

 

ちなみに自分の場合、

時差ボケ対策のために、

飛行機に乗った瞬間から、

現地の時間通りに寝るようにしています。

 

トイレでスクワットする

長時間飛行機で座っていると、

体が強張って、ダルくなってきたり、

最悪の場合、エコノミー症候群になったりします。

しかし、体を動かすスペースはあまり無いので、

自分の場合、定期的に飛行機のトイレで

スクワットをしています。

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

パラダイス山元の飛行機のある暮らし―――年間最多搭乗1022回「ヒコーキの中の人」が贈る空の過ごし方

パラダイス山元の飛行機のある暮らし―――年間最多搭乗1022回「ヒコーキの中の人」が贈る空の過ごし方

シンプルなMPC最適化モデリングの数式導出とPythonサンプルコード

Model Predictive Control: Classical, Robust and Stochastic (Advanced Textbooks in Control and Signal Processing)

Model Predictive Control: Classical, Robust and Stochastic (Advanced Textbooks in Control and Signal Processing)

目次

はじめに

先日、

凸最適化のモデリングツールCVXやCVXPY, CVXGENを紹介しましたが、

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

シンプルなモデル予測制御MPCの最適化問題であれば、

モデリングツールを使わずに、MPCの最適化式を

二次凸最適の標準形に変形して、

ソルバーのみで解くことができます。

myenigma.hatenablog.com

 

モデリングツールは、

任意の凸最適化問題を解くことができるので、

最適化問題を設計している段階では便利ですが、

モデリングコードを解析するために計算時間がかかってしまうので、

最適化システムを高速で回したい場合は、

モデリングの部分を自作実装したほうが良いです。

 

そこで、今回は非常にシンプルなMPC最適化問題を、

下記の資料をベースに

数式導出し、Pythonで実装してみましたので紹介したいと思います。

 

シンプルなMPC問題

今回の例では、

下記の標準的なMPC問題を解きたいとします。

f:id:meison_amsl:20170205052109p:plain

x0は初期状態であるとします。

 

また今回の例題では、

下記の安定なシステムである

MPCパラメータを使うことにします。

A = np.matrix([[0.8, 1.0], [0, 0.9]])
B = np.matrix([[-1.0], [2.0]])
(nx, nu) = B.shape
 
N = 10  # number of horizon
Q = np.eye(nx)
R = np.eye(nu)

 

入力制約のみのMPCモデリングの数式導出とPythonサンプルコード

数式導出

まず初めに、

入力制約条件のみのMPCモデリングの数式を導出したいと思います。

入力プロファイルを一つのベクトルとすると、

下記のように、線形方程式で状態ベクトルのプロファイルを計算することができます。

f:id:meison_amsl:20170207063125p:plain

 

そこで、上記の線形方程式の考えから、

上記のMPCの問題は、

下記のような二次凸問題に変形することができます。

f:id:meison_amsl:20170207063349p:plain

 

上記の式において、入力制約を追加しない場合は、

MPCの制約無しの解となります。

入力の最大値や最小値などの、

入力制約を入れたい場合は、

下記の式のようなに二次凸最適化の不等式制約を入れることで、

入力値の制約条件を追加することができます。

f:id:meison_amsl:20170207072420p:plain

f:id:meison_amsl:20170207072427p:plain

f:id:meison_amsl:20170207072438p:plain

 

Pythonコード

上記の数式をPythonのnumpyで実装したのが、

下記のコードになります。

それぞれのMPCパラメータと、

入力制約(オプション)を入力すると、

最適なMPCの状態プロファイルと、

入力プロファイルを計算してくれます。

 

def opt_mpc_with_input_const(A, B, N, Q, R, P, x0, umax=None, umin=None):
    """
    optimize model predictive control only with input constraints
    (if you want to solve a problem with state constraints, you can use opt_mpc_with_state_const())
    return
        x: state
        u: input
    """
    (nx, nu) = B.shape

    # calc AA
    Ai = A
    AA = Ai
    for i in range(2, N + 1):
        Ai = A * Ai
        AA = np.vstack((AA, Ai))

    # calc BB
    AiB = B
    BB = np.kron(np.eye(N), AiB)
    for i in range(1, N):
        AiB = A * AiB
        BB += np.kron(np.diag(np.ones(N - i), -i), AiB)

    RR = np.kron(np.eye(N), R)
    QQ = scipy.linalg.block_diag(np.kron(np.eye(N - 1), Q), P)

    H = (BB.T * QQ * BB + RR)

    gx0 = BB.T * QQ * AA * x0
    P = matrix(H)
    q = matrix(gx0)

    if umax is None and umin is None:
        sol = cvxopt.solvers.qp(P, q)
    else:
        G = np.zeros((0, nu * N))
        h = np.zeros((0, 1))

        if umax is not None:
            tG = np.eye(N * nu)
            th = np.kron(np.ones((N * nu, 1)), umax)
            G = np.vstack([G, tG])
            h = np.vstack([h, th])

        if umin is not None:
            tG = np.eye(N * nu) * -1.0
            th = np.kron(np.ones((N * nu, 1)), umin * -1.0)
            G = np.vstack([G, tG])
            h = np.vstack([h, th])

        G = matrix(G)
        h = matrix(h)

        sol = cvxopt.solvers.qp(P, q, G, h)

    u = np.matrix(sol["x"])

    # recover x
    xx = AA * x0 + BB * u
    x = np.vstack((x0.T, xx.reshape(N, nx)))

    return x, u

 

凸二次ソルバーとしては、

以前紹介したcvxoptを使いました。

myenigma.hatenablog.com

 

下記はモデリングツールcvxpyを使って、

MPCを解いた場合と、

上記の関数を使った場合の結果です。

 

下記のシミュレーションでは、

下記の入力制約を追加しています。

f:id:meison_amsl:20170207081844p:plain

 

両方の出力で同じ結果が得られていることがわかります。

f:id:meison_amsl:20170207073417p:plain

 

状態・入力制約を含んだMPCモデリングの数式導出とPythonサンプルコード

数式導出

前述の方法だと、入力制約は追加できますが、

状態制約が付けられないので、

別の方法で数式導出する必要があります。

 

まず初めに入力プロファイルと状態プロファイルを一つのベクトルZにまとめた場合、

下記の二次凸最適化の標準形に変形することができます。

f:id:meison_amsl:20170207075337p:plain

f:id:meison_amsl:20170207075705p:plain

 

入力制約と、状態制約は先程と同じように、

不等式制約とすることで、計算することができます。

f:id:meison_amsl:20170207072420p:plain

f:id:meison_amsl:20170207080959p:plain

f:id:meison_amsl:20170207080207p:plain

 

Pythonコード

上記の数式をPythonのnumpyで実装したのが、

下記のコードにです。

 

それぞれのMPCパラメータと、

入力制約, 状態制約(それぞれオプション)を入力すると、

最適なMPCの状態プロファイルと、

入力プロファイルを計算してくれます。

 

def opt_mpc_with_state_constr(A, B, N, Q, R, P, x0, xmin=None, xmax=None, umax=None, umin=None):
    """
    optimize MPC problem with state and (or) input constraints

    return
        x: state
        u: input
    """
    (nx, nu) = B.shape

    H = scipy.linalg.block_diag(np.kron(np.eye(N), R), np.kron(np.eye(N - 1), Q), np.eye(P.shape[0]))

    # calc Ae
    Aeu = np.kron(np.eye(N), -B)
    Aex = scipy.linalg.block_diag(np.eye((N - 1) * nx), P)
    Aex -= np.kron(np.diag([1.0] * (N - 1), k=-1), A)
    Ae = np.hstack((Aeu, Aex))

    # calc be
    be = np.vstack((A, np.zeros(((N - 1) * nx, nx)))) * x0
    #  print(be)

    # === optimization ===
    P = matrix(H)
    q = matrix(np.zeros((N * nx + N * nu, 1)))
    A = matrix(Ae)
    b = matrix(be)

    if umax is None and umin is None:
        sol = cvxopt.solvers.qp(P, q, A=A, b=b)
    else:
        G, h = generate_inequalities_constraints_mat(N, nx, nu, xmin, xmax, umin, umax)

        G = matrix(G)
        h = matrix(h)

        sol = cvxopt.solvers.qp(P, q, G, h, A=A, b=b)

    fx = np.matrix(sol["x"])

    u = fx[0:N * nu].reshape(N, nu).T

    x = fx[-N * nx:].reshape(N, nx).T
    x = np.hstack((x0, x))

    return x, u

def generate_inequalities_constraints_mat(N, nx, nu, xmin, xmax, umin, umax):
    """
    generate matrices of inequalities constrints

    return G, h
    """
    G = np.zeros((0, (nx + nu) * N))
    h = np.zeros((0, 1))
    if umax is not None:
        tG = np.hstack([np.eye(N * nu), np.zeros((N * nu, nx * N))])
        th = np.kron(np.ones((N * nu, 1)), umax)
        G = np.vstack([G, tG])
        h = np.vstack([h, th])

    if umin is not None:
        tG = np.hstack([np.eye(N * nu) * -1.0, np.zeros((N * nu, nx * N))])
        th = np.kron(np.ones((N, 1)), umin * -1.0)
        G = np.vstack([G, tG])
        h = np.vstack([h, th])

    if xmax is not None:
        tG = np.hstack([np.zeros((N * nx, nu * N)), np.eye(N * nx)])
        th = np.kron(np.ones((N, 1)), xmax)
        G = np.vstack([G, tG])
        h = np.vstack([h, th])

    if xmin is not None:
        tG = np.hstack([np.zeros((N * nx, nu * N)), np.eye(N * nx) * -1.0])
        th = np.kron(np.ones((N, 1)), xmin * -1.0)
        G = np.vstack([G, tG])
        h = np.vstack([h, th])

    return G, h

 

下記のグラフはcvxpyを使って、計算した結果の比較です。

下記のシミュレーションでは、

下記の入力制約と状態制約を追加しています。

f:id:meison_amsl:20170207081754p:plain

 

モデリングツールを使った結果と

同じ結果が得られていることがわかります。

f:id:meison_amsl:20170207081429p:plain

 

MPCのTrackingサンプルコード

上記のMPCのコードは、

すべて状態を0にするように制御する

いわゆるRegulatorのコードでしたが、

ある目標値に状態値を近づけるように制御する

いわゆるMPC Tracking の方法を紹介したいと思います。

数式の導出

Trackingの数式の導出に関しては、

下記の書籍の2章と3章を参考にして下さい。

モデル予測制御―制約のもとでの最適制御

モデル予測制御―制約のもとでの最適制御

Predictive Control with Constraints

Predictive Control with Constraints

Python サンプルコード

"""
MPC tracking sample code
author: Atsushi Sakai
"""

import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg

from cvxopt import matrix
import cvxopt

DEBUG_ = False

def get_mat_psi(A, N):
    psi = np.matrix(np.zeros((0, A.shape[1])))

    for i in range(1, N + 1):
        psi = np.vstack((psi, A ** i))

    return psi


def get_mat_gamma(A, B, N):
    (nx, nu) = B.shape
    gamma = np.zeros((nx, nu)) + B

    for i in range(1, N):
        tmat = (A ** i) * B + gamma[-nx:, :]
        gamma = np.vstack((gamma, tmat))

    return gamma


def get_mat_theta(A, B, N):
    AiB = B
    (nx, nu) = B.shape
    theta = np.kron(np.eye(N), AiB)

    tmat = np.zeros((nx, 0))

    for i in range(1, N):
        t = np.zeros((nx, nu)) + B
        for ii in range(1, i + 1):
            t += (A ** ii) * B
        tmat = np.hstack((t, tmat))

    for i in range(1, N):
        theta[i * nx:(i + 1) * nx, :i] += tmat[:, -i:]

    return theta


def model_predictive_control(A, B, N, Q, R, T, x0, u0):

    (nx, nu) = B.shape

    du = np.matrix([0.0] * N).T

    psi = get_mat_psi(A, N)
    gamma = get_mat_gamma(A, B, N)
    theta = get_mat_theta(A, B, N)

    QQ = scipy.linalg.block_diag(np.kron(np.eye(N), Q))
    RR = scipy.linalg.block_diag(np.kron(np.eye(N), R))

    H = theta.T * QQ * theta + RR
    g = - theta.T * QQ * (T - psi * x0 - gamma * u0)

    P = matrix(H)
    q = matrix(g)
    sol = cvxopt.solvers.qp(P, q)
    du = np.matrix(sol["x"])

    fx = psi * x0 + gamma * u0 + theta * du

    ffx = fx.reshape(N, nx)
    ffx = np.vstack((x0.T, ffx))

    u = np.cumsum(du).T + u0

    return ffx, u

def test1():
    print("start!!")
    A = np.matrix([[0.8, 1.0], [0, 0.9]])
    B = np.matrix([[-1.0], [2.0]])
    (nx, nu) = B.shape

    N = 50  # number of horizon
    Q = np.diag([1.0, 1.0])
    R = np.eye(nu)

    x0 = np.matrix([2.0, 1.0]).T
    u0 = np.matrix([0.0])

    T = np.matrix([1.0, 0.25] * N).T

    x, u = model_predictive_control(A, B, N, Q, R, T, x0, u0)

    # test
    tx = x0
    rx = x0
    for iu in u[:, 0]:
        tx = A * tx + B * iu
        rx = np.hstack((rx, tx))

    if DEBUG_:
        plt.plot(x[:, 0])
        plt.plot(x[:, 1])
        plt.plot(u[:, 0])
        plt.grid(True)
        plt.plot(rx[0, :].T, "xr")
        plt.plot(rx[1, :].T, "xb")

        plt.show()

    for ii in range(len(x[0, :]) + 1):
        for (i, j) in zip(rx[ii, :].T, x[:, ii]):
            assert (i - j) <= 0.0001, "Error" + str(i) + "," + str(j)

    target = T.reshape(N, nx)
    for ii in range(len(x[0, :]) + 1):
        assert abs(x[-1, ii] - target[-1, ii]) <= 0.3, "Error"

 

上記のコードを実行すると、

下記のようにMPCで目標値に追従させる

入力プロファイルを計算することができます。

f:id:meison_amsl:20170217055226p:plain

 

すべてのPythonサンプルコード

すべてのPythonサンプルコードは、

下記のGithubリポジトリで公開しています

github.com

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

Model Predictive Control: Classical, Robust and Stochastic (Advanced Textbooks in Control and Signal Processing)

Model Predictive Control: Classical, Robust and Stochastic (Advanced Textbooks in Control and Signal Processing)

『優秀なエンジニアになる方法』を読んで

Book

スゴ腕のエンジニアになる ─ITエンジニアの“今

スゴ腕のエンジニアになる ─ITエンジニアの“今"と“これから"

目次

はじめに

下記の『優秀なエンジニアになる方法』という

電気情報通信学会誌を読みました。

 

色々、自分が出来ていなくて

心が痛い内容が多かったので、

戒めのためにいくつかの内容を参照しておきます。

 

優秀な自発性

  • 本来の仕事をこなしながら、それ以外の仕事を探すこと (ライがパソコンのソフトウェアをインストールしたように)

  • 同僚のためまたは組織のため、手伝ってあげること (ライが同僚のプログラムを直すために手伝ってあげたように)

  • だれの責任でもない、突然表れる仕事を自分から進んでやり、立派にこなすこと

  • プロジェクトを始めたら、実現するまで根気よくやり遂げること (ライが残業をして新しいソフトウェアをイン ストールしたように)

 

適切な発表

普通のエンジニアがよく失敗するのは、

ただ情報を伝える ことだけを考えて、

その後の反響をまでも考えた

メッセージ の伝達ができないことである。

聴衆が変っても発表のやり方 を変えようとしない

 

正しい仕え方

優秀なエンジニアはもっと有利な仕え方を早くから覚える。

それは良い右腕になることである。

つまり自分で点を取 るより、アシストをする方が重要なことが多い。

組織やリーダーを成功させるように手伝いながら、

やるべき事とそのや り方について独立的で、かつ鋭い判断力を発揮する。

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

最適化問題における双対問題とは?

optimization

経済数学入門16:線型計画法(双対性)

経済数学入門16:線型計画法(双対性)

目次

  • 目次
  • はじめに
  • 双対問題とは?
  • なぜ双対問題を考えるのか?
  • 非線形問題における双対問題
  • 線形計画法の双対問題を導出してくれるPythonライブラリdual
  • 様々な最適化問題の双対問題まとめ
    • 線形計画問題 1
    • 線形計画問題 2
    • Second-order cone programming
  • より詳しい双対問題の説明を知りたい人は
  • 参考資料

はじめに

最適化問題を勉強していると、

双対問題[そうついもんだい] (dual problem)という言葉が出てきます。

 

実は最近の凸最適化のソルバーなどが、

高速に問題を解けるようになったのは、

この双対問題という考えを使うことにより、

実現していることが多いようです。

 

今回、この双対問題の概要と、

なぜ双対問題を考えると、

最適化を解きやすくなるのかを

紹介したいと思います。

 

続きを読む