長方形に円の詰め込み最適化について

Question

前提・実現したいこと

過去にいくつかの質問があるようですが、解決までに至りません。
やりたいことは、以下の図９です。
http://www.orsj.or.jp/archive2/or63-12/or63_12_762.pdf
そのための前哨戦として以下の最適化問題をやろうとしています。線形ではなく、非線形でないと,,,,という話もありますが、

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

エラーメッセージ

該当のソースコード

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

#n:並べたい円の数,r:半径,W:横幅,UB高さ
n=15
r=0.2
W=2.2
UB=2.6

x=[LpVariable('x%d' %i, lowBound=0) for i in range(n)]
y=[LpVariable('y%d' %j, lowBound=0) for j 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
m += H
for i in range(n):
    for j in range(n):
        m += x[i]+r*2 <= x[j]+W*(1-u[i][j])
        m += y[i]+r*2 <= 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*2
    m += y[i] <= H-r*2
    m += r <= x[i]
    m += r <= y[i]
      
m.solve()

#
#for i in range(n):
#    print("x[i] value:", x[i].value())
#    print("y[i] value:", y[i].value()) 


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

for i in range(n):
      circle_n = plt.Circle((x[i].value(),y[i].value()),r)
      ax.add_patch(circle_n)
        
#グラフのタイトル
plt.title("Sectional View")

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

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

fig.savefig("test1.png")
コード

試したこと

円の作画は中心点の設定、長方形のゼロ点は左下であるので、r*2で半径ではなく直径で円が重ならないようにしています。
この結果は、以下のようになります。ところが、n=15とするとプログラムが終了しません。5列3段になると思っているのですが（最終的には半径の異なる円を最適に詰め込みたいのですが）。

Comments

[y_waiwai]

このままではコードが読めないので、質問を編集し、<code>ボタンを押し、出てくる’’’の枠の中にコードを貼り付けてください

Answer 1

過去の類似質問に回答していた者です。
https://teratail.com/questions/328071#reply-454147
https://teratail.com/questions/328308#reply-454344

・n=15とするとプログラムが終了しないのは、単に計算時間が掛かっているだけだと思われます。
・以下、回答例を提示します。円の重複判定には上下左右の関係しか見ていないため、隙間が生まれます。

from pulp import *
import matplotlib.pyplot as plt

W = 20
UB = 30

n = 6
r = [1, 2, 3, 4, 5, 6]

H = LpVariable("H")
x = [LpVariable(f'x{i}', lowBound=0) for i in range(n)]
y = [LpVariable(f'y{i}', lowBound=0) for i in range(n)]
u = [[LpVariable(f'u{i}{j}', cat=LpBinary) for j in range(n)] for i in range(n)]
v = [[LpVariable(f'v{i}{j}', cat=LpBinary) for j in range(n)] for i in range(n)]

m = LpProblem(sense=LpMinimize)
m += H

for i in range(n):
    for j in range(n):
        m += x[i] + r[i] <= x[j] - r[j] + W * (1 - u[i][j])
        m += y[i] + r[i] <= y[j] - r[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[i]
    m += y[i] <= H - r[i]
    m += r[i] <= x[i]
    m += r[i] <= y[i]

m.solve()

### 描画 ###

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

for i in range(n):
    print(f'(x, y) = ({x[i].value()}, {y[i].value()})   r = {r[i]}')
    circ = plt.Circle((x[i].value() ,y[i].value()), r[i])
    ax.add_patch(circ)

plt.show()

(x, y) = (7.0, 21.0)   r = 1
(x, y) = (2.0, 18.0)   r = 2
(x, y) = (3.0, 13.0)   r = 3
(x, y) = (16.0, 4.0)   r = 4
(x, y) = (5.0, 5.0)   r = 5
(x, y) = (14.0, 16.0)   r = 6

--- 追記（コメントへの回答） ---

枠で囲んでみると、これ以上Hを最小化できないことがわかる。

Comments

[KAZENOMACHI]

ありがとうございました。ｒをすべて２にして、n=10にすると、なかなか終了しません。n=7までは同じくらいの速度で終了するのですが。また、重複は上下で判定していますが、r=1は一番下でも配置できそうですし、r=6はもう少し下でも配置できそうですが、、、。完全に詰め込むことはできるのでしょうか？また、0.1の余裕を持たせ詰め込むことはできるのでしょうか？

[nanoseeing]

>> ｒをすべて２にして、n=10にすると、なかなか終了しません。
計算が終了するまで待つしかありません。

>> r=1は一番下でも配置できそうですし、r=6はもう少し下でも配置できそう。
回答欄に円を枠で囲った画像を追記しました。Hは最小化されています。

>> 0.1の余裕を持たせ詰め込むことはできるのでしょうか？
元の半径から0.05を足せば良いと思います。

[KAZENOMACHI]

参考にしたpdfのように正方形の枠線での詰込み配置ではなく、円周の重複を避けた詰込みができると思ったのですが、そういう意味ではなかったのですね。これは長方形の断面積を計算して、円の断面積を最大にして、かつ、Hを最小にするような最適化にするとよいのでしょうか？

[KAZENOMACHI]

n = 16
r = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
にしても、同じ円が5列x3段+1個になると思ったのですが、以下のエラーとなります。どうしてなのか？正方形の枠でも十分に入ると思うのですが。

File "C:\Users\05896\Desktop\ダクト配列NEW.py", line 34, in <module>
m.solve()

File "C:\Users\05896\Anaconda3\lib\site-packages\pulp\pulp.py", line 1737, in solve
status = solver.actualSolve(self, **kwargs)

File "C:\Users\05896\Anaconda3\lib\site-packages\pulp\apis\coin_api.py", line 101, in actualSolve
return self.solve_CBC(lp, **kwargs)

File "C:\Users\05896\Anaconda3\lib\site-packages\pulp\apis\coin_api.py", line 159, in solve_CBC
raise PulpSolverError("Pulp: Error while executing "+self.path)

PulpSolverError: Pulp: Error while executing C:\Users\05896\Anaconda3\lib\site-packages\pulp\apis..\solverdir\cbc\win\64\cbc.exe

[nanoseeing]

>> n = 16でエラーが起きる

変数名の重複が原因のようです。u{i}_{j} のように変数名を変更してください。
https://ja.stackoverflow.com/questions/55563/pulp%E3%81%AE%E3%82%A8%E3%83%A9%E3%83%BC%E3%81%AB%E3%81%A4%E3%81%84%E3%81%A6

>> 参考にしたpdfのように正方形の枠線での詰込み配置ではなく、円周の重複を避けた詰込みができると思ったのですが、そういう意味ではなかったのですね。

そのとおりです。結局のところ、正方形の詰込み問題を解いているに過ぎません。

>> これは長方形の断面積を計算して、円の断面積を最大にして、かつ、Hを最小にするような最適化にするとよいのでしょうか？

断面積が何を指すのかわかりませんが、いずれにしろ線形計画問題に落とし込めないと、pulpライブラリでは解けないかと思います。

私の思いついた解法としては、以下があります。
A. 非線形計画問題に対応したソルバーを使用する。
B. 円を正多角形で近似して、線形計画問題に落とし込む。厳密解は得られない。
C. 参考pdfのように、BL法による貪欲的な探索をする。厳密解は得られない。

私自身は数理最適化に詳しいわけではないので、参考までに。

[KAZENOMACHI]

確かに正方形を詰め込んだことと同じですね。正方形の面積：その中にはいる円の面積は78.5%なので、余裕率として21.5%は見ておくという考えをすると実運用ではいいのかもしれませんが、ぎりぎりまで詰め込みたいときの可能性を知りたい場合には、使えませんね。非常に勉強になりました。もう一度、数理最適化から調べて再挑戦したいと思います。