N圆最密堆积、最小外接正方形的matlab求解（二维、三维等圆Packing 问题）

惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

0 前言

本文最开始的想法，是想到一个我有若干个球，究竟需要多大的箱子可以把这些球装下。

后来大概查了一下资料，发现这其实是一个最优化的问题，而且似乎还非常复杂。对于圆来说，还有比较简单的优化方法，但是对于复杂图形来说，则非常困难。

自己也不是专门搞最优化的，就不编自己的函数了，直接用matlab现成的函数进行求解。下面举两个例子作为说明。

更好的例子可参见这篇2019年的论文：
王正理《等圆Packing问题高效求解算法研究》
20以下的最密堆积可以直接在wiki百科上查表得到，词条为：Circle packing in a square
10000以下的最密堆积可以在这个网站查到：
http://hydra.nat.uni-magdeburg.de/packing/csq/csq.html

1 N个圆的最小外接正方形求解

这个先用2维平面问题作为引入。

优化函数用的是fmincon()函数。
这里之前也试过PSO算法，但是估计维度太多，有些空穴小球移动完全不影响最终外接矩形，所以收敛速度异常的慢。经常出现某一些圆挤到一块，但某些圆离得很远。所以还是用传统的最优化方法。

目标优化函数，也就是让其最小的函数，就是正方形的变长。
边界大致设置了左右边界，防止圆离开边界太远。
非线性约束，就是计算每个圆之间的距离，防止发生重叠。

在计算过程中，最容易出现的现象就是圆和圆之间距离过大，就提前结束计算了，比如下面这个图：

为了让其收敛，需要手动的将各个圆之间的距离缩小。我这里用的方法比较简单，就是将所有坐标除以1.5或者其它大于1的数，让其整体向左下角压缩。通常在计算一定步骤之后，再进行压缩，然后继续计算，反复多次，可以让其收敛到较好的结果。

最终代码如下：

clear
clc
close all
%计算圆的最小外接矩形（任意数量N）
%利用fmincon算法
%求解函数最小值

N=13;%圆的数量

%设置optimoptions的输入
fun = @MinValue;
lb = -2*ones(1,N*2);%[xi,yi,xi,yi,...]
ub = (4*N-2)*ones(1,N*2);

A = [];
b = [];
Aeq = [];
beq = [];

%按照网格布置初始圆
NMesh=ceil(sqrt(N));
[XMesh,YMesh]=meshgrid(0:4:(4*NMesh-1),0:4:(4*NMesh-1));
XYMesh=[XMesh(:),YMesh(:)]';
x0=XYMesh(1:2*N);
x0=x0+rand(1,2*N)-0.5;%添加随机性

%非线性约束，目的是让圆和圆之间不重叠
nonlcon = @NolinearCondition;%设定非线性约束条件

options = optimoptions('fmincon','Display','iter','Algorithm','interior-point','MaxFunctionEvaluations',inf);
x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);
%放大圆形，然后再计算一次。相当于让圆膨胀一下，然后利用约束条件再迭代一下
x=x/1.2;
options = optimoptions('fmincon','Display','iter','Algorithm','sqp',&#