面白い最適化問題！

問題へのリンク

問題概要

二次元平面上で、原点からの距離が $r_{1}, r_{2}, \dots, r_{N}$ であるような $N$ 点の凸包の面積として考えられる最大値を求めよ。

制約

$3 \le N \le 8$
$1 \le r_{i} \le 1000$

考えたこと

凸包というところが面倒だが、要は

$N$ 本のうちの何本かを選んで
それを適切な順序に並び替えたときの半径を改めて $r_{0}, r_{1}, \dots, r_{K-1}$ として、 $r_{i}, r_{i+1}$ のなす角を $\theta_{i}$ としたとき
$r_{0}r_{1}\sin \theta_{0} + r_{1}r_{2} \sin \theta_{1} + \dots + r_{K-1}r_{0} \sin \theta_{K-1}$ (の半分) の最大値を求めよ

という問題を考えれば良い。ちゃんと最適化問題として書き下すと、次のようになる。

Maximize:

$r_{0}r_{1}\sin \theta_{0} + r_{1}r_{2} \sin \theta_{1} + \dots + r_{K-1}r_{0} \sin \theta_{K-1}$

s.t.

$\theta_{0} + \theta_{1} + \dots + \theta_{K-1} = 2\pi$
$0 \le \theta_{i} \le \pi$

ここで、 $\theta_{i} \gt \pi$ のケース (原点が凸包に含まれなくなるケース) は考える必要はない。なぜならば、その場合には、面積を小さくさせることなく、凸包が原点を含むように点を反転することができるためだ。

ラグランジュの未定乗数法

この手の最大化問題はラグランジュの未定乗数法で解ける。本当は $\theta_{i} \ge 0$ や $-1 \le \cos \theta _{i} \le 1$ などといった不等式制約があるので KKT 条件を考える方が正確なのだが、大雑把にラグランジュの未定乗数法を使えば大体うまくいくイメージ。厳密には、KKT 条件が加わることで「最適解が不等式制約の端点に張り付く」といった状態が表現できる (参考)。

$L(\theta_{i}, \lambda) = (r_{0}r_{1}\sin \theta_{0} + r_{1}r_{2} \sin \theta_{1} + \dots + r_{K-1}r_{0} \sin \theta_{K-1}) - \lambda(\theta_{0} + \theta_{1} + \dots + \theta_{K-1} - 2\pi)$

として、偏微分 = 0 として

$\lambda = r_{0}r_{1} \cos \theta_{0}$
$\lambda = r_{1}r_{2} \cos \theta_{1}$
$\dots$
$\lambda = r_{K-1}r_{0} \cos \theta_{K-1}$
$\theta_{0} + \theta_{1} + \dots + \theta_{K-1} = 2\pi$

となる。よって、 $\lambda$ の値を適切に定めることで、

$\theta_{0} + \theta_{1} + \dots + \theta_{K-1} = 2\pi$

が成り立つようにすればよい。 $\lambda$ の値によって、 $\theta_{0} + \theta_{1} + \dots + \theta_{K-1}$ の値は単調減少となることに着目して、二分探索法によって求められることがわかる。

なお、適当にやると $\cos \theta_{i} \gt 1$ や $\cos \theta_{i} \lt -1$ になってしまうような場合もある。そのようなときは $\cos \theta_{i} = 1$ や $\cos \theta_{i} = -1$ となるように $\theta_{i}$ を端点に貼り付ければよい。このことは、厳密には KKT 条件を用いた考察によって正当化できる。

コード

#include <bits/stdc++.h>
using namespace std;

template<class T> inline bool chmax(T& a, T b) { if (a < b) { a = b; return 1; } return 0; }
template<class T> inline bool chmin(T& a, T b) { if (a > b) { a = b; return 1; } return 0; }
const double EPS = 1e-10;  // to be set appropriately
const double INF = 1<<29;
const double PI = 3.141592653589793238462643383279502884L;

int main() {
    auto solve = [](const vector<double> &r) -> double {
        int K = (int)r.size();
        if (K < 3) return 0;
        
        // binary search
        double low = -INF, high = INF, sum_theta = 0;
        vector<double> theta(K);
        for (int iter = 0; iter < 100; ++iter) {
            double lambda = (low + high) / 2;
            sum_theta = 0;
            for (int i = 0; i < K; ++i) {
                double c = lambda / r[i] / r[(i+1)%K];
                chmin(c, 1.0 - EPS), chmax(c, -1.0 + EPS);
                theta[i] = acos(c);
                sum_theta += theta[i];
            }
            if (sum_theta < PI*2) high = lambda;
            else low = lambda;
        }

        // calc area
        double res = 0;
        for (int i = 0; i < K; ++i) res += r[i] * r[(i+1)%K] * sin(theta[i]) / 2;
        return res;
    };
    
    int N;
    cin >> N;
    vector<double> R(N);
    for (int i = 0; i < N; ++i) cin >> R[i];
    sort(R.begin(), R.end());
    
    double res = 0;
    for (int bit = 0; bit < (1<<N); ++bit) {
        vector<double> r;
        for (int i = 0; i < N; ++i) if (bit & (1 << i)) r.push_back(R[i]);
        do {
            chmax(res, solve(r));
        } while (next_permutation(r.begin(), r.end()));
    }
    cout << fixed << setprecision(10) << res << endl;
}

けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

AOJ 2373 HullMarathon (JAG 冬合宿 2010 day3-E) (700 点)

問題概要

制約

考えたこと

Maximize:

s.t.

ラグランジュの未定乗数法

コード