はじめに

rhinoceros/grasshopper内で、サクッと最適化ループが回せると便利です。

grasshopperには、"Galapagos"という最適化コンポーネントが付属していますが、これは大域最適化(遺伝的アルゴリズム、焼き鈍し法)を志向したもので、それなりに計算コストがかかる、GUI操作を必要とする、といったことがネックでした。
　　

問題設定

以下のようなサーフェスデータを用意します(見やすくするために、等高線を黒線で示しています)。このサーフェス上に任意の初期座標(赤×印)を定め、そこから最も近い最小値(局所最少解)を求めます。

コンポーネント

GH-Pythonコンポーネントの"Surf"の端子に先ほどのサーフェスのデータ、"Point"の端子に初期座標のデータを入れます。
最適化計算の履歴が"history"端子から出力されます。

コード

コードを以下に示します。

gh-python

import rhinoscriptsyntax as rs
import Rhino
import Rhino.Geometry as geom
import time

# SurfをBrepに変換する
B = Surf.ToBrep()

# PointをBrep上に投影する
pp = geom.Intersect.Intersection.ProjectPointsToBreps(
     [B],
     [Point],
     geom.Vector3d(0, 0, -1),
     0.000001  #torelance
     );

# 近傍点のu,vを返す
_,u,v = Surf.ClosestPoint(pp[0]);

# 勾配を返す関数
def get_gradient(S, x):
    # u,vの座標とベクトルを返す
    _,pt,vecs = S.Evaluate(x[0], x[1], 1)  #(u,v,derivative)
    # u,v方向の勾配を取得
    grad = [vecs[0][2],vecs[1][2]]
    return pt, grad

# 最適化の計算設定
x = [u,v]  # 初期値
alpha = 0.05 # 係数
hist = [] #履歴
pt,_ = get_gradient(Surf, x)

# 勾配降下法
dz = 1
step = 0
st = time.time() #処理開始時間
prev_dz = 0
while (dz > 0.000001)&(step<1000)&(abs(pt[0])<10)&(abs(pt[1])<10):
    pt,grad = get_gradient(Surf, x)
    x = [i - alpha*j for i,j in zip(x, grad)]
    hist.append(pt)  # 履歴に追加
    
    # 収束判定値
    if step!=0:
        dz = abs(prev_z - pt[2])
    prev_z = pt[2]
    
    #次のステップへ
    step += 1
    

# 収束値
pt,_ = get_gradient(Surf, x)
hist.append(pt)

print("elapsed step:", step, "time:", time.time()-st)

# ghコンポーネントの出力
history = hist

結果

初期座標を赤×印、収束した座標を赤○印で示します。

いい感じに収束しました。

初期値を変えて、いくつか実験してみます。

・ケース1

・ケース2

・ケース3

・ケース4

・ケース5

・ケース6

ケース5、6は最適化の設計空間(このサーフェスだと-10〜10)からはみ出したタイミングで、反復計算を打ち切っています。
計算速度は初期値にもよりますが、最大でも0.002s程度でした。

Grasshopper上の操作

インタラクティブに動かすことができます。

まとめ

Grasshopper上で単純な勾配降下法を実装した。