経緯度情報を基にした最大被覆問題を解きたい

Question

実現したいこと

数理最適化の一つである最大被覆問題において，198個の経緯度情報(X座標，Y座標)を基にしたポ集合(点)をカバー(円に含む)する入力した値の数だけの集合とそれを中心とした円(需要範囲)が描画される実行結果を導出したい．

発生している問題・分からないこと

経緯度情報を持つ198個の点と入力した値の数だけの点とそれを中心とした円は描画される状態にあるが、ただ入力した値の数だけの点がランダムに配置されてしまう。全ての集合をカバーするような実行結果を算出するためにどうすればいいか知りたい。

該当のソースコード

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import random, math, time
import pulp


#1．入力データ生成：住居の設置点を読み込む
def make_data_ad():
    #1-1．住居の位置情報を格納したcsvファイルを読み込む
    df = pd.read_csv("X_Y_data(rep_housing).csv", encoding="utf-8")
    #1-2．変数[X]に設置する住居のX座標値(経度)を代入
    X = df["X_CODE"]
    #1-3．変数[Y]に設置する住居のY座標値(緯度)を代入
    Y = df["Y_CODE"]
    #1-4．変数[X_ad]・[Y_ad],を値として返す
    return X,Y

#2．入力データ生成：資源回収物拠点の設置可能点を生成
def make_data_si(sigen, upper_x, upper_y):
    random.seed(10)
    #2-1．変数[si_X]に資源回収物拠点の設置可能点を代入
    si_X = [random.uniform(137.689, upper_x) for i in range(sigen)]
    #2-2．変数[si_X]に資源回収物拠点の設置可能点を代入
    si_Y =[random.uniform(34.682, upper_y) for i in range(sigen)]
    #2-3．[si_X]・[si_Y]を値として返す
    return [si_X, si_Y]

#3．住居と資源回収物拠点間の直線距離を測定する
def distance(si):
    ad1,ad2 = make_data_ad()
    for X_ad, Y_ad in zip(ad1, ad2):
        unit_ad = [X_ad, Y_ad]
    #3-1．変数[dx]に住居のX座標値から資源回収物拠点のX座標値を引いたx軸方向の距離を代入
    dx = unit_ad[0] - si[0]
    #3-2．変数[dy]に住居のY座標値から資源回収物拠点のY座標値を引いたy軸方向の距離を代入
    dy = unit_ad[1] - si[1]
    #3-3．x軸方向の距離の二乗とy軸方向の距離の二乗の和の平方根を求める
    return math.sqrt(dx **2 + dy **2)

#4．需要変数(二値変数)行列を生成
def make_matrix(ps, r):
    #5．許容領域内に分析対象(住居)が収まっているかを判断する
    def check(a):
        #5-1．「もし施設間の距離が許容領域の半径以下なら，要素(資源回収物拠点)を1に、それ以外は0にする」
        if a <= r:
           return 1
        else:
           return 0
    #4-1．フィールド上にプロットするためにn次正方行列に整形する用のリスト[unit_list]を作成
    unit_list = []
    #4-2．変数[n,m]に[X]と[Y]の座標値を一つにまとめたものを代入する
    for n, m in zip(ps[0], ps[1]):
    #4-3．変数[unit]に変数[n,m]を代入
        unit = (n, m)
    #4-4．変数[unit_list]に[unit]を追加
        unit_list.append(unit)
    #4‐5．[unit_list]から変数p1,p2にそれぞれX・Y座標値を代入し、変数[distance(p1, p2)]に入っている距離が許容領域の半径に収まるかどうか判断する
    return [check(distance(p1)) for p1 in unit_list]

