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

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

　おおむね好評だったようで嬉しい限りです。コメントをくださった方、ありがとうございます。
　もっと詳しく説明しろとのコメントがいくつか寄せられていたので、解説します。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」はこんな感じです。

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」はこうなりました。

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年間あるので、いろいろと考えて見ることにします。