外点罚函数法
还没搞定…
function [A,B]=Exterior(f,p,eps)
format long;
d=1;
e1=1e-4;
k=1;
A(k)=0;B(k)=0;
r(k)=1;a=2;
while d>eps
x1=A(k);x2=B(k);
if vpa(p(x1,x2))>0
u1=1;
else
u1=0;
end
F=f+r(k)*u1*p^2;
Fx1=diff(F,'x1');
Fx2=diff(F,'x2');
Fx11=diff(Fx1,'x1');
Fx12=diff(Fx1,'x2');
Fx21=diff(Fx2,'x1');
Fx22=diff(Fx2,'x2');
for n=1:100
F1=subs(Fx1); %求解梯度值和海森矩阵
F2=subs(Fx2);
F11=subs(Fx11);
F12=subs(Fx12);
F21=subs(Fx21);
F22=subs(Fx22);
if(double(sqrt(F1^2+F2^2))<=e1) %梯度法最优值收敛条件
A(k+1)=double(x1);B(k+1)=double(x2);
break;
else
X=[x1 x2]'-([F11 F12;F21 F22])\[F1 F2]';
x1=X(1,1);x2=X(2,1);
end
end
d=norm((A(k+1)-A(k)),(B(k+1)-B(k)));
r(k+1)=a*r(k);
k=k+1;
end
A(k)
B(k)
double(subs(f))
clear all;
clc
syms x1 x2 x;
gama=1.1;
f=x1^2+(1/3)*x2^2;
p=x1+x2-1;
eps=1e-4;
[A,B]=Exterior(f,p,eps)
function f=myfun(x)
f=x(1)^2+(1/3)*x(2)^2;
function p=myfun2(Aeq,beq,x)
p=Aeq*[x(1);x(2)]-beq;