Python实现Powell方法

题目：

代码：（思路以后再写）

import numpy as np
from sympy import symbols, diff,solve

count = 0  # 迭代次数
maxtimes = 40  # 最大迭代次数
e = 0.01

X0 = [0,0]  # 初始点
# 初始方向
d1 = [1,0]
d2 = [0,1]

alpha = symbols('alpha')

def fun(x1,x2):
    return 10*(x1+x2-5)**2 + (x1 - x2)**2

while count < maxtimes:
    count = count + 1
    print('-------第{}次迭代-----'.format(count))
    # (1)
    F0 = fun(X0[0],X0[1])
    print('F0 = ',F0)
    
    X1 = list(np.array(X0) + alpha*np.array(d1))
    res1 = fun(X1[0],X1[1])
    print(res1)
    a = sympy.solve(res1,alpha)   # 写入需要解的方程体
    a1 = complex(a[0]).real  # 取出复数的实部
    a1 = round(a1,4)
    print('a1 = ',a1)
    
    # X1 = [a1,X1[1]]
    X1 = list(np.array(X0) + a1*np.array(d1))
    print('X1 = ',X1)
    
    F1 = fun(X1[0],X1[1])
    print('f1=',F1)
    
    delta1 = abs(F0-F1)
    print('delta1 = ',delta1)
    
    # (2)
    X2 = list(np.array(X1) + alpha*np.array(d2))
#     print('X2=',X2)
    
    res2 = fun(X2[0],X2[1])
    print(res2)
    b = sympy.solve(res2,alpha)   # 写入需要解的方程体
    a2 = complex(b[0]).real  # 取出复数的实部
    a2 = round(a2,4)
    print('a2 = ',a2)
    
    # X2 = [X2[0],a2]
    X2 = list(np.array(X1) + a2*np.array(d2))
    print('X2 = ',X2)
    
    F2 = fun(X2[0],X2[1])
    print('f2=',F2)
    
    delta2 = abs(F2-F1)
    print('delta2 = ',delta2)
    # 取沿 d1 d2 搜索方向增量的最大者
    d_delta = max(delta1,delta2)
    print('d_delta = ',d_delta)
    
    # 终点 x2 的反射点和函数值为
    X3 = list(2*np.array(X2)-np.array(X0))
    print('X3 = ',X3)
    F3 = fun(X3[0],X3[1])
    print('f3 = ',F3)
    
    # (3) 鲍威尔条件判断
    if F3 > F0:
        print('---不满足判断条件，继续迭代...')
        if F2 < F3:
            X0 = X2
    else:
        d3 = list(np.array(X2)-np.array(X0))
        print('d3 = ',d3)
        X0 = list(np.array(X2) + alpha*np.array(d3))
        
        res3 = fun(X0[0],X0[1])
        print(res3)
        c = sympy.solve(res3,alpha)   # 写入需要解的方程体
        a3 = complex(c[0]).real  # 取出复数的实部
        a3 = round(a3,4)
        print('a3 = ',a3)
        
        X0 = list(np.array(X2) + a3*np.array(d3))
        
    if F3 < e:  # 达到了精度，结束迭代
        print('满足Powell条件，结束迭代')
        print('共迭代了{}次'.format(count))
        print('此时的坐标点：',X0)
        print('此时的极小值：',round(fun(X0[0],X0[1]),7))
        break

结果（太长，只截取了部分）：