外点、内点和混合罚函数法（最优化3）

作者: 小火伴 | 来源:发表于2018-01-17 10:53 被阅读1706次

外点、内点和混合罚函数法（最优化3）
SEO如何做好网站布局？(快速排名事半功倍)
4个注意点教你避免网站过度优化
NeoWorld 2.0 大版本上线 3 个月，你必须知道的几件
如何做已有功能的优化？
2018-12-18服务器性能优化基础
内归因和外归因的最佳点
4. 功能点的优化
Mysql多表查询
2，肩袖部常见6个肌肉起止点：

最优化实验报告

——外点、内点和混合罚函数法

实验目的

之前我们已经实验过无约束最优化问题，这次我们将实验一下，在有约束条件下，优化算法应该怎么做。

由于处理约束条件的办法不同，约束优化法可以分为直接间接两类。间接法的思想是将问题转化成无约束的，然后无约束方法求解，虽然算法复杂，但可以求解高维问题，效率高、鲁棒性好。直接法构造迭代过程，使迭代点在可行域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()

外点、内点和混合罚函数法（最优化3）
最优化实验报告 ——外点、内点和混合罚函数法实验目的之前我们已经实验过无约束最优化问题，这次我们将实验一下，在...
SEO如何做好网站布局？(快速排名事半功倍)
子爵SEO认为，现在做SEO最重要的有四点：网站布局、外链建设、内链优化、文章撰写。关于外链建设和内链优化，子爵S...
4个注意点教你避免网站过度优化
要做SEO优化，首先一点是需要优化我们自身的网站，我们都知道，优化网站一般就是优化关键词、外内链等。所以事情都有一...
NeoWorld 2.0 大版本上线 3 个月，你必须知道的几件
转眼间，NeoWorld 2.0 大版本上线 3 个月了。 3 个月以来，游戏内持续优化和完善，游戏外紧跟政策和趋...
如何做已有功能的优化？
功能点的优化是最基础的工作 1、对功能点的不断优化就是迭代 2、不要期望用新加功能来解决老的功能问题 3、功能优化...
2018-12-18服务器性能优化基础
第1个优化点: DNS第2个优化点:TCP --> I/O第3个优化点:请求,响应第4个优化点:Processin...
内归因和外归因的最佳点
《爱的五种能力》节选，讲了一个例子，主要是丈夫和婆妹过于亲密，冷落了妻子，妻子不应该外归因，应该内归因，要有...
4. 功能点的优化
功能点的优化是最基础的工作对功能点的不断优化就是迭代不要期望用新功能来解决老功能的问题功能点的优化和新功能设...
Mysql多表查询
Mysql多表查询，知识点包括：笛卡尔积查询[1] 内连接[^2] 外连接（左外连接、右外连接）[^3] uni...
2，肩袖部常见6个肌肉起止点：
1三角肌前束起点：锁骨外1/3前面近固点：肩屈、肩内旋、肩水平内收、辅助肩外展中束起点：肩峰外侧缘近固点...