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の法則
定常問題を解くならガウスザイデル法で収束するまで反復計算を行うので構いませんが、非定常問題を解くのにガウスザイデル法を用いてしまうと、正確ではない間違った温度傾斜が生じてしまいます。それは、少し困ります（しかし、現在出版されている”表計算＋差分化 熱伝導”書籍は、多くがガウスザイデル法を用いています）。
というわけで、この計算シートではグラフの下に（動画の）上のメッシュ部分をコピーするだけのセルがあり、上のメッシュ部分は「自分の周囲に相当する下のメッシュを用いてメッシュ間の計算を行う」という仕組みになっています。こうするとなぜ、ガウスザイデル法でなく、ヤコビ法の計算になるかは、「エクセルの反復計算時のセル間の計算順序」を念頭に置いて考えてみると、「なるほど！」と納得できると思います。「道具を理解し使いこなす」ということの意味を、心から実感させられます。
だから、この二次元非定常熱伝導問題を解くエクセル・シートは、とても「単純明快な計算がされているだけ」なのですが、実に巧妙に考え抜かれた上で作られた「素晴らしき一品」です。しかも、セルを少しづつ書き変えていくことで、かなり実用的な計算もできるのです。このシートを見たときには、本当に驚かされました。
「どんな風に書きかえると何がわかるか」とか「さらにこんなこともできる」といったことは、明日にでも書いてみようと思います。（というわけで、この続きとして、温度が時間を追って上昇していく機能も付けてみました）
2009-12-11[n年前へ]
■エクセルでシミュレーション Vol.9 二次元非定常熱伝導問題の温度変化グラフも作ってみよう
この記事は、「エクセルでシミュレーション Vol.8 二次元非定常熱伝導問題を解こう」の続きです。前回の記事では、比較的薄い円形ベルトが回転するモデルで、ベルトの厚みがベルト内の熱拡散に対して十分に薄いという仮定の下での近似を行ったと考えると）三次元非定常熱伝導問題を解いているとも言えなくない・・・というような二次元非定常熱伝導問題を解く、エクセルの表計算シートを作り眺めてみました。もちろん、VBAはもちろん関数も使わず、ただセル間の加減乗除演算のみで、熱伝導方程式を差分化したものです。
今日は、計算領域の任意の部分の温度変化をグラフにする部分を前回のエクセルシートに追加してみようと思います。つまり、ベルト状に温度センサでも貼り付けたら、時々刻々と一体どんな温度が測定されるのかをシミュレーション計算時に表示させたい、というわけです。
このような需要は非常に多くあると思いますが、一体それをどのようにすれば作ることができるでしょうか？エクセルのシート上は、非定常問題を解いているので、その瞬間その瞬間の温度分布しかデータが残っていないように思われ、スタート時からの温度変化をグラフにするのは難しいように思えてしまいます。
しかし、「ガウスザイデル法でなく、ヤコビ法を用いて反復計算を行うために使ったテクニック」、すなわち、
エクセルは（セル間の循環参照がある場合に用いる）反復（収束）計算時には、表の左上をスタート地点として「Zの法則 」の順序にしたがって（つまり、左上→右、そして一個下の行をまた左→右という順番で）逐次的に解くという計算順序になります。ということを利用すると、「計算領域の任意の部分の温度変化をグラフにする」ということが、いとも簡単に実現できてしまうのです。種明かしは後回しにして、まずはそのようなことをしているエクセルシート動画を下に張り付けてみます。動画の上部分は前回と同じく「回転するベルトを平面上の等高線色付きグラフとして表示したものです。そして、その下にあるのが、「等高線色付きグラフ右横部にある灰色円部分」の温度変化を横軸：時間・縦軸：温度で時系列グラフにしたものです。
下の折れ線グラフを見れば、「等高線色付きグラフ右横部にある灰色円部分」の温度が時間を追って上昇していくようすをよく理解することができるのではないでしょうか。
種明かしはこうです。エクセルのシートのずっと下方（下の行のセルに）"=温度変化を知りたいセル"という式を入れます。ここでいう温度変化を知りたいセルというのは、たとえば、X24といったセルを示す番号です。ですから、実際に（一番下のセルに）入れるのは、"=X24"といったものになります。
そして、そのセルの上にあるセルをクリックし、=と入力した上で、その下のセルをクリックし"Enter"を押します。つまり、一個下の（さきほど"=X24"という式を入力した）セルを参照するようにするのです。そして、そのセルを選択した上で、上の方までズルズルズル～とコピーしてしまいます。
エクセルの（セル間の循環参照がある場合に用いる）反復（収束）計算時の計算順序を考えてみれば、今「ズルズルズル～とコピー」を行った列には、「古い計算結果が上、一番下が最新の計算結果」という順序のデータが格納されていくことになります。（上から下に計算が逐次行われていることを考えると）計算ステップごとに、一個値が上に移動していくので、結果としてそのように時系列的なデータを保持できる、というわけです（結局のところ、ヤコビ法を用いるために、バッファエリアを下に設けたのと全く同じパターンです）。
あとは、計算スタート時点からの（シミュレーション上の）経過時間なども同じように作り、散布折れ線図でも挿入すると、「上記のような「任意の部分の温度モニター機能付きの疑似三次元非定常熱伝導問題を解くエクセルシートのできあがり」となるわけです。もちろん、温度モニターは（上段で行ったことを他のセルに対してもしてやるだけで）いくつでも設置することができます。実に単純・簡単な（けれど巧妙な）実装ですがとても便利な機能です。
さて、ベルトの温度モニター機能も付きましたので、今度はセンシングしたベルトの温度を用いて、ヒーターの温度を（まずは簡単なPID制御でも使って）制御しベルトの温度分布を適正に調節する機能例でも、簡単に実装し（ハードウェアをいじっている）気分にでも浸ってみることにしましょうか。
2009-12-14[n年前へ]
■続 エクセルの計算ワークシートをRubyでC言語に変換してみよう
前回、「エクセルの計算ワークシートをRubyでC言語に変換してみよう」でサンプルに作ったエクセルのシートが反復収束計算をする必要がない（自己参照型のシートになっていなかったこと）に気づいたので、今日はエクセル・シートを作り直し、サーバ内に置いておきましたここに置いておきます。置いたエクセル・ファイル"ex2.xls"は、5×5のセルの周囲と中央のセルが固定条件で、それ以外の（5×5内の）セルが自分の上下左右のセルの平均の値をとるというラプラス方程式を差分化したエクセルシートです。そして、xls2c.rbがエクセル・ファイルをC言語ソースに変換するRubyスクリプトで、ex2.xlsをC言語に変換した結果ファイルが、ex2.cになります。
元々の記事は、「とても無意味なこと・けれど技術発想的には何だか（的外れ具合が）面白いこと」をしよう、という考えのもとにやってみたのですが、非定常熱伝導計算の復習をしているうちに、非定常問題を解くシートに対しては、意外に使えるかもしれない、という気がしてきました。
プログラムというものをまったく知らない人が、熱伝導方程式の物理的な意味を理解上で、表計算でシミュレーションをしてみた（習ってみた）あとに、（エクセルのセルと直感的に対応した）自動出力したＣ言語ソースをもとに「printfの使い方だけを習い」色んな計算をしてみれば、そして注目するセル＝箇所の時間変化・経過を出力させてみたりすれば、案外面白く感じるかもしれない・・・という風にも思ったのです。
もちろん、こんな単純自動変換Cソースでは一瞬のうちに物足りなくなるはずなので、上のような時間を経ることで、プログラム言語を学ぼうとするまでの小さな「足がかり」「キッカケ」になるのではないか、という間違っていそうなそんなことをふと思ったのです。いわば、エクセルを入口に差分化シミュレーション・プログラムの世界に"Hello World"する、というわけです。ゼロから学ぶのと、何か叩き台があって、それをほんの少しづついじりながら学ぶのとでは、気楽さが全然違うのではないか、といったことを「後付けで」考えてみたりしたのです。
ちなみに、エクセルシート→C言語 変換などで、配列を使わずにグローバル変数を使ったのは、まずは何より簡単な変換スクリプトにしたいということ、そして、変換後のC言語ソースをコンパイルしたときに、少しは早く動作すると良いなと考えて、そんな作りにしてみました。
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制御のような古典的な制御をもう少し復習してみようと思います。
そして、せっかく熱伝導方程式を計算してみたりしているわけですから、そういった微分方程式（つまり現在から未来を示す式）を使うことで実現できる、最適化制御についても考えてみることにします。
