【拉格朗日乘数法】[NOI2012]骑行川藏

题目描述

分析

首先，我们来考虑一下部分分。

N=1

直接算即可，v1=Ek1×s1−−−−−√+v′1T=s1v1

N=2

v1=E1k1×s1−−−−−√+v′1v2=E−E1k2×s2−−−−−√+v′2T=s1v1+s2v2
我们猜测T关于E1的函数是一个单谷函数，对E1三分一下就可以了

N≤100&N≤1000

我不会QAQ

N≤10000

首先，我们将题目要求转化为数学式子

E≥∑i=1Nkisi(vi−v′i)2T=∑i=1Nsivi

求Tmin
基于贪心的思想，显然E=∑Ni=1kisi(vi−v′i)2时，答案最优。

拉格朗日乘数法

具体证明我不会。。。。
例题

在g(x1,x2,x3,…,xn)=c的约束条件下，求f(x1,x2,x3,…,xn)的最值。

我们先来讲一下结论。

结论

当f(x1,x2,x3,…,xn)取到最值时，f(x1,x2,x3,…,xn)和g(x1,x2,x3,…,xn)=c的法向量（也就是梯度向量）平行。

我们将例题具体化，比如f(x,y)=x2+y2g(x,y)=xy=3，求f(x,y)min
先来肉眼观察，很容易发现f(x,y)min=6
然后我们用拉格朗日乘数法试试。
设

∇f(x,y)=λ∇g(x,y)

我们就可以列出方程

⎧⎩⎨2x2yxy=λy=λx=3

解得

λ=±2

所以{xy=3√=3√or{xy=−3√=−3√
此时f(x,y)=6，和我们观察出的答案一致。

怎么求多元函数梯度
我也不太会，但是如果这个多元函数可以分解为多个一元函数相加的话，我们求每个一元函数的导数，然后组成的向量就是多元函数的梯度了。比如∇f(x,y)=(2x,2y)，这对于这道题来说已经够用了。

使用拉格朗日乘数法解决这道题

我们再来看那两个式子

E=∑i=1Nkisi(vi−v′i)2T=∑i=1Nsivi

就是要求在E=∑Ni=1kisi(vi−v′i)2的约数条件下，Tmin
套用拉格朗日乘数法，列出方程。

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪−s1v21−s2v22⋯−snv2n=2λk1s1(v1−v′1)=2λk2s2(v2−v′2)=2λknsn(vn−v′n)

移一下项，可以得到2λkisi(vi−v′i)v2i=−1
由于{vi>0vi>v′i，显然可以得到λ<0。
2kisi(vi−v′i)v2i=−1λ，vi随λ的增大而增大。
令h(vi)=2λkisi(vi−v′i)v2i，h(vi)随vi的增大而减小。

我们可以先二分λ的大小，再根据h(vi)=−1二分求出vi的值，再根据此时的E来调整λ的大小。

代码

#include<cstdio>
#include<cmath>
#include<algorithm>
#define MAXN 10000
#define EPS 1e-7
using namespace std;
int n;
double s[MAXN+10],k[MAXN+10],v[MAXN+10],vv[MAXN+10],E,ans;
inline double sqr(double x){
    return x*x;
}
inline double cal(int i,double vi,double lambda){
    return vi*vi*2*lambda*k[i]*(vi-v[i]);
}
double Divide_Conqure(int i,double l,double r,double lambda){
    double mid;
    int t=0;
    while(++t<60){
        mid=(l+r)/2;
        if(cal(i,mid,lambda)<-1)
            r=mid;
        else
            l=mid;
    }
    return l;
}
double check(double lambda){
    int i;
    double En=0;
    for(i=1;i<=n;i++){
        vv[i]=Divide_Conqure(i,max(v[i],0.0),1e3,lambda);
        En+=k[i]*s[i]*sqr(vv[i]-v[i]);
    }
    return En;
}
double Divide_Conqure(double l,double r){
    double mid;
    int t=0;
    while(++t<70){
        mid=(l+r)/2;
        if(check(mid)<E)
            l=mid;
        else
            r=mid;
    }
    return l;
}
//　f(v1, v2, v3, ..., vn) = ∑ si / vi 
//　　g(v1, v2, v3, ..., vn) = ∑ kisi(vi-vi')^2
// (si'vi-sivi')/vi^2=-si/vi^2=2*lambda*kisi(vi-vi')
//vi^2*2*lambda*ki(vi-vi')和vi成反比&&=-1 
//vi^3*2*lambda*ki-vi^2*lambda*vi'*2*ki=-1/lambda vi和lambda成正比
int main()
{
    scanf("%d%lf",&n,&E);
    int i;
    for(i=1;i<=n;i++)
        scanf("%lf%lf%lf",s+i,k+i,v+i);
    Divide_Conqure(-1e3,0);
    for(i=1;i<=n;i++)
        ans+=s[i]/vv[i];
    printf("%.8lf\n",ans);
}