⛳

Metaheuristics.jl を使ってジオフェンスを計算する (数理最適化 Advent Calendar 2023)

2023/12/15に公開

数理最適化

ジオフェンス

tech

本記事は数理最適化 Advent Calendar 2023 の15日目の記事です。

Juliaの Metaheuristics.jl パッケージとGPSデータを利用して、ジオフェンスの構築をやってみた結果を紹介します。

はじめに

まずは聞いたことがない人もいるかもしれないので、タイトルにある「ジオフェンス」というものと、この記事でやったことについて概要を説明します。

ジオフェンスとは

ジオフェンスというのは地理上の仮想的な境界線内のことです。

似たような概念はいろいろあると思いますが、ある種のテリトリーみたいなものです。

例えば市街地にジオフェンスを考え

ジオフェンスを跨いだ回数で重量課金を設計する
ジオフェンスを跨ぐことを仮想的に禁止して罰金を課す

といった、ちょっとしたサービスの料金体系などをつくることができます。

他にもお店や観光スポットにジオフェンスを考え、ジオフェンスを跨いだ端末にPush通知でなにかの情報を送信するなどのアプローチが考えられます。

きっかけ

今回の記事のネタを考えていたとき、こちらの発表を聞く機会がありました（ACMの論文のページはこちらです）。

発表が面白かったので、これを再現しようと思いました。

論文ではジオフェンス最適化については遺伝的アルゴリズムGAを使っていたのですが、この記事ではGAの実装まではやっておらず、タイトルのパッケージ (Metaheuristics.jl くん) に計算を丸投げしました。

やったこと

具体的にやったことは以下のとおりです。

論文と同じデータ (GeoLife GPSデータ) を前処理 (座標変換) し、ジオフェンスなしで似たような図を再現する。
最適化 (論文中の fitness 関数) を定義し、最適化してもらって、結果を見比べる。

データ処理

GPSデータの変換

GeoLife GPSデータは Microsoft Research が公開している有名なGPSデータです。

位置情報データや軌跡データを扱う人は必ず使ったことがあると思います。GPS的な情報のログデータがファイルに保存されていて、読み込むと (lat, lon) の座標値が得られます。もしデータを処理したい方がいたら、このあたりのページを見ていただくと自分でプログラムが書けると思います。

lat,lon,flag,alt,date,datestr,timestr,prefix
39.984071,116.320147,0,492,39744.2914236111,2008-10-23,06:59:39.0,20081023065939
39.98414,116.320122,0,492,39744.2914814815,2008-10-23,06:59:44.0,20081023065939
39.984167,116.320114,0,492,39744.2915393518,2008-10-23,06:59:49.0,20081023065939

座標データの付与

データを読み込むまではすぐに終わります。次に論文を確認すると、座標変換を行っています。

All coordinate values were converted from a geographic coordinate system (EPSG:4326) to a Cartesian coordinate system (EPSG:32650) by GeoPandas [19], a Python library.

これは何かというと、地球をモデル化するときの球の設定や、もしくは局所的に平面を想定する場合の基準点などでバリエーションがあるため、処理をした場合には処理を合わせておく必要があるということです。詳しくはこちらの本を読んでください。

https://amzn.to/3t5j7jX

ここで躓きポイントが1つあるのですが、座標変換の使いやすいライブラリが Julia ですぐに手に入りませんでした (おそらく真面目に探せばあると思いますが…)。仕方ないので座標変換だけ python の pyproj を使って処理し、先に全部のGPSデータを Trajectory として置いておくことにしました。

処理したCSVは、このように(x, y)の値を持っています。

lat,lon,flag,alt,date,datestr,timestr,prefix,x,y
39.906488,116.388619,0,107,39747.5890277778,2008-10-26,14:08:12.0,20081026134407,447741.40541719145,4417557.259469705
39.906493,116.38867,0,107,39747.5890856481,2008-10-26,14:08:17.0,20081026134407,447745.7685376378,4417557.784575468
39.906478,116.38872,0,106,39747.5891435185,2008-10-26,14:08:22.0,20081026134407,447750.03098671953,4417556.090477691

可視化

論文に登場する2つのPOIのうち、片方の座標を計算したのもを例題に使います。

Workers’ Cultural Palace: A temple located next to the Gate of Heavenly Peace. The position of the POI is set to latitude = 39.90874, longitude = 116.39368.

# WCP (Worker's Cultural Palace)
const WCP_x = 448175.70093339786
const WCP_y = 4417804.258708443

あまり遠いデータを全部扱う必要はないため、論文に従って (?)、(x, y) の座標で500離れたデータを無視するフィルタリングをかけておきました。

ここまでくると、論文のジオフェンス計算の入力に使われたデータが模倣できているはずです。見てみます。論文の図（ジオフェンス付き）は、上記論文の Fig.2(i) から抜粋させて引用させていただきました。私の再現の図の方は地図背景 (OpenStreetMap) がないことと、ユーザのprefixで色を変えていることが違いますが、概ね同じデータを用意できたと思います。

