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