题目描述
分析
首先,我们来考虑一下部分分。
N=1
直接算即可, v1=Ek1×s1−−−−−√+v′1T=s1v1
N=2
v1=E1k1×s1−−−−−√+v′1v2=E−E1k2×s2−−−−−√+v′2T=s1v1+s2v2
我们猜测
T
关于
N≤100 & N≤1000
我不会QAQ
N≤10000
首先,我们将题目要求转化为数学式子
求 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
然后我们用拉格朗日乘数法试试。
设
我们就可以列出方程
解得
所以 {xy=3√=3√or{xy=−3√=−3√
此时 f(x,y)=6 ,和我们观察出的答案一致。
怎么求多元函数梯度
我也不太会,但是如果这个多元函数可以分解为多个一元函数相加的话,我们求每个一元函数的导数,然后组成的向量就是多元函数的梯度了。比如
∇f(x,y)=(2x,2y)
,这对于这道题来说已经够用了。
使用拉格朗日乘数法解决这道题
我们再来看那两个式子
就是要求在 E=∑Ni=1kisi(vi−v′i)2 的约数条件下, Tmin
套用拉格朗日乘数法,列出方程。
移一下项,可以得到 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);
}