「NOI2012」骑行川藏

题目

点这里看题目。

---

题目描述

蛋蛋非常热衷于挑战自我，今年暑假他准备沿川藏线骑着自行车从成都前往拉萨。

川藏线的沿途有着非常美丽的风景，但在这一路上也有着很多的艰难险阻，路况变化多端，而蛋蛋的体力十分有限,因此在每天的骑行前设定好目的地，同时合理分配好自己的体力是一件非常重要的事情。

由于蛋蛋装备了一辆非常好的自行车，因此在骑行过程中可以认为他仅在克服风阻做功（不受自行车本身摩擦力以及自行车与地面的摩擦力影响）。

某一天他打算骑 \(n\) 段路，每一段内的路况可视为相同：对于第 \(i\) 段路,我们给出有关这段路况的 \(3\) 个参数 \(s_i,k_i,v_i'\)，其中 \(s_i\) 表示这段路的长度，\(k_i\) 表示这段路的风阻系数，\(v_i'\) 表示这段路上的风速（\(v_i'\gt 0\) 表示在这段路上他遇到了顺风，反之则意味着他将受逆风影响）。

若某一时刻在这段路上骑车速度为 \(v\)，则他受到的风阻大小为 \(F=k_i(v-v_i')^2\)（这样若在长度为 \(s\) 的路程内保持骑行速度 \(v\) 不变，则他消耗能量（做功）\(E=k_i(v-v_i')^2s\) ）。

设蛋蛋在这天开始时的体能值是 \(E_U\)，请帮助他设计一种行车方案,使他在有限的体力内用最短的时间到达目的地。请告诉他最短的时间 \(T\) 是多少。

数据规模与约定

所有数据满足 \(1\le n\le 10^4,E_U\le 10^8,s_i\in (0,10^5],k_i\in (0,15],v_i'\in (-100,100)\)。

提示

必然存在一种最优的体力方案满足：蛋蛋在每段路上都采用匀速骑行的方式。

分析

在“提示”的提示下，该问题等价于以下最优化问题：

\[\begin{matrix} \min&\displaystyle\sum_{i=1}^ns_iv_i^{-1}\\ \text{s.t.}&\displaystyle\sum_{i=1}^nk_is_i(v_i-v_i')^2=E_U\\ &v_i\in \mathbb R_+&i=1,2,3,\dots,n \end{matrix} \]

第一个条件原本应该是 \(\le E_U\)，但是显然取到 \(\min\) 时有 \(=E_U\) 成立。

先忽略 \(v_i\in \mathbb R_+\) 的条件（这实际上是不等式条件），则剩下的部分可以用拉格朗日乘子法解决。记拉格朗日函数：

\[L(v,\lambda)=\sum_{i=1}^ns_iv_i^{-1}+\lambda\left[\sum_{i=1}^nk_is_i(v_i-v_i')^2-E_U\right] \]

现在解决 \(L\) 的无约束最值问题。令其各维偏导为 \(0\)：

\[\begin{cases} \frac{\partial}{\partial v_i} L=-s_iv_i^{-2}+2\lambda k_is_i(v_i-v_i')=0&\forall i=1,2,3\dots,n&\\\ \frac{\partial}{\partial \lambda} L=\sum_{i=1}^nk_is_i(v_i-v_i')^2-E_U=0 \end{cases} \]

注意到第一类方程只涉及 \(\lambda\) 和 \(v_i\) 两个变量，我们尝试将 \(v_i\) 解出来，便于代入第二类方程求解。因为 \(v_i\neq 0\)，我们可以将其整理成一个三次方程的形式：

\[2\lambda k_i\cdot v_i^2(v_i-v_i')-1=0 \]

注意到在最值点处一定有 \(v_i\ge v_i'\)，所以一定有 \(\lambda>0\)。于是可以分离变量：

\[v_i^2(v_i-v_i')=\frac{1}{2\lambda k_i} \]

令 \(f_i(x)=x^2(x-v_i')\)，讨论 \(f_i\) 的形态：

• 如果 \(v_i'\le 0\)，则 \(x\in (-\infty,0]\) 时 \(f_i(x)\le 0\)；\(x\in (0,+\infty)\) 时 \(f_i(x)>0\)。

  因为 \(\frac{1}{2\lambda k_i}>0\)，所以 \(f_i(x)=\frac{1}{2\lambda k_i}\) 存在且仅存在正根。更进一步地，根据 \(f_i\) 在 \((0,+\infty)\) 上的单调性可知存在且仅存在唯一正根，这个正根就应该是我们想要的 \(v_i\)。

• 如果 \(v_i'>0\)，则 \(x\in (-\infty,v_i']\) 时 \(f_i(x)\le 0\)；\(x\in(v_i',+\infty)\) 时 \(f_i(x)>0\)。

  类似讨论可知，方程存在且仅存在唯一正根，我们取它为 \(v_i\) 即可。

所以，已知 \(\lambda\) 时算出一个 \(v_i\) 可以直接二分或牛顿迭代。令 \(g(\lambda)=\sum_{i=1}^nk_is_i[v_i(\lambda)-v_i']^2-E_U\)。现在我们已经可以求出一个 \(g(\lambda)\)，但怎么找出 \(g(\lambda)\) 的根？

简单地分析一下，可以发现 \(\lambda\) 变大时，\(\frac{1}{2\lambda k_i}\) 变小，于是 \(v_i(\lambda)\) 统一变小且保持 \(v_i(\lambda)>v_i'\)，于是 \(g(\lambda)\) 也会变小。这就暗示我们 \(g(\lambda)\) 在 \((0,+\infty)\) 上就是单调的。又因为 \(\lambda\) 足够接近 \(0\) 时 \(g(\lambda)>0\)，充分大时 \(g(\lambda)<0\)，就可以用二分找出根了。

---

因为我高兴，所以下面是无关紧要的运算：

• 第一个是 \(f_i(x)=x^2(x-v'_i)\)：

  \[f_i'(x)=3x^2-2v_i'x \]

  可见在 \(x>\max\{v_i',0\}\) 时总有 \(f_i'(x)>0\)。

• 第二个是 \(v_i\)。方程两边同时对 \(\lambda\) 求导：

  \[2k_iv_i^2(v_i-v_i')+6\lambda k_iv_i^2\cdot \frac{\text d v_i}{\text d\lambda}-4\lambda k_iv_i'v_i\cdot \frac{\text d v_i}{\text d\lambda}=0 \]

  可以得到：

  \[\frac{\text dv_i}{\text d\lambda}=\frac{v_i(v_i'-v_i)}{\lambda(3v_i-2v_i')} \]

  因为 \(v_i>\max\{v_i',0\},\lambda>0\)，故这个导数始终是 \(<0\) 的。

• 第三个是 \(g(\lambda)\)：

  \[g'(\lambda)=2\sum_{i=1}^nk_is_i(v_i-v_i')\cdot \frac{\text dv_i}{\text d\lambda} \]

  容易发现 \(g'(\lambda)\) 在 \(\lambda>0\) 时始终 \(<0\)。

---

最后说一说“提示”。不清楚证明它可以不可以先划分成 \(n\) 段（要求每段内部 \(v\) 相同），再证明最值点处全局 \(v\) 相同，再令 \(n\rightarrow +\infty\) 从而产生“极限”。

如果有懂哥的话请教授一下 /kel。

代码

#include <bits/stdc++.h>

#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )

const int NEWTON_ROUND = 30;
const int SEARCH_ROUND = 50;
const int MAXN = 1e4 + 5;

template<typename _T>
inline void Read( _T &x ) {
    x = 0; char s = getchar(); bool f = false;
    while( s < '0' || '9' < s ) { f = s == '-', s = getchar(); }
    while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
    if( f ) x = -x;
}

template<typename _T>
inline void Write( _T x ) {
    if( x < 0 ) putchar( '-' ), x = -x;
    if( 9 < x ) Write( x / 10 );
    putchar( x % 10 + '0' );
}

double resis[MAXN], dist[MAXN], velo[MAXN];

int N; double E;

inline double FindRoot( const double &a, const double &b ) {
    double u = std :: max( 1., b / a );
    rep( stp, 1, NEWTON_ROUND )
        u = u * ( 1 - ( a * u - b ) / ( 3 * a * u - 2 * b ) ) + 1. / u / ( 3 * a * u - 2 * b );
    return u;
}

inline bool Chk( const double &lamb ) {
    double su = 0;
    rep( i, 1, N ) {
        double c = 2 * lamb * resis[i];
        double v = FindRoot( c, c * velo[i] );
        su += dist[i] * resis[i] * ( v - velo[i] ) * ( v - velo[i] );
    }
    return su >= E;
}

int main() {
    Read( N ); scanf( "%lf", &E );
    rep( i, 1, N ) scanf( "%lf %lf %lf", &dist[i], &resis[i], &velo[i] );
    double x = 0;
    while( true ) {
        if( Chk( x + 1 ) )
            x = x + 1;
        else {
            double l = x, r = x + 1, mid;
            rep( stp, 1, SEARCH_ROUND ) {
                mid = ( l + r ) * 0.5;
                if( Chk( mid ) ) l = mid;
                else r = mid;
            }
            x = l;
            break;
        }
    }
    double ans = 0;
    rep( i, 1, N ) {
        double c = 2 * x * resis[i];
        ans += dist[i] / FindRoot( c, c * velo[i] );
    }
    printf( "%.10f\n", ans );
    return 0;
}