一個半徑為R的空心球，如果用半徑為r的小球體填充，至少需要多少個小球體？

01-11

1.R/r為正整數，2.請給出詳細的數學建模過程。

注意：為方便起見，文中示意圖全部為二維情形——使用半徑為 $r$ 的圓盤填充半徑為 $R$ 的圓，題中所求小球填充空心球的三維情形示意圖請自行類比和腦補。

--------------------------------------------------------------
常規的建模思路：
一、問題分析

要求使用數量儘可能少的半徑為r的小球體來填充半徑為R的空心球，本文將問題理解為：

當空心球內的各個小球不發生相對位移（相對於空心球）的前提下，若無論如何也不能再塞入一個小球，此時空心球內的小球個數 $N$ (正整數)最少為多少？

這裡的「填充」一詞需要加以解釋和限制，由於是尋求小球的最小數量，也就是使用最少數量的小球有效地佔據空心球內部空間，而「填充」一詞是否需要滿足各小球之間、一部分小球與空心球內表面（記為 $S_{in}$ ）之間存在接觸的限制是模糊的（經過查詢現代漢語字典，我不及格的語文水平仍然沒能幫我確定下來），下面分兩種情況分別考慮：

（它們的區別在於對「填充」一詞的限制不同）

Q1：小球之間、小球與空心球內表面 $S_{in}$ 之間沒有必須接觸的限制，下文稱作「離散情形」，如下圖1：

（圖1：大圓半徑為 $R$ ，小圓盤半徑為 $r$ ，深色小圓已填入，淺色為試圖填入的小圓）

此示意圖中 $2r< R< 3r$ ，易知空心球內最少填充1個小球，可保證第二個淺色的小球無法再塞入。

這種情形與生活中常用的「填充」一詞含義略有不同，但是為了尋求填入小球的最少數目，在不要求小球與空心球相互接觸的前提下，這種「填充」會以最大的效率佔據空心球的內部空間。

Q2：任一小球至少與其他小球存在一個接觸點、且空心球內表面至少與某一個小球存在接觸點，從而保證空心球內填入的小球在空間上的物質分布具有連續性，以滿足「填充」的要求，下文稱作「連續情形」，如下圖2：

（圖2：圓內最少填充2個小圓，第三個淺色的小圓無法再塞入）

下面我們會發現，這兩類情形的求解實際上只相差一個限制條件，且連續情形包含在離散情形中，前者為後者的一種極限。

先從離散情形進行建模。

二、模型假設

1.空心球的內徑為 $R$ ；

2.半徑為 $r$ 的小球體均為剛體。

三、主要符號假設

$N$ ——空心球內填入的小球個數；

$i$ ——空心球內小球的編號( $i=1,2,...,N$ )；

$Omega _0$ ——空心球內全部小球編號的集合；

$S_{in}$ ——空心球的內表面。

四、建立模型

1.目標分析

填充空心球的小球數最少，即：

$min N$

2.約束分析

我們可以以空心球的球心作為原點，建立空間直角坐標系，則各小球的球心位置可記為位矢 $vec{r_i}$ ( $i=1,2,...,N$ )。

（1）顯然，各小球必須位於空心球內，則有：

$left| vec{r_i} ight| leq R-r$

（2）任意兩個小球在空間上不存在重疊的區域，即一個球不能進入另一個球內部，也不能發生形變。若記空心球內全部小球的編號為集合 $Omega _0={1,2,...,N}$ ，則有:

$left| vec{r_i}-vec{r_j} ight| geq 2r~~(i e j,~i,jinOmega _0 )$

（2）根據前面的分析，空心球被「填充」完畢意味著無法再塞入一個半徑為 $r$ 的小球，但此時仍然存在許多縫隙，每一個縫隙都意味著可以再填入一個半徑較小的小球，記這些小球中半徑最大者為 $r_{max}$ 。那麼，我們保證 $N$ 個小球的確可以「填充」這個空心球的標準就是：

$r_{max}< r$

（圖3：四個黑色圓盤示意實際填入的小球，6個藍色圓盤示意小球之間和小球與 $S_{in}$ 之間產生的各個縫隙可能填入的最大小球）

