最优化实验报告
——外点、内点和混合罚函数法
实验目的
之前我们已经实验过无约束最优化问题，这次我们将实验一下，在有约束条件下，优化算法应该怎么做。
由于处理约束条件的办法不同，约束优化法可以分为直接间接两类。间接法的思想是将问题转化成无约束的，然后无约束方法求解，虽然算法复杂，但可以求解高维问题，效率高、鲁棒性好。直接法构造迭代过程，使迭代点在可行域D中，一步一步降低目标函数值。算法简单，但是计算量大，一般只用于求解不等式约束问题。
罚函数法又称乘子法，是指将有约束最优化问题转化为求解无约束最优化问题：其中M为足够大的正数， 起"惩罚"作用， 称之为罚因子，F(x, M )称为罚函数。内部罚函数法也称为障碍罚函数法。
这种方法是在可行域内部进行搜索，约束边界起到类似围墙的作用，如果当前解远离约束边界时，则罚函数值是非常小的，否则罚函数值接近无穷大的方法。在进化计算中，研究者选择外部罚函数法的原因主要是该方法不需要提供初始可行解。其中B(x)是优化过程中新的目标函数，Gi和Hj分别是约束条件gi(x)和hj(x)的函数，ri和cj是常数，称为罚因子。
本次实验将验证间接法中的内点、外点、混合罚函数法。
实验环境
本次实验，使用的程序语言为Python3.5.2，使用符号系统解析法推导计算，依赖的工具包为：
import math
import sympy
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
要想在IDE（集成开发环境）中实现画图的交互式效果，还得写
plt.ion()
实验内容及步骤
- 定义函数和符号
- 求解最优X
- 画出构造函数变化图
1. 外点罚函数法基本原理与步骤
问题：
minf(X)=−x1+x2
s.t.{g1(X)=lnx2≥0h1(X)=x1+x2−1≥0
构造一个增广函数：
F(X,Mk)=f(X)+Mkα(X)
其中惩罚函数：
α(X)=∑j=1m[hj(X)]2+∑i=1l[gi(X)]2u(gi(X))
α(X){0,当X∈D时,0,当X∉D时.
u(gi(X))={0,当gi(X)≥0时,1,当gi(X)<0时.i=1,2,⋯,l.
Mk
是一个逐渐增大的参数，称为惩罚因子。
当X∈D
时
F(X,Mk)=f(X)
算法流程图：
- 内点罚函数法基本原理与步骤
问题：
x(1+1)3+x213min
s.t.{g1(X)=x1−1≥0,g2(X)=x2≥0.
构造一个增广函数：
F(X,rk)=f(X)+rk∑i=1l1gi(X)
其中rk
称为障碍因子。
算法流程图：
- 混合罚函数法基本原理与步骤
问题：
minf(X)=−x1+x2
s.t.{g1(X)=lnx2≥0h1(X)=x1+x2−1≥0
构造一个增广函数：
F(X,rk)=f(X)+rk∑i=1l1gi(X)+1rk∑j=1m[hj(X)]2
算法流程图：
实验结果及分析
- 外点罚函数法：
起初，实验发生了错误，使用数值计算的方法做最速下降收敛增广函数，在接近理论最优解时，梯度爆炸，变得很大，经过多次参数调整，也只能求得不精确（+-0.02）的解。
然后使用了 符号解析法 ，配合 最优步长法 ，取得了极好的效果。
起始点：
迭代中逐渐增大非可行域的值（Z轴数值F明显变化，最优解也在变）：
经过几十次迭代，就取得最后结果：
- 内点罚函数法
吸取调参经验教训，结合间接求解有约束问题能转化成无约束问题，并且函数显然在定义域内连续可导。所以，这次直接用了解非线性方程组，求极值取极限的方法。
很快取得了结果：
- 混合罚函数法
复用上一题程序，只需改进解的可行范围检查即可得：
程序
失败的数值计算方法——会爆炸.py
'''
失败的数值方法
'''
import sympy
import math
import numpy as np
C = 10 # 放大系数
M = 1 # 惩罚因子
epsilon = 1e-5 # 终止限
f=lambda x:-x[0] + x[1]
h=lambda x:x[0] + x[1] - 1
g=lambda x:max(-math.log(x[1]),0)
alpha=lambda x:h(x)**2+g(x)**2
F=lambda x:f(x)+M*alpha(x)
# 梯度下降来最小化F
def GD(x):
delta_x = 1e-11 # 数值求导
t = 0.0001 # 步长
# t=sympy.symbols('t')
e = 0.01 # 极限
# my_print(e)
np.array(x)
for i in range(500):
grad = np.asarray([(F([x[0] + delta_x, x[1]]) - F(x)) / delta_x, (F([x[0], x[1] + delta_x]) - F(x)) / delta_x])
# print('x',x)
print('g',grad)
# t=sympy.solve(sympy.diff(F(x-t*grad),t),t)
x = x - t * grad
# my_print(np.linalg.norm(grad))
if (np.linalg.norm(grad) < e):
# print(np.linalg.norm(grad))
# print('g', grad)
break
return list(x)
pass
x = [-0.5, 0.2]
for i in range(20):
x = GD(x)
print('-------', i, x)
if (M * alpha(x) <= epsilon):
print(alpha(x))
break
else:
M = C * M
print(x, f(x))
外点罚.py
'''
这是有用的
外点罚函数
'''
import math
import sympy
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
plt.ion()
fig = plt.figure()
ax = Axes3D(fig)
def draw(x,index,M):
# F = f + MM * alpha
# FF = sympy.lambdify((x1, x2), F, 'numpy')
Z = FF(*(X, Y,M))
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow',alpha=0.5)
ax.scatter(x[0], x[1], FF(*(x[0],x[1],M)), c='r',s=80)
ax.text(x[0], x[1], FF(*(x[0],x[1],M)), 'here:(%0.3f,%0.3f)' % (x[0], x[1]))
ax.set_zlabel('F') # 坐标轴
ax.set_ylabel('X2')
ax.set_xlabel('X1')
plt.pause(0.1)
# plt.show()
# plt.savefig('./image/%03d' % index)
plt.cla()
C = 10 # 放大系数
M = 1 # 惩罚因子
epsilon = 1e-5 # 终止限
x1, x2 = sympy.symbols('x1:3')
MM=sympy.symbols('MM')
f = -x1 + x2
h = x1 + x2 - 1
# g=sympy.log(x2) if sympy.log(x2)<0 else 0
g = sympy.Piecewise((x2-1, x2 < 1), (0, x2 >= 1))
# u=lambda x:
alpha = h ** 2 + g ** 2
F = f + MM * alpha
# 梯度下降来最小化F
def GD(x,M,n):
# F = f + M * alpha
# delta_x = 1e-11 # 数值求导
# t = 0.0001 # 步长
e = 0.001 # 极限
# my_print(e)
np.array(x)
for i in range(15):
t = sympy.symbols('t')
grad = np.asarray(
[sympy.diff(F, x1).subs([(x1, x[0]), (x2, x[1]),(MM,M)]), sympy.diff(F, x2).subs([(x1, x[0]), (x2, x[1]),(MM,M)])])
# print('g',grad)
# print((x-t*grad))
# print(F.subs([(x1,(x-t*grad)[0]),(x2,(x-t*grad)[1])]))
t = sympy.solve(sympy.diff(F.subs([(x1, (x - t * grad)[0]), (x2, (x - t * grad)[1]),(MM,M)]), t), t)
print('t',t)
x = x - t * grad
print('x', x)
# print('mmm',M)
draw(x,n*10+i,M)
# my_print(np.linalg.norm(grad))
# print(type(grad))
if (abs(grad[0]) < e and abs(grad[1]) < e):
# print(np.linalg.norm(grad))
print('g', grad)
break
return list(x)
pass
x = [-0.5, 0.2]
X = np.arange(0, 4, 0.25)
Y = np.arange(0, 4, 0.25)
X, Y = np.meshgrid(X, Y)
FF = sympy.lambdify((x1, x2,MM), F, 'numpy')
# ff=sympy.lambdify((x1, x2), f, 'numpy')
for n in range(5):
x = GD(x,M,n)
# print(FF(*(x[0],x[1],M)))
print('-------', n, x)
if (M * alpha.subs([(x1, x[0]), (x2, x[1])]) <= epsilon):
print("aaaaaaaaaaaa", M*alpha.subs([(x1, x[0]), (x2, x[1])]))
break
else:
M = C * M
# print(M)
print("实际最优值",x,FF(*(x[0],x[1],M)))
print("理论最优值",FF(*(0,1,0)))
内点罚.py
'''
内点罚函数
'''
import sympy
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
plt.ion()
fig = plt.figure()
ax = Axes3D(fig)
def draw(ff, x, name):
X = np.arange(0, 4, 0.25)
Y = np.arange(0, 4, 0.25)
X, Y = np.meshgrid(X, Y)
Z = ff(*(X, Y))
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow', alpha=0.5)
ax.scatter(x[0], x[1], ff(*(x)), c='r', s=80)
ax.text(x[0], x[1], ff(*(x)), 'here:(%0.3f,%0.3f)' % (x[0], x[1]))
ax.set_zlabel('f') # 坐标轴
ax.set_ylabel('X2')
ax.set_xlabel('X1')
plt.pause(0.2)
plt.savefig(name)
plt.cla()
# 极值法最小化F
def OPT(F, x1, x2, r):
x = [[]] * 2
# 1 F(x1,x2,r)
# 2 偏导
F_x1 = sympy.diff(F, x1)
F_x2 = sympy.diff(F, x2)
# print(F_x1)
# print(F_x2)
# 3 求解
# try:
# sol=[sympy.solve(F_x1,x1,real=True),sympy.solve(F_x2,x2,real=True)]
# except NotImplementedError:
sol = sympy.nonlinsolve([F_x1, F_x2, sympy.Eq(F_x1, F_x2)], [x1, x2])
# print('sol_set:', sol)
# exit(0)
# sol=[sol[0][0],sol[1][0]]
# print(sol)
# 4 取极限判断
x_set = []
for s in sol:
x_temp = [[]] * 2
for i, x_ in enumerate(s):
x_ = sympy.limit(x_, r, 0)
x_temp[i] = x_
# print(x_temp)
x_set.append(x_temp)
return x_set
# s_set=[[x01,x02,x03],[]]
pass
if __name__ == '__main__':
####################### start
# 定义
x1, x2 = sympy.symbols('x1:3')
r = sympy.symbols('r')
f = 1 / 3 * (x1 + 1) ** 3 + x2
g1 = x1 - 1
g2 = x2
alpha = 1 / g1 + 1 / g2
F = f + r * alpha
# c = 0.1 # 缩小系数
# r = 10 # 障碍因子
# epsilon = 1e-5 # 终止限
# x = [3, 4]
ff = sympy.lambdify((x1, x2), f, 'numpy')
x_set = OPT(F, x1, x2, r)
# print(x_set)
# 判断得到符合条件的解
account=0
for x in x_set:
if (g1.subs(x1, x[0]) >= 0 and g2.subs(x2, x[1]) >= 0):
print(x)
draw(ff, x, './image1/6_'+str(account))
print("实际"+str(account), x, ff(*(x)))
account+=1
print("理论最优值", ff(*(1, 0)))
混合罚.py
'''
混合罚函数法
'''
import sympy
from OP6_2 import OPT, draw
if __name__ == '__main__':
# x = [0.00000, 1.00000]
####################### start
# 定义
x1, x2 = sympy.symbols('x1:3')
r = sympy.symbols('r')
f = -x1 + x2
h = x1 + x2 - 1
# g = sympy.log(x2)
g = x2 - 1
F = f + r * 1 / g + 1 / sympy.sqrt(r) * h ** 2
# epsilon = 1e-5 # 终止限
ff = sympy.lambdify((x1, x2), f, 'numpy')
x_set = OPT(F, x1, x2,r)
# 判断得到符合条件的解
account = 0
for x in x_set:
if (g.subs([(x1,x[0]),(x2,x[1])])>=0 and h.subs([(x1,x[0]),(x2,x[1])])==0):
print(x)
draw(ff, x, './image1/7_' + str(account))
print("实际" + str(account), x, ff(*(x)))
account += 1
print("理论最优值",ff(*(1,0)))
im2gif.py
'''
转成动图
'''
import os
import imageio
def create_gif(image_list, gif_name):
frames = []
for image_name in image_list:
frames.append(imageio.imread(image_name))
# Save them as frames into a gif
imageio.mimsave(gif_name, frames, 'GIF', duration=0.2)
return
def main():
os.chdir('./image/')
image_list = os.listdir('./')
gif_name = 'created_gif.gif'
create_gif(image_list, gif_name)
if __name__ == "__main__":
main()
网友评论