「最適化数学で新年の挨拶」の解説

 みなさん、あけましておめでとうございます。
 新年開始とともに、ニコニコに動画を投稿しました。

 おおむね好評だったようで嬉しい限りです。コメントをくださった方、ありがとうございます。
 もっと詳しく説明しろとのコメントがいくつか寄せられていたので、解説します。2と9の文字は、両方とも4つの点とその周りを公転する4点、更にそれを回る4点を元に作られます。以下、それぞれを恒星、惑星、衛星と呼ぶことにします。4衛星を2組にわけ2本の直線を作り、その交点に文字を描かせます。
 それぞれの点の座標は、時刻をtとして次のように表しました。

p1=(a1*cos(2(2t-e1)),a1*sin(2(2t-e1)))+p5
p2=(a2*cos(2(2t-e2)),a2*sin(2(2t-e2)))+p6
p3=(a3*cos(2(2t-e3)),a3*sin(2(2t-e3)))+p7
p4=(a4*cos(2(2t-e4)),a4*sin(2(2t-e4)))+p8
p5=(a5*cos(3(t-e5)),a5*sin(3(t-e5)))+(c5,d5)
p6=(a6*cos(3(t-e6)),a6*sin(3(t-e6)))+(c6,d6)
p7=(a7*cos(t-e7),a7*sin(t-e7))+(c7,d7)
p8=(a8*cos(t-e8),a8*sin(t-e8))+(c8,d8)

(c5,d5)~(c8,d8)が4恒星の座標、p5~p8が4惑星の座標、p1~p4が4
衛星の座標となります。a1~a8は公転半径の大きさ、e1~e8は公転での位相を決定します。またp1とp2を通る直線、p3とp4を通る直線を結びその交点を求め、これを点pとします。pの座標は、複雑な形になりますが方程式を解けば式として表すことが可能です。a*,c*,d*,e*の合計24変数を調整し、pが文字を描くようにすれば良いわけです。
 なおそれぞれの式のtには1,2,3の係数が掛かっていますが、これは各点が回転する速度を表します。t=2Piで全ての点が元に戻るようにしたかったので、係数は整数である必要があります。この値は、実は手動で調整しました。手で描くとわかるのですが、「2」も「9」も一度普通になぞって逆方向もなぞり元の位置に戻ると、X軸方向は3往復、Y軸方向は1往復します。このことから縦に動く線とその3倍の速さで横に動く線ができることを期待してp5~p8でのtの係数、すなわちは惑星の公転速度は1および3と決定しました。これは上手く働いたようです。p1~p4でのtの係数、これは衛星の公転速度になりますが、いろいろと試行した末に上の形になっています。今見ると変な係数の付け方で計算していますが、どうせ最適化すれば同じ式になる(ことが期待できる)ので気にしなくて良いかと思います。
 さて、pにどのような軌跡を描かせたいのかを決めます。これはフリーハンドで描いてその座標を読みました。まず「2」はこんな感じです。
2.PNG

t=0Pi (-0.437739, 0.64179)
t=1Pi/28 (-0.286784, 0.88061)
t=2Pi/28 (-0.0754772, 0.940305)
t=3Pi/28 (0.166043, 1.)
t=4Pi/28 (0.437739, 0.940305)
t=5Pi/28 (0.588657, 0.910438)
t=6Pi/28 (0.679222, 0.761181)
t=7Pi/28 (0.739611, 0.611923)
t=8Pi/28 (0.709435, 0.432838)
t=9Pi/28 (0.588657, 0.223886)
t=10Pi/28 (0.528304, 0.164191)
t=11Pi/28 (0.407526, -0.0149336)
t=12Pi/28 (0.226395, -0.164191)
t=13Pi/28 (0.0150882, -0.283581)
t=14Pi/28 (-0.166043, -0.432838)
t=15Pi/28 (-0.316997, -0.582095)
t=16Pi/28 (-0.467915, -0.731353)
t=17Pi/28 (-0.649046, -0.850743)
t=18Pi/28 (-0.8, -0.970133)
t=19Pi/28 (-0.679258, -1.)
t=20Pi/28 (-0.467915, -1.)
t=21Pi/28 (-0.256608, -0.970133)
t=22Pi/28 (-0.0754772, -0.970133)
t=23Pi/28 (0.13583, -0.940305)
t=24Pi/28 (0.286784, -0.940305)
t=25Pi/28 (0.437739, -0.970133)
t=26Pi/28 (0.618869, -0.970133)
t=27Pi/28 (0.739611, -0.940305)
t=28Pi/28 (0.8, -0.910438)
t=29Pi/28 (0.739611, -0.940305)
t=30Pi/28 (0.618869, -0.970133)
t=31Pi/28 (0.437739, -0.970133)
t=32Pi/28 (0.286784, -0.940305)
t=33Pi/28 (0.13583, -0.940305)
t=34Pi/28 (-0.0754772, -0.970133)
t=35Pi/28 (-0.256608, -0.970133)
t=36Pi/28 (-0.467915, -1.)
t=37Pi/28 (-0.679258, -1.)
t=38Pi/28 (-0.8, -0.970133)
t=39Pi/28 (-0.649046, -0.850743)
t=40Pi/28 (-0.467915, -0.731353)
t=41Pi/28 (-0.316997, -0.582095)
t=42Pi/28 (-0.166043, -0.432838)
t=43Pi/28 (0.0150882, -0.283581)
t=44Pi/28 (0.226395, -0.164191)
t=45Pi/28 (0.407526, -0.0149336)
t=46Pi/28 (0.528304, 0.164191)
t=47Pi/28 (0.588657, 0.223886)
t=48Pi/28 (0.709435, 0.432838)
t=49Pi/28 (0.739611, 0.611923)
t=50Pi/28 (0.679222, 0.761181)
t=51Pi/28 (0.588657, 0.910438)
t=52Pi/28 (0.437739, 0.940305)
t=53Pi/28 (0.166043, 1.)
t=54Pi/28 (-0.0754772, 0.940305)
t=55Pi/28 (-0.286784, 0.88061)

