[BZOJ2876][Noi2012]骑行川藏（二分答案+拉格朗日乘数法）

Address

https://www.lydsy.com/JudgeOnline/problem.php?id=2876

Solution

为了方便某些萌新，讲下拉格朗日乘数法。
作用：求多元函数 f(x1,x2,...,xn)$f(x_1,x_2,...,x_n)$ 在约束 ϕ(x1,x2,...,xn)=0$\varphi (x_1,x_2,...,x_n)=0$ 下的极值。引入变量 λ$\lambda$ ，令拉格朗日函数：

F(x1,x2,...,xn,λ)=f(x1,x2,...,xn)+λϕ(x1,x2,...xn)$$F(x_1,x_2,...,x_n,\lambda )=f(x_1,x_2,...,x_n)+\lambda \varphi (x_1,x_2,...x_n)$$

对 F$F$ 的每个变量（包括 λ$\lambda$ ）分别求偏导。
（事实上对 λ$\lambda$ 求偏导后的结果就是 ϕ(x1,x2,...,xn)$\varphi (x_1,x_2,...,x_n)$ ）
（多元函数 F$F$ 对一个自变量 x$x$ 的偏导数，定义为将 F$F$ 除 x$x$ 之外的自变量全部视为常数后对 F$F$ 求导）
当 F$F$ 对每个自变量的偏导数都等于 0$0$ 时，
f(x1,x2,...,xn)$f(x_1,x_2,...,x_n)$ 在约束 ϕ(x1,x2,...,xn)=0$\varphi (x_1,x_2,...,x_n)=0$ 下取得极值。
回到原问题。
发现体力不是一个标准的约束。
但如果体力没用完，则一定可以在某一天速度更快。
所以，存在一种最优方案使得体力恰好用完。
同样我们也能得到，存在一种最优方案使得对于每个 i$i$ 有 vi≥v′i$v_i\ge v_i^{\prime }$ 。
（ vi$v_i$ 表示第 i$i$ 天实际速度）
目标函数：

f(v1,v2,...,vn)=∑i=1nsivi$$f(v_1,v_2,...,v_n)=\sum _{i=1}^n\frac{s_i}{v_i}$$

约束条件：

ϕ(v1,v2,...,vn)=−EU∑i=1nki(vi−v′i)2si=0$$\varphi (v_1,v_2,...,v_n)=-E_U\sum _{i=1}^n{k}_i(v_i-v_i^{\prime }{)}^2{s}_i=0$$

引入系数 λ$\lambda$ ，构造拉格朗日函数：

∑i=1n(sivi+λki(vi−v′i)2si)−EU$$\sum _{i=1}^n(\frac{s_i}{v_i}+\lambda k_i(v_i-v_i^{\prime }{)}^2{s}_i)-E_U$$

对 vi$v_i$ 求偏导得到：

−siv2i+λki(2vi−2v′i+v′2i)si$$-\frac{s_i}{v_i^2}+\lambda k_i(2v_i-2v_i^{\prime }+v_i^{\prime 2})s_i$$

要求上式为 0$0$ 。

−1v2i+λki(2vi−2v′i)=0$$-\frac{1}{v_i^2}+\lambda k_i(2v_i-2v_i^{\prime })=0$$

2λki(vi−v′i)v2i=1$$2\lambda k_i(v_i-v_i^{\prime })v_i^2=1$$

易得上式左边在 vi≥v′i$v_i\ge v_i^{\prime }$ 即上式左边的值大于等于 0$0$ 时单调递增。
故在 λ$\lambda$ 确定时二分答案可以算出 vi$v_i$ 。
同时， λ$\lambda$ 越大， vi$v_i$ 越小。
所以 λ$\lambda$ 的值也可以二分，通过 ϕ(v1,v2,...,vn)$\varphi (v_1,v_2,...,v_n)$ 的值来调整 λ$\lambda$ 。
最后得出了 λ$\lambda$ 的值之后就能计算取到极值时的 vi$v_i$ 。

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define For(i, a, b) for (i = a; i <= b; i++)
using namespace std;
const int N = 1e4 + 5;
const double eps = 1e-12;
int n;
double eu, k[N], s[N], vp[N], ans;
double vsol(double lambda, int i) {
    double l = max(vp[i], 0.0), r = 1e5;
    while (r - l > eps) {
        double mid = (l + r) / 2, val;
        val = 2.0 * lambda * k[i] * (mid - vp[i]) * mid * mid;
        if (val < 1) l = mid;
        else r = mid;
    }
    return (l + r) / 2;
}
double calc(double lambda) {
    int i; double res = 0, tmp;
    For (i, 1, n) tmp = vsol(lambda, i),
        res += k[i] * s[i] * (tmp - vp[i]) * (tmp - vp[i]);
    return res;
}
int main() {
    int i;
    double x, y, lambda, l = 0, r = 1e5;
    cin >> n >> eu;
    For (i, 1, n) scanf("%lf%lf%lf", &s[i], &k[i], &vp[i]);
    while (r - l > eps) {
        double mid = (l + r) / 2;
        if (calc(mid) > eu) l = mid;
        else r = mid;
    }
    For (i, 1, n) ans += s[i] / vsol((l + r) / 2, i);
    printf("%.10lf\n", ans);
    return 0;
}