坐标变换法求解求解无约束优化问题的例题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() # 显示模块中的所有绘图对象