题解-洛谷P2571 [SCOI2010]传送带

本文主要讲述了某蒟蒻历经一天一夜尝试 $O(1)$ 做法的一点点思考，以及对梯度下降法及各种优化算法较为详细的介绍

如有不足之处，或有更好的想法，望不吝赐教～

题目传送门

「洛谷P2571 [SCOI2010]传送带」

$O(1)$ 做法の尝试

推柿子

容易得到，最优路径一定是从 $A$ 移动到 $AB$ 上的某点 $M$，再由 $M$ 经平面移动到 $CD$ 上的某点 $N$，最后由 $N$ 移动到 $D$。

所花时间：

$$t=\frac{AM}{P}+\frac{MN}{R}+\frac{ND}{Q}$$

找到最优的 $M$、$N$ 即可。

直接设坐标比较麻烦，可以设线段的比值「取 $0$ ~ $1$」

设$AM=xAB,ND=yCD(x,y\in [0,1])$

则：

$$\begin{gathered} x_M-x_A=x(x_B-x_A),y_M-y_A=x(y_B-y_A) \\ {x}_M=x_A+x(x_B-x_A),y_M=y_A+x(y_B-y_A) \end{gathered}$$

同理：

$$\begin{gathered} x_N-x_D=y(x_C-x_D),y_N-y_D=y(y_C-y_D) \\ {x}_N=x_D+y(x_C-x_D),y_N=y_D+y(y_C-y_D) \end{gathered}$$

$$t=\frac{x\cdot AB}{P}+\frac{MN}{R}+\frac{y\cdot CD}{Q}$$

$$\begin{array}{rl}MN& =\sqrt{(x_M-x_N{)}^2+(y_M-y_N{)}^2}\\ & =\sqrt{(x_A+x(x_B-x_A)-x_D-y(x_C-x_D){)}^2+(y_A+x(y_B-y_A)-y_D-y(y_C-y_D){)}^2)}\end{array}$$

这一大串式子看起来很难受，不妨简化一下。

设：

$$c_1=\frac{AB}{P}{c}_2=\frac{CD}{Q}{c}_3=\frac{1}{R}$$

则

$$t=c_1\cdot x+c_2\cdot y+c_3\cdot MN$$

设：

$$\begin{array}{rlrl}& c_4=x_B-x_A& c_5=x_D-x_C& c_6=x_A-x_D\\ & c_7=y_B-y_A& c_8=y_D-y_C& c_9=y_A-y_D\end{array}$$

则

$$MN=\sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}$$

$$t=f(x,y)=c_1\cdot x+c_2\cdot y+c_3\cdot \sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}$$

c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9))

答案即为求解 $f(x,y)「x,y\in [0,1]」$ 的最小值。

对 $x$、$y$ 求偏导：

$$\begin{array}{rl}\frac{\partial f}{\partial x}& =c_1+c_3\cdot \frac{2\cdot (c_4\cdot x+c_5\cdot y+c_6)\cdot c_4+2\cdot (c_7\cdot x+c_8\cdot y+c_9)\cdot c_7}{2\sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}}\\ & =c_1+c_3\cdot \frac{(c_4\cdot x+c_5\cdot y+c_6)\cdot c_4+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_7}{\sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}}\\ \frac{\partial f}{\partial y}& =c_2+c_3\cdot \frac{2\cdot (c_4\cdot x+c_5\cdot y+c_6)\cdot c_5+2\cdot (c_7\cdot x+c_8\cdot y+c_9)\cdot c_8}{2\sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}}\\ & =c_2+c_3\cdot \frac{(c_4\cdot x+c_5\cdot y+c_6)\cdot c_5+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_8}{\sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}}\end{array}$$

令 $\frac{\partial f}{\partial x}=\frac{\partial f}{\partial y}=0$，解出所有驻点：

$$\{\begin{array}{l}\begin{array}{r}{c}_1+c_3\cdot \frac{(c_4\cdot x+c_5\cdot y+c_6)\cdot c_4+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_7}{\sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}}=0\\ c_2+c_3\cdot \frac{(c_4\cdot x+c_5\cdot y+c_6)\cdot c_5+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_8}{\sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}}=0\end{array}\end{array}$$

将 $c_1$、$c_2$ 移到等号右边，将 $\sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}$ 乘过去，两边开平方：

