【Python】最速上升法求函数最大值

运用回归设计和响应曲面分析的方法来解决实验工作者通过实验寻找最佳的工艺条件，在实验者所选定的一个小区域上运用回归设计拟合线性回归方程，用最速上升法向最佳区域逼近因子区域，通常，将此法重复使用若干次，就可达到最佳区域。

利用最速上升法寻找最佳区域的步骤为：
(1)在变量 的某区域内，通过正交设计，拟合一个线性回归方程；
(2)由拟合的线性回归方程，找出最速上升路线；
(3)沿着这条上升路线进行一系列实验，直到y的值不再明显增大为止；
(4)在(3)中y值不再明显增大的那一点的邻域内重复(1)、(2)、(3)；
(5)当拟合的线性回归方程不再显著，y的曲面具有明显的弯曲时，拟合二次回归，进
而找出最佳条件。

# Method of Steepest Ascent 範例7.5, page 690


import math as m

def f(x,y):                                #定義函數f(x,y)
    return 4*x*y-m.pow(x,4)-m.pow(y,4)+4


def Gf(x,y):                               #定義f(x,y)的一次偏導數
    x1=4*y-4*m.pow(x,3)
    x2=4*x-4*m.pow(y,3)
    return [x1,x2]

def GGf(x,y):                              #定義f(x,y)的二次偏導數
    x1=-12*m.pow(x,2)
    x2=4
    x3=-12*m.pow(y,2)
    return [x1,x2,x3]

def h(t,x,y):              #定義f(x,y)沿著梯度方向的單變數函數h(t)
    u=Gf(x,y)
    return f(x+t*u[0],y+t*u[1])

def dh(t,x,y):             #定義h(t)的一次導函數
    u=Gf(x,y)
    v= Gf(x+t*u[0],y+t*u[1])
    return u[0]*v[0]+u[1]*v[1]

def ddh(t,x,y):            #定義h(t)的二次導函數
    u=Gf(x,y)
    x=x+t*u[0]
    y=y+t*u[1]
    v=GGf(x,y)
    return v[0]*m.pow(u[0],2)+2*v[1]*u[0]*u[1]+v[2]*m.pow(u[1],2)
    

x=4        #隨意起始點
y=3
max=-100.7    #隨意函數初值
m1=f(x,y)     #起始點含數值
while abs(max-m1)>0.0000001:      #Steepest Ascent method
    max=m1
    u=Gf(x,y)
    t=0.3
    s=t-dh(t,x,y)/ddh(t,x,y)
    while abs(s-t)>0.0000001:    #尋找dh的零根 (牛頓法)
       t=s
       s=t-dh(t,x,y)/ddh(t,x,y)
    
    x=x+s*u[0]   #新的起始點
    y=y+s*u[1]
    m1=f(x,y)

print(x)      #印出最大值發生的位置
print(y) 
print(max)

https://baike.baidu.com/item/%E6%9C%80%E9%80%9F%E4%B8%8A%E5%8D%87%E6%B3%95/19190718?fr=aladdin#4