题解：P17737 [ICPC 2019 Xuzhou R] Factory

· · 题解

先把“每种饮料都指定一台机器，整升交给它生产”称为一个方案。记方案 s 中机器 j 的总工作时间除以 N 为 a_{js}。除以 N 是为了让线性规划中的系数保持在较小的范围内。

虽然一个方案不拆分饮料，但按比例混合多个方案，就可以表示任意连续生产安排。具体地，设第 i 种饮料由机器 j 生产的比例为 p_{ij}。对每种饮料独立地按这些比例选择机器，再把得到的各个方案按照出现的概率混合，就能还原所有 p_{ij}。因此，使用这样的方案不会损失最优解。

现在换一个角度：在 N 分钟内，最多能生产多少份完整订单？设按照方案 s 生产 z_s 份订单，可以写出线性规划：

\begin{aligned} \text{最大化}\quad &F=\sum_s z_s,\\ \text{满足}\quad &\sum_s a_{js}z_s\le 1 &&(1\le j\le K),\\ &z_s\ge 0. \end{aligned}

如果完整问题的最优值为 F^*，原题答案就是 N/F^*，因为生产量和时间可以按相同比例缩放。这里虽然有 K^N 个方案，却只有 K 条机器容量约束。这正适合列生成：开始只加入 K 个方案，第 j 个方案把所有饮料都交给机器 j；之后每次只加入一个有用的新方案。

假设已经用单纯形法求出了当前这些方案的最优值 F，那么当前方案确实能完成订单，所以 U=N/F 是答案的上界。为了寻找下一个方案，考虑这个有限线性规划的对偶：给机器赋予非负权值 y_j，最小化 \sum_j y_j，并要求每个已加入方案满足 \sum_j a_{js}y_j\ge1。这些权值可以从单纯形表中直接读出：原始松弛变量的检验数取负，就是对应的 y_j。

接下来寻找加权耗时最小的方案。各饮料的选择相互独立，第 i 种饮料只需要交给使 t_{ij}y_j 最小的机器。因此一次 O(NK) 扫描，就能求出最小值，并构造对应方案：

v=\min_s\sum_j a_{js}y_j =\frac1N\sum_{i=1}^N\min_{1\le j\le K}(t_{ij}y_j).

若 v<1，这个方案违反了当前的对偶约束，也就是一个有机会改善答案的新变量。把它加入单纯形表，重新优化即可。所有以前加入的方案都保留，避免反复丢弃已有信息。

更关键的是，这次扫描还给出了答案的下界，而且不必依赖对偶问题已经被精确求解。对于任意非负权值 y_j，生产第 i 种饮料的加权耗时至少为 \min_j(t_{ij}y_j)；如果一份订单在 T 分钟内完成，每台机器的工作时间又至多为 T。将所有饮料的加权耗时相加，就得到

T\sum_j y_j\ge\sum_i\min_j(t_{ij}y_j)=Nv.

所以 L=Nv/\sum_jy_j 是答案的下界。代码会把浮点误差产生的微小负权值截成 0，再用这个公式计算下界。

为减轻大量相同耗时引起的浮点退化，代码把机器约束的右端项从 1 轻微调整为 h_j=1+10^{-5}j/K。这个调整也计入误差分析：当前解对应的实际机器工作时间至多为 Nh_j/F，因此把上界改成 U=(1+10^{-5})N/F，仍然是原题的合法上界。下界按原题计算，无须修改。扰动后的对偶目标相应变成 \sum_jh_jy_j，寻找新方案的公式保持不变。

当 L\ge0.982U 时，输出 A=(L+U)/2。设真实答案为 T^*，由 L\le T^*\le U 可得

\frac{|A-T^*|}{T^*} \le\frac{U-L}{2L} \le\frac{0.018}{2\times0.982} <0.01.

这就满足题目的相对误差要求，并留有一定的数值误差余量。

设总共额外加入 R 个方案，单纯形进行了 P 次换基。每次寻找新方案需要 O(NK)，每次换基需要 O(K(K+R))，总时间为 O((R+1)NK+PK(K+R))，空间为 O(NK+K(K+R))。这里不能把 R,P 当作有严格常数上界的参数。

