外点牛顿罚函数优化matlab代码

close all
clear all
clc
syms x1 x2 M;              %e为罚因子。
m（1）=1;                    %c为递增系数。赋初值。
c=8;
a（1）=20;
b（1）=20;                              
f=x1^2+x2^2+M*（（20-x1）^2+（2*x1-x2-30）^2）;
f0（1）=500;

%求偏导、海森元素
fx1=diff（f‘x1‘）;
fx2=diff（f‘x2‘）;
fx1x1=diff（fx1‘x1‘）;
fx1x2=diff（fx1‘x2‘）;
fx2x1=diff（fx2‘x1‘）;
fx2x2=diff（fx2‘x2‘）;

%外点法M迭代循环
for k=1:100                                                          
    x1=a（k）;x2=b（k）;M=m（k）;
    %牛顿法求最优值
    for n=1:100                                                      
        f1=subs（fx1）; %求解梯度值和海森矩阵
        f2=subs（fx2）;
        f11=subs（fx1x1）;
        f12=subs（fx1x2）;
        f21=subs（fx2x1）;
        f22=subs（fx2x2）;
        if（double（sqrt（f1^2+f2^2））<=1e-6）                                  %最优值收敛条件
            a（k+1）=double（x1）;b（k+1）=double（x2）;f0（k+1）=double（subs（f））;
            break;
        else
            X=[x1 x2]‘-inv（[f11 f12;f21 f22]）*[f1 f2