坐标变换法求解求解无约束优化问题的例题python实现

题目描述：

设目标函数为：

取初始点为：

用坐标轮换法求解最优点（极值点）

解：

使用坐标轮换法进行求解无约束优化问题时，需要求解最优步长α，而求解最优步长首先确定函数搜索区间，在这里选择的时进退法进行求解搜索区间，然后用黄金分割法求解α 的近似最优解。

以下为python代码，我是根据这个函数只有x1和x2两个变量为前提自己编写的代码，不适用于三个及三个以上变量的优化问题。（我发现得出的结果并不是最优解，但是与最优解还是较为接近的，但是误差还是比较大（捂脸））

# 坐标轮换法（第一题）  编译环境：JupyterNotebook
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# 初始值的设置
x_0 = np.array([[-2.0,2.2]])  # 初始矩阵x_0
epoch = 100 # 算法的迭代次数


# 目标函数
def function(x):
    fx = (4 + (9.0/2.0)*(x[0][0]) - 4*(x[0][1]) + (x[0][0])**2 + 2*(x[0][1])**2 - 2*(x[0][0])*(x[0][1]) + (x[0][0])**4 - 
    2*((x[0][0])**2)*(x[0][1]))
    return fx

# 用进退法找到搜索区间中用到的函数，用于计算φ(α)
def fAlpha(x, a, judge):
    if(judge == 0):
        return (function(x + a*np.array([[1, 0]])))
    if(judge == 1):
        return (function(x + a*np.array([[0, 1]])))
    pass

#进退法求搜索区间
def SearchRegion(x, judge):  # x为变量的矩阵，judge为判断盖茨迭代为变量x1还是变量x2
    a_0 = 0
    h = 0.1
    a_1 = a_0
    a_2 = a_0 + h
    
    while(1):
        f1 = fAlpha(x, a_1, judge)
        f2 = fAlpha(x, a_2, judge)

        # 判断前进还是后退
        if(f2 < f1):
            a_3 = a_2 + h
            f3 = fAlpha(x, a_3, judge)
            
            #判断搜索区间
            if(f2 <= f3): # 满足高低高条件，直接输出搜索区间
                return np.array([a_1, a_3])
            if(f2 > f3): # 不满足高低高条件，继续搜索
                h = 2*h
                a_1 = a_2
                a_2 = a_3
        
        # 判断前进还是后退
        if(f2 >= f1):
            h = -h
            # 对调a_1和a_2
            t = a_1
            a_1 = a_2
            a_2 = t
            # 对调f1和f2
            t = f1
            f1 = f2
            f2 = t
            
            a_3 = a_2 + h
            f3 = fAlpha(x, a_3, judge)
            
            # 判断搜索区间
            if(f3 >= f2): # 满足高低高条件，直接输出搜索区间
                return np.array([a_3, a_1])
            if(f3 < f2): # 不满足高低高条件，继续搜索
                h = -2*h
                a_1 = a_2
                a_2 = a_3
        pass
    pass

# 用黄金分割法求最有步长
def GoldenSection(x, search_region, judge):
    a = search_region[0]
    b = search_region[1]
    
    a_1 = b - 0.618*(b - a)
    a_2 = a + 0.618*(b - a)
    
    E_golden_section = 0.01 # 黄金分割法收敛精度预设为0.01
    
    f1 = fAlpha(x, a_1, judge)
    f2 = fAlpha(x, a_2, judge)
    
    # 循环搜索最优步长
    while((abs((b-a)/b) > E_golden_section) & (abs((f2-f1)/f2) > E_golden_section)):
        if(f1 <= f2):
            b = a_2
            a_2 = a_1
            a_1 = b - 0.618*(b - a)
            
        if(f1 > f2):
            a = a_1
            a_1 = a_2
            a_2 = a + 0.618*(b - a)
        pass
    
    f1 = fAlpha(x, a_1, judge)
    f2 = fAlpha(x, a_2, judge)
    return ((a + b)/2)

# 坐标轮换法求最终结果
E_univariate_search = 0.01  # 坐标轮换法的精度预设为0.01
x = x_0
for i in range(epoch):
    x0 = x
    search_region = SearchRegion(x, 0)
    a = GoldenSection(x, search_region, 0)
    x = x + a*np.array([[1, 0]])
    search_region = SearchRegion(x, 1)
    a = GoldenSection(x, search_region, 1)
    x = x + a*np.array([[0, 1]])
    
    fanshu = x - x0
    if ((((fanshu[0][0])**2 + (fanshu[0][1])**2))**0.5 < E_univariate_search):
        break
    pass

print("最优解为：", x[0][0], "," ,x[0][1])
print("极值点为：" , function(x))

运行结果如下：
（不是很正确，因为我通过一个点一个点的试发现还有更小的解，不过相差不大）

以下是目标函数的图像，以及相应的python代码：

import matplotlib.pyplot as plt #绘图用的模块  
from mpl_toolkits.mplot3d import Axes3D #绘制3D坐标的函数  
import numpy as np  

# 目标函数
def fun(x,y):
    fx = 4 + (9.0/2.0)*x - 4*y + (x)**2 + 2*(y)**2 - 2*x*y + (x)**4 - 2*((x)**2)*y
    return fx
  
fig1 = plt.figure() # 创建一个绘图对象  
ax = Axes3D(fig1) # 用这个绘图对象创建一个Axes对象(有3D坐标)  
X = np.arange(-10, 10, 1)  
Y = np.arange(-20, 50, 1) # 创建了从-10到25，步长为0.1的arange对象  
# 至此X,Y分别表示了取样点的横纵坐标的可能取值  
# 用这两个arange对象中的可能取值一一映射去扩充为所有可能的取样点  

X,Y = np.meshgrid(X,Y)  
Z = fun(X,Y) # 用取样点横纵坐标去求取样点Z坐标  

plt.title("Image for the Function") # 总标题  
ax.plot_surface(X, Y, Z, rstride = 1, cstride=1, cmap=plt.cm.coolwarm) # 用取样点(x,y,z)去构建曲面  
ax.set_xlabel('x label', color = 'r') 
ax.set_ylabel('y label', color = 'g')  
ax.set_zlabel('z label', color = 'b') # 给三个坐标轴注明

plt.show() # 显示模块中的所有绘图对象