最優(yōu)化實(shí)驗(yàn)報(bào)告
——外點(diǎn)趾娃、內(nèi)點(diǎn)和混合罰函數(shù)法
實(shí)驗(yàn)?zāi)康?/h2>
之前我們已經(jīng)實(shí)驗(yàn)過無約束最優(yōu)化問題，這次我們將實(shí)驗(yàn)一下缔御，在有約束條件下抬闷，優(yōu)化算法應(yīng)該怎么做。
由于處理約束條件的辦法不同，約束優(yōu)化法可以分為直接間接兩類笤成。間接法的思想是將問題轉(zhuǎn)化成無約束的评架，然后無約束方法求解，雖然算法復(fù)雜炕泳，但可以求解高維問題纵诞，效率高、魯棒性好培遵。直接法構(gòu)造迭代過程浙芙，使迭代點(diǎn)在可行域D中，一步一步降低目標(biāo)函數(shù)值籽腕。算法簡單嗡呼，但是計(jì)算量大，一般只用于求解不等式約束問題皇耗。
罰函數(shù)法又稱乘子法南窗，是指將有約束最優(yōu)化問題轉(zhuǎn)化為求解無約束最優(yōu)化問題：其中M為足夠大的正數(shù)， 起"懲罰"作用郎楼， 稱之為罰因子万伤，F(xiàn)(x, M )稱為罰函數(shù)。內(nèi)部罰函數(shù)法也稱為障礙罰函數(shù)法呜袁。
這種方法是在可行域內(nèi)部進(jìn)行搜索壕翩，約束邊界起到類似圍墻的作用，如果當(dāng)前解遠(yuǎn)離約束邊界時(shí)傅寡，則罰函數(shù)值是非常小的放妈，否則罰函數(shù)值接近無窮大的方法。在進(jìn)化計(jì)算中荐操，研究者選擇外部罰函數(shù)法的原因主要是該方法不需要提供初始可行解芜抒。其中B(x)是優(yōu)化過程中新的目標(biāo)函數(shù)，Gi和Hj分別是約束條件gi(x)和hj(x)的函數(shù)托启，ri和cj是常數(shù)宅倒，稱為罰因子。
本次實(shí)驗(yàn)將驗(yàn)證間接法中的內(nèi)點(diǎn)屯耸、外點(diǎn)拐迁、混合罰函數(shù)法。
實(shí)驗(yàn)環(huán)境
本次實(shí)驗(yàn)疗绣，使用的程序語言為Python3.5.2线召，使用符號系統(tǒng)解析法推導(dǎo)計(jì)算，依賴的工具包為：
import math
import sympy
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
要想在IDE（集成開發(fā)環(huán)境）中實(shí)現(xiàn)畫圖的交互式效果多矮，還得寫
plt.ion()
實(shí)驗(yàn)內(nèi)容及步驟
- 定義函數(shù)和符號
- 求解最優(yōu)X
- 畫出構(gòu)造函數(shù)變化圖
1. 外點(diǎn)罰函數(shù)法基本原理與步驟
問題：
minf(X)=?x1+x2
s.t.{g1(X)=lnx2≥0h1(X)=x1+x2?1≥0
構(gòu)造一個(gè)增廣函數(shù)：
F(X,Mk)=f(X)+Mkα(X)
其中懲罰函數(shù)：
α(X)=∑j=1m[hj(X)]2+∑i=1l[gi(X)]2u(gi(X))
α(X){0,當(dāng)X∈D時(shí),0,當(dāng)X?D時(shí).
u(gi(X))={0,當(dāng)gi(X)≥0時(shí),1,當(dāng)gi(X)<0時(shí).i=1,2,?,l.
Mk
是一個(gè)逐漸增大的參數(shù)缓淹，稱為懲罰因子哈打。
當(dāng)X∈D
時(shí)
F(X,Mk)=f(X)
算法流程圖：
- 內(nèi)點(diǎn)罰函數(shù)法基本原理與步驟
問題：
x(1+1)3+x213min
s.t.{g1(X)=x1?1≥0,g2(X)=x2≥0.
構(gòu)造一個(gè)增廣函數(shù)：
F(X,rk)=f(X)+rk∑i=1l1gi(X)
其中rk
稱為障礙因子。
算法流程圖：
- 混合罰函數(shù)法基本原理與步驟
問題：
minf(X)=?x1+x2
s.t.{g1(X)=lnx2≥0h1(X)=x1+x2?1≥0
構(gòu)造一個(gè)增廣函數(shù)：
F(X,rk)=f(X)+rk∑i=1l1gi(X)+1rk∑j=1m[hj(X)]2
算法流程圖：
實(shí)驗(yàn)結(jié)果及分析
- 外點(diǎn)罰函數(shù)法：
起初讯壶，實(shí)驗(yàn)發(fā)生了錯(cuò)誤料仗，使用數(shù)值計(jì)算的方法做最速下降收斂增廣函數(shù)，在接近理論最優(yōu)解時(shí)伏蚊，梯度爆炸立轧，變得很大，經(jīng)過多次參數(shù)調(diào)整躏吊，也只能求得不精確（+-0.02）的解氛改。
然后使用了 符號解析法 ，配合 最優(yōu)步長法 颜阐，取得了極好的效果。
起始點(diǎn)：
迭代中逐漸增大非可行域的值（Z軸數(shù)值F明顯變化吓肋，最優(yōu)解也在變）：
經(jīng)過幾十次迭代凳怨，就取得最后結(jié)果：
- 內(nèi)點(diǎn)罰函數(shù)法
吸取調(diào)參經(jīng)驗(yàn)教訓(xùn)，結(jié)合間接求解有約束問題能轉(zhuǎn)化成無約束問題是鬼，并且函數(shù)顯然在定義域內(nèi)連續(xù)可導(dǎo)肤舞。所以，這次直接用了解非線性方程組均蜜，求極值取極限的方法李剖。
很快取得了結(jié)果：
- 混合罰函數(shù)法
復(fù)用上一題程序，只需改進(jìn)解的可行范圍檢查即可得：
程序
失敗的數(shù)值計(jì)算方法——會爆炸.py
'''
失敗的數(shù)值方法
'''
import sympy
import math
import numpy as np
C = 10 # 放大系數(shù)
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 # 數(shù)值求導(dǎo)
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))
外點(diǎn)罰.py
'''
這是有用的
外點(diǎn)罰函數(shù)
'''
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') # 坐標(biāo)軸
ax.set_ylabel('X2')
ax.set_xlabel('X1')
plt.pause(0.1)
# plt.show()
# plt.savefig('./image/%03d' % index)
plt.cla()
C = 10 # 放大系數(shù)
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 # 數(shù)值求導(dǎo)
# 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("實(shí)際最優(yōu)值",x,FF(*(x[0],x[1],M)))
print("理論最優(yōu)值",FF(*(0,1,0)))
內(nèi)點(diǎn)罰.py
'''
內(nèi)點(diǎn)罰函數(shù)
'''
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') # 坐標(biāo)軸
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 偏導(dǎo)
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 # 縮小系數(shù)
# 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("實(shí)際"+str(account), x, ff(*(x)))
account+=1
print("理論最優(yōu)值", ff(*(1, 0)))
混合罰.py
'''
混合罰函數(shù)法
'''
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("實(shí)際" + str(account), x, ff(*(x)))
account += 1
print("理論最優(yōu)值",ff(*(1,0)))
im2gif.py
'''
轉(zhuǎn)成動(dòng)圖
'''
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()