#6．グラフ上に結果を描写する
def draw(ps, fs, r):
    #6-0．make_data_adを実行
    ad1,ad2 = make_data_ad()
    #6-1．変数[fig]にx軸、y軸における描画領域の大きさを代入する
    plt.figure(figsize=(5, 5))
    #6-2．変数[ax]にグラフの設定を保存するための枠を作成
    ax = plt.axes()
    #6-3．変数[x]にグラフ軸の数値幅と増幅単位の設定を保存
    x = np.arange(137.69, 137.76, 0.1)
    #6-4．住居の位置をグラフ上に描画する
    plt.scatter(ad1, ad2)
    #6-5．for文を用いて資源回収物拠点のX・Y座標値を変数[(i, j)］に，変数[k]に住宅の代表点から資源回収物拠点までの距離を代入する
    for k, (i, j) in zip(fs, zip(ps[0], ps[1])):
        #6-6．もし資源回収物拠点を設置する場合の描画設定を行う
        if k.value() == 1:
           #6-6-1．資源回収物拠点の色は赤色
           plt.scatter(i, j, s=2, color='red')
           #6-6-2．許容領域の描画設定を行う
           c = patches.Circle(xy=(i, j), radius=r, fill=False)
          #6-6-3．許容領域を描画する
           ax.add_patch(c)
        #6-7．資源回収物拠点を設置しない場合の設定を行う
        else:
          #6-7-1．passする
            pass
    #6-8．縦横比が同じ長さになるように調整。
    plt.axis('scaled')
    ax.set_aspect('equal')
    #6-9．グリッド線を表示
    plt.grid()
    #6-10．ⅹ軸のラベルを設定
    ax.set_xlabel('経度', fontname="MS Gothic")
    #6-11．y軸のラベルを設定
    ax.set_ylabel('緯度', fontname="MS Gothic")
    #6-12．x軸の下限値と上限値とその間隔を設定
    ax.set_xticks(np.linspace(137, 138, 10))
    #6-13．y軸の下限値と上限値とその間隔を設定
    ax.set_yticks(np.linspace(34, 35, 10))
    #6-14．グラフを描画
    plt.show()

#7．住居と資源回収物拠点における集合被覆問題を解く
def solver_set(ps, r):
    #7-1．変数[size]に資源回収物拠点のx座標値の個数を代入
    size = len(ps[0])
    #7-2．変数[cs]に「#3需要変数(二値変数)を作成」の関数を代入
    cs = make_matrix(ps, r)
    #7-3．変数[prob]に集合被覆問題を解く命令を代入
    prob = pulp.LpProblem('max-cover', sense = pulp.LpMaximize)
    #7-4．変数[fs]に変数[size]に格納されている値を変数[i]に移行し，この問題で解く際に必要となる関数を作成
    fs = [pulp.LpVariable('f{}'.format(i), lowBound=0, cat = 'Binary') for i in range(size)]
    #7-5．変数[prob]にも目的関数を追加
    prob += pulp.lpSum(fs)  
    #7-6．変数[cs]から変数[c]に制約条件を追加
    for c in cs:
        prob += pulp.lpDot(c, fs) >= fs[c]
    #7-7．変数[s]に時間計測の命令を代入
    s = time.time()
    #7-8．変数[status]に集合被覆問題が解けたかどうかの確認を行う命令を代入
    status = prob.solve()
    #7-9．変数[status]の結果を表示
    print('Status', pulp.LpStatus[status])
    #7-10．設置された資源回収物拠点の数を表示
    print('z 施設数:', prob.objective.value())
    #7-11．集合被覆問題の回答処理に掛かった時間を表示
    print('処理時間', time.time() - s)
    #7-12．関数[draw(ps, fs, r)]を実行
    draw(ps, fs, r)

#許容領域の半径
r = 0.03
#資源回収物拠点の個数
sigen = 10
#乱数上限
upper_x = 137.758
upper_y = 34.778
#実行
solver_set(make_data_si(sigen, upper_x, upper_y), r)

試したこと・調べたこと

• teratailやGoogle等で検索した
• ソースコードを自分なりに変更した
• 知人に聞いた
• その他

上記の詳細・結果

