现代设计方法 张大可版 梯度法 python实现

采用梯度法求解最优化问题，示例基于书籍范例编写，并不适合所有情况，可以按需求进行修改补充。

代码注解已在程序注释标出，新手编写的程序较为简单，欢迎大佬批评指导

# 开发人员：盛夏的原谅色
# 开发时间：2023/10/27 9:09
# 本文件仅适用于学习分享，不可用作商业用途侵权必究
# 梯度法python实现

from sympy import symbols,diff,solve,Function,Eq,N
# 设定对应函数
def equation(x1,x2):
    # 现代设计方法 张大可 第四章4.2函数，初始点（0,0），精度要求0.1
    fx=(x1-1)**2+(x2-2)**2
    # 现代设计方法 张大可 第四章例4-3函数，初始点（2,2），精度要求0.005
    # fx=x1**2+25*x2**2
    return fx
# 梯度计算函数
def gradient(x_1,x_2):
    m = symbols('x_1')
    n = symbols('x_2')
    # 分别计算关于x1和x2的偏导数方程
    derivative1 = diff(equation(m, n), m)
    derivative2 = diff(equation(m, n), n)
    # 将初始点带入偏导方程，得到梯度
    x_g1 = N(derivative1.subs(m, x_1))
    x_g2 = N(derivative2.subs(n, x_2))
    # 以元组形式返回梯度数据
    return x_g1,x_g2
# 最优步长计算函数
def step(x1,x2):
    a = symbols('a')
    # 调用梯度函数计算梯度
    x = gradient(x1,x2)
    # 将梯度函数返回的梯度数据取出
    x_g1,x_g2=x[0],x[1]
    # 带入a值得到带步长参数的迭代点
    x_1,x_2 = x1-x_g1*a,x2-x_g2*a
    # 将带a的迭代点带入方程求导，得到关于a的导数方程
    derivative = diff(equation(x_1,x_2), a)
    # 避免出现a=0导致的程序报错
    try:
        # 令关于a的方程为0（存在极值的必要条件）,求解得到最佳步长a的数值
        a_value = solve(derivative, a)[0]
    except:
        # 如果出现a为0的情况，将a赋值为0，避免求解函数遍历不到a而报错
        a_value = 0
    # 将步长带入得到迭代点
    m = x1-a_value*x_g1
    n = x2-a_value*x_g2
    # 以元组形式返回迭代点，步长，梯度数据
    return m,n,a,x_g1,x_g2
# 主函数输入初始点，精度要求，循环计算
def main():
    # 输入初始点，精度要求
    x_1 = round(float(input('请输入初始点x1:')),7)
    x_2 = round(float(input('请输入初始点x2:')),7)
    acc0 = round(float(input('请输入精度要求:')),7)
    i = 0
    # while循环，i<20确保不会陷入死循环，可以根据情况修改
    while i<20:
        # 将初始值暂存，方便计算精度时调用
        x0 = [x_1,x_2]
        # 调用step函数获得迭代点，步长，梯度数据
        x = step(x_1,x_2)
        # 从返回的元组中取出相应的数据
        x_1 = x[0]
        x_2 = x[1]
        a = x[2]
        x_g1=x[3]
        x_g2=x[4]
        # 进行精度的计算
        acc = (x_g1**2+x_g2**2)**0.5
        i+=1
        # 比较精度是否满足精度要求，若满足则输出步长，最优点，最优解等数据
        if acc<acc0:
            # 将最优点代入函数得到最优解
            value=equation(x0[0],x0[1])
            print('满足精度要求，此时精度为:',acc,'步长为:',a)
            print('最优点为：',x0,'对应最优解为:',value)
            break
        else:
            # 不满足精度继续进行迭代计算
            print('此时精度为:',acc,'不满足要求继续')
            print('初始点为：', x0, '迭代点为:', x_1, x_2)
            continue

# 执行主函数
if __name__ == "__main__":
    main()