Grasshopperで単純な最適化ループを回す

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

grasshopperには、"Galapagos"という最適化コンポーネントが付属していますが、これは大域最適化(遺伝的アルゴリズム、焼き鈍し法)を志向したもので、それなりに計算コストがかかる、GUI操作を必要とする、といったことがネックでした。
　　
# 問題設定
以下のようなサーフェスデータを用意します(見やすくするために、等高線を黒線で示しています)。このサーフェス上に任意の初期座標(赤×印)を定め、そこから最も近い最小値(局所最少解)を求めます。
![ViewCapture20221002_173453.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/4954716a-0087-8116-40f2-1b7b93d3ef2c.png)


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

![スクリーンショット 2022-10-02 17.36.20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/eeec3a78-8ba6-36a3-11d7-276f2d34878e.png)


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

```python: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

```


# 結果
初期座標を赤×印、収束した座標を赤○印で示します。
![ViewCapture20221002_173436.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/4521f98f-33c4-474a-7ecd-02674a96433f.png)
いい感じに収束しました。


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

・ケース1
![ViewCapture20221002_171820.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/5cd08f51-42fa-e202-a480-29793bd2411a.png)
・ケース2
![ViewCapture20221002_171840.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/5e406b47-a032-5003-053e-4dbf2a92ea14.png)
・ケース3
![ViewCapture20221002_171911.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/a31de79d-3ec8-1de8-1e5a-53ec38f56faa.png)
・ケース4
![ViewCapture20221002_171943.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/2719339e-edbd-6463-5aae-67a4b0a085df.png)
・ケース5
![ViewCapture20221002_172030.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/8a4fe774-ceae-5a85-379b-3d01c7f493f0.png)
・ケース6
![ViewCapture20221002_172051.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/f16a042a-89fc-5a1c-4a95-edfc13b35fe0.png)


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

# Grasshopper上の操作
インタラクティブに動かすことができます。
![動画.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/d42401a2-616f-f3ea-8991-5fb078291846.gif)

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