问题解答
给定点集 \(\mathcal{P}=\{P_1,P_2,\cdots,P_n\}\) ，求包围该点集的最小椭圆 \(\mathcal{E}\) 。
椭圆大小
给定椭圆的半轴 \(a\) 、\(b\) ，如何定义椭圆的大小？
- 面积：\(\mathrm{Area}(\mathcal{E})=\pi ab\)
- 周长：\(\mathrm{Perimeter}(\mathcal{E})=\int_0^{2\pi} \sqrt{(a\sin t)^2+(b\cos t)^2}~\mathrm{d}t\)
对于周长 \(\mathrm{Perimeter}(\mathcal{E})\) 没有一个简单的初等函数表达形式，故采用面积 \(\mathrm{Area}(\mathcal{E})\) 来定义椭圆大小。
问题建模
数学抽象
将二次曲线形式的椭圆方程 \[ Ax^2+Bxy+Cy^2+Dx+Ey+F=0 \] 变换为中心形式，即 \[ A(x-x_0)^2+B(x-x_0)(y-y_0)+C(y-y_0)^2=f \] 其中 \[ x_0=\frac{2CD-BE}{B^2-4AC}\qquad y_0=\frac{2AE-BD}{B^2-4AC} \]
\[ f=-F+Ax_0^2+Bx_0y_0+Cy_0^2 \]
将中心形式方程记为二次型，即 \[ \underbrace{\begin{bmatrix} x-x_0 \\ y-y_0 \end{bmatrix}^{\mathrm{T}}}_{\vec{\boldsymbol{x}}^{\mathrm{T}}}\cdot \underbrace{\begin{bmatrix} A/f&B/2f \\ B/2f&C/f \end{bmatrix}}_{M}\cdot\underbrace{\begin{bmatrix} x-x_0 \\ y-y_0 \end{bmatrix}}_{\vec{\boldsymbol{x}}}=1 \] 对于实对称矩阵 \(M\) ，存在正交矩阵 \(Q\) ，使得 \[ Q\cdot M \cdot Q^{\mathrm{T}}=\begin{bmatrix} \lambda_1&0 \\ 0&\lambda_2 \end{bmatrix} \] 若 \(\det(Q)=1\) （逆时针旋转为正），则有 \[ a^2=1/\lambda_1 \qquad b^2=1/\lambda_2 \] 于是有 \[ \mathrm{Area}(\mathcal{E})=\pi ab=\pi\sqrt{\frac{1}{\lambda_1\lambda_2}}=\pi\sqrt{\frac{1}{\det(M)}} \] 即当 \(\det(M)\) 取最大值时，椭圆面积最小。
优化建模
中心形式方程中有六个未知数，可将其中的 \(f\) 作为归一化量，则设计变量为 \((A,B,C,x_0,y_0)\) 。
\(A(x-x_0)^2+B(x-x_0)(y-y_0)+C(y-y_0)^2=f\) ，可令 \(f\) 为常数
显然，矩阵 \(M\) 的两个特征值均需大于零，考虑到 \[ \lambda_1\lambda_2=\frac{4AC-B^2}{4f^2}\qquad \lambda_1+\lambda_2=\frac{A+C}{f} \] 于是有 \[ 4AC-B^2>0\qquad (A+C)f>0 \] 显然，\(A\) 、\(C\) 、\(f\) 同号。
考虑到点集 \(\mathcal{P}=\{P_1,P_2,\cdots,P_n\}\) 在椭圆 \(\mathcal{E}\) 内，则有 \[ A(x_i-x_0)^2+B(x_i-x_0)(y_i-y_0)+C(y_i-y_0)^2\le f \] 其中，\((x_i,y_i)\) 为第 \(i\) 个点的坐标。
综上，该优化问题为 \[ \begin{array}{l} \min.\quad f(A,B,C,x_0,y_0)=\left.\displaystyle\frac{4f^2}{4AC-B^2}\right|_{f=\mathrm{positive~constant}}\\[8pt] \mathrm{st}.\quad \left\{\begin{array}{l} 4AC-B^2>0 \\[5pt] A>0 \\[5pt] C>0 \\[5pt] A(x_i-x_0)^2+B(x_i-x_0)(y_i-y_0)+C(y_i-y_0)^2\le f\quad(i=1,2,\cdots,n) \end{array}\right. \end{array} \]
计算程序
计算步骤如下：
- 使用 MATLAB 内置函数
convhull求给定点集 \(\mathcal{P}=\{P_1,P_2,\cdots,P_n\}\) 的凸包，并记为 \(\mathcal{P}'\) ； - 按照优化思想，求凸包 \(\mathcal{P}'\) 的最小外包圆，并记其圆心为 \(P_0\) 、半径为 \(r_0\) ；
- 对优化问题 \(\min.f\) ，补充约束条件 \(\sqrt{4f^2/(4AC-B^2)} \le r_0^2\) ，即凸包 \(\mathcal{P}'\) 的最小外包椭圆面积不超过最小外包圆面积；
- 将 \(P_0\) 和 \(r_0\) 作为优化问题 \(\min.f\) 的初值，调用函数
fmincon即可。
由 \(4AC-B^2>0\) 和 \(\sqrt{4f^2/(4AC-B^2)} \le r_0^2\) 可得出 \(4AC-B^2 \ge 4f^2/r_0^4\)
MATLAB 代码如下：
1 | clc;clear;close all; |