また「9」はこうなりました。
9.PNG

t=0Pi (0.65, 0.830183)
t=Pi/21 (0.273674, 0.924517)
t=2Pi/21 (0.034217, 1.)
t=3Pi/21 (-0.273674, 1.)
t=4Pi/21 (-0.444739, 0.924517)
t=5Pi/21 (-0.65, 0.716971)
t=6Pi/21 (-0.615804, 0.547154)
t=7Pi/21 (-0.513173, 0.396214)
t=8Pi/21 (-0.205282, 0.377363)
t=9Pi/21 (0.034217, 0.433969)
t=10Pi/21 (0.239478, 0.547154)
t=11Pi/21 (0.444739, 0.667552)
t=12Pi/21 (0.513132, 0.735849)
t=13Pi/21 (0.615762, 0.84906)
t=14Pi/21 (0.513132, 0.547154)
t=15Pi/21 (0.478935, 0.45282)
t=16Pi/21 (0.478935, 0.30188)
t=17Pi/21 (0.410501, 0.0188773)
t=18Pi/21 (0.273674, -0.226423)
t=19Pi/21 (0.136848, -0.509426)
t=20Pi/21 (-0.0000205261, -0.811306)
t=21Pi/21 (-0.136848, -1.)
t=22Pi/21 (-0.136848, -1.)
t=23Pi/21 (-0.0000205261, -0.811306)
t=24Pi/21 (0.136848, -0.509426)
t=25Pi/21 (0.273674, -0.226423)
t=26Pi/21 (0.410501, 0.0188773)
t=27Pi/21 (0.478935, 0.30188)
t=28Pi/21 (0.478935, 0.45282)
t=29Pi/21 (0.513132, 0.547154)
t=30Pi/21 (0.615762, 0.84906)
t=31Pi/21 (0.513132, 0.735849)
t=32Pi/21 (0.444739, 0.667552)
t=33Pi/21 (0.239478, 0.547154)
t=34Pi/21 (0.034217, 0.433969)
t=35Pi/21 (-0.205282, 0.377363)
t=36Pi/21 (-0.513173, 0.396214)
t=37Pi/21 (-0.615804, 0.547154)
t=38Pi/21 (-0.65, 0.716971)
t=39Pi/21 (-0.444739, 0.924517)
t=40Pi/21 (-0.273674, 1.)
t=41Pi/21 (0.034217, 1.)

 改めてみると字が下手すぎですね。フリーハンドで描いたのは、別に味のある文字を使いたかったとかいうのではなく、単に時間が無かったからに他なりません。
 時刻tにおけるpの座標をp(t)とします。「2」の場合はt=0~55Pi/28、「9」の場合はt=0~41Pi/21まで変化させ、上で求めた各点とp(t)との距離の2乗をそれぞれの時刻で計算し、合計します。この式は24変数を含むスカラーになることがわかります。この巨大な式を評価式として最小化すれば良いわけです。ただ評価式が複雑すぎて解析的な解は求まらないかと思います。そこで最小化には遺伝的アルゴリズム滑降シンプレックス法を組み合わせた自作のプログラムを使用しました。なお滑降シンプレックス法は、名前こそ似ていますが線形計画法で登場するシンプレックス法とは直接の関係は無いので注意が必要です。
 結果はこのようになりました。
「2」の場合

a1=-0.281248,a2=-0.217597,a3=-0.133595,a4=0.135594,
a5=0.233587,a6=-0.533821,a7=-0.865466,a8=-0.774918,c5=0.296162,
c6=0.128995,c7=-0.0170449,c8=-0.310357,d5=1.06265,d6=-0.197444,
d7=-0.199307,d8=-0.308207,e1=0.311494,e2=0.218138,e3=0.0167031,
e4=1.52846,e5=0.952566,e6=0.0411405,e7=1.25651,e8=1.3886

「9」の場合

