Pythonによる円（2次元）の詰め込み問題

Question

前提・実現したいこと

pythonのmatpoltlibを用いて、容器の中に円をしたから詰め込み、高さが最も低くなるように詰め込んでいきたいと考えております。

プログラムはエラーが起こらずにコンパイルできましたが、円の数を変更しても1つしか出力されずに困っております。
また、その円も下から詰め込まれずに宙に浮いております。
何かアドバイスを頂けると助かります。

下記の文献を参考に取り組んでおります。
http://www.orsj.or.jp/archive2/or63-12/or63_12_762.pd

発生している問題・エラーメッセージ

エラーメッセージはでません。

エラーメッセージ

該当のソースコード

from pulp import *
m = LpProblem(sense=LpMinimize)
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

#n：個数
#r=円の半径
#W=容器の幅
#UB=容器の高さ
n=10
r=0.1
W=2.2
UB=2.6

fig = plt.figure()
ax = fig.add_subplot(111,aspect='equal')
ax.set_xlim([0,W])
ax.set_ylim([0,UB])

x=[LpVariable("x%d" %i, lowBound=0) for i in range(n)]
y=[LpVariable("y%d" %i, lowBound=0) for i in range(n)]

u=[[LpVariable("u%d%d" %(i,j),cat=LpBinary) for j in range(n)] for i in range(n)]
v=[[LpVariable("v%d%d" %(i,j),cat=LpBinary) for j in range(n)] for i in range(n)]

H=LpVariable("H")

m += H
for i in range(n):
	for j in range(n):
		m += x[i]+r <= x[j]+W*(1-u[i][j])
		m += y[i]+r <= y[j]+UB*(1-v[i][j])
		if i < j:
			m += u[i][j]+u[j][i]+v[i][j]+v[j][i] >= 1 

for i in range(n):
	m += x[i] <= W-r
	m += y[i] <= H-r

for x in range(n):
	for y in range(n):
		x=r*n
		y=r*n
		circle_n=plt.Circle((x,y),r)
		ax.add_patch(circle_n)

#グラフのタイトル
plt.title("Sectional View")

#x軸ラベル
plt.xlabel("Hole diameter")

#y軸ラベル
plt.ylabel("depth")

fig.savefig("test1.png")

試したこと

パラメーターを少しずつ変更しながら原因を探しておりますが、まだわかっておりません。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Answer 1

原因がわからない場合は、printを使ったデバッグを覚えましょう。

円が1つしか表示されない →　座標(x,y)や半径rがおかしいのでは？
と考え、x,y,rなどの変数を表示してみましょう。

下記のようにprintを追記してみました。

from pulp import *
m = LpProblem(sense=LpMinimize)
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

#n：個数
#r=円の半径
#W=容器の幅
#UB=容器の高さ
n=10
r=0.1
W=2.2
UB=2.6

fig = plt.figure()
ax = fig.add_subplot(111,aspect='equal')
ax.set_xlim([0,W])
ax.set_ylim([0,UB])

x=[LpVariable("x%d" %i, lowBound=0) for i in range(n)]
y=[LpVariable("y%d" %i, lowBound=0) for i in range(n)]

u=[[LpVariable("u%d%d" %(i,j),cat=LpBinary) for j in range(n)] for i in range(n)]
v=[[LpVariable("v%d%d" %(i,j),cat=LpBinary) for j in range(n)] for i in range(n)]

H=LpVariable("H")

m += H
for i in range(n):
    for j in range(n):
        m += x[i]+r <= x[j]+W*(1-u[i][j])
        m += y[i]+r <= y[j]+UB*(1-v[i][j])
        if i < j:
            m += u[i][j]+u[j][i]+v[i][j]+v[j][i] >= 1 

for i in range(n):
    m += x[i] <= W-r
    m += y[i] <= H-r

### 追記 ###
print("x         :", x)
print("x type    :", type(x))
print("x[0]      :", x[0])
print("x[0] type :", type(x[0]))
print("x[0] value:", x[0].value())
###########

for x in range(n):
    for y in range(n):
        x=r*n
        y=r*n
        circle_n=plt.Circle((x,y),r)
        ax.add_patch(circle_n)

#グラフのタイトル
plt.title("Sectional View")

#x軸ラベル
plt.xlabel("Hole diameter")

#y軸ラベル
plt.ylabel("depth")

fig.savefig("test1.png")

出力結果

x         : [x0, x1, x2, x3, x4, x5, x6, x7, x8, x9]
x type    : <class 'list'>
x[0]      : x0
x[0] type : <class 'pulp.pulp.LpVariable'>
x[0] value: None

これを見れば、xはリストで、
xの中身は「LpVariable」であることがわかります。

また、value()を使ってxの中身を取り出すと、
「None」になってしまっていることがわかります。

なぜ、「None」になっているか考えれば、
m.solve()を書き忘れていることにも気付けるはずです。

Comments

[yoshih]

printでデバックするテクニックは非常に勉強になりました。
数式を解き終わったと宣言してあげないといけないという理解でよろしいでしょうか。

また、m.solve()を追加してみましたが変わらず1つしか円が生成されませんでした。
forの変数が間違っているのでしょうか。

お忙しいところ恐縮ですが何卒宜しくお願い致します。

[nanoseeing]

m.solve()を実行することで、制約条件に従って問題を解いてくれます。

forの変数が間違っているかというご質問に関しては、その通りなのですが、ご自身で確認してみましたか？　forの文の中でx、yをprintしてみて想定どおりの数値になっているか確認してみてください。

いきなり複雑な問題に挑戦するのは難しいです。とりあえず、シンプルな問題を考え、想定どおりの答えになっているか確認してみてください。

例えば、幅が2の長方形に、半径1の円を2つ、詰め込むことを考えます。答えは、中心座標が(1, 1)と(1, 3)になるはずですよね。

[nanoseeing]

なお、pulpは線形計画問題を解くためのライブラリだと思うのですが、円の詰め込みは、線形な制約条件で表すことができるのでしょうか。円の重複判定には中心間の距離が必要で、2乗が出てきませんか？

そのあたりの数学に関しては詳しくないので間違っていたら申し訳ないのですが。

[yoshih]

printで値を見て、僕も円の詰め込みは線形な制約条件では表せないと思いました。そもそもuとvは上下左右の位置しかみておらず、最密充填に詰め込むことができないので、もう一度作り直す必要がありそうです。。。
重複判定は使用せずに、円と接する距離でプログラムしてみようと考えてます。
次はprintを活かしながら実行していきたいと思います。また質問することがあると思いますので、そのときは宜しくお願い致します。

[KAZENOMACHI]

https://teratail.com/questions/328071
で横から参加しました。同じようなことを考えている方がいるのは勇気づけられます。長方形に大きさの違う円をHを超えないように、最大数配列する最適化、もちろん、配列できなかった円の番号は把握しておきたいのですが。
http://www.orsj.or.jp/archive2/or63-12/or63_12_762.pdf
の図９のような感じです。