下面是完整的 C++17 实现。单纯形表的前 K 列始终保留原始松弛变量，它们既用于读取机器权值，也用于把新方案转换到当前基下。

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

struct Simplex {
    static constexpr double EPS = 1e-7;
    int k;
    vector<vector<double>> a;
    vector<double> b, c;
    double value = 0;

    // 前 k 列是原始松弛变量，后面是已经加入的生产方案。
    Simplex(int n, const vector<double>& sum) : k(sum.size()),
        a(k, vector<double>(2 * k)), b(k), c(2 * k) {
        for (int i = 0; i < k; ++i) {
            b[i] = n / sum[i];
            a[i][i] = b[i];

            // 对机器容量加入微小扰动；最终计算上界时会补偿。
            b[i] *= 1 + 1e-5 * (i + 1) / k;

            a[i][k + i] = 1;
            c[i] = -a[i][i];
            value += b[i];
        }
    }

    void pivot(int r, int e) {
        int m = c.size();
        double v = a[r][e];

        b[r] /= v;
        for (double& x : a[r]) x /= v;

        for (int i = 0; i < k; ++i) if (i != r) {
            v = a[i][e];
            b[i] -= v * b[r];
            for (int j = 0; j < m; ++j) {
                a[i][j] -= v * a[r][j];
            }
            a[i][e] = 0;
        }

        v = c[e];
        value += v * b[r];
        for (int j = 0; j < m; ++j) {
            c[j] -= v * a[r][j];
        }
        c[e] = 0;
    }

    // 加入一个方案，重新优化当前有限个方案的线性规划。
    void add(const vector<double>& col) {
        // 前 k 列保存了基变换，用它们转换新列。
        for (int i = 0; i < k; ++i) {
            double v = 0;
            for (int j = 0; j < k; ++j) {
                v += a[i][j] * col[j];
            }
            a[i].push_back(v);
        }

        double reduced = 1;
        for (int j = 0; j < k; ++j) {
            reduced += c[j] * col[j];
        }
        c.push_back(reduced);

        while (true) {
            int e = max_element(c.begin(), c.end()) - c.begin();
            if (c[e] <= EPS) return;

            int r = -1;
            for (int i = 0; i < k; ++i) if (a[i][e] > 1e-9) {
                if (r == -1) {
                    r = i;
                    continue;
                }
                double u = b[i] / a[i][e];
                double v = b[r] / a[r][e];
                if (u < v) r = i;
            }
            pivot(r, e);
        }
    }
};

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int n, k;
    cin >> n >> k;

    vector<array<unsigned char, 50>> t(n);
    vector<double> sum(k, 0);

    for (auto& row : t) {
        for (int j = 0; j < k; ++j) {
            int v;
            cin >> v;
            row[j] = v;
            sum[j] += v;
        }
    }

    Simplex lp(n, sum);
    vector<double> y(k), col(k);

    while (true) {
        double sy = 0;
        for (int j = 0; j < k; ++j) {
            y[j] = max(0.0, -lp.c[j]);
            sy += y[j];
        }

        fill(col.begin(), col.end(), 0);

        // 每种饮料独立选择加权耗时最小的机器。
        for (const auto& row : t) {
            int p = 0;
            double best = row[0] * y[0];

            for (int j = 1; j < k; ++j) {
                double v = row[j] * y[j];
                if (v < best) {
                    best = v;
                    p = j;
                }
            }
            col[p] += row[p];
        }

        double v = 0;
        for (int j = 0; j < k; ++j) {
            col[j] /= n;
            v += col[j] * y[j];
        }

        // 上界补偿容量扰动，下界始终针对原问题计算。
        double upper = (1 + 1e-5) * n / lp.value;
        double lower = min(upper, n * v / sy);

        if (lower >= 0.982 * upper) {
            cout << fixed << setprecision(10)
                 << (lower + upper) / 2 << '\n';
            return 0;
        }

        lp.add(col);
    }
}