google上で「数理最適化　最大被覆問題　Python」と検索し、最大被覆問題を取り扱うPythonのプログラムを参照し、変更を加えた。
下記が変更したコードの部分
def solver_set(ps, r, sigen):
#7-1．変数[size]に資源回収物拠点のx座標値の個数を代入
size = len(ps[0])
#7-2．変数[cs]に「#3需要変数(二値変数)を作成」の関数を代入
cs = make_matrix(ps, r)
#7-3．変数[prob]に集合被覆問題を解く命令を代入
prob = pulp.LpProblem('max-cover', sense = pulp.LpMaximize)
#7-4．変数[fs]に変数[size]に格納されている値を変数[i]に移行し，この問題で解く際に必要となる関数を作成
fs = [pulp.LpVariable('f{}'.format(i), lowBound=0, cat = 'Binary') for i in range(size)]
xs = [pulp.LpVariable('x{}'.format(i), lowBound=0, cat = 'Binary') for i in range(size)]
prob += pulp.lpSum(xs) == sigen
#7-6．変数[cs]から変数[c]に制約条件を追加
for i in range(size):
prob += pulp.lpDot(cs[i], fs) >= xs[i]
#7-7．変数[s]に時間計測の命令を代入
s = time.time()
#7-8．変数[status]に集合被覆問題が解けたかどうかの確認を行う命令を代入
status = prob.solve()
#7-9．変数[status]の結果を表示
print('Status', pulp.LpStatus[status])
#7-10．設置された資源回収物拠点の数を表示
print('z 施設数:', prob.objective.value())
#7-11．集合被覆問題の回答処理に掛かった時間を表示
print('処理時間', time.time() - s)
#7-12．関数[draw(ps, fs, r)]を実行
draw(ps, fs, r)

補足

PythonのVersionは2023.22.1であり、ソルバーにはpulpを使用しています。

Comments

[ikedas]

参考までに聞きたいのですが、この問題を解こうと思った動機は何かありますか。

[退会済みユーザー]

卒業研究で施設の最適配置に関連するテーマに設定していて、その最適配置を導出する方法として最大被覆問題を選んだのが動機です。

[ikedas]

そうですか。
大学には学術研究と教育の機能がありますから、卒業研究をサポートすることもその機能の一部なんですよね。
卒業研究も含め大学でやる学習・研究の相談は、大学で然るべき対応をしているはずですから、そちらを頼ってください。
たとえば研究室の教員や院生などの先輩に相談するとか、同級生とディスカッションするとかです。そう言った支援を受けて当然なだけのお金をあなたは大学に払っているはずです。

大学で学んでいる方なら、当サイトに質問しにこないでください。そうではなく、大学で学んだ成果を活かして近い将来に当サイトで質問に回答できるようになってください。

[fana]

質問するにしても，もう少し話を整理した方が良いのでは．
例えば現状が

> 入力した値の数だけの点がランダムに配置されてしまう

という結果になるのは何故なのか？　という話とか，

あとは，この問題の何が難しい箇所なのか？　みたいな話とか．
（緯度と経度　という話だとこれを単純な２次元平面上の点としては扱えない，みたいな部分が難しいのかな？）

[ikedas]

いや、質問の内容が云々以前に、卒研の「正解」を質問しに来ているのがそもそもおかしい。卒業研究がなんなのかわかってないんじゃないですか。

[bsdfan]

コード全体をみていませんが、何を最大化したいのか、そのときの制約がなにか、あたりがちゃんと設定できていないように見えます。

[fana]

マルチポストっぽい：
https://qiita.com/yamamonxd0828/questions/687e160af188a4851ce4

[ikedas]

> 何を最大化したいのか
単に「最大被覆問題」という問題の類型があるというだけでしょう。ググると、「COVID-19のワクチン接種施設 (複数) を自治体内のどこに配置すれば、一定の距離内で利用可能な世帯の数を最大化できるか」みたいな例が出てきます。

質問の日本語も突っ込みどころ満点ですが、大学の卒業研究を丸投げしてくるとはteratail史上最大級の丸投げ質問ではないでしょうか。

[bsdfan]

今のコードではfsの合計を最大化しようとしていますが、おそらく最大化したいものはそれじゃないはず。それに続く制約条件も何を意味しているのかちょっと不明です。
(さすがに、このコードを完成させるだけの内容が卒業研究のすべてだとは思いたくないのでコメントしています。が、マルチポストっぽいのでこれ以上は控えます)

[ikedas]

元にしたと見られるネット上のコードを見ると、距離の合計を最大化するという意図で正しいようです (距離を直線で測っているとかみんな砂漠か草原にテントでも張って住んでるのかという単純化されたモデルなので、あまり深く考えてもしかたがないのかと。資源回収拠点の稼働率の指標としてはありかもしれない)。あと、制約条件のコードはなくてもいいみたいな説明もある。
マルチポストのもう一方の投稿を見ると「いつも同じ結果しか出ない」というのが質問者さんが抱えている問題のようですが、そうなるように改変してあるコードなのだから当たり前ですね。

現状の延長に、質問者さんの得たい結果が出るコードはないと思います。というかそれを考えるのが研究でしょ。

