hirax.net::inside out::2009年12月

最新記事(inside out)へ  |   年と月を指定して記事を読む(クリック!)

2009年11月 を読む << 2009年12月 を読む >> 2010年1月 を読む

2009-12-01[n年前へ]

NEWS今昔物語「女と男と旅に出る」編 (2004年09月02日) 

5年前のNEWS(未来)を振り返ってみて思うこと

 洋式トイレで男性が小便をする時、一体どういう風に「する」のが一番汚れが少なくてすむのかは、昔から調べてみたいと思っています。「座ってするなら、どの下向き角度で放水を行い・どういう面に・どういう角度で水流をあてるのが良いのか」あるいは「立ってするなら、一体便器のどの部分にどういう角度で当てるのが一番良いのか」を調べてみたいと思っているのです。はたして、ライオン株式会社の研究報告のように、本当に「水たまりの中央を狙うのが一番飛沫が飛び散らない」のかが、どうしても納得ができないのです。他の条件に比べて、どういう風に「飛沫飛散が変化するから、汚れが改善する」という「過程」と「結果」を実際に眺めてみたい、というわけです。そういうわけで、"ToDo"リストに「洋式便所での小便シミュレータ作成」を、追加しておくことにします。

 「美人になる湯」のニュースを読むと、ふと、こんなことを考えます。温泉は好きですが、自分が男性であるせいか、「美人になる湯」よりも、ただ浸かっているだけで、「ハードウェア作りが得意になる湯」とか「語学力が向上する湯」とか「絵を描くテクニックが上達する湯」といった温泉があった方がうれしいように思います。そういう、快適かつ実践的な温泉はないものでしょうか。

 とはいえ、「美人になる湯」ならお湯に浸かっているだけでも効果があることもあるかもしれませんが、私の欲しい「温泉」は「怠け者の夢」としか言いようがないものかもしれません。

(記事を書いた時の)ひとこと

 あなたが男性だとして、洋式便所で小便をするときには、座ってする派だろうか、それとも立ってする派だろうか。「小便の飛び散り」による汚れを気にして、座ってする人も多いことだろう。一体、「座る派 v.s. 座らない派」の比率はどの程度なのだろうか?今回は、「男」と「女」という観点から気になるニュースをいくつか集めてみた。

「男の立ちション」は「的外れ」

 8月3日に開催された日本家政学会の研究発表会で、ライオン株式会社が「洋式トイレのニオイの原因で落としにくい"ふち裏汚れ"」の大きな原因が「洋式トイレで男性が立ちションを行い、便器内の水たまりの"手前"や"奥"を狙っている」であるという研究報告を行った。小便の飛沫が飛び散らないようにするためには、なるべく中央を「狙う」べきらしい。これは、洋式トイレで立って小便をする派の男性にとっては、少し覚えておくと良い豆知識だろう。

 ドイツでは、「座って小便をする奴="Sitzpinkler"」という単語は「弱虫」を意味する。しかし、そのドイツでも「立ちションをすると警告を発する装置」が大いに売れていて(Oliver@「スラッシュドットジャパン」さんの記事参照)、掃除の手間を考えるともはやドイツ男も座りションベンをしなければいけないのかもしれない、というニュースが8月18日に流れた。男らしさの象徴たる立ちションも、「過去の遺物として流されてしまう」のかもしれない。

 ところで、母親が幼児にトイレ・トレーニングをすることが実は多く、「男らしさの象徴たる立ちション」という考えがそもそも「的外れ」だという鋭い指摘もされている…。