a1=-0.197302,a2=0.145299,a3=-0.197725,a4=-0.200862,
a5=-0.524991,a6=0.779246,a7=-0.472475,a8=1.29185,c5=0.117048,
c6=0.560647,c7=0.611663,c8=-0.875045,d5=0.415859,d6=-0.461184,
d7=1.32743,d8=-0.340324,e1=1.35564,e2=-0.552802,e3=-1.25281,
e4=0.787604,e5=0.866171,e6=-0.150847,e7=-0.496076,e8=4.18745

 円の半径a*の値が負になったりしていますが、p1~p8の式のまま素直に描画すれば良いので特に問題はありません。最適化の限界で、またそもそも変数の数の限界から、元の点とぴったり同じ軌跡を描くことはできません。できあがった動画を見ても、往復で軌跡がずれているのがわかります。なお、円の数を増やせばもっと綺麗な文字の形になるのかとのコメントがありましたが、確かに綺麗になる可能性はありますが変数の数が増えるので探索空間が爆発的に増大し最適化が困難になるかもしれません。この辺りのさじ加減にはなかなか難しいものがあります。
 以上が「2」と「9」の描き方です。「0」については楕円で代用しています。あのように棒を動かすとその棒上の点が楕円を描くことは簡単に証明可能です。
 好評だったので来年の元旦にもまた何か投稿しようと思います。今度は遺伝的プログラミングなどを使って、どの円とどの円を組にするかなども含めて最適化すると面白いかな、と思っています。あるいは全く違う生物的なネタも考えています。ただ、「2010」という文字が簡単すぎて今ひとつ創作意欲をかき立てられないですね。まだ1年間あるので、いろいろと考えて見ることにします。

8 個のコメント

1 個のping

Skip to comment form

    • ワタボン on 2009年1月5日 at 1:40 PM
    • 返信

    ニコニコ動画から来ました情報系の大学生です。
    すばらしい動画ですね。
    ですが何をやっているのかさっぱりです(汗)
    で質問なんですが、こういうのをやろうと思ったらどういう本などが参考になるのでしょうか?

    • めろん on 2009年1月5日 at 6:36 PM
    • 返信

    ニコニコ動画から来ました工業系の高校生です。
    すごかったんですが全く分かりません・・・orz
    この記事見たら余計に混乱しましたわわわわわ。
    まだ高1なので卒業する頃には理解してやります!!

    • gotcha@ぴっちぶれんど on 2009年1月5日 at 10:05 PM
    • 返信

    >ワタボンさん
    私が愛用してるのは、この本です。
    http://www.amazon.co.jp/dp/4874085601
    滑降シンプレックス法は、ここに載ってるのをほぼ丸写しして使っています。実用的な数値計算の本としてはバイブル的な存在だと思います。
    遺伝的アルゴリズムに関しては、どの本でも良いから読んで自分で簡単なプログラムを作ってみるのが一番良いですね。この本なんかが評判良いです。
    http://www.amazon.co.jp/dp/4274066274
    今回のような動画を作るには、私の場合は過去からの蓄積したライブラリがあったのですが、それでも時間が無くてスパゲッティプログラムになってしまいました。一から作ろうとすると大変だと思います。
    >めろんさん
    いや、いろいろ端折った説明で申し訳ないです。私も高校生の頃だったら理解できなかったと思います。ただ、分からない単語があってもググればいろいろ出る今の時代はうらやましい。
    要するにですね、CPUをフルに使って力任せにさまざまな円の配置を試したら、たまたま上手く数字が描けるような組み合わせが見つかったということなんです。さまざまな配置を試す際に、ランダムに配置を作るよりは多少効率が良い方法を使ったといったところでしょうか。ちなみに、なぜあの配置でちゃんと文字が描けるのかは、up主でも説明不可能です。

  1. 衛星地図

    衛星地図はズームダウンすることにより、航空写真よりも、もっと大きなエリアで見ることができます。もちろん、衛星地図はかなり細かいところまで見ることもできます。

  2. 古い記事にコメントしてすみません。
    動画、楽しませていただきました。:-)
    この動画にPhunのタグが(タグロックされずに)付けられているのですが、実際にPhunで作成されているのであれば、ぜひシーンを公開していただけないでしょうか?
    動画の説明文で、
    > Phunなどを使ってピタゴラ的な感じの年賀状が作れないだろうかと思って始めたのですが、日頃から最適化のお世話になっている身分としてはこういうものが良いかなと方針転換しました。
    と書かれていることから、「Phunではないかも」と思いつつも、もしPhunで実装されていたら詳しく見てみたいと思いコメントさせていただきました。
    よろしくお願いいたします。

    • gotcha@ぴっちぶれんど on 2010年10月27日 at 11:29 PM
    • 返信

    >Tatt様
    まぎらわしい書き方をして申し訳ありません。Phunではないのです。連続した静止画を作製するところまでは、全て自作のライブラリを使いました。
    円の位置や大きさは本文に書いたので、頑張ればPhunでも作れるかもしれませんが…

  3. なるほど!ありがとうございます!
    Phunで言うところの「トレーサー」の色の重ね方がPhunにそっくりで驚きました。 🙂

  4. いいまちがい と とらえたのでしょうか。

  1. […] この方のブログへ行き解説を読んでみました。(http://kogarashi.net/archives/12) これによると、 […]

コメントを残す

メールアドレスは公開されません