现代设计方法 张大可版 内点法与外点法的python实现

       注解已在程序注释给出，本方法只适合较为简单的惩罚函数，若设置的惩罚函数比较复杂，会出现卡住一直出不了结果的情况。本程序经示例验证无误，欢迎各位大佬批评指正。

# 开发人员：盛夏的原谅色
# 开发时间：2023/11/13 13:01
# 本文件仅适用于学习分享，不可用作商业用途侵权必究

# 内/外点惩罚法（解析法求解）
from sympy import symbols,diff,solve,Function,Eq,N
# 内/外点惩罚函数，可替换为自己需要的
def fun(x1,x2,r):
    # y=4*x1-x2**2-12-r*(x1**2+x2**2-25)**-1+r**-1*(x1**2+x2**2-10*x1-10*x2+34)**2  # 习题5.3混合惩罚函数（过于复杂，无法用解析法解出）
    y=x1**2+x2**2+r*(1-x1)**2  # 现代设计方法 张大可版 例5-2示例外点惩罚函数
    return y
# 计算迭代点
def dif(r0):
    x1 = symbols("x1")
    x2 = symbols("x2")
    r = symbols("r")
    # 分别对x1/x2求偏导，得到导函数
    derivative1 = diff(fun(x1, x2, r), x1)
    derivative2 = diff(fun(x1, x2, r), x2)
    # 令两式等0构建方程
    eq1=Eq(derivative1,0)
    eq2=Eq(derivative2,0)
    # 联立两方程求解出含r的解
    x_r = solve((eq1,eq2),(x1,x2))
    # 打印输出
    print(x_r)
    # 分别取出x1/x2含未知量r的解
    x1_val = x_r[x1]
    x2_val = x_r[x2]
    # 将未知量r的值带入计算得到结果
    result1 = x1_val.subs(r, r0)
    result2 = x2_val.subs(r, r0)
    # 返回x1,x2
    return result1,result2
def main():
    # 输入初始点
    x0 = input("请输入初始点(x1,x2),用 间隔：").split(" ")
    # 改变初始点的数据类型
    x = [float(i) for i in x0]
    r0 = float(input("请输入初始惩罚因子："))  # 一般选取1
    C = float(input("请输入惩罚因子递减/增系数：")) # 系数<1可用于内点法，系数>1可用于外点法
    acc0 = float(input("请输入精度要求1："))  # 输入精度要求
    # 初始点暂存
    X0 = [x[0], x[1]]
    k = 0
    # 默认迭代100次
    while k<100:
        # 计算迭代后点位
        X=dif(r0)
        # 计算精度值
        acc=abs((abs(X[0])**2+abs(X[1])**2)**0.5-(abs(X0[0])**2+abs(X0[1])**2)**0.5)
        # 打印输出
        print(acc)
        # 满足精度要求则打印输出结果，终止迭代
        if acc<=acc0:
            print("满足精度要求输出X({0},{1})".format(X[0],X[1]))
            break
        # 不满足精度要求则继续迭代，更新初始点，惩罚因子，打印输出当前结果
        else:
            X0[0],X0[1]=X[0],X[1]
            print(X0)
            r0=C*r0
            k+=1
            print("不满足精度要求继续，此时X({0},{1})".format(round(X[0],3),X[1]))
            continue
# 调用主函数
if __name__ == "__main__":
    main()