基于scipy优化的填充函数法的非期望行为

Question

大家好，对于一个uni项目，我必须编写一个代码来实现一个算法，该算法使用填充函数来确定全局最小值

• 第0步：在可行域中生成一个随机点X0，并启动k=0
• 第1步：从这一点局部最小化，找到x*\k
• 第2步：利用x*定义一个填充函数，并对该函数进行局部最小化，直到找到一个点（x\u k+1）！=(例如：f（x+1）<；f（x\k）
• 第3步：设置k=k+1并返回步骤1

教授提出的填充函数如下：

其中x颚化符是可行空间中的随机点

以下是我对实施的看法：

# OBJECTIVE woods function (feasible reason -10<= xi <=10
def f(x):
    x1 = x[0]
    x2 = x[1]
    x3 = x[2]
    x4 = x[3]
    f1 = 100 * (x1 ** 2 - x2) ** 2 + (x1 - 1) ** 2 + (x3 - 1) ** 2 + \
         90 * (x3 ** 2 - x4) ** 2 + 10.1 * ((x2 - 1) ** 2 + (x4 - 1) ** 2) + 19.8 * (x2 - 1) * (x4 - 1)
    return f1  

# DEF FILLED
def filled(x_star,x_k):
    ro = 10 ** -6
    c=(10*(len(x_k)**(1/2))) + 1
    eps=0.05
    '''  tau=c*((math.exp(eps**2))/(math.exp(eps**2)-1))
    eta = np.linalg.norm(x_k - np.random.uniform(-10, 10, len(x_k))) ** 2
    phi = tau * (math.exp(min(0, f(x_k) - f(x_star)+ro)) ** 3)'''
    tau = 100
    eta = np.linalg.norm(x_k - np.random.uniform(-10, 10, len(x_k))) ** 2
    phi = min(0, (tau *(f(x_k) - f(x_star) + ro)) ** 3)
    ff = eta + phi
    return ff

def cb(x):
    if(x != x_star1).all():
        if f(x_star1) > f(x):
            return True  

def mainf(max_it):
    global it
    it=0
    nfunc = 0
    niter = 0

    while it < max_it:
        step = 0
        if it == 0:
            x0 = np.random.uniform(-10, 10, 4) # controllare sempre la regione di interesse data dal problema ed aggiustare x0 di conseguenza, se non migliora ad ogni step restringere la regione???
        min1 = optimize.minimize(f, x0, bounds=region)
        '''nfunc += min1['nfev']
        niter += min1['nit']'''
        x_star = min1['x']
        global x_star1
        x_star1 = x_star
        step = step + 1

        if step > 0:
            min2 = optimize.minimize(filled,x_star,args=x0,callback=cb, bounds=region)
            nfunc += min2['nfev']
            '''b =min2.hess_inv.todense() #if i want to print matrix b for bfgs'''
            niter += min2['nit']
            x0_prec = x_star1  # mi serve esclusivamente per stampare il valore della f.o precedente
            x0 = min2['x']
            step = 0

        it = it +1

        if it % 100 == 0:
            print('\nFILLED ITERATION: ', it,
                  '\n', min2,
                  '\nf(x*) = \t', f(x0),)
            ''' print('\n##################### UPDATING ##################### ',
                  '\nFILLED ITERATION: ', it,
                  '\nNUMBER OF ITERATION F = \t', nfunc,
                  '\nNUMBER OF ITERATION SUB_ROUTINE = \t', niter,
                  '\nNEW X_0 (X*) = \t', x0,
                  '\nF(x*) = \t', f(x0),
                  '\nF(x*-1) = \t', f(x0_prec),
                  '\n#################################################### ')'''  

# LAUNCH
region = ((-10,10),(-10,10),(-10,10),(-10,10))
mainf(1000)

结果如下（仅每100次迭代打印一次）：

FILLED ITERATION:  100 
       fun: 577.5439225647822
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([ 3.80521846e+10, -1.54833810e+10, -4.01513113e+10, -1.19568671e+09])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 105
      nit: 3
   status: 0
  success: True
        x: array([1.00000086, 1.00000107, 1.00000083, 1.00000159]) 
f(x*) =      1.1568608740354235e-10

FILLED ITERATION:  200 
       fun: -1.2231896380796914e+25
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([7.00795898e+24, 3.57388094e+23, 6.30719727e+24, 3.22349966e+23])
  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 10
      nit: 1
   status: 0
  success: True
        x: array([-10., -10., -10., -10.]) 
f(x*) =      2304082.0

FILLED ITERATION:  300 
       fun: 457.02318754757766
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([-1.10431945e+10,  1.21949498e+10,  5.53201735e+09,  2.71537503e+10])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 145
      nit: 4
   status: 0
  success: True
        x: array([0.99999833, 0.9999989 , 0.99999895, 0.99999918]) 
f(x*) =      6.860892674573475e-10

FILLED ITERATION:  400 
       fun: 486.9572281371511
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([ 2.45998492e+10, -1.13562475e+10, -5.95760544e+09,  1.24355462e+10])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 165
      nit: 3
   status: 0
  success: True
        x: array([1.00000164, 0.99999653, 0.99999717, 1.00000206]) 
f(x*) =      9.943858291905142e-09

FILLED ITERATION:  500 
       fun: 665.3558594771831
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([-2.23804012e+10, -1.75026323e+10, -5.86394467e+10,  3.99424992e+09])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 160
      nit: 3
   status: 0
  success: True
        x: array([0.99999841, 0.99999842, 0.99999884, 1.00000119]) 
f(x*) =      1.3785899179572299e-09

FILLED ITERATION:  600 
       fun: -1.2230622314507554e+25
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([-7.00740944e+24,  3.57363398e+23, -6.30669626e+24,  3.22327632e+23])
  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 10
      nit: 1
   status: 0
  success: True
        x: array([ 10., -10.,  10., -10.]) 
f(x*) =      2304002.0

FILLED ITERATION:  700 
       fun: -1.223062231450772e+25
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([-7.00740901e+24,  3.57363828e+23, -6.30669583e+24,  3.22327632e+23])
  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 15
      nit: 2
   status: 0
  success: True
        x: array([ 10., -10.,  10., -10.]) 
f(x*) =      2304002.0

FILLED ITERATION:  800 
       fun: 591.9197854392103
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([-1.75122200e+10, -1.22348657e+10,  4.27795139e+10, -2.66557905e+10])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 130
      nit: 3
   status: 0
  success: True
        x: array([1.00000198, 1.00000166, 1.00000205, 1.00000216]) 
f(x*) =      1.0189519266743117e-09

FILLED ITERATION:  900 
       fun: 32.3560340558523
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([3.71659764e+10, 3.80271187e+10, 7.81573243e+10, 1.06411141e+11])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 185
      nit: 4
   status: 0
  success: True
        x: array([0.99999737, 0.99999752, 1.00000158, 1.00000202]) 
f(x*) =      9.043566134285354e-10

FILLED ITERATION:  1000 
       fun: -511999980101.2271
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([8.44794391e+15, 4.30829929e+14, 7.60321360e+15, 3.88618035e+14])
  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 10
      nit: 1
   status: 0
  success: True
        x: array([-10., -10., -10., -10.]) 
f(x*) =      2304082.0

解应该是大约x=[1]和f（x*）=0
这里主要有两个问题，我正在努力解决

1. 在一些迭代中，例如迭代100，结果是非常令人满意的，但是向前看f（x）应该总是改进并收敛到全局最小值（0），但是很明显迭代300比迭代100稍微差一些
2. 更大的问题，有时极小化以消息b'收敛结束：投影的范数\u梯度\lt=_结果X*向量总是[10]，函数值完全不存在

说了这么多，你们能帮我理解问题出在哪里吗？非常感谢
*****编辑******
我想我设置填充函数的方式可能有问题，我试着按照图片中的公式，但我不确定我做得对不对