Answer 1

これは回答ではありませんが、後述する理由からあえて投稿しています。

---

全ての集合をカバーするような実行結果を算出するためにどうすればいいか知りたい。

それは質問者さんの卒業研究のテーマそのものではないですか。それを考えるのはあなたです。

ところで、ご質問の問題は卒業研究の対象なのですよね。　あなたが行うのは「研究」であって、その成果をもってあなたは卒業資格を得るのですよね。

「研究」とは、「物事を深く考えたり、詳しく調べたりして、真理、理論、事実などを明らかにすること。」(日本国語大辞典) です。特に学術研究の場合、調査や実験の結果を解釈し考察した結果、今までになかった独自の知見をもたらすような営為を指します。

「研究」は今までの勉強と違って、正解があることが最初からわかっていてそれを当てれば合格、なのではありません。あなた自身が、他の人がまだ発見していない新たな真理、理論、事実を提示できなければ「研究」とは言えません。

逆にいうと、他の人の発見を自分の研究として発表したらそれは「盗用」となります。

どの大学・学術機関にも「研究者倫理規程」のようなものがあります。あなたも「研究」をするのですから、当然このような倫理規程の対象となります。その規定の中で、盗用は最も厳しく禁じられている行為の一つです。

さて、今後この質問に回答がついたとします。そしてあなたがその回答の内容を使って卒業研究の論文を書いたとします。論文が受理された後に学校関係者の誰かがこちらの回答を見つけ、内容が一致すると判断されれば、それは「盗用」と認定される可能性があります。——たとえ回答を投稿した人が「自分の回答は自由に使ってもらって構わない」と思っていたとしても関係ありません。内容が一致するのなら客観的に見てそれは「盗用」です。

盗用と認定されれば、基本的にはその論文は撤回されます。卒業研究という場合、大抵はその研究の成果が学位授与をはじめとした卒業資格の条件になりますから、卒業も取り消されるかもしれないということです。

ですので、もしもこの質問に回答がついたら、質問者さんはその内容を自分の卒業研究に取り入れないように細心の注意を払うべきでしょう。このことを言いたいためにこの投稿をしました。

---

当サイトでも、学校の授業での課題を質問しにくる学生が残念ながらなくならないのですが (勉強は学校でやるべきですし、個人的にはカンニングとみなされるおそれが高いと考えるので、本人自身のためにもやめてほしいです)、研究となるとそれより格段に慎重さが要求されます。これから「研究」を始める以上、あなたも研究者の仲間入りをするのであり、研究倫理を守らなければなりません。

研究倫理について詳しいことは、あなたが指導を受けている教員の方に聞いてください。

あと、もう一つ指摘しておきます。文章でもソースコードでも、他人の書いたものを参考にしたのなら出所を明記するのが約束です。学校でそう指導しているはずですよ。

Comments

[fana]

まぁ訊き方ですよね．

質問内容が研究としての 独自性/新規性/etc 的な主張部分そのものに関わるようだとＮＧでしょうけど，
そうでない（ 例えば，結局のところ質問内容が pulp とかいう何かの使い方でしかない，みたいな）場合であれば問題ない，みたいな．

本件の問題を pulp とかいう何かを用いて解くことが研究と言えるのか？　すなわちそこに　独自性/新規性/etc　が存在するのか否か，私には分かりませんが……
本件で示されている問題は「最大被覆問題」とかいうジャンルの単なるシンプルな例題の１つでしかなくて，研究内容はもっと先のところにあるのではないだろうか…？　という気もします．

いずれにせよ，教員etc に質問することを考えるべきであろうと思いますが．

[ikedas]

> 研究内容はもっと先のところにあるのではないだろうか
そうだと思いますよ。「PuLPを使って問題を解きました」だけでは研究にならないので。では親身になって相談に乗ってやればいいのかというと、それで味を占められては困る。
ネットで得た半端な知識をもとに適当に書き飛ばした卒業研究なんて、通してもらえるわけがない。だからさっさと学校に戻ってほしいです。というか、まず学校で勉強のしかたを勉強してくれと思う。
でも、「なんか叱られたみたいだけどアカウント消して行方をくらませばいいや」と思ってるんなら、見通しは暗いかもしれませんね。まあ私の知ったことではないですが。