2019自作CNCマシン・レーザーカッターについてMac+Arduino:自作CNCマシンの記録。CNCシールドV3.5+bCNCを使用中。BluetoothモジュールおよびbCNCのPendant機能でスマホからもワイヤレス操作可能。その他、電子工作・プログラミング、機械学習などもやっています。MacとUbuntuを使用。レーザーシールド制作・改良*CNCマシンの制作記録は2016/04/10〜の投稿に書いてあります。scipy.optimize.linprog今回はscipyを使ってTSP(巡回セールスマン問題)を解いてみました。前回までのPulpのような最適化ソルバー/モデラーであれば0-1整数計画法(混合整数計画法)として計算可能でしたが、scipyの線形計画法の関数であるlinprogは整数計画法まではカバーされていないため、離散値(整数値)ではなく連続値として計算しなければいけません。Pulpを使ったほうがコーディングも楽だし処理も早いと思いますが、今回は敢えてscipyで可能かどうかという試みです。アルゴリズム的には基本的に前回と同じですが、今回は連続値を扱うのとscipyのフォーマットで書き直すため微妙に異なる点があります。データは前回同様TSPLIBのアメリカ48都市を用いることにします。目的関数とベースの制約:最小化する式は以下のようになります。2とりあえず各都市において最低2本の経路が接続されていなければならないという条件を加えておきます。j)は経路として採用された場合は1、不採用なら0になりますが、scipyの場合はこの値が連続値になるので0.0〜1.0の小数になります。上記をscipy.optimize.linprogのフォーマットで書き直すと以下のようになります。j)のマトリクスになりますが1次元配列として用意します。48*47/2=1128個の値を含んだ配列になります。j)の値の範囲(今回は0.0〜1.0の連続値)この制約で一度答えを出してみます。simplexを用いることにします(そうしないとWarningが出る)。j)の値がそれぞれどのくらいになっているか見てみます。j)は0か0.5か1.0(縦軸)になっています。j)をx[n]のインデックス番号に変換した値です(n=0〜1127)。実際これを地図に落とし込むと以下のようになります。j)の値が0.0〜1.0の連続値になっているため、この段階では0.5と1.0の二種類の経路が混在しています(もしバイナリ設定が可能なPulpなどのライブラリを使えば0か1の表現になるのでもう少し簡単になるのですが)。部分巡回路除去:次の段階として、この状態から部分巡回路を選択して、その部分巡回路に追加の制約を与えます。目につくのが、右上にある7都市からなる部分巡回路と下部右寄りにある五角形の部分巡回路です。j)=1.0で構成されているもの(実際は0.5の経路群を含めても可)。複数ある場合は都市数が少ない方をターゲットとします。40}として、以下の制約を追加して再度解きます。40)]であり、st_flatというのは経路群を都市群に変換したものです。simplex")閉じている部分巡回路に対して部分巡回路外の都市と接続するように仕向けるための制約です。j)の組み合わせが2個以上という制約です。2scipy.optimize.linprogでは基本的に「<=」が使用されるため「>=」の場合はをマイナス反転するようです。j)が条件に従う場合は-1でそうでない場合は0という配列にしておきます(この部分が多少面倒)。b_ubのほうにもマイナスをかけて「-2」にしておきます。この結果、五角形の部分巡回路は以下のようにつながります。次のターゲットは右上にある7都市の部分巡回路です。同様に制約を追加して解いていきます。部分巡回路がない場合(接続数1):部分巡回路除去の手順を繰り返すと徐々に部分巡回路はつながっていき、以下のような状態になります。この状態では分離した部分巡回路はありません。このようなひとつながりでありつつ巡回路としては不完全な状態のときには、条件を変えて次のターゲットを選ぶことにします。23}で構成されています。{41}の1箇所で全体に対してつながっているので、この都市群に最低でも二本の経路でつながるような制約を追加します。42}も{35}の1箇所で辛うじてつながっています。このような1箇所でしか接続されていない都市群を探し出して二経路接続の制約を追加するというのが今回の手段になります。これらの都市群の数は全体の都市数の半分未満の数でなければならないという条件と、もし複数ある場合は合計都市数の少ない方を選択するということにしておきます。j)=1.0の部分巡回路/部分経路を足し合わせる・足し合わせた都市群は全体都市数の半数未満でなければならないという感じです。この部分はかなり面倒な手続きになりますが、そのような都市群を探し出す関数を用意することにしました(あとで登場)。42}に制約を与えると以下のようになります。j)=0.5の赤い経路は消えてくれました。以後は、これまでと同様に・分離した部分巡回路(全体の半分未満の都市数からなる)・分離した部分巡回路がない場合、1箇所でつながっている都市群という条件で制約を追加していきます。部分巡回路がない場合(接続数3以上):前述の手順を繰り返していると、以下のようなパターンもでてきます。全体としてはつながっており、1箇所だけでつながっている部分もありません。j)=0.5で構成される赤い部分巡回路で他の都市群に接続している部分が3箇所以上あるものを選び、他の都市群への接続数を2箇所以下に制約します。j)=1.0の部分経路は赤い巡回路に足し合わせます。40}を今回のターゲットとします。この都市群は{13}と{22}と{28}の3箇所において他の都市群とつながっているので、この接続数3を2以下に制約します。以下がこの部分のコード。e])b_ub.append(2)この制約で一旦解いてみると、以下のようになります。j)=1.0の部分巡回路だけになりました。部分巡回路がなくなってひとまとまりになります。これで最適解が得られました。制約を与える手順:今回は0.5の経路と1.0の経路(実際は0の経路もあるけれども表示されない)があるので、制約を与えるターゲットをどのように選ぶかが複雑でしたが、要約すれば以下のような感じになります。・ターゲットとなる部分巡回路もしくは足し合わせた都市群の数は全体都市数の半数未満とする・分離した部分巡回路あるいは都市群がある場合は、その都市群がターゲット・分離した部分巡回路がない場合は、0.5の部分巡回路とそれに付随する1.0の部分経路を足し合わせた都市群をターゲットとする・足し合わせた都市群が他の都市群と接続している箇所が、2未満の場合は2以上の制約をあたえる・足し合わせた都市群が他の都市群と接続している箇所が、3以上の場合は2以下の制約をあたえるこの条件分けで繰り返しループ処理させるコードが以下になります。simplex")最初のst_elim()というのが、ターゲットとなる都市群を探し出す関数です。続くwhileループ内でターゲットとなる都市群に条件(他都市との接続数)に応じて制約を追加していきます。今回の場合1ループで約10秒、48都市を解決するのに約2分かかりました。他のサンプルでは試していないので、まだ完璧ではないかもしれませんが(これから検証)、一応scipyを用いて試してみたという感じです。constraints)を試してみました。前回同様Pulpを使ってLPで解くことにします。よるTSPの解法がある程度理解できたのでアルゴリズムを整理し直してみました。事前準備:まず、ダウンロードした「att48.tsp」をPandasで読み込みます。Y座標」となっており、主にX座標とY座標だけ抜き出します。File)」があるので、その部分を含めないようにします。厳密解データ(解答)も読み込んでおきます。厳密解は巡回順の整数値が並んでいるだけです。[0]draw_plot(T_opt)表示用draw_plot()関数を用意しておき厳密解を表示させます。実際このような巡回路となれば正解というわけです。48都市あるので各都市のインデックス番号を0〜47としておきます。TSPLIBでは1〜48の番号になっていますが、とりあえず今回は0〜47で計算することにします。LPアルゴリズム:48都市あるので、2都市間経路の組み合わせは48*47/2=1128となります。0)は同じものとみなします(順列ではなく組み合わせのほうが数が減るので)。ED)前回同様Pulpを使って最小化する式を用意します。ここで注意すべきことは、i<jという大小関係になっていることと、変数xがバイナリLpBinaryになっているということです。j)]はその組み合わせが経路として採用されれば1で不採用なら0となります。j)]を積算しその値が最小値になるようにLPで解いていきます。この制約の部分を前回よりも簡潔に書き直してみました。j)が二つ必要ということになります。以下のような感じです。二つ採用(都市47において)一旦この制約でsolve()し、varValue=1となる経路(採用された経路)を表示させます。するとこのような複数の部分巡回路が出来上がります。この段階では各辺がバラバラなので部分巡回路ごとにグループ分けする関数を用意して仕分けします。一応、辺数の少ない順に並べ替えています。constraintsを追加し部分巡回路除去を行います。35)]に制約を与えることにしました(前回は複数の部分巡回路に制約を与えていましたが)。部分巡回路は完結した閉路なので、最低でも二本の経路を外につなげるように仕向けます。それが以下の制約になります。j)の組み合わせのうち、iかjのどちらか一方がその部分巡回路に含まれて、他方が含まれていない組み合わせが二つ以上必要という条件にします。j)のどちらか一方が部分巡回路に含まれているかどうかという条件です。することで、部分巡回路内にある最低でも二つの都市が他の部分巡回路の都市へつながろうとします。この条件を与えたあと一旦solve()して結果を見ます。部分巡回路に多少変化が見られるので、何度か繰り返すと部分巡回路がまとまって一つの大きな巡回路になります。while文で最終的にひとつの巡回路になるまで繰り返すコードが以下です。ST[0])今回は15ループ(3秒ほど)で48都市を最適化できました。上記の状態だとまだバラバラの辺の集合なので、以下の処理で最終的な都市順路に並び替えます。0]という順路リストが得られます。逆順になっている場合もあるので必要に応じて反転させます。TSPLIBにあるデータは都市番号1から始まっているので、それぞれに+1すれば同じ値になります。アニメーション:また変化をアニメーション表示するには以下(前半の設定は省略)。gifplt.close()HTML(ani.to_jshtml())以下が作成したアニメーション。赤い部分は制約を与える部分巡回路。Reflect以下が全体のコード。j)の組み合わせが倍の2256通りあるパターンについても試しています。組み合わせが多いのでその分時間がかかりますが、同じアルゴリズムで解くことができます。掲載されているGrbl1.1対応ソフトの一つでもあります。2017年現在では、このページの内容はやや古いものとなっています(Grbl0.9以前のもの)。