論文の図 (ジオフェンス付き)	処理したデータとWCP

ここから、図 (Fig.2(i)) に出てくるジオフェンス (赤色の領域) を求めるのが最適化パートです。

最適化

最適化パートを説明します（Advent Calendarはここから書けという説もあります）。

ジオフェンス計算の気持ち

設計者の気持ち (?) を想像します。

GPSデータの人にあるPOI (ここではWCPのこと) について何か宣伝をしたいと思います。すると、あるPOIの情報をフェンス内のユーザに届けたいです。しかし、あまり遠くの人に情報を送っても効果が薄そうです。そのため

対象とするPOIに関係ある場所の近くに起く
データの多くの人にささる (多くの人が円の内側を通る)

を満たしたい気持ちになります。これを関数に落とせば、何らかの形で最適化できるということになります。

定式化 (論文より)

論文ではジオフェンスを $i=(x, y, r)$ として (つまり中心 $(x, y)$ で半径が $r$ の部分が仮想領域) 表現します。ここで $i$ を評価するための fitness 関数 $F(i) := f(i) + g(i)$ として、以下の2つの項を導入しています。

上に書いた通り、 $f(i)$ の方でPOIとの関連度を、 $g(i)$ の方で入力データに対するカバー率を持ちました。

式の図は $f(i)$ の方です:

f(i) := \frac{1}{2}\{\mathit{dist}(POI, i.xy) + i.r + |\mathit{dist}(POI, i.xy) - i.r|\}

こちらの式の図は $g(i)$ の方です:

g(i) := \mu \max(0, \mathit{cr}_\mathit{limit} - cr(i))

カバー率 $cr(i)$ ですが、入力データのうち、今作ったジオフェンス $i$ をかすっている割合とします (たぶん、そんな感じの量です)。

目的関数の実装

上で説明したものを実装しました。

const WCP_x = 448175.70093339786
const WCP_y = 4417804.258708443

const μ = 2000
const crlimit = 0.5

function f(x, y, r)
    term1 = sqrt((x - WCP_x) ^ 2 + (y - WCP_y) ^2)
    term2 = r
    term3 = abs(term1 - r)
    return (term1 + term2 + term3) / 2
end

# カバー率と思っている量
function cr(x, y, r, data)
    counter = 0
    for df in data
        if sum(((df1.x .- x) .^ 2 .+ (df1.y .- y) .^ 2) .<= r^2) > 0
            counter += 1
        end
    end
    counter / length(data)
end

# ここに元データ (GPSのDataFrames) が入っている。
# お行儀がまったく良くなさそうな実装!!!!
list_dfs = load("list_dfs.jld2")["data"];
function g(x, y, r; data=list_dfs)
    μ * max(0, crlimit - cr(x, y, r, data))
end


# i = (x, y, r) を3次元ベクトルvとして思う
fitness(v) = f(v[1], v[2], v[3]) + g(v[1], v[2], v[3])

Metaheuristics.jl による最適化

あとはfitnessを最適化すると問題が解けそうです。上で定義した fitness をライブラリに丸投げして終了です！ (本当はGAを実装しようと思ったのですが間に合いませんでした)。

結果だけ述べると、最適化パートの実装はこんな感じになっています。楽ちんですね。

using Metaheuristics

# (x, y, r) のなんとなくの範囲を設定
bounds = [WCP_x - 1000 WCP_x + 1000; WCP_y - 1000 WCP_y + 1000; 100 500]

# 何も設定していないECA()でとりあえずfitnessを最適化する
result = optimize(fitness, bounds, ECA())

# 結果を取得する
X, Y, R = minimizer(result)

計算したジオフェンスの可視化

最後に可視化します。

using Colors
using Plots
gr()

# Plots.jl で円を描く関数
function circleShape(h, k, r)
    theta = LinRange(0, 2*pi, 500)
    h .+ r * sin.(theta), k .+ r * cos.(theta)
end

# 描画
fig = plot(size=(500, 500))

colors = distinguishable_colors(length(list_dfs))

for (i, df) in enumerate(list_dfs)
    scatter!(fig, df.x, df.y, ms=2, c=colors[i], alpha=0.25, label=nothing)
end

# ジオフェンス
plot!(fig, [X], [Y], color=:red, marker=:+, markersize=10, label=nothing)
plot!(fig, circleShape(X, Y, R), color=:black, seriestype=:shape, fillalpha=0.2, label=nothing, aspect_ratio=:equal)

scatter!(fig, [WCP_x], [WCP_y], marker=:x, c=:red, label=nothing)

論文の図と、今回可視化した図を比較してみます。論文の図は上で例に使ったFig.2(i)のことです。

論文の図	今回計算したジオフェンス

似たような結果が得られました！！！！

(1つしか再現の図がないのは、他のパラメータのジオフェンスがそこまできれいに見つけられなくて、そっ閉じしたからです)。

まとめ

今年は Metaheuristics.jl を使ってジオフェンス計算の再現をやってみました。

ソースコードとデータはそのうち真面目に整理しておきます (来年までには…)。