システムとモデリング

SysML, Matlab/Simulink,他モデリング全般について。たまにプロジェクトマネジメント。

微分方程式を行列で差分化して解く?

ストラング先生の計算理工学に、行列で微分方程式を差分化して解く方法が載っていました。 浅学のためこういった解法は初めてであまり理解もできていませんが、練習してみたいと思います。

世界標準MIT教科書 ストラング:計算理工学

世界標準MIT教科書 ストラング:計算理工学

例として下記のような微分方程式を解いてみます。

$$ \begin{equation} \frac{dy}{dx}=3x^{2}, y(0) = 0 \end{equation} $$

この方程式の解は y = x^{3}ですが、これを行列で差分化して確認してみます。

yを前進差分で差分化する行列は例えば以下のようなものになります。 $$ A= \begin{bmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{bmatrix} $$

これにyをかけると以下のようになります。このように前進差分を取ることができます。 $$ \begin{bmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} y(0) \\ y(1) \\ y(2) \\ y(3) \end{bmatrix}  = \begin{bmatrix} y(0) \\ y(1) - y(0) \\ y(2) - y(1) \\ y(3) - y(2) \end{bmatrix} $$

今回は右辺が 3x^{2}なので $$ \boldsymbol{y} = \begin{bmatrix} y(0) \\ y(1) \\ y(2) \\ y(3) \end{bmatrix} $$ とすると 行列を使った方程式は下記のようになります。 $$ \begin{bmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{bmatrix} \boldsymbol{y}  = \begin{bmatrix} 0 \\ 3 \\ 12 \\ 27 \end{bmatrix} $$

これを解くと $$ \boldsymbol{y}  = \begin{bmatrix} 0 \\ 3 \\ 15 \\ 42 \end{bmatrix} $$ が得られます。ただしこの値は y=x^{3}とはまったく一致しません。 おそらくステップ数が4しかないために誤差が大きかったものと思われます。

次はステップ数を増やして解法を確認してみます。

今回はここまでになります。

書評:PROJECT MANAGEMENT FOR THE UNOFFICIAL PROJECT MANAGER

プロジェクトマネジメント洋書の書評になります。

最近出版される技術系和書は入門書が多く、洋書と内容面で隔たりを感じています。プロジェクトマネジメント分野の洋書を読むのは今回が初めてですがAmazonでこの本が人気だったのでKindleで購入してみました。Kindleではわからない単語をその場で和訳できるので比較的サクサク読めました。もともとの英語も平易な単語が多く全体的に読み易かったです。

 

Project Management for the Unofficial Project Manager: A FranklinCovey Title

Project Management for the Unofficial Project Manager: A FranklinCovey Title

 

 内容

この本の内容は次の言葉に集約されています。

People + Process = Sccess 

 噛み砕いて言うと、リーダーシップ(People)だけでも、技法(Process)だけでもプロジェクトを成功に導くことはできない、その両方が必要、という内容です。

 リーダーシップとしては次の4つの態度を取ることが必要と説いていて、本書で繰り返し繰り返し説明されています。

  • Demonstrate respect
  • Listen first
  • Clarify expectation
  • Practice accountability

プロセスでは各フェーズごとに必要なスキルセットとツールとその例が物語形式で紹介されています。内容はPMBOKに準拠しています。

  • Initiate
  • Plan
  • Execute
  • Monitor and Control
  • Close

ツールはステークスホルダー分析、要件定義、WBSリスク管理など標準的なものですがドキュメントの例が載っているので実際に再現することが可能なようになっています。

最後に、タイトルの「Unofficial Project Manager」とは「プロジェクトマネジメントの教育を受けないままプロジェクトマネジメントを任された人」のことです。そういう意味では日本のプロジェクトマネージャーの大半が「Unofficial Project Manager」なのではないでしょうか。

和書と比べてどうか?

 私の勝手なプロジェクトマネジメント系和書の印象ですが、だいたい「ただのPMBOKの解説書」と「著者固有の経験談」の両極端に大別されると思います。前者は実施例に乏しく全く実用的でない。後者は著者の経験が技術として抽象化されておらずただの場当たり的なノウハウが貯まるだけ、加えてIT企業の例が大半で他業種に活用し辛いといった印象です。プロジェクトマネジメント本として求められるのはこの両者の中間と思うのですが、この本はまさにそれに該当する本だと思います。このような位置付けの本は和書ではなかなかお目にかかれないと思います。あえて和書で挙げるなら、下記の本はプロセスに重点を置いていますが、かなり近いと思います。

 

 ただし、わざわざ英語でこの本を読むべきかどうかと言われると微妙です。日本人著者でない翻訳書には優秀な本も多いので、まずそちらを読まれたほうが良いと思います。一例ですが下記の本などはかなりおすすめです。

 

アート・オブ・プロジェクトマネジメント ―マイクロソフトで培われた実践手法 (THEORY/IN/PRACTICE)

アート・オブ・プロジェクトマネジメント ―マイクロソフトで培われた実践手法 (THEORY/IN/PRACTICE)

 

 

今回はここまでにします。

 

 

Juliaで最適化計算

今回はJuliaのパッケージ"JuMP.jl"を使って最適化問題を解いてみます。

JuMP.jlについて

Juliaの最適化計算パッケージです。 www.juliaopt.org

特徴は下記のようですね。

  • 本来の数式に近い記述が可能
  • 速い
  • 記述がsolverに依存しない

例題を実行する。

下記ドキュメントにある例題を実行してみます。

JuMP — Julia for Mathematical Optimization — JuMP -- Julia for Mathematical Optimization 0.18 documentation

using JuMP
using Clp

m = Model(solver = ClpSolver())
@variable(m, 0 <= x <= 2 )
@variable(m, 0 <= y <= 30 )

@objective(m, Max, 5x + 3*y )
@constraint(m, 1x + 5y <= 3.0 )

print(m)

status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("y = ", getvalue(y))

これを実行すると下記が出力されます。

Max x1 + 2 x2 + 5 x3
Subject to
 -x1 + x2 + 3 x3 <= -5
 x1 + 3 x2 - 7 x3 <= 10
 0 <= x1 <= 10
 x2 >= 0
 x3 >= 0
Optimal Solutions:
x1 = 10.0
x2 = 2.1875
x3 = 0.9374999999999999

問題なく最適化計算できたようです。

練習問題を解く

筆者お気に入りのこの教科書の問題を解いてみます。

数理モデリング入門 ―ファイブ・ステップ法― 原著第4版

数理モデリング入門 ―ファイブ・ステップ法― 原著第4版

例題3.4

植え付け可能な625エーカーの土地を持つ家族農業経営を考えよう。植え付けを考慮されている作物はトウモロコシ、小麦、オート麦である。1,000エーカーフートの水が灌漑用に必要で、家族は1週あたり300時間仕事ができると見積もられる。その他のデータは下表にまとめられている。最大利益を得るために各作物の植え付け面積を決定しよう。

必要量 (エーカーあたり) トウモロコシ 小麦 オート麦
灌漑用水(エーカーフート) 3.0 1.0 1.5
労働力(人時/週) 0.8 0.2 0.3
収益(ドル) 400 200 250
変数 説明
x1 トウモロコシの植え付け面積 (>=0)
x2 小麦植え付け面積 (>=0)
x3 オート麦植え付け面積 (>=0)

これをJuliaで実装した場合下記のようなコードになります。

using JuMP
using Clp

m = Model(solver=ClpSolver())

@variable(m, x1 >=0)
@variable(m, x2 >=0)
@variable(m, x3 >=0)


@objective(m, Max, 400x1 + 200x2 + 250x3)

@constraint(m, 3.0x1 + 1.0x2 + 1.5x3 <=1000)
@constraint(m,  0.8x1 + 0.2x2 + 0.3x3 <=300)
@constraint(m,  x1 + x2 + x3 <= 625)

print(m)
status = solve(m)

println("Optimal Solutions:")
println("x1 = ", getvalue(x1))
println("x2 = ", getvalue(x2))
println("x3 = ", getvalue(x3)

このとき出力は以下のようになります。

Max 400 x1 + 200 x2 + 250 x3
Subject to
 3 x1 + x2 + 1.5 x3 <= 1000
 0.8 x1 + 0.2 x2 + 0.3 x3 <= 300
 x1 + x2 + x3 <= 625
 x1 >= 0
 x2 >= 0
 x3 >= 0
Optimal Solutions:
x1 = 187.49999999999997
x2 = 437.5
x3 = 0.0

教科書の解答では x1=187.5, x2=437.5, x3=0.0 なので、正解ですね。

今日はここまでにします

JuliaでModelicaライクなシミュレーション

今回はJuliaのパッケージ"Modia.jl"を使用してModelicaのような非因果モデリングでシミュレーションしてみます。

Modelicaについては以前の記事を参照してください。

otepipi.hatenablog.com

Modia.jlについて

Modia.jlはJulia言語による方程式ベースのシステムモデリング環境です。作者はModelicaを使用した経験をもとにこのライブラリを制作しており、Modelicaと非常に似た作りになっています。

github.com

https://modiasim.github.io/Modia.jl/slides/Systems-Modeling-and-Programming-Slides.pdf

https://modiasim.github.io/Modia.jl/slides/Modia-JuliaCon-2017.pdf

ModelicaではなくあえてModiaを使う理由として制作者は以下の点を挙げています。

  • Modelicaの言語としての発展速度が衰えてきている。
  • Modelicaの言語仕様は膨大で理解が難しい
  • Modelicaはデータの取り扱いをはじめとしたアルゴリズム部分が弱い

Juliaは元から科学計算用に開発された高速な言語で今後も発展性が期待できるためこちらに乗り換えたほうが良いということでしょうか。

方程式ベースでモデリングできる点はモデリング時の式変形の負荷が少なく素晴らしいので試してみます。

実装

Modia.jlで捕食者ー非捕食者モデルのシミュレーションをしてみます。

$$\frac{dx}{dt}=x+y$$ $$\frac{dy}{dt}=-x+y$$ $$x………捕食者の個体数$$ $$y………非捕食者の個体数$$ $$初期値を下記のように仮定する。$$ $$x(0)=y(0)=1000$$

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

using Modia
using Plots
@model FirstOrder begin
   x = Variable(start=1000)   # start means x(0)
   y = Variable(start=1000)   # start means y(0)
@equations begin
   der(x) = x + y        # der() means time derivative
   der(y) = -x + y
   end
end;
result = simulate(FirstOrder, 1) # 1秒でシミュレーション停止

t=result["time"]
x=result["x"]
y=result["y"]

plot(t,x,label="x")
plot!(t,y,label="y")

グラフは以下のようになります。 f:id:Otepipi:20190210125211p:plain

本来なら描画にはModiaMath.plotを使うのですが、これにはPyPlotが必要でかつ今回の環境ではPyPlotがうまくインストールできなかったので代用としてPlotsを使用しました。 方程式ベースでシミュレーションできています。プログラムの書き方もModelicaとそっくりに見えますね。

なお、Modelicaで実装したコードは下記になります。

model TEST
  Real x (start = 1000);
  Real y (start = 1000);
 equation
 der (x) = x + y ;
 der(y) = -x + y ;
end TEST ;

さらにこのプログラムをコンパイルして実行するには下記の命令が必要です。

simulate(TEST, stopTime=1)

最後に結果をプロットします。

plot({x,y})

結果は下図のようになります。 f:id:Otepipi:20180221194825p:plain

計算結果は同じですね。

今回はここまでにします。

R Markdownでプロジェクトマネジメント(ガントチャート編

前回に引き続いてR Markdwnを使ってプロジェクト計画書を作成していきます。今回はWBSテーブルからガントチャートの生成を行います。

前回の記事は↓

otepipi.hatenablog.com

ガントチャートライブラリ

ガントチャートを生成するにはライブラリを使用します。今回の案件に使えそうなライブラリはいくつかありますが、ガントチャートのレイアウトが綺麗だと思った「lares」ライブラリを使用します。CRANには登録されていないようです。

github.com

ガントチャートの使い方は↓にあります。この記事にあったガントチャートの例が綺麗だったのでlaresを使うことにしました。

datascienceplus.com

https://datascienceplus.com/wp-content/uploads/2018/11/Captura-de-pantalla-2018-11-09-a-las-2.35.56-p.-m..png

 

今回は使用しませんでしたが、他のガントチャートライブラリには例えば「Plotly」があります。

plot.ly

Plotlyのガントチャートはマウスオーバーでタスクの詳細がポップアップされるのは良いですがレイアウトが寂しいと感じたので使用を見送りました。パラメーター設定で綺麗にできるかもしれませんが……

 

ガントチャートの実装

前回に引き続き↓のwbs.xlsxからガントチャートを作成します。

f:id:Otepipi:20190127175934p:plain

wbs.xlsx

まずはlaresをインストールしなくてはいけませんので下のように宣言します。CRANに登録されていないのでgithubからインストールする必要があります。

install.packages("devtools")
library(devtools)
install_github("laresbernardo/lares")

そしてガントチャート生成のコードは以下のようになります。


library(readxl) ##エクセルデータ読み込みのためのパッケージ
dat_wbs <- read_excel("WBS.xlsx",sheet = 1) #エクセルファイルの読み込み。
library(lares)

plot_timeline(event = dat_wbs$タスク名称, #ガントチャートの作成
              start = dat_wbs$開始日, 
              end = dat_wbs$終了日, 
              label = dat_wbs$担当者, 
              group = dat_wbs$大項目,
              save = FALSE,#チャートを出力するかどうか
              title = "xxプロジェクト",
              subtitle = "ガントチャート")

 

この結果生成されるガントチャートは↓のようになります。タスクの数が少なくスカスカのチャートですが、簡単なガントチャートを作成することができました。

f:id:Otepipi:20190203095825p:plain

生成されたガントチャート

今回はここまでにします。
 

R Markdownでプロジェクトマネジメント(WBS編

今回はプロジェクトマネジメント関連のお話になります。

プロジェクトマネジメントに必須のドキュメントと言えばプロジェクト計画書ですが、本来はプロジェクト期間中常にアップデートされていかなければなりません。

ただ、忙しい日々の業務の中でプロジェクト計画書を更新していくのは並大抵の苦労ではありません。実際は殆ど更新されず、無用の長物になってしまっているのではないでしょうか。

今回は楽にプロジェクト計画書を更新するための試みとしてR Markdownを使ったプロジェクトマネジメントを考えていきたいと思います。

R Markdownとは

R Markdownは基本的にはMarkdownですが、RチャンクによりRのコードを実行できることに特長があります。Jupyter notebookのRバージョンと言った方がわかりやすいかもしれません。このR MarkdownはRStudioを使用してword,html,pdf形式で発行することができますので、他人との共有も簡単です。

これにより各エクセルシートで更新されている

をR Markdownで書かれた計画書に取り込むことで、常に最新の情報にアップデートされたプロジェクト計画書が手に入ると考えています。

(各エクセルシートは常に更新されているという前提です。これはこれでやや無理のあるというか難しい仮定なのですが……)

 

今回R Markdownでプロジェクト計画書を作成するにあたっては下記を参考にしました。

qiita.com

gihyo.jp

mrunadon.github.io

 

RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−

RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−

 

エクセルのWBSを取り込む

ここでWBSをエクセルで管理しているとして、それをR Markdownに取り込んでみます。

例として下のようにWBS.xlsxを作成しています。

f:id:Otepipi:20190127175934p:plain

WBS.xslx

これをR Markdownに取り込みます。プロジェクト計画書に取り込み、htmlで発行した結果下記のようになります。エクセルの内容が取り込まれており、エクセルの中身を変更した場合でも、それが反映されるようになります。

f:id:Otepipi:20190127184605p:plain

R Markdownのコードは下のようになります。Rのコードの解説はコメントアウトされています。

 

---
title: "プロジェクト計画書"
author: "プロジェクトマネージャー"
date: "2019年1月27日"
output: 
  html_document:
    
    df_print: "paged" #テーブルの表示設定
    toc: yes #目次あり。見出しレベル2まで表示
    toc_float: true #ページ左に見出し一覧表示
  word_document:
    toc: yes
---



```{r setup, include=FALSE} 
## 使用するライブラリの読み込み
library(tidyverse) ##データ操作、グラフ化
library(readxl) ##エクセルデータ読み込み
```


## WBS
ここにプロジェクトの作業項目一覧を記す

```{r WBStable} 
dat_wbs <- read_excel("WBS.xlsx",sheet = 1) #エクセルファイルの読み込み。
knitr::kable(dat_wbs) # テーブルの表示
```

 

 また、R Markdwon側で表を加工することができます。例えば進捗にバーチャートをつけることができます。

f:id:Otepipi:20190127191900p:plain

 

 

これは下記のようなコードで表現できます。

library(formattable) ##キレイな表のライブラリ
formattable(dat_wbs,list(進捗 = color_bar("tomato"))) #wordには対応していません

 

これによって進捗状況などの視覚化はR Markdwon側で行うことができます。

これによってggplotをはじめとしたRの強力な視覚化ライブラリを使うことができ、より見やすい計画書になります。

 

今回はここまでにします。

 

 

 

 

遅れて新年のご挨拶

新年を迎えて半月経とうとしていますが、遅れての新年の挨拶になります。

今年もよろしくお願いいたします。

 

お正月休み期間は、今までほとんど手を付けていなかった制御工学について勉強しておりました。

使用したテキストを書いておきたいと思います。

読み終わった本 

はじめての制御工学 改訂第2版 (KS理工学専門書)

はじめての制御工学 改訂第2版 (KS理工学専門書)

 

古典制御について書かれています。自分が読んだはじめての制御工学についての本でしたので書評できるほどではありませんが、一応最後まで読み通すことができました。

 

 

演習で学ぶ基礎制御工学

演習で学ぶ基礎制御工学

 

 「はじめての制御工学」のあとはこの本でひたすら演習しました。ラプラス変換をゴリゴリ手計算したのは学生の時以来で、手がとても疲れました。本当に演習しか無いので、初見の問題の解き方については答えをちら見するか、「はじめての制御工学」に戻ることになります。「はじめての制御工学」の守備範囲ではない問題も多々ありました。

はじめての現代制御理論 (KS理工学専門書)

はじめての現代制御理論 (KS理工学専門書)

 

 現代制御理論についてはこの本で勉強しました。本当に基礎的な内容に絞らているというのは読んでいてなんとなくわかりました。現代制御に触れたのもこの本が初めてで、状態空間表現など見たこともありませんでした。

 

演習で学ぶ現代制御理論

演習で学ぶ現代制御理論

 

 次は案の定このシリーズで演習です。今回ゴリゴリ計算するのはラプラス変換ではなく行列演算です。行列式逆行列の求め方は大学のはじめに習ったきりでほとんど忘れていましたがこの本で思い出すことができました。しかしこの本を手計算で走破するのは非常に骨が折れました。

 

途中の本

Control Systems Engineering

Control Systems Engineering

 

 制御工学の分野も洋書のほうが進んでいるだろうか?と思って購入した本。やはりこちらのほうが内容は深く、また実際の応用例についての言及が多いなという印象です。ただ、(化学工学と比べて)和書とそこまで差はないかなと思います。

 

 

制御のためのMATLAB

制御のためのMATLAB

 

 制御理論をプログラムで実装する際のヒントとして購入しました。ツールボックスの使用を前提としていたのでツールボックスも購入しました。他のMATLAB関係の本だとSimulinkでの実装になりコードでの実装例があまり無いのではと懸念して、この本を購入しました。この本はガッツリMATLABコードでの実装なので非常に参考になります。ただツールボックスの部分でブラックボックス化してしまっているのが残念です。

 

今回はここまでです。もう少し詳しくなれば制御工学についての記事も書いていきたいと思います。

では、本年もよろしくお願いいたします。