注意：這六個藍色小球為各個縫隙中可能填入的最大小球，因而存在如球1和2的重疊情形，這並不矛盾。考慮它們的目的在於找到其中的最大者，在圖中也就是球6。只要保證球6的半徑小於 $r$ ，就意味著再也無法塞入一個半徑為 $r$ 的小球，表明此時使用4個小球即可填充該空心球。

因此，題目可以轉化為一個非線性的單目標最小規劃：

$min ~N~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\ s.t.left{egin{matrix} left | vec{r_i} ight |leq R-r~(iin Omega _0)~~~~~~~~~~~~~~\ left| vec{r_i}-vec{r_j} ight| geq 2r~~(i e j,~i,jinOmega _0 )\ r_{max}< r~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ end{matrix} ight.$

於是，此時的關鍵在於找到一個計算 $r_{max}$ 的演算法。

3.計算 $r_{max}$ 的演算法

觀察前面的圖3容易發現，我們需要找到所有可能存在的藍色小球，而它們又分為兩類：

（1）一類由三個黑色小球與空心球內表面 $S_{in}$ 圍成，如圖3中的藍色球3、4、5、6 所示，這裡稱為邊界縫隙球；

（2）另一類由四個黑色小球圍成，如圖3中的藍色球1、2 所示，這裡稱為內部縫隙球。

先考慮內部縫隙球。

Ⅰ、計算由等徑球構成的某縫隙內可塞入的最大球半徑
事實上，我們從圖3中不難發現，以球1為例，在二維情形中，尋求該縫隙（球Ⅰ、Ⅱ、Ⅲ圍起來的縫隙）中的可填入最大球問題實質是一個三點求定圓問題。

（圖4：圓1與定圓Ⅰ、Ⅱ、Ⅲ均相切時達到最大）

二維情形下，為了求解圓Ⅰ、Ⅱ、Ⅲ圍起來的縫隙中所能填入的最大圓，我們先設其圓心的位矢為 $vec{r}~$ ，顯然當該圓與三球均相切時達到最大，半徑為 $ilde{r}$ ，則有：

$ilde{r} =left| vec{r_1}-vec{r}~$

類似地，在三維情形下，我們同樣可以一個由四個球圍起來的縫隙中所能塞入的最大球半徑為：

$ilde{r} =left| vec{r_i}-vec{r}~$

於是，我們只需要找出這些小球圍出來的全部縫隙，並一一求解能填入的最大球，它們其中的最大者即為內部縫隙球集合的 $r_{max}^{in}$ 。下面就根據這個想法來設計一個遍歷內部縫隙的演算法。

Ⅱ、基於圖論描述小球的空間關係

這裡採用圖論中的有向圖模型來描述全體小球。

任取一個小球作為樹根，記其編號為1，則該小球球心的位矢為 $vec{r_1}$ ，現尋找包裹球1的第一層小球構成的集合 $S_1$ 。

(圖5:球2~7為包裹小球1的第一層小球)

這裡採用的閾值 $4r$ 是基於一個考慮：當與球1有關的某縫隙恰好能放入一個半徑為 $r$ 的小球時，與這個內部縫隙球相切的小球與球1的最遠距離為 $4r$ 。

當兩個小球的距離超過 $4r$ 時，存在兩種情況：存在一個內部縫隙球可以與這兩個球均接觸，則該縫隙球的半徑一定大於 $r$ ；沒有任何內部縫隙球同時與兩者產生接觸。前者不滿足題目要求的填充要求，後者可認為兩者是被其它小球分隔開的。

於是，對於球1，以其球心為原點， $4r$ 為半徑，所有球心處於該範圍內的小球均被定義包裹球1的第一層小球，它們的編號構成集合 $S_1$ :

$S_1= { i|left| vec{r_i}-vec{r_1} ight|< 4r ,iinOmega _0-{1} }$ ，此時記 $Omega _1=Omega _0-{1}$ 。

因而，以球1為父親， $S_1$ 中的元素均為它的兒子；再以 $iin S_1$ 為父親，可以得到它的兒子集合為

$S_2^i= { j|left| vec{r_j}-vec{r_1} ight|< 4r ,jinOmega _1-{i} }$ ，此時記 $Omega_2=Omega_1-S_1$

對 $S_1$ 中元素的所有兒子集合求並，可得包裹球1的第二層小球的編號集合為：

$S_2=igcup_{iinOmega_1}S_2^i$

如此遞歸，至 $Omega_M=S_M$ 時所有小球已被遍歷，遞歸結束。在這種分層標準下，全部小球會以球1為核心，包裹共 $M$ 層。

同時，可得到表示空心球中全部小球的有向圖。

例如，由該演算法可將圖5描述為一個有向圖：

（圖6：基於分層迭代關係由圖5所得的有向圖）

但是，這樣的有向圖並不能直觀地展示出其中包含的數學性質，如果定義頂點可重複使用，我們可以將其轉化為一棵有向樹 $T$ ，如下：

（圖7：基於分層迭代關係由圖5所得的有向樹）

圖中標紅的點是指在上一層或本層中存在重複的分支點，不需要往下繼續分支。

容易發現：

（1）這個有向樹中第 $m$ 代（記球1為第0代）後代中黑色點代表的編號組成的集合正是包裹球1的第 $m$ 層小球的編號集合 $S_m$ ；

（2）每一個黑色分支點 $i$ 的所有兒子構成的集合就是前面所求的包裹該球的第一層小球（已去掉遍歷過的小球）的集合 $S_m^i$ ，有 $S_m=igcup_{iin S_{m-1}}S_m^i$ 。

到這裡，我們利用一個特殊的有向樹描述了空心球內小球的空間分布關係。

當我們找到了第一層有哪些小球時，也就確定了第一層中包含哪些縫隙，確定第二層小球的同時就確定的第一層與第二層之間的縫隙，依次遞歸至全部小球，便遍歷了所有的內部縫隙。

接下來就基於這個結構來建立尋找 $r_{max}^{in}$ 的演算法。

Ⅲ、計算 $r_{max}^{in}$

基本想法是：從球1出發，先從集合 $S_1$ 中取出3個小球 $p,q,r$ （互異），根據Ⅰ中的結論，它們與球1一起可以確定一個被圍在其中的最大球半徑為： $ilde{r}_{p,q,r} =left| vec{r_i}-vec{r}~$

其中， $vec{r}~$ 為該縫隙球球心的位矢， $ilde{r}_{p,q,r}$ 為其半徑。

但是，這樣的隨機組合產生的縫隙球存在與其它小球有空間重疊部分的可能性，我們需要進行相應的排除。顯然，只需要保證這些縫隙球不與 $S_1-{p,q,r}$ 的小球不重疊即可（因為它們隔得最近）：

$left| vec{r}~$

同樣地，對第 $m$ 層（ $mleq M-1$ ）小球進行夾層縫隙的檢索，只不過分成了 $left| S_{m-1} ight|$ （表示集合 $S_{m-1}$ 的元素個數）個子集—— $S_m^i(iin S_{m-1})$ ，我們分別完成檢索即可：

$ilde{r}_{p,q,r} =left| vec{r_j}-vec{r}~$

$left| vec{r}~$

據此，我們可以得到內部縫隙球的最大半徑為：

$r_{max}^{in}=max{ ilde{r}_{p,q,r}}$

接下來，我們再考慮邊界縫隙球。

Ι、確定位於邊界的小球集合

與前面確定包裹球1的第一層小球的思路相同，這裡認為：當一個小球的球心距離空心球內表面 $S_{in}$ 小於 $3r$ 時，則認為它屬於邊界球集合 $Omega _{bou}$ :

$Omega _{bou}={i|~vec{r_i}>R-3r,iin Omega_0}$

Ⅱ、判斷邊界縫隙無法再填入半徑為 $r$ 小球的模型

與前面對內部縫隙的處理不同，這裡採用取半徑為 $r$ 的小球與三個相鄰的填充小球相切，考察該球是否與 $S_{in}$ 存在交點來判斷是否可填入。

如下圖，篩選出的邊界球為紅線涉及到的全部小球：

（圖8：邊界小球示意圖）

需要注意兩點性質：

（1）如果把所有邊界小球的球心連接起來，會構成一個凹/凸多面體（兩者均有可能），這些球心即為該多面體的頂點，每一個面均為三角形；

（2）當一個小球與相鄰三個球相切時，直觀來看就是將它擱在一個面的三個頂點上(如圖8的紅色小球)，此時它的球心與三定點的距離均為 $2r$ 。

基於以上兩點，我們先取三個相鄰的填充小球 $p,q,rin Omega _{bou}$ ，對它們與 $S_{in}$ 構成的邊界進行試探的小球球心為 $vec{r}~$ ，半徑為 $r$ ，則試探小球無法塞入的條件是：

$left| vec{r_i}-vec{r} ~$ R-r ~~~~~~~~~~~~~~~~~" eeimg="1">

Ⅲ、遍歷全部的邊界縫隙

當我們拿出一個邊界縫隙時，由Ⅱ中的約束可以保證該縫隙不能再塞入一個半徑為 $r$ 的小球，那麼我們接下來只需要遍歷全部的邊界縫隙即可，也等價於遍歷前述多面體的全部外表面。

與前面遍歷內部縫隙的演算法類似，同樣先採用一個有向圖結構來描述全部的邊界小球，再找到所有的p,q,r組合即可，這裡不再贅述。

綜上，我們可以得到求解本問題（離散情形）的模型為：

$min ~N~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\ s.t.left{egin{matrix} left | vec{r_i} ight |leq R-r~~(iin Omega _0)~~~~~~~~~~~~~\ left| vec{r_i}-vec{r_j} ight| geq 2r~~(i e j,~i,jinOmega _0 )\ r_{max}^{in}=max{ ilde{r}_{p,q,r}}<r~~~~~~~~~\ left| vec{r_i}-vec{r} ~$ R-r ~~~~~~~~~~~~~~~~~~~
end{matrix}
ight." eeimg="1">

對於連續情形，需要增加的唯一約束為：對於任意一個小球 $i$ ，至少有一個填充小球與其接觸。（事實上邊界不需要考慮）

我們定義一個0-1變數 $omega (vec{r_i})$ :

$omega (vec{r_i},vec{r_j})=left{egin{matrix} 1,~~left| vec{r_i}-vec{r_j} ight|=2r \ 0,~~left| vec{r_i}-vec{r_j} ight|> 2r<br /> end{matrix}<br /> ight.<br /> \i<br /> e j,~i,jin Omega _0~~~~~~~~~~~~~~~~~~~~$

則該約束可描述為：

$sum_{jinOmega _0}{omega (vec{r_i},vec{r_j})} geq 1~,~~~i=1,2,...,N$

於是連續情形的模型為：

$min ~N~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\ s.t.left{egin{matrix} left | vec{r_i} ight |leq R-r~~(iin Omega _0)~~~~~~~~~~~~~~~~~~~~\ left| vec{r_i}-vec{r_j} ight| geq 2r~~(i e j,~i,jinOmega _0 )~~~~~~~\ r_{max}^{in}=max{ ilde{r}_{p,q,r}}<r~~~~~~~~~~~~~~~~~\ left| vec{r_i}-vec{r} ~$ R-r ~~~~~~~~~~~~~~~~~~~~~~~~~~~
\
sum_{jinOmega _0}{omega (vec{r_i},vec{r_j})} geq 1~,~~~i=1,2,...,N
end{matrix}
ight." eeimg="1">

五、模型演算法的複雜性問題

上面建立的模型為全局搜索模型，結合一些分析（主要是球的密堆問題），可知模型的演算法複雜度是 $O(n!)$ 。至於本問題是否為NP完全問題，我還沒有仔細考慮（事實上我也不會）。

常規的全局最優解建模思路就是這樣。

如果問題規模太大 $Rgg r$ ,而且計算時間無法接受，同時可以接受局部最優解，那麼可以採取啟發式演算法進行妥協。這裡就不談了。

--------------------------------

全文是邊做題邊寫的，省略了部分證明和說理的過程，國慶有時間的話就改改，沒時間就算了。

：）

有問題可以私信或評論。

是否也可以這樣考慮：考慮到太陽等天體大概都是球體形狀，並且這種匯聚力來源於萬有引力，由於縫隙的存在，因此能給出一個上限 Nmax = Vb/Vs，這樣從上限數目開始用匯聚演算法，直到找到一個數目與Vb最接近為止。