「温泉の科学」と「美人の湯」

 今年の夏は、「"温泉の素"を入れた温泉」や「不当表表示温泉」とか「温泉でないといいつつ実は正真正銘の温泉(だけど無許可…)」など、各種パターの「温泉」についてのニュースが巷を騒がせた。そこで、色んな温泉の科学を読んで、温泉に詳しくなってみるのも「エンジニア」的に面白いと思う。また、女性には「お肌が滑らかになるという "美人の湯"の謎 」もきっと興味が惹かれることだろう。

 ちなみに、「美人の湯」度のトップは「奥熊野温泉([体験レポート:http://www2.sen-shu.ne.jp/yatakarasu/mikumano/mikumano_004.htm)」と梅香丘温泉で、いずれも和歌山県である。美人になりたいと思ったら、和歌山へ温泉旅行をするに限る?

「高速道路上を飛行するハト」と「長くて短いカメの家出」

 7月27日のSCIENTIFIC AMERICAN.COMチューリッヒ大学のPeter LippらがGPS機器をハトに負わせて帰巣するようすを確認したところ、旅慣れたハトほど高速道路上空を飛び、しかも長距離になればなるほど高速道路の上を飛びたがる、というという研究結果が掲載された(日本語訳)。それが楽で安全だかららしいが、GPSを載せて高速道路上を移動するだなんて、あんまり賢く無さそうなハトも実は人間とたいして変わらない知恵を持っていることに驚く。

 高速道路上を滑空するハトの記事とは全く逆に、9月1日に英国で記事になった亀の家出の話も面白い。30歳代の亀が飼い主の元から家出したが、長く3ヶ月かかって、それでもわずか5kmしか進まなかったという。この話を「自分の歩み」や「遠くへ行きたいけれど…」といった自分の気持ちに重ねて、ふとしみじみする人も多いかもしれない。

2009-12-02[n年前へ]

Ruby版 Thinkpad 加速度センサ類取得クラスを書きました 

 Lenovo(旧IBM)のノートPC ThinkPad にはハードディスク保護を目的として、加速度センサが付いています。SSDを搭載している機種でも、(その他のハードウェアは共通化されているからだろうと思いますが)やはり加速度センサが付いているということです。

 以前、C++でThinkpad加速度センサ類取得クラスを作りましたが(このC++コードを眺めてみると、不要な部分が残っていたので、その点については後で修正します)、今回はRubyでそのクラスを組んでみました。作成したRubyソースは、ここに置きました(accelerometer.rb)。ソース最下部には、使用方法もコメントアウトした上で付けておきました。また、このバージョンはC++バージョン同様、最近のThinkpadでも動く実装になっています。

 読みやすさのために改行を適当に入れたソースは、下記のようになります。

# jun hirabayashi jun@hirax.net 2009.12.02

class Accelerometer
  require 'Win32API'

  attr_accessor :x
  attr_accessor :y
  attr_accessor :temperature
  attr_accessor :presentState

  def initialize
    begin
      @sensorFunctionAPI=Win32API.new("Sensor",
        "ShockproofGetAccelerometerData", ['P'], 'V')
      @buffer=' '*17
      if @sensorFunctionAPI.Call(@buffer)
        parse(@buffer)
        @offsetX=-@latestRawAccelDataX.to_i;
        @offsetY=-@latestRawAccelDataY.to_i; 
      end
    rescue
      'We have some problem.'
      exit!
    end
  end
  
  def parse(result)
    @presentState,@latestRawAccelDataX,@latestRawAccelDataY,
    @latestAccelDataX,@atestAccelDataY,@temperature,
    @atestZeroG_X,@atestZeroG_Y=result.unpack("iSSSSCSS")
  end 

  def getAccelerometerData
     @sensorFunctionAPI.Call(@buffer)
     parse(@buffer)
     @x=@latestRawAccelDataX.to_i+@offsetX.to_i;
     @y=@latestRawAccelDataY.to_i+@offsetY.to_i;
     return [@x,@y]
  end
end

 このRubyクラスを使った加速度取得スクリプト例は、下記のようになります。ちなみに、Xが(ノートPCに向かって)左右方向で、Yが奥行方向になります。これで、1秒おきにThinkpadの傾斜をコンソール出力します。

require 'accelerometer'

accelerometer=Accelerometer.new
10.times{ |i|
  puts accelerometer.getAccelerometerData.join(',')
  sleep 1
}

 MA5授賞式で面白いものを作っている開発者を見たせいか、何か色々作ってみたくなります。4年ほど前には、この加速度センサテクニックを使って「「未来の立体ディスプレイ」を作る」とか色々作ったような覚えがありますが、まだまだ色々な応用がありそうな気がしてきました。というわけで、まずは準備のために、今日はRuby版 Thinkpad 加速度センサ類取得クラスを書いてみました。

2009-12-03[n年前へ]

NEWS今昔物語「他人と自分」編 (2004年06月00日) 

5年前のNEWS(未来)を振り返ってみて思うこと

 この時、つまり5年以上前のNEWSは「評価」「比較」ということを主題に集めたように思う。NEWSを書きながら、他者との評価・比較や他者からの評価・比較といったことが念頭に浮かんだ。だから、「他人と自分」編というサブタイトルを付けた。

 最近のNEWSのひとつに「事業仕分け」がある。これもひとつの「他者による評価」「他者との比較」である。5年前の記事中にスーパーコンピュータに関する記事もあるが、今現在、次世代スーパーコンピュータも俎上(そじょう)に上がっている。一体、スーパーコンピュータに関する「事業仕分け」では、どのような「評価」「比較」がされていくのだろうか。

 ところで、「自分への視線を発見・記録するメガネ」だが、自分への評価を眺めることができて、さらに、「自分を対象とした事業仕分け」まで解析・表示することができる「新カメラ」すら、いつの日か生まれるかもしれない。まるで、SFのようだが、これまでの科学技術の歴史はSF本の中のものをたくさん具現化してきた。もしも、そんな「個人対象の事業仕分けカメラ」が生み出されてしまったら…考えるだけで恐ろしい。

(記事を書いた時の)ひとこと

 「他人」と「自分」という「人」と、「機械」や「機械ができること」についてのNEWSを集めてみました。

「センター試験」と「2次試験」の数学の得点は無関係?

 東北大学の理学部志願者の「センター試験」と「二次試験」の各教科成績を調べ、外国語などはセンター得点が高ければ2次得点も高いという相関が表れたが(例:基礎統計学から見た入試)、数学では相関が極めて低かった、との 森田康夫氏による分析結果が5月30日に報告された。 その点数のズレの原因は、大学が重視(評価)する二次試験は「論理的思考能力」を重視し、センター試験は「計算能力」を重視(評価)しているためであるという。

 最近は、エンジニアに対する評価をどのように行うか、という論議もよく見かける。 しかし、誰もが納得するような評価を行うということは極めて難しい問題である。少なくとも大学入試には解ける問題しか出題されないが、この評価制度の場合は解くことが未来永劫できない難問だったりするのかもしれない。

自分への視線を発見・記録するメガネ

 「赤外線を投光し、角膜からの反射光が瞳の中心に位置するかどうかを解析することにより、自分への視線を発見・記録する」というメガネの研究が5月4日に視線入力に関するシンポジウムで報告された。このメガネをかければ、自分が何を眺めたかだけでなく、誰が自分を見つめたかまで記録される(MPEGビデオ)。逆に言えば、自分の視線が何処を向いているかを記録されてしまうのだから、サングラスをかけて自分の視線を隠す人が増えるかもしれない。

 そういえば、携帯電話で撮影した顔写真を送信すると、その「顔」がどの芸能人に似ているか評価・採点してくれるサービスなどもある。同じように、他人が自分を眺めながら「どんな評価・採点をしているか」まで視線から知らされてしまったら…? そんな「評価スカウター」があったら、人の視線・評価が気になってたまらなくなってしまうに違いない。

高性能計算用途のWindowsと2テラFLOPSのGRAPE-DR

 先日、Microsoftが「高性能コンピューティング向けのWindows」の計画を開始した。 並列計算をWindowsで行いやすくする、という計画である。とはいえ、私たちの手元にあるようなクライアント用のWindowsではなく、あくまでサーバー用のWindowsの話だ。ExcelからグリッドPCを用いた演算機能を利用できるようにするという製品もあるが、残念ながらそういった環境はまだ一般的にはなっていない。一般的でもお手軽でもないからこそ、「スーパー」コンピュータなのである。

 当初は「20万円で作ったスーパーコンピュータ」という売り文句で有名になったGRAPEだって進化を続け、先日は2ペタフロップスを実現しようとするGRAPE-DR計画発表された。汎用計算ができるわけではなく、あくまで専用計算機であるから、比較するわけにはいかないが、2テラFLOPSというのはスゴイ世界だ。何しろ、一秒間に2千兆(2000000000000000)回もの浮動小数点演算を行ってしまうのだから。

米国防省の検閲語句がフォント解析で解読される

 米国防省のイラクに関する文書から「検閲により消されてしまった文字」を、消された単語の幅を推定し、フォント(字体)の種類を考慮して(その単語の幅になる)文字の組み合わせを調べだし、妥当な単語・名称をしらみつぶしに検索することで、消された文字を復元したという5月13日の Natureで報告された。コンピュータ上の電子辞書と文章解析ソフトだけで、米国防省の検閲を突破できたわけである。

 ところで、こんな風に機密保持のために塗り潰した部分が解読されてしまったのも、そもそも紙資料の画像で機密文章を公開したからである。いつまで紙文化が続いていくのだろうか。

2009-12-04[n年前へ]

Thinkpad 加速度センサでExcel3次元グラフを未来のディスプレイ風にしよう 

 Ruby版 Thinkpad 加速度センサ類取得クラスを作ったので、さっそく何か加速度センサを使って遊んでみることにしましょう。まずは、とても簡単に「未来のディスプレイ」を作る風のことを、Ruby+Microsoft Excelで行ってみることにしましょう。

 というわけで、Thinkpadを傾けると、それに応じてエクセルの三次元グラフを(傾けた方向に応じた)さまざまな方向から眺めることができるという次のようなスクリプトを書いてみました。

 このスクリプトはまずエクセルで(適当な値を代入した)三次元グラフを作成し、(その後の30秒間にあなたがグラフを色々カスタマイズする時間を与えた上で)、30秒経た後にThinkpadの傾きに応じて三次元グラフを好きな方向からリアルタイムにグリグリ動かしながら眺めることができる、というスクリプトです。(動画は「Thinkpad 加速度センサ+RubyによるExcel3次元グラフ動画」で観ることができます)

require 'accelerometer'
require 'win32ole'

excel=WIN32OLE.new("excel.application")
excel['Visible']=TRUE
excel.Workbooks.Add()

excel.Range("a1")['Value']=1
excel.Range("a2")['Value']=1
excel.Range("a3")['Value']=1
excel.Range("a4")['Value']=1
excel.Range("b1")['Value']=1
excel.Range("b2")['Value']=2
excel.Range("b3")['Value']=2
excel.Range("b4")['Value']=1
excel.Range("c1")['Value']=1
excel.Range("c2")['Value']=3
excel.Range("c3")['Value']=2
excel.Range("c4")['Value']=1
excel.Range("d1")['Value']=1
excel.Range("d2")['Value']=1
excel.Range("d3")['Value']=1
excel.Range("d4")['Value']=1
excel.Range("a1:d4").Select()
excelchart=excel.Charts.Add()
excelchart['Type']=-4100
sleep 30 # wait time for custumize
accelerometer=Accelerometer.new
300.to_i.times do
  accelerometer.getAccelerometerData
  x=accelerometer.x
  y=accelerometer.y
  r=Math.sqrt(x*x+y*y)
  elevation=90-r.to_i
  rot=360-(Math.atan2(y,r)+Math::PI)/(2*Math::PI)*360
  excelchart.rotation=rot  # 0 to 360
  excelchart.elevation=elevation # -90 to 90
  sleep(0.2)
end
excel.ActiveWorkbook.Close(0)
excel.Quit()
 このスクリプトを動かせば、Thinkpadを傾ければ、低い角度からグラフを眺めたり、左右に倒せば左右の異なる方向からグラフを眺めたりする…ということができます。「未来風」かはさておき、現実世界で行う動作に適切に対応した動作、つまりは「実感・体感できる」ソフトウェアを作る、というのは何だかとても面白いものです。

 今回作ったスクリプトを動かしているようす、Thinkpad+Ruby+Microsoftエクセルでグラフを未来風にグリグリ動かしている様子を撮影した動画は、「Thinkpad 加速度センサ+RubyによるExcel3次元グラフ動画」として、公開しました。

Thinkpad 加速度センサでExcel3次元グラフを未来のディスプレイ風にしよう






2009-12-05[n年前へ]

Thinkpad 加速度センサ+RubyによるExcel3次元グラフ動画 

 「Thinkpad 加速度センサでExcel3次元グラフを未来のディスプレイ風にしよう」を操作しているようすを、ケータイのカメラで撮影してみました。それが、下の動画になります。Ruby版 Thinkpad 加速度センサ類取得クラスを使い、Excel 2007のグラフを視点を操作している、という具合の動画です。パースペクティブを適切につけておくと、上下方向の視線移動に関してはかなり自然に眺めることができます。(この動画実行のためのRubyソースは上記記事を参照してください)

 Excelの三次元グラフが、ElevationとRotationという、天頂近くで精度が低下しやすい座標指定方式なのと、Thinkpadの加速度センサが2軸方式で、鉛直線中心の回転を取得することができなかったりするところが(東芝のPCなどでは3軸取得が可能だったりするのですが)、「少しの変さ」をかもしだしていたりしますが、そこはそれ、Microsoft ExcelのグラフをノートPCの方向を変えると、それに応じた色々な方向から眺めることができる、というのは何だか面白い、とは思いませんか?

Thinkpad 加速度センサ+RubyによるExcel3次元グラフ動画






2009-12-06[n年前へ]

Thinkpad加速度センサ取得用C++クラスの手直しをしました 

 Lenovo(IBM) Thinkpad加速度センサ取得用C++クラス(関連記事・Thinkpad加速度センサ取得用C++クラス新しいThinkpad にも対応した加速度センサ値取得プログラム)を少し手直ししておきました。動作は全く変わりませんが前回の修正の際に不要な部分が残っていたので、その点について直しました。

Gravity ball  Thinkpad加速度センサ取得用C++クラスをまとめたヘッダファイルソース(Sensordll.h)、および、使用サンプルソース・バイナリ(sample.cpp・sample.exe)は、ここに置いておきました(古いバージョンは、サブディレクトリに置いてあります)。

サンプル・アプリケーション例では、よく意味がわらないままに、"Temprature"も出力するようにしてあります。

sample.exe 1000 10
という風にコマンドラインからアプリケーションを実行すると、
0=x, 0=y,35=temp.
0=x, 0=y,35=temp.
0=x, 0=y,35=temp.
0=x, 0=y,35=temp.
1=x, 1=y,35=temp.
-16=x, -3=y,35=temp.
-1=x, -1=y,35=temp.
とった具合に左右方向の傾斜と奥行き方向の傾斜(とtemprature)を出力します。

2009-12-07[n年前へ]

エクセルの計算ワークシートをRuby計算スクリプトに変換してみよう 

 下記で使ったエクセルファイルは、自己参照部分がなく、つまり反復計算が必要ないものを使ってしまいました。反復収束計算が必要なエクセルファイルを「続 エクセルの計算ワークシートをRubyでC言語に変換してみよう」 でサーバ内に置きましたので、下記記事を読んだ後には、上記記事(さらにその後に続く記事など)を引き続きご覧ください。


 目的に対して達成方法があまりに過剰、あるいは、あまりに的外れで荒唐無稽なものというのは、それはそれで何だか面白いような気がします。

 今日は、そんなものを、ふと、けれど無性に作ってみたくなりました。

 先日、Microsoft Excelのような「表計算ソフト」スプレッドシートを使った「静電界計算や非定常熱伝導シミュレーション計算」をみっちりしてきました。離散化された物理式を表計算ソフトウェアを用いて計算を行うというのは、「セル」というメッシュと空間が感覚的に近く感じられることから、とても自然にシミュレーション計算を行う空間分割と計算領域との対応を感じることができるのが、とても便利で良いと思っています。

 もちろん、セル間の計算を「反復計算」を行うことで、複雑なことを気にせず、計算が自然に終わってしまう、というのも実に「自然」で「お手軽」で良い、と思っています。

 とはいえ、私はエクセルは苦手です。毎年何回か、エクセルでのシミュレーション計算を(悩んでいる)人に教えるという作業をするのですが、それでもエクセルを使うのは苦手です。そこで、エクセルで作った離散化シミュレーション用.xlsシートをRuby Script(プログラム)に自動的に変換する、Rubyスクリプトを作ってみることにしました。そう、「エクセルで物理シミュレーションを学んだ人が、そのエクセルシートをそのまま多言語に移植することができる」スクリプトを作ってみよう、と考えてみたのです。

 作ったコードxls2rb.rbは次のようになります(近日中にサーバに置いて置きます)。これは、「反復計算を使った」エクセルファイルを読み込み、それをRubyスクリプトに変換するスクリプトです。「表計算ソフト」スプレッドシートを使いった「静電界計算や非定常熱伝導シミュレーション計算」は、四則演算と周囲の空間が持つ値の計算だけで行うことができます。だから、そういう「条件」下では、そんなスクリプトを作るのも比較的簡単です。初期化や反復計算をそれぞれ「関数」化しているとはいえ、変数はすべてグローバル変数ですし、そこら辺の「エクセル感」も再現してみたスクリプトです。

require 'win32ole'

def getAbsolutePath(filename)
  fso=WIN32OLE.new('Scripting.FileSystemObject')
  fso.GetAbsolutePathName(filename)
end

def getAlphabet(n)
  val=n.to_f/25.0 
  mod=n%25
  result=''
  result+=('a'..'z').to_a[(n-val)/25] if val>1
  result+=('a'..'z').to_a[mod]
end

def defFunc(state,book)
  src=''
  book.Worksheets.each do |sheet|
    y=1
    sheet.UsedRange.Rows.each do |row|
      x=0
      record=[]
      row.Columns.each do |cell|
        if state==:init &&cell.Value!=''
          record<<'  '+'$'+getAlphabet(x)+y.to_s+
            '='+cell.Value.to_s if cell.Value.to_s!=''
        end
        if state==:calc &&cell.Formula!=''
          t=cell.Formula.sub(/[=$]/,'')
          t=t.downcase.gsub(/([a-zA-Z]+\d+)/,'$\1')
          record<<'  '+'$'+getAlphabet(x)+y.to_s+'='+t
        end
        x+=1
      end
      src+=record.join("\n")+"\n" if record.join('').gsub("\n",'')!=''
      y+=1
    end
  end
  return src
end

filename=getAbsolutePath(ARGV[0])
excel=WIN32OLE.new('Excel.Application')
book=excel.Workbooks.Open(filename)
initsrc=defFunc(:init,book)
calcsrc=defFunc(:calc,book)
book.close
excel.quit
GC.start

puts <<INIT
# autocreated ruby script from excel file
# jun hirabayashi jun@irax.net http://www.hirax.net

def init
#{initsrc}end

def calc
#{calcsrc}end

10.times do
  calc
end
INIT

 ちなみに、このスクリプトを

ruby xls2rb.rb ex.xls
という具合にして(引数は変換したいXLSファイルです)、実行すると、次のようなスクリプトができあがります。ちなみに、ここで使ったex.xlsファイルは、ラプラス方程式を3X3のセルに離散化し解く、実に単純なエクセルシートです。(この辺りのファイル一式を週末にでもサーバ上に置いておきます)
# autocreated ruby script from excel file
# jun hirabayashi jun@irax.net http://www.hirax.net

def init
  $a1=1.0
  $b1=2.0
  $c1=3.0
  $a2=1.0
  $b2=2.0
  $c2=5.0
  $a3=1.0
  $b3=2.0
  $c3=3.0
end

def calc
  $a1=1
  $b1=2
  $c1=3
  $a2=1
  $b2=2
  $c2=$b2+$c1
  $a3=1
  $b3=2
  $c3=3
end

10.times do
  calc
end
 つまり、今日私が作ったスクリプトの動作は、エクセルが最初に行う初期化ルーチンを"init"関数として定義(作成)し、次に行う反復計算を"calc"関数として定義(作成)し、それを(適当に決めた)10回繰り返すスクリプトを作り出す、という具合です。ポイントは初期化時には、cell.Valueを用いることで初期数値を設定し、逐次計算時にはcell.Formulaを使うことで数式を使う、という「使い分け」になります。反復計算自体は、ガウスザイデル法によって行われます。

 あとは適当に、知りたい値を出力する(たとえば、"puts $c2"といった)出力処理文でも書き足せば、はい、シミュレーション・rubyプログラムのできあがり、というわけです。

離散化された物理計算をするために作成したエクセル・シートを、Rubyスクリプトに自動で変換し、好きに加工できるプログラム…って便利なような、そもそもそんな用途ってある訳ないような…というのが正直で的確な感想だと思います。つまり、それはかなり無意味なツールです。

 他の言語に置き換えて、高速化指向の変換プログラムを書いてみるのも(一回くらいはやってみようと思いますが、今ひとつ「魅力」に欠けるように思います。

 今日作った、「離散化された物理現象を計算するために作成したエクセル・シートを、Rubyスクリプトに自動で変換するプログラム」なんていうものを、(何だか馬鹿馬鹿しいけれど)意外に少し面白い笑える、そう思える人が一人でもいたら、幸いに思います。

2009-12-08[n年前へ]

Excel 2007で科学技術計算風の等高線グラフを作る 

 表計算ソフト(スプレッドシート)でシミュレーション計算をした後には、必ず計算結果を図示したくなります。「離散化された空間がセルで分割された表に似ていること」そして、表上で計算された結果をグラフ化することが簡単なこと、がMicrosoft Excelのような表計算ソフト(スプレッドシート)の好ましい点なわけですから、必ずと言って良いくらい「表計算ソフト(スプレッドシート)でシミュレーション計算」をした後には、グラフ作成の作業がグリコのオマケのように、必ずと言って良いほど、ついてきます。

 しかし、残念なことに、Excel 2003以前でも、自分でグラフをカスタマイズしなければ、(たとえば)Excel2007で科学技術計算風の等高線グラフを作ることはできません。Excel 2007になっても、実に残念なことに、いまだに同じなのです。一見陰影がついた綺麗なグラフに見えたとしても、それはまるで、何だか昔ながらのエクセルのグラフをただシュールにリアルな筆致で描いたかのようなグラフになってしまいます。

 そこで、今回はExcel 2007用の科学技術計算風の等高線グラフ・テンプレート(右上のチャートがその例になります)を作ってみることにしました。…といっても、私が行ったのは(人から頂いた)カスタマイズされた綺麗なグラフが張り付けられているたExcel 2000形式の.xlsファイルをExcel 2007で開き、その後、新しいエクセルブックを作成した後にグラフも含め全コピーし、グラフをさらに適当にカスタマイズした上で、デザイン-テンプレートとして保存するということだけです。つまりは、ただフォーマット変換作業と、手順解説記述、ということになります。

 まずはこのファイル(rainbow3DContor.crtx)を、C:\Users\Administrator\AppData\Roaming\Microsoft\Templates\Charts といった自分の(ユーザの)Excel用チャートテンプレートを保持するディレクトリに置きます。すると、
挿入-その他のグラフ-すべての種類のグラフ-テンプレート-作ったテンプレート
という手順で、そのカスタマイズされたグラフを一瞬で作ることができるようになります。

 もちろん、一回グラフにした後に、
デザイン-グラフの種類変更-テンプレート-作ったテンプレート
で変更する、という手順でも構いません。どちらの手順でも構いませんが、これで自然な虹色の縦軸(値軸)が0~100の等高線鳥瞰図を作ることができるようになります。

 私自身はExcelは、今でも「時間泥棒」だと思っています。一見、地道に作業を続けているように見えても、それは単に時間を非効率的に売っている作業をしているように思えてしまいます。だから、あまり、Excelが好きというわけではありません。

 とはいえ、エクセルを使わざるを得ない状況というものもやはりあるわけで、そんな場合のために、今日はExcel 2007で科学技術計算風の等高線グラフを作るためのテンプレートファイルを作ってみました。

Excel 2007で科学技術計算風の等高線グラフを作るExcel 2007で科学技術計算風の等高線グラフを作る






2009-12-09[n年前へ]

エクセルの計算ワークシートをRubyでC言語に変換してみよう 

 この記事のスクリプトや、(下記のエクセルファイルとは違いますが、より有用そうな)エクセルファイル例、そして、そのエクセルファイルを変換したC言語ソースを保存し、また、つらつら考えことなどを「続 エクセルの計算ワークシートをRubyでC言語に変換してみよう」 に書きましたので、下記記事を読んだ後には、上記記事(さらにその後に続く記事など)を引き続きご覧ください。


 先日、「エクセルの計算ワークシートをRuby計算スクリプトに変換してみよう」でエクセルで作った離散化シミュレーション用.xlsシートをRuby Script(プログラム)に自動的に変換する、Rubyスクリプトを作ってみました。基本的には四則演算で(反復収束計算手法を用いることにより)「この世界を(エクセル上で)シミュレーションするスプレッドシート」を作ったなら、それを他のプログラミング言語に「変換する」アプリケーション例として、エクセルのシートをRuby言語に変換するソフトを作ってみたわけです。目標は「何だかすげー!、けど、役に立たねー!」です。

 というわけで、引き続き、今日はエクセルのシートをC言語に変換するRubyスクリプトを書いてみました。もちろん、今回も、基本的には四則演算と数値だけで作られ、反復収束計算手法を用いることで最終的な計算結果を得るような「シート」を前提にしています。

 というわけで、書いたコード"xls2c.rb"はこんな具合です。

require 'win32ole'

def getAbsolutePath(filename)
  fso=WIN32OLE.new('Scripting.FileSystemObject')
  fso.GetAbsolutePathName(filename)
end

def getAlphabet(n)
  val=n.to_f/25.0 
  mod=n%25
  result=''
  result+=('A'..'Z').to_a[(n-val)/25] if val>1
  result+=('A'..'Z').to_a[mod]
end

def defFunc(state,book)
  src=''
  book.Worksheets.each do |sheet|
    y=1
    sheet.UsedRange.Rows.each do |row|
      x=0
      record=[]
      row.Columns.each do |cell|
        if state==:def &&cell.Value.to_s!=''
          record<<'  float '+getAlphabet(x)+y.to_s+'=0.;'
        end
        if state==:init &&cell.Value!=''
record<<'  '+getAlphabet(x)+y.to_s+
            '='+cell.Value.to_s+';' if cell.Value.to_s!=''
        end
        if state==:calc &&cell.Formula!=''
t=cell.Formula.sub(/[=$]/,'')
record<<'  '+getAlphabet(x)+y.to_s+'='+t+';'
        end
        x+=1
      end
      if record.join('').gsub("\n",'')!=''
        src+=record.join("\n")+"\n" 
      end
      y+=1
    end
  end
  return src
end

filename=getAbsolutePath(ARGV[0])
excel=WIN32OLE.new('Excel.Application')
book=excel.Workbooks.Open(filename)
defsrc=defFunc(:def,book)
initsrc=defFunc(:init,book)
calcsrc=defFunc(:calc,book)
book.close
excel.quit
GC.start

puts <<INIT
// autocreated C source from excel file
// jun hirabayashi jun@irax.net http://www.hirax.net

#include "stdio.h"

// gloval variables
#{defsrc}
void init(void){
#{initsrc}}

void calc(void){
#{calcsrc}}

int main(){
  int i;
  init();
  for(i=0;i<10;i++){ calc(); }
  printf("C2=%f",C2);
  return 0;
}
INIT
 このスクリプトを、
ruby xls2c.rb ex.xls > sample.c
という具合に実行すると、下のようなC言語ソースができあがります。
// autocreated C source from excel file
// jun hirabayashi jun@irax.net http://www.hirax.net

#include "stdio.h"

// gloval variables
  float A1=0.;
  float B1=0.;
  float C1=0.;
  float A2=0.;
  float B2=0.;
  float C2=0.;
  float A3=0.;
  float B3=0.;
  float C3=0.;

void init(void){
  A1=1.0;
  B1=2.0;
  C1=3.0;
  A2=1.0;
  B2=2.0;
  C2=5.0;
  A3=1.0;
  B3=2.0;
  C3=3.0;
}

void calc(void){
  A1=1;
  B1=2;
  C1=3;
  A2=1;
  B2=2;
  C2=B2+C1;
  A3=1;
  B3=2;
  C3=3;
}

int main(){
  int i;
  init();
  for(i=0;i<10;i++){ calc(); }
  printf("C2=%f",C2);
  return 0;
}
 もちろん、今回もスプレッドシートの「セル」はすべてグローバル変数として取り扱い、エクセルが最初に行う初期化ルーチンを"init"関数として定義(作成)し、次に行う反復計算を"calc"関数として定義(作成)し、それを(適当に決めた)for文で10回繰り返す(計算を収束させるためには、実際にはもっと多く繰り返し計算をさせることが必要でしょう)、というソースです。

 さて、このC言語ソースは、(たとえば、Borland C++ Compilerを使うなら)

bcc32 sample.c
という具合で、sample.exeというバイナリができあがります。

 いつか、エクセルのシートを変換し・作成したC言語プログラムをコンパイルすると、どれだけ遅くなるのか・どれだけ早くなるのかを、色々な環境で確かめてみたい、と思っています。

2009-12-10[n年前へ]

エクセルでシミュレーション Vol.8 二次元非定常熱伝導問題を解こう 

 表計算ソフト(スプレッドシート)プログラムであるMIcrosoft Excelを使い、(VBAはもちろん関数も使わず)ただセル間の加減乗除演算のみで、二次元非定常熱伝導問題を解きで遊んでみることにしてみました。伝熱方程式を差分化し、二次元の非定常解析を行ってみる、というものです。二次元といっても、比較的薄い円形ベルト(シート)上のものが回転するというモデルで、ベルト(シート)の厚みが(ベルト内の熱拡散に対して)十分に薄いという仮定の下での近似を行った、と考えれば、(そういう条件下の)三次元非定常熱伝導問題を解いているとも言えるシミュレーション計算モデルです。

 計算領域は右のようなものになり、2本のローラにより回転駆動させている「下から5点のヒータで温められている茶色のベルト」という計算するモデルです。その概略を大雑把に描いてみたのが、上図になります。

 というわけで、下の動画が、「ベルトが回転し動き」「途中にあるヒーターでベルトを温め」「ベルトの平面内では熱伝導が生じ」「周りの空気とベルトの境界から熱が放熱され」・・・といった2.2秒の間に生じることをすべて計算し・結果刻々と表示するシートで非定常伝熱計算を行っているようすです。(熱伝導率等はテキトーな値です)

ちなみに、このシートは、毎年行われている「スプレッドシートを用いた電子写真シミュレーション実習」テキストの「有限差分法による定着プロセスの熱伝導計算(非定常問題)」で使う教材用シート(関連テキスト )にもとづいています。

 動画の上の部分が、10×20の領域に差分化されたメッシュの温度分布を数値で示しており、動画の下部部分には温度分布をグラフで示しています。このグラフを眺めていれば、20℃の室温下で5点のヒータにより温度を上げられていく(回転する)ベルトが、次第に温度分布を持ちながらやがて定常的な温度分布になっていくさまを見てとることができます。

 この計算シートは、熱伝導方程式

をCFL条件(Courant-Friedrichs-Lewy Condition)下で陽的に差分化し、エクセルの反復計算を用いることで、解いています。しかも、メッシュ間の計算は、ヤコビ法を用いた計算を行っています。

 「メッシュ間の計算は、ヤコビ法を用いた計算を行っている」と、ひとことで書きましたが、これは実は凄いことです。なぜかというと、エクセルは(セル間の循環参照がある場合に用いる)反復(収束)計算時には、表の左上をスタート地点として「Zの法則 」の順序にしたがって(つまり、左上→右、そして一個下の行をまた左→右という順番で)逐次的に解くということをします。つまり、何も考えずにただ自分のセルの周囲との演算を行うと、「あるセルの計算を行うときには、自分の上と左のセルはすでに値が更新され新しい値が入っているが、右と下のセルはまだ古い値が入っている」という状態になります。つまり、ガウス・ザイデル法で計算されることになります(ヤコビ法は全セルの古い値を用いて計算を行った後に、全セルの値を新しい計算結果で置き換えますが、ガウス・ザイデル法は各メッシュ(セル・格子点)の計算を行ったら、その計算結果を、次のメッシュの計算時に即時に反映させて計算を行います)。

 定常問題を解くならガウスザイデル法で収束するまで反復計算を行うので構いませんが、非定常問題を解くのにガウスザイデル法を用いてしまうと、正確ではない間違った温度傾斜が生じてしまいます。それは、少し困ります(しかし、現在出版されている”表計算+差分化 熱伝導”書籍は、多くがガウスザイデル法を用いています)。

 というわけで、この計算シートではグラフの下に(動画の)上のメッシュ部分をコピーするだけのセルがあり、上のメッシュ部分は「自分の周囲に相当する下のメッシュを用いてメッシュ間の計算を行う」という仕組みになっています。こうするとなぜ、ガウスザイデル法でなく、ヤコビ法の計算になるかは、「エクセルの反復計算時のセル間の計算順序」を念頭に置いて考えてみると、「なるほど!」と納得できると思います。「道具を理解し使いこなす」ということの意味を、心から実感させられます。

 だから、この二次元非定常熱伝導問題を解くエクセル・シートは、とても「単純明快な計算がされているだけ」なのですが、実に巧妙に考え抜かれた上で作られた「素晴らしき一品」です。しかも、セルを少しづつ書き変えていくことで、かなり実用的な計算もできるのです。このシートを見たときには、本当に驚かされました。

 「どんな風に書きかえると何がわかるか」とか「さらにこんなこともできる」といったことは、明日にでも書いてみようと思います。(というわけで、この続きとして、温度が時間を追って上昇していく機能も付けてみました

エクセルでシミュレーション 二次元非定常熱伝導問題を解いてみようエクセルでシミュレーション Vol.8 二次元非定常熱伝導問題を解こう






2009-12-11[n年前へ]

エクセルでシミュレーション Vol.9 二次元非定常熱伝導問題の温度変化グラフも作ってみよう 

 この記事は、「エクセルでシミュレーション Vol.8 二次元非定常熱伝導問題を解こう」の続きです。前回の記事では、比較的薄い円形ベルトが回転するモデルで、ベルトの厚みがベルト内の熱拡散に対して十分に薄いという仮定の下での近似を行ったと考えると)三次元非定常熱伝導問題を解いているとも言えなくない・・・というような二次元非定常熱伝導問題を解く、エクセルの表計算シートを作り眺めてみました。もちろん、VBAはもちろん関数も使わず、ただセル間の加減乗除演算のみで、熱伝導方程式を差分化したものです。

 今日は、計算領域の任意の部分の温度変化をグラフにする部分を前回のエクセルシートに追加してみようと思います。つまり、ベルト状に温度センサでも貼り付けたら、時々刻々と一体どんな温度が測定されるのかをシミュレーション計算時に表示させたい、というわけです。

 このような需要は非常に多くあると思いますが、一体それをどのようにすれば作ることができるでしょうか?エクセルのシート上は、非定常問題を解いているので、その瞬間その瞬間の温度分布しかデータが残っていないように思われ、スタート時からの温度変化をグラフにするのは難しいように思えてしまいます。

 しかし、「ガウスザイデル法でなく、ヤコビ法を用いて反復計算を行うために使ったテクニック」、すなわち、

 エクセルは(セル間の循環参照がある場合に用いる)反復(収束)計算時には、表の左上をスタート地点として「Zの法則 」の順序にしたがって(つまり、左上→右、そして一個下の行をまた左→右という順番で)逐次的に解くという計算順序になります。
ということを利用すると、「計算領域の任意の部分の温度変化をグラフにする」ということが、いとも簡単に実現できてしまうのです。種明かしは後回しにして、まずはそのようなことをしているエクセルシート動画を下に張り付けてみます。動画の上部分は前回と同じく「回転するベルトを平面上の等高線色付きグラフとして表示したものです。そして、その下にあるのが、「等高線色付きグラフ右横部にある灰色円部分」の温度変化を横軸:時間・縦軸:温度で時系列グラフにしたものです。

 下の折れ線グラフを見れば、「等高線色付きグラフ右横部にある灰色円部分」の温度が時間を追って上昇していくようすをよく理解することができるのではないでしょうか。

 種明かしはこうです。エクセルのシートのずっと下方(下の行のセルに)"=温度変化を知りたいセル"という式を入れます。ここでいう温度変化を知りたいセルというのは、たとえば、X24といったセルを示す番号です。ですから、実際に(一番下のセルに)入れるのは、"=X24"といったものになります。

 そして、そのセルの上にあるセルをクリックし、=と入力した上で、その下のセルをクリックし"Enter"を押します。つまり、一個下の(さきほど"=X24"という式を入力した)セルを参照するようにするのです。そして、そのセルを選択した上で、上の方までズルズルズル~とコピーしてしまいます。

 エクセルの(セル間の循環参照がある場合に用いる)反復(収束)計算時の計算順序を考えてみれば、今「ズルズルズル~とコピー」を行った列には、「古い計算結果が上、一番下が最新の計算結果」という順序のデータが格納されていくことになります。(上から下に計算が逐次行われていることを考えると)計算ステップごとに、一個値が上に移動していくので、結果としてそのように時系列的なデータを保持できる、というわけです(結局のところ、ヤコビ法を用いるために、バッファエリアを下に設けたのと全く同じパターンです)。

 あとは、計算スタート時点からの(シミュレーション上の)経過時間なども同じように作り、散布折れ線図でも挿入すると、「上記のような「任意の部分の温度モニター機能付きの疑似三次元非定常熱伝導問題を解くエクセルシートのできあがり」となるわけです。もちろん、温度モニターは(上段で行ったことを他のセルに対してもしてやるだけで)いくつでも設置することができます。実に単純・簡単な(けれど巧妙な)実装ですがとても便利な機能です。

 さて、ベルトの温度モニター機能も付きましたので、今度はセンシングしたベルトの温度を用いて、ヒーターの温度を(まずは簡単なPID制御でも使って)制御しベルトの温度分布を適正に調節する機能例でも、簡単に実装し(ハードウェアをいじっている)気分にでも浸ってみることにしましょうか。

2009-12-12[n年前へ]

ケツメイシ「トレイン」 

 ホームに入ってくる電車を待つ人たちの中に混じり、電車に乗り込む。空いている席に座り、向かいの窓の向こうに見える夜の街並みを眺めながら、ノートPCを取り出してメールの返事をいくつか書いているうちに、気づくといつの間にか寝ている。

ただ悠々と走ってく その積み荷の重さ
誰が知ってるだろう でも構わず
音立て 飛び立て 夢の空へ。
 毎日走り続ける列車には、たくさんの人が乗っていて、そんな人たちを眺め、その人たちの毎日は一体どんなものなのだろうか、と考えてみる。そんな時、ケツメイシ「トレイン」(下の動画)が頭に流れることがある。そして、そんなBGMを心の中でパワープレイしながら、色々なことを考えているうちに、いつの間にかいつも寝てしまう。

2009-12-13[n年前へ]

各種言語からのThinkpadの加速度センサ値取得方法 

 ThinkpadなどノートPCでは、ハードディスクの衝撃回避のために、加速度センサを搭載していたりします。私を含め、Thinkpadの加速度センサ値を取得し、色々遊んでいる人は多いので、今回は「各種言語からのThinkpadの加速度センサ値取得方法」へのリンクを作ってみました。私の知っているものを並べただけですので、「おいおい、この言語版もあるぜ」という情報などをjun@hirax.netまで頂ければ、リストを更新しておこうと思います。



2009-12-14[n年前へ]

続 エクセルの計算ワークシートをRubyでC言語に変換してみよう 

 前回、「エクセルの計算ワークシートをRubyでC言語に変換してみよう」でサンプルに作ったエクセルのシートが反復収束計算をする必要がない(自己参照型のシートになっていなかったこと)に気づいたので、今日はエクセル・シートを作り直し、サーバ内に置いておきましたここに置いておきます。置いたエクセル・ファイル"ex2.xls"は、5×5のセルの周囲と中央のセルが固定条件で、それ以外の(5×5内の)セルが自分の上下左右のセルの平均の値をとるというラプラス方程式を差分化したエクセルシートです。そして、xls2c.rbがエクセル・ファイルをC言語ソースに変換するRubyスクリプトで、ex2.xlsをC言語に変換した結果ファイルが、ex2.cになります。

 元々の記事は、「とても無意味なこと・けれど技術発想的には何だか(的外れ具合が)面白いこと」をしよう、という考えのもとにやってみたのですが、非定常熱伝導計算の復習をしているうちに、非定常問題を解くシートに対しては、意外に使えるかもしれない、という気がしてきました。

 プログラムというものをまったく知らない人が、熱伝導方程式の物理的な意味を理解上で、表計算でシミュレーションをしてみた(習ってみた)あとに、(エクセルのセルと直感的に対応した)自動出力したC言語ソースをもとに「printfの使い方だけを習い」色んな計算をしてみれば、そして注目するセル=箇所の時間変化・経過を出力させてみたりすれば、案外面白く感じるかもしれない・・・という風にも思ったのです。

 もちろん、こんな単純自動変換Cソースでは一瞬のうちに物足りなくなるはずなので、上のような時間を経ることで、プログラム言語を学ぼうとするまでの小さな「足がかり」「キッカケ」になるのではないか、という間違っていそうなそんなことをふと思ったのです。いわば、エクセルを入口に差分化シミュレーション・プログラムの世界に"Hello World"する、というわけです。ゼロから学ぶのと、何か叩き台があって、それをほんの少しづついじりながら学ぶのとでは、気楽さが全然違うのではないか、といったことを「後付けで」考えてみたりしたのです。

 ちなみに、エクセルシート→C言語 変換などで、配列を使わずにグローバル変数を使ったのは、まずは何より簡単な変換スクリプトにしたいということ、そして、変換後のC言語ソースをコンパイルしたときに、少しは早く動作すると良いなと考えて、そんな作りにしてみました。

続 エクセルの計算ワークシートをRubyでC言語に変換してみよう






2009-12-15[n年前へ]

ComThreadを使った「制御プログラムの作り方」図解編 

 以前、Rubyでテキストや(比較的単純低速な)バイナリデータ送受信のためのクラスであるComThreadを作りました(情報ページ)。それでは、どんな風に使うかということも、これまた以前、ComThreadを使った「制御プログラムの作り方」に書いたことがあります。

 ComThreadという名前が、実にありがちな名前であるので、他のライブラリともしかしたら混同されているのかもしれませんが、「送受信に2つのポートを必要とするのか?」という質問があったので、「そんなことはないですよ」とここに書いておきます(もしかしたら、最初に作成した際の紹介記事の2番目の例が、機器間のシリアル通信をモニタリングする例だったから、そういう風に勘違いされたかもしれません)。

 最初の記事の2番目の例、すなわち、機器間のシリアル通信をモニタリングする例は、AとBという機器間で行われている通信をモニタしたい、AとBは通常通り通信しあう関係にしたままで、その互いのやり取りを知りたい、という場合に、シリアル通信の内容を調べる、というプログラムです。つまり、簡略図にすると、A<->Bという関係をA<->PC<->Bというにしてやり、PCはAから受信した内容をBに送り、Bから受信した内容をAに送るわけです。だから、PCからみると、機器Aへの接続を行うCOMポートと、機器Bへの接続を行うCOMポートの2つが必要だったというわけです。

 もちろん、それは2機器間の通信をモニタリングする例だったからで、(最初の紹介記事の)一番最初に挙げた例のように、1つの機器相手に送受信を行うだけであれば、当たり前のはなしですが、COMポートは1つしか必要としません。

 さて、これだけではつまらないので、今日は、ComThreadを使った「制御プログラムの作り方」で書いたことを図にしてみました。つまり、前回、文章とコードで書いた内容を「図」にしてみました。それが下の図になります。

 前回書いたように、(センサーやモーターといった)各ハードウェア機器に対するモデル・クラスをConThreadクラスを継承させて作り、また、制御を行うコントローラをThreadクラスを継承して作り、そのコンローラインスタンスに、各モデルのインスタンスを渡してやる、というのがComThreadを使い、複数機器を制御するのであれば、一番整理しやすいような気がします。

 つまり、(各ハードウェア)モデル・インスタンス(スレッド)を作成し、コントローラ・インスタンス(スレッド)を作成して、(各ハードウェア)モデル・インスタンスをコントローラ・インスタンスに渡した上で、必要な期間だけそれらのスレッドを動かしておく、というようなスクリプトを作る、という具合です。

ComThreadを使った「制御プログラムの作り方」図解編






2009-12-16[n年前へ]

エクセルでシミュレーション Vol.10 二次元非定常熱伝導問題シミュレーション+P(ID)制御エクセルシートを作ってみよう 

 「エクセルでシミュレーション Vol.9 二次元非定常熱伝導問題の温度変化グラフも作ってみよう」で、二次元非定常熱伝導問題を解くシートを使いながら(実用的範囲では三次元問題と大差がない)、ヒーターで回転するベルトを高温にした場合の、その回転ベルトの温度分布がどのように変化していくかを計算し、グラフ化し、さらにベルト中の一部分の温度を時系列的にモニターする機能を付けた、シートで遊んでみました。

 そこで、今回は、センシングしたベルトの温度を用いて、ヒーターの温度を簡単なPID制御を使ってフィードバック制御するエクセル・シートを作ってみることにしましょう。

PID制御とは、調整量を(現在の)出力値と目標値との差に比例=Proportionalした量、その(過去からの)積分=Integralに応じた量、および(その瞬間に次にどのように変化するかという)微分=Differentialに応じた量にしたがって変える制御です。言いかえれば、現在(Proportional)・過去(Integral)・未来(Differential)の挙動に応じて、制御調整量を変えてやろうという制御です。比較的、古典的な制御手法ですが、現在でも、もっとも一般的な制御手法です。

 まずは、前回のエクセル・シートでヒーター部分を単純に(計算をさせ始めてから=回転ベルトに対するヒーターを動かし始めてから)100℃にし続けた場合の計算過程を示してみます。つまり、何の制御もしない場合です。その計算結果が、下の動画です。上の動画で上に示したグラフが、回転ベルトの温度分布を示したグラフです(回転ベルトを切り開いたように温度分布を示しています)。動画中の下の折れ線グラフは(上のコンター図で灰色丸部の温度の時間変化を示したグラフです)。ちなみに、縦軸は0℃から150℃までで、横軸が時間軸です。

 下の折れ線グラフを見れば、上のコンター図で灰色丸で示した部分の、温度が時間に応じて上昇し、やがて100℃近くになっていることがわかります。

 さて、次にPID制御を行ってみることにしましょう。…とはいえ、最初は簡単のために、「灰色丸で示した部分の温度と目標調整温度である100℃との差をヒーターに足す(ネガティブ・日―ドばっく)」という「P=比例成分=現在の違い」だけを利用したP(ID)制御を行ってみます。つまり、現在の出力と目標出力の差異のみをヒーター出力に(適当な比率で)足し合わせてみるのです。言いかえれば、ベルトの温度が0℃なら、ヒーターを200℃くらいにすることで、急激にベルトを温め温度を調整し、目標温度100度と現在のベルト温度の差が小さくなってきたら、ヒーターの温度を110℃くらいに変えてやる・・・というような制御をしてみます。そんなシミュレーション計算を行ってみた結果が、下の計算結果になります。ちなみに、ヒーターの制御温度はエクセルのシートでB32セルで計算されています。また、灰色丸のセンサ取り付け部分(を示した)の温度を示す下の折れ線グラフは、先ほどと同じく、縦軸は0℃から150℃までの温度を示し、横軸が時間軸となっています。

 この結果動画を見るとわかるように、「P=比例成分=現在の違い」だけを利用したP(ID)制御(=積分成分と微分成分を使わない)制御では、ベルトの温度は早く立ち上がりますが、ベルトの温度は振動して、なかなかすぐに安定してくれません。しかも、早く温度を立ち上げようとすると、温度振動は大きくなってしまい、温度振動を防ごうとすると、なかなか温度が早く立ち上がってくれない、という相反する関係になっています。いわば、(ある程度、減衰がある)強いバネと弱いバネの振動と同じような現象が起きてしまいます。

 さて、それでは、一体どうしたら良い制御ができるのでしょうか。・・・せっかく、簡易にエクセルのような(表計算=スプレッドシート・ソフトウェアを使った)熱伝導方程式を使ったシミュレーション計算を行うことができる環境があるわけですから、今回扱ったPID制御のような古典的な制御をもう少し復習してみようと思います。

 そして、せっかく熱伝導方程式を計算してみたりしているわけですから、そういった微分方程式(つまり現在から未来を示す式)を使うことで実現できる、最適化制御についても考えてみることにします。

2009-12-17[n年前へ]

とても気になる文芸雑誌の「yom yom」 

 文芸雑誌の「yom yom (ヨムヨム) 」には、私が好きでとても読みたい作家たちの名前が並んでいる。だから、いつも買おうかどうか悩む。そして、手に取り頁を少しめくってみる。

  カバンに一冊入れておけば、いつでもどこでも読書を楽しめる文庫本。新しい雑誌「yom yom」は、手軽な文庫で読書に親しんでいる読者に向けた、「読む」楽しみをもっと拡げるための〈ヨムヨム雑誌〉です。
 バックナンバーも気になるものばかりで、一通り揃えて読んでみたくなる。
 「yom yom」はエッセイや紀行文、ブックガイド、インタビューはもちろん、小説も読み切りですので、最初から最後までまるごと一冊、毎号すべて楽しむことが出来ます。お好きなところから頁を開いてみて下さい。
 けれど、お気に入りばかり集めてあるものなのに、これが不思議と本屋でなかなか買えないことが多い。

 好きなものばかりが集まっていると、逆に自分が読みたいものが見えなく・見つけにくくなってしまったりするものなのかもしれない。あるいは、豪華なパフェを前にすると、自分が何を食べたいのがわからなくなってしまうのと似ているのかもしれない。

 文芸雑誌「yom yom」は、私にとっては(かつての)読売ジャイアンツみたいなものなのかもしれない。豪華すぎて、ひとりひとりが見えなくなってしまっているのだろうか。

2009-12-18[n年前へ]

論理の先にある非論理的なサンタの正体 

 時間をかけ、目に見えて役立つように見えるわけでもない文章をただ書き続ける、ということがなかなかできません。以前は、「できるかな?」と題して、そういうことができていたような気がしますが、最近は「技術の話、けれど、技術だけでもない話」というものを書くことができていないように思います。けれど、そういうことを、また再開していこう、と考えています。

 「できるかな?」では、(特に初期以外の時期は)「技術的な計算をしていくと、その計算結果を通して(何らかの)「オチ」にたどり着き、話が終わる」ということを意識して書いていました。それはなぜかというと、「○×が良ければ良いほど、△□は良い」といったような、何かの特異点や境目がないような話、つまりオチ(境目)のないは、何だかつまらない当たり前の技術問題に過ぎないように思えたからです。

 逆に言うと、その「技術が解き明かすオチ」を思いついてから(ある程度そうなるだろうと確信してから)、初めて計算を始め・文章にまとめていました。それは、少し研究・論文書き、といったことと似ているのかもしれせん。

 だから、時に、(たとえば月刊化学でされたような)こんな質問をされても、

Q: 調べていて,意外な結論が導きだされたりすることはあるのですか?
 いつも、こんな感じの答えをしていました。
A: いいえ、結論(オチ)まで考えてから実際の計算や 実験に取りかかるので,意外な答えに至ることはありません。ただ、一回だけ、「意外な答え」に気づき驚いたことがあります。それは、「サンタが街にやってくる」という話です.

 クリスマスシーズンが近づいてきました。だから、もう一度この「サンタが街にやってくる」という話を、ここに挙げておこうと思います。「できるかな?」の内容は3度ほど書籍化されていますが、この話は、一番最初に出版された幻の本「できるかな?―うれしはずかし無敵の科学 (ホームページ・ブックス) 」にしか収録されていません。けれど、すべての「できるかな?」の中でも、私の一番のお気に入りの話です。

 「サンタはいないのだろうか」と自問自答しながら文章を書いている中で、ふとサンタというものの正体に気づかされたとき、自分が考えもしなかった「本当の答え」を見つけたような気がしました。

 最初は少し数学的な話から入りますが、その論理的な計算から導き出される答えは、実に非論理的な「答え」です。けれど、それが私が一番心で納得できた「答え」です。だから、数学部分を乗り越えて、最後まで読んで欲しいと願います。

 ちなみに、、「サンタが街にやってくる」では珍しくBGMをつけてあります。静かに読みたい方は、消音ボタンを押してから、リンクをクリックして頂ければ、幸いです。

2009-12-19[n年前へ]

名前付きセルのあるエクセルのワークシートをC++言語プログラムに変換してみよう 

これまで「続 エクセルの計算ワークシートをRubyでC言語に変換してみよう」や「エクセルの計算ワークシートをRuby計算スクリプトに変換してみよう」で、エクセルの表計算シートを他の言語のプログラムに変換する、なんていうことをしてきました。ここで想定している「エクセルの表計算シート」というのは、「エクセルで作った離散化シミュレーション用.xlsシート」を主眼に置いています。つまり、何らかの方程式を差分化し、エクセルシート上でその差分化された空間を表現し、反復収束計算ことで求めたい結果を得る、そんなためのエクセルシートを主眼に置いています。

 ところで、そんな風にエクセルシートを使うときには、セルに名前をつけたくなります。たとえば、真空の透磁率をA1セルに入れたら、そのセルを"A1"セルでなく、"eps"なんていう名前で参照したくなります。そこで、たとえばこんな風にセルに名前をつけます。すると、そのセルの値を、"A1"というようにも参照できるのに加えて、"eps"なんていう名前でもアクセスすることができるようになります。そうすれば、他のセルでは"=eps*4"なんていう風に数式を入力することができるわけです(参考ビデオ)。

 そこで、今回は名前付きセルを使ったエクセルのワークシート(簡単のために、一枚目のシートだけを使ったエクセルブックを前提にしています)をC++言語プログラムに変換するRubyスクリプトを書いてみました。変換用のRubyスクリプトや、名前付きセルを使った(ラプラス方程式を解く)エクセルシートや、エクセルシートから変換・作成されたC++言語ソースファイルは、ここに置いておきました。このサンプルのエクセルファイルでは、中心のセルの値を変化させるために入れるセルに"centerVal"という名前を付けています。そのため、変換されたC++ソースには、

  C3=centerVal;
といった記述が現われています。エクセルでは、「名前」はセルにではなくブックの直下で管理されているので、Ruby変換スクリプトでは、
names={}
book.names.each do |name|
  cell=$1.gsub('$','') if /!([^!]+)$/=~name.RefersTo
  names[cell]=name.name
end
といったようなコードにしてあります。

 名前付きセルを使ったエクセルシートを使い物理定数や中間変数などをわかりやすく表現した上で、離散化シミュレーションをシート上の「いかにもわかりやすく「差分化された空間で」実感した後に、(その自分自身が作ったスプレッド・シートを元に)変数の名前が(少なくともA1とかD4とかいう名前よりは)わかりやすいシミュレーションプログラムが自動生成されるとしたら、つまり、自分自身がスプレッド・シート上で作ったものが、C++言語の(もちろん他の言語でも)プログラム・ソースとして眺めることができるとしたら、…これって結構面白いと思いません?

2009-12-20[n年前へ]

年賀状作りをモノクロ版画で簡単にしよう 

 年賀状を書く時期になってきましたので、そんな時に(少しは便利な)消しゴム版画風画像作成ソフト「ナンシー"小"関風パッチもん版画」作成ソフト」の使い方でも、動画として(下に)張り付けておくことにします。行っている作業は、適当な画像をソフトウェアのアイコンにドラッグ&ドロップし、画像を作成・(Sキーを押すことで)保存し、(ESCキーを押すことで)ソフトウェアを終了させ、後は適当に”ペイント”ソフトウェアで文字を手描きしている、という具合です。

 カラー化したバージョンや、高画質化したものも作成し、さらにそれをWEBアプリケーションにしたものも作成したのですが、(少し前にサーバのクラスタ構成を変更した際に)動きがおかしくなっているようなので(その問題修正は後回しにして)、とりあえず白黒版画のアプリケーション版の紹介動画を今回は張っておきます。

2009-12-21[n年前へ]

続 「サンタが街にやってくる」 

 自分が間違えて使っていたことを、「それは間違っているよ」と指摘して頂くのは本当にありがたいことです。

 サンタのオチは秀逸だと思います。でも、「犯」とは呼ばないと思います。
 新明解辞書で、「犯」という言葉をひきました。すると、こんなような言葉が出てきました。
 犯す=そうすれば当然罰せられることをする 
 あぁ、本当に仰られる通りなんですね。「犯」ではないですね。言葉を考え直します。心から、感謝させて頂きます。

2009-12-22[n年前へ]

「メディア異人列伝」というエッセイ集 

 噂の真相に12年にわたり、連載されていた、「メディア異人列伝」が好きだった。単なるインタビューというのではなく、永江朗が綴る地の文中に「インタビューイ」の言葉が挿入される。何だか、永江朗という擦り硝子越しにインタビューイが描かれているような、そんな連載だったように思う。

  その単行本「メディア異人列伝 」を読んだ。すると、やはり、描かれているのは確かにインタビューイだが、それを描き出す永江朗の視点の印象の方が遥かに強く残った。つまり、これは、素の「伝記」というよりは、永江朗が書くインタビューイを見ながら書いたデッサン、というか、まるでエッセイ集のように思えた。

 そして、だからこそ、400頁を超える本だが、一気に読みふけってしまった。こういう本は、ぜひ文庫本にして欲しい。そして、本屋の棚の中にいつでもあると良いな、と思う。

2009-12-23[n年前へ]

続 rubyscript2exe.rb のエラー 

 Rubyスクリプトを実行形式(Windowsで言えば.EXE形式)にするrubyscript2exe.rbというスクリプトがある。このスクリプトは、ただ一つのファイルを(変換したいRubyスクリプトと同じディレクトリに置くだけで使うことができるので)実に便利で重宝する。

 しかし、このrubyscript2exe.rbが動作しないことがある。そんなシチュエーションのひとつは以前書いた、rubyscript2exe.rbが"Frozen String error"を出力し動かない場合である。

 そして、他にも、

uninitialized constant
Gem::RubyGemsVersion (NameError)
とか
Couldn't execute this command (rc=256):
というメッセージが吐かれ、Rubyスクリプトを実行形式に変換することができないことがある。

 そんな時は、DOSプロンプトから、

set RUYBOPT=
とか、
set RUBYOPT=-Ke -rkconv
とか、適当にRUBYOPTを設定してやると、rubyscript2exe.rbが動くようになる。つまり、
ruby rubyscript2exe.rb hoge.rb
という具合で、hoge.rbをhoge.exeに変換することができるようになる。

2009-12-24[n年前へ]

現実の物理現象速度を再現してしまうと、違和感のあるアニメになってしまう 

 画像電子学会研究会予稿 09-03-08 の栗田・青木「リアルバーチャリティー -映像の中のウソー」の中の「アニメーションの中では、現実の物理現象の速度をきちんと再現した映像は現実的に見えなくなってしまう」という話が面白かった。この予稿の内容は、コマーシャル用のアニメーション作成中に、アニメーションを作るプロでない人が感じた興味深かったことをまとめたものである。

 木に積もった雪がドスンと落ちる様子や、温泉の湯気が立ち上るようすを実際に撮影し、それをコンピュータ・グラフィクスでアニメーションにしてみると、今一つ感覚と違うという。現実の物理現象速度を再現してしまうと、違和感のあるアニメになってしまうというのである。

 現実の速度と同じ湯気では「湯気の動きが早すぎて、寒く感じてしまう」し、現実と同じ速度でドスンと落ちる雪は「遅すぎて、(登場人物が)驚くような感じにはできない」というのだ。つまり、現実とアニメの中での「速度感覚」というものは、何かしら異なっている、という。(自然とは異なる)非等時的なデフォルメがなければ、「自然なアニメーション」にはならない、というのである。

 逆に、アニメーション世界の時間感覚を現実に持ち込んだらどう感じるものなのだろうか。もしかしたら、それは「現実よりもさらに現実的に見える世界」になっていたりするのだろうか。

2009-12-25[n年前へ]

エクセルでのセルの「名前」管理 

名前付きセルのあるエクセルのワークシートをC++言語プログラムに変換してみよう」でこう書きました。

 エクセルでは、「名前」はセルにではなくブックの直下で管理されている

 エクセルのワークシートをC++言語プログラムに変換するプログラムを書くとき、最初は、Cellオブジェクトが名前を保持している、たとえばNameなんていうプロパティがあるのではないかと思いこんでいました。けれど、プログラムを書いている途中で、「Cellオブジェクトは自分に付けられている名前情報を保持していない」ということに気づき、どこに名前情報が保持されているかを調べてみたら、エクセルのワークブック直下に、つまり"book.names"というような形で、Nameの情報を示す配列が格納されていたのでした。(セルはちなみに、”ブック”の下の”ワークシート”の下の”行”の下の”列”の下にいます)

 そのことに気付いた瞬間は、(エクセルのワークシートをC++言語プログラムに変換するプログラム)の「コードの量が増えそうだな」と思い暗い気分になりました。しかし、少し考え、エクセルを作る側から考えてみれば、「名前はワークブックが保持する作りにするよなぁ」と思ったのです。

 もしも、セルが名前を保持する作りにしていたら、エクセルの動き上、ちょっと面倒になってしまうからです。それは、たとえば、どこかのセルに"velocity"という名前を付けていたとして、(それと異なる)またどこかのセルで"=0.1*velocity"なんていう式を使っていたとします。すると、その"=0.1*velocity"という計算をするときに、(セルが名前を保持する作りにしていたら)その名前を保持するセルが見つかるまでセル検索を行わなければならないからです。

 もちろん、その検索自体は、「エクセルにおける循環参照時の計算順序」で書いたような「エクセルの計算プロセスの詳細」を考えてみれば、エクセルで「(複数ブック・複数シート内のすべての)計算」を一回するのにつき一回だけですむはずだと思いますが、それでも検索作業はかなり面倒でしょう。

 大体、セルが「自分がどういう名前で他のセルから呼ばれているか」なんて知る必要がないわけです。他のセルの名前を知りたいことはあっても、自分の名前なんかどうでも良いのです(今回、私が書いたようなプログラムを書くのでなければ、ですね)。だから、エクセルでのセルの「名前」オブジェクトをブック直下に置くようにするようにするのは、当然の作りです。

 と、エクセルの作りに納得しつつ、(もちろん、自分のプログラムを書くのには、少し面倒だなぁと思いながらも)前回書いたような

names={}
book.names.each do |name|
  cell=$1.gsub('$','') if /!([^!]+)$/=~name.RefersTo
  names[cell]=name.name
end
といったコードを書いたわけです。

 ちなみに、蛇足ですが、上記コードはbookがnameオブジェクトの配列を保持していて、nameオブジェクトはnameという自分の名前と、どこのセルを意味するかを指し示すRefersToというプロパティを持っているので、名前で参照することができるHashを作っているという具合になります。

2009-12-26[n年前へ]

「テレビニュース」=「私たち自身」 

 ふと、今日は、以前書いた文章を見直し、清書してみることにしました。

 私は新聞記者が「起きた出来事を一人でも多くの人に一秒でも速く届けるという役割」を重視する姿勢はあまり価値があるとは思っていません。「起きた出来事」を明らかにするには時間がかかる以上、「起きた(と思われるが、その確からしさは決して高くない)ことを多くの人に速く届けてしまうという自らの判断」の必要性・危険性を「その書き手」はいつも疑い・考え続けるべきだと思います。そして、その必要性・危険性を考え続けながら、それでも、長い時間をかけて事実を明らかにしていく、というやり方が大事だと思っています。どれだけ、長い時間をかけても、です。

 新聞記者たちが、そんな期待に応えてくれているかどうかどうかは疑問に思っています。少なくとも、長い時間をかけてみても、コストパフォーマンスが低ければ、仕事にはならないからです。

 とりあえず、私たちの中にある(少なくとも私の中にある)刹那的な野次馬根性に対して、記者たちは私たちの野次馬根性・期待に対して100%以上の仕事をしていると思っています。

 さて、次は新聞ではなく、テレビの話です。ただ、基本的には多少の差はあれど、とても似ている話です。

 以前、「テレビニュース」はあなたの人生にとってなくてはならないものですか?というアンケート結果を流すTVを見ました。その時、「テレビニュース」というものに対する(私たちの)イメージは、すべて自分たちの上に降りかかってくるに違いない、と確信しました。私たちのイメージの中にある「テレビニュース」という言葉はすべて「ぼくら自身」と全て置き換えられるに違いない、と思いながら「テレビニュース」はあなたの人生にとってなくてはならないものですか?というアンケートへの答えが、刻々とTV画面に流れるようすを知人宅で眺めました。

「テレビニュース」が世の中に与える影響は? -> 悪い
「テレビニュース」は信用できる -> No
「テレビニュース」はあなたの人生にとって
  なくてはならないものですか? -> Yes

 最後の質問まで眺め終わったとき、「テレビニュース」という言葉は、すべて「私たち自身」という言葉で全て置き換えられる、と確かに判断したのです。

2009-12-27[n年前へ]

「ドナドナ」と「コンドルは飛んでいく」 

 細見和之の「ポップミュージックで社会科 (理想の教室) 」が、とても奥深く面白く、まだまだはまっている。今日はまっていた部分は、「ドナドナ」という歌と、その歌詞が歌われた時代だ。「ある晴れた昼下がり 市場へつづく道・・・」と悲しげに歌われるあの、「ドナドナ」である。

 ドナドナは、1940年にが米国に移住してきたユダヤ人、アーロン・ツァイトリンが(東欧ユダヤ人の日常語であった)イディッシュ語で作詞し、ニューヨークの劇場で発表されたものだという。細見和之氏がイディッシュ語版を翻訳した一部を抜粋すると、このような歌詞になる。

荷車のうえに子牛が一頭
縄に縛られて横たわっている。
子牛がうめくと農夫が言う
いったい誰が子牛であれとお前に命じたのか。
お前だって鳥であることができたろうに。
ひとびとは哀れな子牛を縛りあげ
そして引きずっていって殺す。
翼を持つものなら空高く舞い上がり
誰の奴隷にもなりはしない。
 この子牛は、どう読んでも「人」である。日本語の歌詞より、ずっと辛く・苦しい歌だ。

 この歌が、(日本でほど知名度が高いわけではないが)一般的に知られ・歌われるようになったのは、1960年にこの曲が英訳されたジョーン・バエズのアルバムに収録されてからだという(そして、1960年代後半に日本で大ヒットした)。

 ジョーン・バエズの歌を聴いていたその頃の、日本の私たちは一体どのようなことを感じながら、どんなことをしながら、この歌詞を聴いていたのだろう。前回、

 その時代、あの時代、さまざまな時に流れ・歌われた曲を聴きながら、その時折を知る人たちの話を訊いててみたい、と思う。
と書いたが、この曲に関しても、やはりそう思う。

 私がこの歌詞を読んで、ふと思い出したのがサイモン&ガーファンクル が歌った「コンドルは飛んでいく」である。アンデスのフォルクローレの名曲に、ポール・サイモンが英語の歌詞をつけて、よく歌われるようになったあの「コンドルは飛んでいく」である。

地を這うカタツムリになるよりは、スズメになりたい。
もしなれるなら、心からそうなりたいと願う。
人は地面にしばりつけられ、
一番悲しい声を響かせている。
 この曲も同じ時代に歌われたからなのか、ポール・サイモンがやはり同じユダヤ系だからなのか、あるいは、単に悲しげなメロディとモチーフからに過ぎないのか、それとも普遍的な何かがこういった曲の底に流れているからなのか・・・理由はわからない。けれど、「ドナドナ」のメロディと(イディッシュ語版もしくは英語版)歌詞を眺めると、「コンドルは飛んでいく」を、私はふと連想してしまった。

 「ドナドナ」は、日本の小学校で音楽教育を受けた人ならば、多くの人が知っていると思う。人は、たとえば、あなたは、一体「ドナドナ」を聴いて、どんな風景を思い浮かべるだろう。

2009-12-28[n年前へ]

VISTAで使うUSB接続シリアルポート+wincom.rb 

 「TEXCELLのRubyシリアル通信ライブラリ wincom.rb で、(WindowsXP上では受信できるのに)Windows VISTA上では、レガシーなCOM1ではデータ受信できるが、(Prolificのチップ・ドライバーを使った)USB I/F接続シリアルポート接続では受信できなかった」という記事を読んだ。原因は、

 Vista環境にて、USB接続シリアルポートでは、(Windows APIの)ReadFileを実行した際、「読み取ったバイト数」として常にゼロが帰ってくるため。
ということで、USBシリアルI/Fドライバのバグらしいが、考えてみれば、「VISTAマシンでシリアル接続をしたことがないことに気づいたこと」「ATEN製UC-232A・シグマAPO製URS232GFと同タイプを使う可能性もあること」から、リンク先の記事にあった対処法を自分用に書き写しておく。機会があれば、自分でも確認してみよう、と思う。
 対処法:ReadFileの「読み取ったバイト数」は使わず、ClearCommErrorを実行した際に得られたCOMSTAT構造体の「受信バッファにあるデータ・バイト数」を使う(常に0バイトよりはマシだ)。
def receive
# rcvchar=@wcrecv.unpack("a#{irlen[0]}")[0]
rcvchar=@wcrecv.unpack("a#{ilen}")[0]

2009-12-29[n年前へ]

「カルマンフィルタ」と「エクセルで解く2次元非定常熱伝導問題」 

 正月に、(自分用の)汎用「カルマンフィルタ」ライブラリをRubyとCとExcelで書いてみることにした。たとえば、さまざまなデータ、たとえば、信頼性が低く、誤差の大きなセンサデータや、安定性に欠ける実験データから、現実に近い状態量を推定するツールを作ってみることにした。そして、何か(解析式による)モデル計算や各種シミュレーション計算と比較をしてみたり、それらの計算改善へのフィードバック例を作ってみよう、と考えた。

 そこで、扱う題材を考えつつ、実際に上記のようなことを行っている例を探してみた。すると、たとえば、

  1. 簡単なカルマンフィルタの実験
  2. カルマンフィルターと有限要素法を用いた一次元非定常熱伝導問題に対する推定
  3. カルマンフィルタFEMを使った熱の擬似状態フィードバック制御
といったものがある。これらの記事が(下に張り付けた動画でその一端がわかると思うが)実にわかりやすく・面白くて楽しく・役に立ちそうに見える。何というか、つまるところ、魅力を持つに必要な三拍子がすべて備わっている。

 先日、「2次元非定常熱伝導問題を解く」エクセル・シート、しかも、そのシートに、センサ機能/フィードバック機能なども付けてみた。そんな素材・材料が揃ってきたこともあるので、まずは、PID制御で(疑似三次元空間における)温度制御を行う例をいくつか作り、その後は、上記記事を参考にしつつ「(誤差を付加した)センサ→カルマンフィルタ→制御量最適化」という例でも作ってみることにしよう。

2009-12-30[n年前へ]

イディッシュ(Yiddish)語版の「ドナドナ」 

 イディッシュ(Yiddish)語版の「ドナドナ」を聴こうと思い、2001年に'Jewish Entertainment Hour'で歌われたものを、(下に張り付けた)動画で聴いてみた。ドイツ語をもう少し真面目に勉強していたら、なんとなくでも歌詞を聴くことができたのだろうか。

2009-12-31[n年前へ]

Rubyで単純なカルマンフィルタを書いてみた 

 An Introduction to the Kalman Filterの一定出力プラント例に対するカルマンフィルタをPythonで実装した Cookbook / KalmanFilteringを参考に、Rubyで単純なカルマンフィルタを適当に書いてみました。

 動かし方は、

ruby kalman.rb 100
というようにでも起動すると、「出力」→「観測」→「予測」のカルマン・フィルタの動作を99回ほど繰り返す、という具合です。出力結果例は、また後ほど示すころにしますが、開始後およそ5回程度からは、「真の出力値」に近い予測を行うことができています。

 (最後に配列を使ってグラフ化などをするわけでもないことと、一回前のデータ以外は不必要だということがわかりにくくなってしまいそうだったので)本来は不要な配列を使いたくありませんでした。ただ、変数が(一見)増えたように見えるのもどうかと思い、Python例と同じように配列を使った実装にしてしまいました。

 一応、実際に動作する際に「発生すること」「処理すること」の「手順・ループ」を意識しながら、コメント文を書いてみました(正式な用語とは異なるものも多いかもしれませんが)。

 ただし、今回の「一定出力」の例では、「状態方程式」「出力方程式」が登場しないので、そういう拡張をすることを考えた場合には、変数名が適切でないような気がします。そういった辺りのことは、また次の機会に整理し直すことにします。さて、これで、今年を終わりにしたいと思います。

 というわけで、下記がRubyソースになります。

# 簡単なカルマンフィルタRuby実装例

def gaussian(n) # 正規分布作成用適当関数
   sum = 0.0
   ((n*12).to_i).times{ sum += rand() }
   return sum-(12/2*n)
end

operationNum=ARGV[0].to_i   # Operation Num

x=[]               # 出力の真の値
z=[]               # 測定値
xhat=[]           # 事後出力推定値
xhatMinus=[]   # 事前出力推定値
p=[]               # 事後誤差推定値
pMinus=[]       # 事前誤差推定値
kG=[]             # カルマンゲイン

xVall=55.5       # 出力の真の値
q=1.0              # 出力自身が持つノイズの分散
r=2.0              # 観測における誤差の分散
xhat<<20          # 初期出力推定値
p<<0               # 初期誤差推定値

puts 'k,x,z,xhat,kG'
1.upto(operationNum) do |k|
  # DoProcess
  x[k]=xVall+gaussian(q) # 出力値生成
  z[k]=x[k]+gaussian(r)           # 測定値生成
  
  # Time Updade ("Predict")
  xhatMinus[k]=xhat[k-1] # 状態予測を進める
  pMinus[k]=p[k-1]+q   # 誤差共分散計算を進める
  
  # Measurement Update("Correct")
  kG[k] = pMinus[k]/( pMinus[k]+r) # カルマンゲイン算出
  xhat[k] =xhatMinus[k]+ #測定値を用い事後出力推定値算出
       kG[k]*(z[k]-xhatMinus[k]) 
  p[k] = (1-kG[k])*pMinus[k]  # 事後誤差算出

  # この回のオペレーションでの結果を表示
  puts [k,x[k],z[k],xhat[k],kG[k]].join(',')  if k>1
end