$$\{\begin{array}{l}\begin{array}{}\text{(1)}& c_3^2\cdot [(c_4\cdot x+c_5\cdot y+c_6)\cdot c_4+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_7{]}^2=c_1^2\cdot [(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2]\text{(2)}& c_3^2\cdot [(c_4\cdot x+c_5\cdot y+c_6)\cdot c_5+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_8{]}^2=c_2^2\cdot [(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2]\end{array}\end{array}$$

这样一来，方程组的形式就优美了许多，我们可以轻松求解。

用正则表达式稍作整理：

$$c{3}^2\ast ((c4\ast x+c5\ast y+c6)\ast c4+(c7\ast x+c8\ast y+c9)\ast c7{)}^2=c{1}^2\ast ((c4\ast x+c5\ast y+c6{)}^2+(c7\ast x+c8\ast y+c9{)}^2),c{3}^2\ast ((c4\ast x+c5\ast y+c6)\ast c5+(c7\ast x+c8\ast y+c9)\ast c8{)}^2=c{2}^2\ast ((c4\ast x+c5\ast y+c6{)}^2+(c7\ast x+c8\ast y+c9{)}^2)$$

塞到求解方程组的计算器中，即可求出 $x$、$y$ 的值：

$$\{\begin{array}{l}\begin{array}{rl}x& =\frac{c5\cdot c9-c6\cdot c8}{c4\cdot c8-c5\cdot c7}\\ y& =-\frac{c4\cdot c9-c6\cdot c7}{c4\cdot c8-c5\cdot c7}=\frac{c6\cdot c7-c4\cdot c9}{c4\cdot c8-c5\cdot c7}\end{array}\end{array}$$

x=(c5*c9-c6*c8)/(c4*c8-c5*c7),y=(c6*c7-c4*c9)/(c4*c8-c5*c7)

将 $x$、$y$ 代入原方程，可以利用表达式化简工具验证一下：

$$(c3{)}^2\ast (((c4)\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+(c5)\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+c6)\ast (c4)+((c7)\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+(c8)\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+(c9))\ast (c7){)}^2-((c1{)}^2\ast (((c4)\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+(c5)\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+(c6){)}^2+((c7)\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+(c8)\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+(c9){)}^2))$$

$$(c3{)}^2\ast (((c4)\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+(c5)\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+(c6))\ast (c5)+((c7)\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+(c8)\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+(c9))\ast (c8){)}^2-((c2{)}^2\ast (((c4)\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+(c5)\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+(c6){)}^2+((c7)\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+(c8)\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+(c9){)}^2))$$

上面两个式子的值都为 $0$，答案正确。

此时，$x$、$y$ 的取值即为 $f(x,y)$ 的极值点。

将 $x$、$y$ 代入 $f(x,y)$，得：

$$c1\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+c2\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+c3\ast sqrt((c4\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+c5\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+c6)\ast (c4\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+c5\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+c6)+(c7\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+c8\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+c9)\ast (c7\ast (c5\ast c9-c6\ast c8)/(c4\ast c8-c5\ast c7)+c8\ast (c6\ast c7-c4\ast c9)/(c4\ast c8-c5\ast c7)+c9))$$

化简，得：

$$((c1\ast c5-c2\ast c4)\ast c9-c1\ast c6\ast c8+c2\ast c6\ast c7)/(c4\ast c8-c5\ast c7)$$

即

$$f(x,y)=\frac{(c_1\cdot c_5-c_2\cdot c_4)\cdot c_9-c_1\cdot c_6\cdot c_8+c_2\cdot c_6\cdot c_7}{c_4\cdot c_8-c_5\cdot c_7}$$

((c1*c5-c2*c4)*c9-c1*c6*c8+c2*c6*c7)/(c4*c8-c5*c7)

这就完了？

当然没有。。。

首先，你不知道解出来的 $f(x,y)$ 是极大值还是极小值。

其次，就算是极小值，你也别忘了 $x,y\in [0,1]$ 的条件。

再者，边界值也有可能是合法区间内的最值。

因此，可以先判断 $x$，$y$ 的范围，若均在 $[0,1]$ 内则用 $f(x,y)$ 更新答案；然后用四个边界值 $f(0,0),f(0,1),f(1,0),f(1,1)$ 更新答案。

code:

#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))

using namespace std;

long double P,Q,R,ans=INT_MAX;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;

struct node{
    long double x,y;
}A,B,C,D,M,N;

inline long double dis(node a,node b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline long double f(long double x,long double y)
{
    return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}

signed main()
{
    scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
    c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
    c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
    x=(c5*c9-c6*c8)/(c4*c8-c5*c7),y=(c6*c7-c4*c9)/(c4*c8-c5*c7);
    if(x>=0&&x<=1&&y>=0&&y<=1) ans=min(ans,f(x,y));
    ans=Min(f(0,0),f(0,1),f(1,0),f(1,1),ans);
    printf("%.2Lf\n",ans);
    return 0;
}

特判——「一元」

若只想到这里，你只能获得 $60pts$。

还有一点需要处理：分母为 $0$ 的情况。

考虑 $x、y、f(x,y)$ 在 $c_4\cdot c_8-c_5\cdot c_7=0$ 时没有定义。

展开发现，此时 $(x_B-x_A)\cdot (y_D-y_C)=(x_D-x_C)\cdot (y_B-y_A)$。

注意 $AB$、$CD$ 可能为 $0$。

特判这两种情况：

设 $CD=0$，那么容易想到，最优路径一定是从 $A$ 移动到 $AB$ 上的某点 $M$，再由 $M$ 经平面移动到 $D$。

显然这是原先的一种特殊情况：$N$ 与 $C、D$ 重合

此时只需对原式略作修改即可。

所花时间：

$$t=\frac{AM}{P}+\frac{MD}{R}=\frac{x\cdot AB}{P}+\frac{MD}{R}$$

$$\begin{array}{rl}MD=MN& =\sqrt{(x_M-x_N{)}^2+(y_M-y_N{)}^2}\\ & =\sqrt{(x_A+x(x_B-x_A)-x_D-y(x_C-x_D){)}^2+(y_A+x(y_B-y_A)-y_D-y(y_C-y_D){)}^2)}\\ & =\sqrt{(x_A+x(x_B-x_A)-x_D{)}^2+(y_A+x(y_B-y_A)-y_D{)}^2)}\end{array}$$

容易发现，此时 $y$ 的取值不再影响结果，$t$ 变成了关于 $x$ 的一元函数。

重新用一下前面定义的几个参数：

$$c_1=\frac{AB}{P}{c}_2=\frac{CD}{Q}=0c_3=\frac{1}{R}$$

$$\begin{array}{rlrl}& c_4=x_B-x_A& c_5=x_D-x_C=0& c_6=x_A-x_D\\ & c_7=y_B-y_A& c_8=y_D-y_C=0& c_9=y_A-y_D\end{array}$$

$$MD=\sqrt{(c_4\cdot x+c_6{)}^2+(c_7\cdot x+c_9{)}^2}$$

$$t=f(x)=c_1\cdot x+c_3\cdot \sqrt{(c_4\cdot x+c_6{)}^2+(c_7\cdot x+c_9{)}^2}$$

答案即为求解 $f(x)「x\in [0,1]」$的最小值

直接求导：

$$\begin{array}{rl}\frac{\mathrm{d}f}{\mathrm{d}x}& =c_1+c_3\cdot \frac{2\cdot (c_4\cdot x+c_6)\cdot c_4+2\cdot (c_7\cdot x+c_9)\cdot c_7}{2\sqrt{(c_4\cdot x+c_6{)}^2+(c_7\cdot x+c_9{)}^2}}\\ & =c_1+c_3\cdot \frac{(c_4\cdot x+c_6)\cdot c_4+(c_7\cdot x+c_9)\cdot c_7}{\sqrt{(c_4\cdot x+c_6{)}^2+(c_7\cdot x+c_9{)}^2}}\end{array}$$

令 $\frac{\mathrm{d}f}{\mathrm{d}x}=0$，解出所有驻点：

$$c_1+c_3\cdot \frac{(c_4\cdot x+c_6)\cdot c_4+(c_7\cdot x+c_9)\cdot c_7}{\sqrt{(c_4\cdot x+c_6{)}^2+(c_7\cdot x+c_9{)}^2}}=0$$

稍作化简：

$$c_3^2\cdot [(c_4\cdot x+c_6)\cdot c_4+(c_7\cdot x+c_9)\cdot c_7{]}^2=c_1^2\cdot [(c_4\cdot x+c_6{)}^2+(c_7\cdot x+c_9{)}^2]$$

不就是个一元二次方程嘛，手解也能解出来「就是懒得解」

再次用正则表达式整理：

$$c{3}^2\ast ((c4\ast x+c6)\ast c4+(c7\ast x+c9)\ast c7{)}^2=c{1}^2\ast ((c4\ast x+c6{)}^2+(c7\ast x+c9{)}^2)$$

塞进那个工具里，然后就得到了两个解：

$$\begin{array}{rl}{x}_1& =\frac{\sqrt{c_3^2{c}_7^2+c_3^2{c}_4^2-c_1^2}(c_1{c}_4{c}_9-c_1{c}_6{c}_7)+((c_1^2-c_3^2{c}_4^2)c_7-c_3^2{c}_7^3)c_9-c_3^2{c}_4{c}_6{c}_7^2+(c_1^2{c}_4-c_3^2{c}_4^3)c_6}{c_3^2{c}_7^4+(2c_3^2{c}_4^2-c_1^2)c_7^2+c_3^2{c}_4^4-c_1^2{c}_4^2}\\ x_2& =-\frac{\sqrt{c_3^2{c}_7^2+c_3^2{c}_4^2-c_1^2}(c_1{c}_4{c}_9-c_1{c}_6{c}_7)+(c_3^2{c}_7^3+(c_3^2{c}_4^2-c_1^2)c_7)c_9+c_3^2{c}_4{c}_6{c}_7^2+(c_3^2{c}_4^3-c_1^2{c}_4)c_6}{c_3^2{c}_7^4+(2c_3^2{c}_4^2-c_1^2)c_7^2+c_3^2{c}_4^4-c_1^2{c}_4^2}\end{array}$$

x1=(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+((c1*c1-c3*c3*c4*c4)*c7-c3*c3*c7*c7*c7)*c9-c3*c3*c4*c6*c7*c7+(c1*c1*c4-c3*c3*c4*c4*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4)

x2=-(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+(c3*c3*c7*c7*c7+(c3*c3*c4*c4-c1*c1)*c7)*c9+c3*c3*c4*c6*c7*c7+(c3*c3*c4*c4*c4-c1*c1*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4)

同理，当 $AB=0$ 时，$M$ 与 $A、B$ 重合：

$$c_1=\frac{AB}{P}=0c_2=\frac{CD}{Q}{c}_3=\frac{1}{R}$$

$$\begin{array}{rlrl}& c_4=x_B-x_A=0& c_5=x_D-x_C& c_6=x_A-x_D\\ & c_7=y_B-y_A=0& c_8=y_D-y_C& c_9=y_A-y_D\end{array}$$

$$t=\frac{AN}{R}+\frac{ND}{Q}=\frac{AN}{R}+\frac{y\cdot CD}{Q}=c_2\cdot y+c_3\cdot AN$$

$$\begin{array}{rl}AN=MN& =\sqrt{(x_M-x_N{)}^2+(y_M-y_N{)}^2}\\ & =\sqrt{(x_A+x(x_B-x_A)-x_D-y(x_C-x_D){)}^2+(y_A+x(y_B-y_A)-y_D-y(y_C-y_D){)}^2)}\\ & =\sqrt{(x_A-x_D-y(x_C-x_D){)}^2+(y_A-y_D-y(y_C-y_D){)}^2)}\\ & =\sqrt{(c_5\cdot y+c_6{)}^2+(c_8\cdot y+c_9{)}^2)}\end{array}$$

$$t=f(y)=c_2\cdot y+c_3\cdot \sqrt{(c_5\cdot y+c_6{)}^2+(c_8\cdot y+c_9{)}^2)}$$

后面步骤相似，最后化简为：

$$c_3^2\cdot [(c_5\cdot y+c_6)\cdot c_5+(c_8\cdot y+c_9)\cdot c_8{]}^2=c_2^2\cdot [(c_5\cdot y+c_6{)}^2+(c_8\cdot y+c_9{)}^2]$$

$$c{3}^2\ast ((c5\ast y+c6)\ast c5+(c8\ast y+c9)\ast c8{)}^2=c{2}^2\ast ((c5\ast y+c6{)}^2+(c8\ast y+c9{)}^2)$$

$$\begin{array}{rl}{y}_1& =\frac{\sqrt{c_3^2{c}_8^2+c_3^2{c}_5^2-c_2^2}(c_2{c}_5{c}_9-c_2{c}_6{c}_8)+((c_2^2-c_3^2{c}_5^2)c_8-c_3^2{c}_8^3)c_9-c_3^2{c}_5{c}_6{c}_8^2+(c_2^2{c}_5-c_3^2{c}_5^3)c_6}{c_3^2{c}_8^4+(2c_3^2{c}_5^2-c_2^2)c_8^2+c_3^2{c}_5^4-c_2^2{c}_5^2}\\ y_2& =-\frac{\sqrt{c_3^2{c}_8^2+c_3^2{c}_5^2-c_2^2}(c_2{c}_5{c}_9-c_2{c}_6{c}_8)+(c_3^2{c}_8^3+(c_3^2{c}_5^2-c_2^2)c_8)c_9+c_3^2{c}_5{c}_6{c}_8^2+(c_3^2{c}_5^3-c_2^2{c}_5)c_6}{c_3^2{c}_8^4+(2c_3^2{c}_5^2-c_2^2)c_8^2+c_3^2{c}_5^4-c_2^2{c}_5^2}\end{array}$$

y1=(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+((c2*c2-c3*c3*c5*c5)*c8-c3*c3*c8*c8*c8)*c9-c3*c3*c5*c6*c8*c8+(c2*c2*c5-c3*c3*c5*c5*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5)

y2=-(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+(c3*c3*c8*c8*c8+(c3*c3*c5*c5-c2*c2)*c8)*c9+c3*c3*c5*c6*c8*c8+(c3*c3*c5*c5*c5-c2*c2*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5)

至于 $AB$、$CD$ 同时为 $0$ 的情况，只需计算 $A、D$ 间距离即可。

code:

#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))

using namespace std;

long double P,Q,R,ans=INT_MAX;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;

struct node{
    long double x,y;
}A,B,C,D,M,N,t;

inline long double dis(node a,node b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline long double f(long double x,long double y)
{
    return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}

signed main()
{
    scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
    c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
    c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
    if(c4*c8-c5*c7)
    {
        x=(c5*c9-c6*c8)/(c4*c8-c5*c7),y=(c6*c7-c4*c9)/(c4*c8-c5*c7);
        if(x>=0&&x<=1&&y>=0&&y<=1) ans=min(ans,f(x,y));
    }
    else
    {
        if(C.x==D.x&&C.y==D.y)
        {
            x=(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+((c1*c1-c3*c3*c4*c4)*c7-c3*c3*c7*c7*c7)*c9-c3*c3*c4*c6*c7*c7+(c1*c1*c4-c3*c3*c4*c4*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4);
            ans=min(ans,f(x,y));
            x=-(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+(c3*c3*c7*c7*c7+(c3*c3*c4*c4-c1*c1)*c7)*c9+c3*c3*c4*c6*c7*c7+(c3*c3*c4*c4*c4-c1*c1*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4);
            ans=min(ans,f(x,y));
        }
        else if(A.x==B.x&&A.y==B.y)
        {
            y=(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+((c2*c2-c3*c3*c5*c5)*c8-c3*c3*c8*c8*c8)*c9-c3*c3*c5*c6*c8*c8+(c2*c2*c5-c3*c3*c5*c5*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5);
            ans=min(ans,f(x,y));
            y=-(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+(c3*c3*c8*c8*c8+(c3*c3*c5*c5-c2*c2)*c8)*c9+c3*c3*c5*c6*c8*c8+(c3*c3*c5*c5*c5-c2*c2*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5);
            ans=min(ans,f(x,y));
        }
        else if(C.x==D.x&&C.y==D.y&&A.x==B.x&&A.y==B.y) ans=min(ans,dis(A,D));
    }
    ans=Min(f(0,0),f(0,1),f(1,0),f(1,1),ans);
    printf("%.2Lf\n",ans);   
    return 0;
}

特判——「平行」

若只想到这里，你还是只能获得$60pts$。。。

还有最后一种情况：即是 $AB、CD$ 不为 $0$，$AB//CD$ 也会使 $c_4\cdot c_8-c_5\cdot c_7=0$「因为$(x_B-x_A)\cdot (y_D-y_C)=(x_D-x_C)\cdot (y_B-y_A)$」

回到一开始的式子：

$$t=f(x,y)=c_1\cdot x+c_2\cdot y+c_3\cdot \sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}$$

似乎很难下手了……

或许我们可以换一种思路「乱搞」：充分利用几何知识

两条直线平行，我们可以直接考虑以下几种方案：

1. $A\to D$
2. 过 $C$ 向 $AB$ 作垂线，若与 $AB$ 有交点 $M$，则：$A\to M\to C\to D$
3. 过 $D$ 向 $AB$ 作垂线，若与 $AB$ 有交点 $M$，则：$A\to M\to D$
4. 过 $B$ 向 $CD$ 作垂线，若与 $CD$ 有交点 $N$，则：$A\to B\to N\to D$
5. 过 $A$ 向 $CD$ 作垂线，若与 $CD$ 有交点 $N$，则：$A\to N\to D$

还有可能两条线段在同一条直线上，需要特判一下。

对这几种方案分别取最小值，得到正确答案的概率就大了许多。

code:

#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))

using namespace std;

long double P,Q,R,ans=INT_MAX;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;

struct node{
    long double x,y;
}A,B,C,D,M,N,t;

inline long double dis(node a,node b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline long double f(long double x,long double y)
{
    return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}

inline bool eq_direction(node A1,node B1,node A2,node B2)//判断共线向量A1B1、A2B2是否同向
{
    if(B1.x-A1.x&&B2.x-A2.x) return (B1.x-A1.x)*(B2.x-A2.x)>0;
    else return (B1.y-A1.y)*(B2.y-A2.y)>0;
}

inline node foot_point(node P,node A,node B)//过P点作直线AB的垂线,返回垂足
{
    node res;
    //计算AB的解析式
    long double a=B.y-A.y,b=A.x-B.x,c=-a*A.x-b*A.y;
    res.x=(b*b*P.x-a*b*P.y-a*c)/(a*a+b*b);
    res.y=(a*a*P.y-a*b*P.x-b*c)/(a*a+b*b);
    return res;
}

inline bool on_line(node P,node A,node B)//判断在直线AB上的点P是否在线段AB上,只需判断P是否在A、B围成的矩形中即可
{
    return min(A.x,B.x)<=P.x&&P.x<=max(A.x,B.x)&&min(A.y,B.y)<=P.y&&P.y<=max(A.y,B.y);
}

signed main()
{
    scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
    c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
    c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
    if(c4*c8-c5*c7)
    {
        x=(c5*c9-c6*c8)/(c4*c8-c5*c7),y=(c6*c7-c4*c9)/(c4*c8-c5*c7);
        if(x>=0&&x<=1&&y>=0&&y<=1) ans=min(ans,f(x,y));
    }
    else
    {
        if(C.x==D.x&&C.y==D.y)
        {
            x=(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+((c1*c1-c3*c3*c4*c4)*c7-c3*c3*c7*c7*c7)*c9-c3*c3*c4*c6*c7*c7+(c1*c1*c4-c3*c3*c4*c4*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4);
            ans=min(ans,f(x,y));
            x=-(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+(c3*c3*c7*c7*c7+(c3*c3*c4*c4-c1*c1)*c7)*c9+c3*c3*c4*c6*c7*c7+(c3*c3*c4*c4*c4-c1*c1*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4);
            ans=min(ans,f(x,y));
        }
        else if(A.x==B.x&&A.y==B.y)
        {
            y=(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+((c2*c2-c3*c3*c5*c5)*c8-c3*c3*c8*c8*c8)*c9-c3*c3*c5*c6*c8*c8+(c2*c2*c5-c3*c3*c5*c5*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5);
            ans=min(ans,f(x,y));
            y=-(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+(c3*c3*c8*c8*c8+(c3*c3*c5*c5-c2*c2)*c8)*c9+c3*c3*c5*c6*c8*c8+(c3*c3*c5*c5*c5-c2*c2*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5);
            ans=min(ans,f(x,y));
        }
        //else if(C.x==D.x&&C.y==D.y&&A.x==B.x&&A.y==B.y) ans=min(ans,dis(A,D));
        else
        {
            ans=min(ans,dis(A,D)/R);
            if(c4*(C.y-B.y)==c7*(C.x-B.x))//已知AB//CD,若AB//BC,则A、B、C、D四点共线
            {
                ans=min(ans,dis(A,B)/P+dis(B,C)/R+dis(C,D)/Q);//ABCD
                ans=min(ans,dis(A,B)/P+dis(B,D)/R);//ABD
                ans=min(ans,dis(A,C)/R+dis(C,D)/Q);//ACD
                if(!eq_direction(A,C,B,C)&&!eq_direction(A,D,B,D))
                {
                    if(eq_direction(A,C,C,D)) ans=min(ans,dis(A,C)/P+dis(C,D)/max(P,Q));//ACDB
                    else ans=min(ans,dis(A,D)/P);//ADCB
                }
                if(!eq_direction(A,C,A,D)&&!eq_direction(B,C,B,D))
                {
                    if(eq_direction(C,A,A,B)) ans=min(ans,dis(A,B)/max(P,Q)+dis(B,D)/Q);//CABD
                    else ans=min(ans,dis(A,D)/Q);//CBAD
                }
            }
            else//不共线
            {
                t=foot_point(C,A,B);
                if(on_line(t,A,B)) ans=min(ans,dis(A,t)/P+dis(t,C)/R+dis(C,D)/Q);//AMCD
                t=foot_point(D,A,B);
                if(on_line(t,A,B)) ans=min(ans,dis(A,t)/P+dis(t,D)/R);//AMD
                t=foot_point(B,C,D);
                if(on_line(t,C,D)) ans=min(ans,dis(A,B)/P+dis(B,t)/R+dis(t,D)/Q);//ABND
                t=foot_point(A,C,D);
                if(on_line(t,C,D)) ans=min(ans,dis(A,t)/R+dis(t,D)/Q);//AND
            }
        }
    }
    ans=Min(f(0,0),f(0,1),f(1,0),f(1,1),ans);
    printf("%.2Lf\n",ans);   
    return 0;
}

瓶颈

emmm...

仍是$60pts...$

乱搞终究不是正道。。。

不难发现，$O(1)$做法的瓶颈卡在了$AB//CD$使得求导后分母为$0$的情况。

若有大佬可以通过调整参数、用更好的方法求最值的方式来绕过这一特殊情况，请不吝赐教～

梯度下降法

$O(1)$：我死了QAQ……

再次回到起点：

$$t=f(x,y)=c_1\cdot x+c_2\cdot y+c_3\cdot \sqrt{(c_4\cdot x+c_5\cdot y+c_6{)}^2+(c_7\cdot x+c_8\cdot y+c_9{)}^2}$$

非 $O(1)$ 复杂度求解多元函数最值的方法有很多：粒子群算法、模拟退火、三分套三分、牛顿迭代法……

在此介绍梯度下降法。

梯度

了解多元微积分的各位大佬们都知道，梯度是一个向量，指向多元函数的值变化最快的方向，大小即为变化率。

「不了解的可以看看 3Blue1Brown 的作者 Grant Sanderson 讲的多元微积分课程，在b站上已有熟肉【链接】」

百度百科对梯度的定义：

设二元函数$z=f(x,y)$ 在平面区域$D$上具有一阶连续偏导数，则对于每一个点$P(x,y)$都可定出一个向量

$$\{\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\}=f_x(x,y)\hat{i}+f_y(x,y)\hat{j}$$

该函数「也不是不对，但称为向量更好理解」就称为函数$z=f(x,y)$在点$P(x,y)$的梯度，记作$\text{grad}f(x,y)$或$\mathrm{\nabla }f(x,y)$。

于是有：

$$\text{grad}f(x,y)=\mathrm{\nabla }f(x,y)=\{\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\}=f_x(x,y)\hat{i}+f_y(x,y)\hat{j}$$

其中

$$\mathrm{\nabla }=\frac{\partial }{\partial x}\hat{i}+\frac{\partial }{\partial y}\hat{j}$$

称为（二维的）向量微分算子或Nabla算子，且有

$$\mathrm{\nabla }f=\frac{\partial f}{\partial x}\hat{i}+\frac{\partial f}{\partial y}\hat{j}$$

梯度下降法的主要思想

再来介绍梯度下降法。

正如我们先前对 $f(x,y)$ 求偏导，令偏导值为 $0$，解出所有驻点的操作，实际上就是试图直接令梯度值为 $0$。

但是也正如我们遇到的问题，有时直接求解比较困难，我们不妨通过迭代，让梯度值逐渐下降，即逐渐接近极值点。

例如，对于函数 $f(x)=x^2$，求导得 $f^{\prime }(x)=2x$。

不妨设 $x_0=1$，则从点 $(1,1)$ 开始，计算其梯度值：$f^{\prime }(1)=2>0$。

梯度值 $>0$，说明在该点沿使 $x$ 增大的方向「即 $x$ 轴正半轴方向」，$f(x)$ 函数值会增大，增大的速率为 $2$。

容易发现，在 $(3,9)$ 处，$f^{\prime }(3)=6$，增大的速率为 $6$，观察图像也可得知，沿 $x$ 轴正半轴方向，$f(x)$ 的图像变陡，上升速率变快。

而在 $(0,0)$ 处，$f^{\prime }(0)=0$，函数图像在这一点的切线斜率是平的，增大的速率无限接近 $0$，此时梯度值为 $0$，$(0,0)$ 是函数 $f(x)$ 的一个驻点，同时 $f(x)$ 也在此取到最小值。

因此，我们有如下策略：

梯度值 $>0$，向左移动一点；梯度值 $<0$，向右移动一点。

即每次让 $x$ 朝与梯度值符号相反的方向移动，使梯度值逐渐下降，最终趋于 $0$。

用数学语言来描述，即：

$$x_{k+1}=x_k-\mathrm{\nabla }f(x_k)$$

其中 $x_k$ 为第 $k$ 次迭代时点的横坐标，$x_{k+1}$ 为第 $k+1$ 次迭代移动到的点的横坐标，$x_0$ 表示初始横坐标。

$\mathrm{\nabla }f(x)$ 表示在 $x$ 处的梯度值，$\mathrm{\nabla }f(x)=\frac{\text{d}f}{\text{d}x}$。

$x_k-\mathrm{\nabla }f(x_k)$ 表示向与梯度值符号相反的方向移动，符合之前的策略

优化——「学习率」

但这样就有一个问题：

例如 $f(x)=x^2,x_0=10$，则 $\mathrm{\nabla }f(x)=2x$：

$$\begin{array}{rl}{x}_1& =x_0-2\cdot x_0=-10\\ x_2& =x_1-2\cdot x_1=10\\ x_3& =x_2-2\cdot x_2=-10\\ ...\end{array}$$

可以看到，虽然横坐标一直在变化，但一直在 $10、-10$ 之间振荡，函数值始终没有降到最低点。

为此，我们需要一个参数控制移动的距离，这个参数被称作学习率，用 $\eta$ 表示：

$$x_{k+1}=x_k-\eta \cdot \mathrm{\nabla }f(x_k)$$

显然，在刚刚的计算中，$\eta =1$。

当 $\eta >1$ 时，函数值不仅不会降到最低点，甚至会越来越大：

$e.g.\eta =1.05$

$$\begin{array}{rl}{x}_1& =x_0-2\cdot x_0=-11\\ x_2& =x_1-2\cdot x_1=12.1\\ x_3& =x_2-2\cdot x_2=-13.31\\ ...\end{array}$$

而若将 $\eta$ 调低到一个较小值，如 $0.02$：

$$\begin{gathered} x_1=x_0-2\cdot x_0=9.6 \\ {x}_2=x_1-2\cdot x_1=9.2 \\ {x}_3=x_2-2\cdot x_2=8.8 \\ ... \end{gathered}$$

降低的速度太慢，容易超时。

当 $\eta =0.2$ 时，迭代 $10$ 次左右即可降入谷底：

$$\begin{array}{rlr}{x}_1& =x_0-2\cdot x_0=6& ||\mathrm{\nabla }f(x_0)||=12\\ x_2& =x_1-2\cdot x_1=3.6& ||\mathrm{\nabla }f(x_1)||=7.2\\ x_3& =x_2-2\cdot x_2=2.16& ||\mathrm{\nabla }f(x_2)||=4.32\\ ...\end{array}$$

可以看到，右边的梯度值在每次迭代后都会下降，故称为梯度下降法「Gradient Descent」

只要选择合适的学习率，梯度就可以下降到任意小：

$$\underset{k\to \mathrm{\infty }}{lim}f(x_k)=minf(x)$$

可以用泰勒公式进行严格证明，此处不再赘述（逃

因此可以通过梯度值的大小作为终止条件「也可以直接用迭代次数控制精度」

对于二元函数，同样可以用梯度下降法求解极值：

$$f(x_{k+1},y_{k+1})=f(x_k,y_k)-\eta \cdot \mathrm{\nabla }f(x_k,y_k)$$

$e.g.f(x,y)=x^2+2y^2,(x_0,y_0)=(-3.5,-3.5),\eta =0.1$，则$\mathrm{\nabla }f(x,y)=(2x,4y)$

$$\begin{gathered} (x_1,y_1)=(x_0,y_0)-\eta \cdot \mathrm{\nabla }f(x_0,y_0)=(-2.8,-2.1) \\ (x_2,y_2)=(x_1,y_1)-\eta \cdot \mathrm{\nabla }f(x_1,y_1)=(-2.24,-1.26) \\ ... \end{gathered}$$

代码——「梯度下降法」

code：

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

inline long double f(long double x)
{
    return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}

inline long double numerical_diff(long double x)//数值微分法估计一阶导数
{
    long double dx=1e-6;
    return (f(x+dx)-f(x-dx))/(dx*2);
}

inline long double gradient_descent(long double x,long double eta)//梯度下降法
{
    long double eps=1e-6;
    //int cnt=0;
    while(abs(numerical_diff(x))>eps)
    {
        x-=eta*numerical_diff(x);
        //cnt++;
    }
    //cout<<cnt<<endl;
    return f(x);
}

signed main()
{
    cout<<gradient_descent(0,0.02)<<endl;
    return 0;
}

$e.g.f(x)=3x^4-x^3+2x^2-9x+5\sqrt{(x+3{)}^2+(5x+6{)}^2}-25$

经过$9$次迭代达到目标精度。

优化——动量梯度下降法「MGD」

需要指出的是，梯度值为 $0$ 的驻点不一定是函数的极值点，如 $f(x)=x^3$ 在 $x=0$ 处梯度值为 $0$，但并不是函数的极值点：

「实际上，$(0,0)$ 是 $f(x)$ 的拐点，此处函数的凹凸性发生改变」

同时，对于非单峰函数来说，梯度下降法的结果易受到初始值的影响，也就是陷入局部最优解。

现实生活中，一个小球从高处落下，大概率会越过比较低的坎继续下降。

我们同样可以引入惯性来优化：

用梯度模拟受力，使之不直接控制移动距离，而是给小球一个速度 $v$；同时可以引入阻力，也就是速度衰减率 $\beta$ 使小球减速。

于是，整个过程就像一个有动量的小球在空间中来回滚动，故称动量梯度下降法「Gradient Descent with Momentum，简称MGD」

根据动量定理

$$F\mathrm{\Delta }t=m\mathrm{\Delta }v$$

移项，得：

$$\mathrm{\Delta }v=\frac{\mathrm{\Delta }t}{m}\cdot F$$

直接令 $\eta =\frac{\mathrm{\Delta }t}{m}$，则：

$$\mathrm{\Delta }v=\eta \cdot F$$

但在实际编程中，不必那么精细地刻画阻力，直接用一个速度衰减率 $\beta$ 代替即可：

$$v_{k+1}=\beta \cdot v_k-\eta \cdot \mathrm{\nabla }f(x_k)$$

code：

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

inline long double f(long double x)
{
    return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}

inline long double numerical_diff(long double x)
{
    long double dx=1e-6;
    return (f(x+dx)-f(x-dx))/(dx*2);
}

inline long double mgd(long double x,long double eta,long double beta)//动量梯度下降法
{
    long double v=0,eps=1e-6;
    int cnt=0;
    while(abs(numerical_diff(x))>eps)
    {
        v=beta*v-eta*numerical_diff(x);
        x+=v;
        cnt++;
    }
    cout<<cnt<<endl;
    return f(x);
}

signed main()
{
    cout<<mgd(0,0.03,0.05)<<endl;
    return 0;
}

迭代次数：$13$。

「多峰函数的优化效果较明显」

优化——自适应梯度算法「AdaGrad」

对比一下两种算法的真实移动距离 $\mathrm{\Delta }{x}_i$：

	梯度下降法	动量梯度下降法
$\mathrm{\Delta }{x}_i$	$\eta \cdot \frac{\partial f}{\partial x_i}$	$\eta \cdot v_i$

动量梯度下降法只优化了 $\frac{\partial f}{\partial x_i}$ 的部分。

很自然的想到，是否可以优化学习率 $\eta$？

答案是肯定的，这种算法称为自适应梯度算法「AdaGrad【Adapt+Gradient】」，基本思想是开始时快速移动接近目标，然后减速提高精度「梯度对其产生的动力效果逐渐减弱」。实现时记录每次梯度下降的梯度平方和 $h$，将 $\eta$ 除以 $\sqrt{h}$ 进行调整「为防止分母出现 $0$，往往加上一个较小的常量」

code：

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

inline long double f(long double x)
{
    return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}

inline long double numerical_diff(long double x)
{
    long double dx=1e-6;
    return (f(x+dx)-f(x-dx))/(dx*2);
}

inline long double AdaGrad(long double x,long double eta,long double beta)//自适应梯度算法
{
    long double v=0,h=0,grad=numerical_diff(x),eps=1e-6;
  	//int cnt=0;
    while(abs(grad)>eps)
    {
        h+=grad*grad;
        v=beta*v-eta*grad/(sqrt(h)+eps);
        x+=v;
      	grad=numerical_diff(x);
        //cnt++;
    }
    //cout<<cnt<<endl;
    return f(x);
}

signed main()
{
    cout<<AdaGrad(0,0.5,0.05)<<endl;
    return 0;
}

迭代次数：$12$。

「多峰函数的优化效果较明显」

优化——自适应动量算法「Adam」

自适应梯度算法的缺陷：若初始点离极值点较远，可能还没到极值点，小球就跑不动了「相应的有 RMSProp 优化算法，通过把 $h$ 乘以一个衰减率 ${\beta }_2$ 来逐步遗忘之前的梯度，专业一点说是“指数移动平均”，与动量梯度下降法有异曲同工之妙，此处不作详细介绍」

我们也可以将动量梯度下降法「MGD」与自适应梯度算法「AdaGrad」融合在一起，得到自适应动量算法「Adaptive Momentum Estimation，简称 Adam」

在梯度下降过程中，Adam 算法用梯度模拟小球受力改变速度 $v$，同时会增大小球的质量而改变小球移动的难易程度，二者的作用最终影响移动距离。

code：

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

inline long double f(long double x)
{
    return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}

inline long double numerical_diff(long double x)
{
    long double dx=1e-6;
    return (f(x+dx)-f(x-dx))/(dx*2);
}

inline long double Adam(long double x,long double eta,long double beta1,long double beta2)
{
    long double v=0,m=0,grad=numerical_diff(x),t,t1=beta1,t2=beta2,eps=1e-6;
    //int cnt=0;
    while(abs(grad)>eps)
    {
        t=eta*sqrt(1-t2)/(1-t1);
        v+=(1-beta1)*(grad-v);
        m+=(1-beta2)*(grad*grad-m);
        x-=t*v/(sqrt(m)+eps);
        t1*=beta1,t2*=beta2;
        grad=numerical_diff(x);
        //cnt++;
    }
    //cout<<cnt<<endl;
    return f(x);
}

signed main()
{
    cout<<Adam(0,0.5,0.6,0.9999)<<endl;
    return 0;
}

迭代次数：$49$。

综合来看，AdaGrad 的算法似乎更具优势，但毕竟各个算法特点不同，适用的函数也不同，用哪种算法还是应根据实际应用而定。

二维写法

code：

struct node{
    long double x,y;
};

inline long double f(long double x,long double y)
{
    return ...;
}

inline long double part_x(long double x,long double y)
{
    long double dx=1e-6;
    return (f(x+dx,y)-f(x-dx,y))/(dx*2);
}

inline long double part_y(long double x,long double y)
{
    long double dy=1e-6;
    return (f(x,y+dy)-f(x,y-dy))/(dy*2);
}

inline node gd(long double x,long double y,long double eta,long double beta)
{
    long double eps=1e-6;
    node res;
    int cnt=0,T=100;
    while(T--&&abs(part_x(x,y))+abs(part_y(x,y))>eps)
    {
        x-=eta*part_x(x,y);
        y-=eta*part_y(x,y);
        //cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
    }
    res.x=x,res.y=y;
    return res;
}

inline node mgd(long double x,long double y,long double eta,long double beta)
{
    long double vx=0,vy=0,eps=1e-6;
    node res;
    int cnt=0,T=100;
    while(T--&&abs(part_x(x,y))+abs(part_y(x,y))>eps)
    {
        vx=beta*vx-eta*part_x(x,y);
        vy=beta*vy-eta*part_y(x,y);
        x+=vx,y+=vy;
        //cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
    }
    res.x=x,res.y=y;
    return res;
}

inline node AdaGrad(long double x,long double y,long double eta,long double beta)
{
    long double vx=0,vy=0,hx=0,hy=0,gradx=part_x(x,y),grady=part_y(x,y),grad=abs(gradx)+abs(grady),eps=1e-6;
    node res;
    int cnt=0,T=100;
    while(T--&&grad>1e-5)
    {
        hx+=gradx*gradx;
        hy+=grady*grady;
        vx=beta*vx-eta*gradx/(sqrt(hx)+eps);
        vy=beta*vy-eta*grady/(sqrt(hy)+eps);
        x+=vx,y+=vy;
        gradx=part_x(x,y);
        grady=part_y(x,y);
        grad=abs(gradx)+abs(grady);
        //cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
    }
    res.x=x,res.y=y;
    return res;
}

inline node Adam(long double x,long double y,long double eta,long double beta1,long double beta2)
{
    long double vx=0,vy=0,mx=0,my=0,gradx=part_x(x,y),grady=part_y(x,y),grad=abs(gradx)+abs(grady),t,t1=beta1,t2=beta2,eps=1e-7;
    node res;
    int cnt=0,T=100;
    while(T--&&grad>1e-5)
    {
        t=eta*sqrt(1-t2)/(1-t1);
        vx+=(1-beta1)*(gradx-vx);
        vy+=(1-beta1)*(grady-vy);
        mx+=(1-beta2)*(gradx*gradx-mx);
        my+=(1-beta2)*(grady*grady-my);
        x-=t*vx/(sqrt(mx)+eps);
        y-=t*vy/(sqrt(my)+eps);
        t1*=beta1,t2*=beta2;
        gradx=part_x(x,y);
        grady=part_y(x,y);
        grad=abs(gradx)+abs(grady);
        //cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
    }
    res.x=x,res.y=y;
    return res;
}

本题代码

回到此题，直接用 Adam 交了一发，AC～「可能运气较好……」

code：

#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#include<cstdlib>
#include<ctime>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))

using namespace std;

long double P,Q,R,ans=INT_MAX;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;

struct node{
    long double x,y;
}A,B,C,D,M,N,t;

inline long double dis(node a,node b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline long double f(long double x,long double y)
{
    return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}

inline long double part_x(long double x,long double y)
{
    long double dx=1e-6;
    return (f(x+dx,y)-f(x-dx,y))/(dx*2);
}

inline long double part_y(long double x,long double y)
{
    long double dy=1e-6;
    return (f(x,y+dy)-f(x,y-dy))/(dy*2);
}

inline node Adam(long double x,long double y,long double eta,long double beta1,long double beta2)
{
    long double vx=0,vy=0,mx=0,my=0,gradx=part_x(x,y),grady=part_y(x,y),grad=abs(gradx)+abs(grady),t,t1=beta1,t2=beta2,eps=1e-6;
    node res;
    int T=20;
    while(T--&&grad>1e-4)
    {
        t=eta*sqrt(1-t2)/(1-t1);
        vx+=(1-beta1)*(gradx-vx);
        vy+=(1-beta1)*(grady-vy);
        mx+=(1-beta2)*(gradx*gradx-mx);
        my+=(1-beta2)*(grady*grady-my);
        x-=t*vx/(sqrt(mx)+eps);
        y-=t*vy/(sqrt(my)+eps);
        t1*=beta1,t2*=beta2;
        gradx=part_x(x,y);
        grady=part_y(x,y);
        grad=abs(gradx)+abs(grady);
    }
    res.x=x,res.y=y;
    return res;
}

inline long double Rand()
{
    return abs(rand()*rand()*rand()%(10000000)/(long double)(10000000));
}

inline void gradient()
{
    for(int T=1;T<=200000;T++)
    {
        t=Adam(Rand(),Rand(),Rand(),Rand(),Rand());
        if(t.x>=0&&t.x<=1&&t.y>=0&&t.y<=1) ans=min(ans,f(t.x,t.y));
    }
}

signed main()
{
    srand(time(0));
    scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
    c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
    c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
    gradient();
    ans=Min(f(0,0),f(0,1),f(1,0),f(1,1),ans);
    printf("%.2Lf\n",ans);
    return 0;
}

三分套三分法

突然发现三分套三分跑的贼快，于是也打算简单实现一下。

模板就不讲了，直接贴一下之前的博客：「二分法&三分法模板」

上代码：

#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))

using namespace std;

long double P,Q,R,ans=INT_MAX,eps=1e-6;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;

struct node{
    long double x,y;
}A,B,C,D,M,N,t;

inline long double dis(node a,node b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline long double f(long double x,long double y)
{
    return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}

inline long double g(long double x)
{
    long double l=0,r=1,mid;
    while(r-l>eps)
    {
        mid=(l+r)/2;
        if(f(x,mid-eps)<f(x,mid+eps)) r=mid;
        else l=mid;
    }
    return f(x,mid);
}

signed main()
{
    scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
    c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
    c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
    long double l=0,r=1,mid;
    while(r-l>eps)
    {
        mid=(l+r)/2;
        if(g(mid-eps)<g(mid+eps)) r=mid;
        else l=mid;
    }
    printf("%.2Lf\n",g(mid));
    return 0;
}

$$Thanksforreading.$$

————THE——END————