内点法外点法matlab代码,懲罰函數法（內點法、外點法）求解約束優化問題最優值 matlab...

1、 用外點法求下列問題的最優解

方法一：外點牛頓法：

clc

m=zeros(1,50);a=zeros(1,50);b=zeros(1,50);f0=zeros(1,50);%a b為最優點坐標，f0為最優點函數值，f1 f2最優點梯度。

syms x1 x2 e;                                            %e為罰因子。

m(1)=1;c=10;a(1)=0;b(1)=0;                               %c為遞增系數。賦初值。

f=x1^2+x2^2+e*(1-x1)^2;f0(1)=1;

fx1=diff(f,'x1');fx2=diff(f,'x2');fx1x1=diff(fx1,'x1');fx1x2=diff(fx1,'x2');fx2x1=diff(fx2,'x1');fx2x2=diff(fx2,'x2');%求偏導、海森元素。

for k=1:100                                                          %外點法e迭代循環.

x1=a(k);x2=b(k);e=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))<=0.001)                                  %最優值收斂條件

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]';

x1=X(1,1);x2=X(2,1);

end

end

if(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))<=0.001)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=0.001)   %罰因子迭代收斂條件

a(k+1)   %輸出最優點坐標，罰因子迭代次數，最優值

b(k+1)

k

f0(k+1)

break;

else

m(k+1)=c*m(k);

end

end

方法二：外點梯度法：

clc

m=zeros(1,50);a=zeros(1,50);b=zeros(1,50);f0=zeros(1,50);

syms d x1 x2 e;

m(1)=1;c=10;a(1)=0;b(1)=0;

f=x1^2+x2^2+e*(1-x1)^2; f0(1)=1;

fx1=diff(f,'x1');

fx2=diff(f,'x2');

for k=1:100

x1=a(k);x2=b(k);e=m(k);

for n=1:100

f1=subs(fx1);

f2=subs(fx2);

if(double(sqrt(f1^2+f2^2))<=0.002)

a(k+1)=double(x1);b(k+1)=double(x2);f0(k+1)=double(subs(f));

break;

else

D=(x1-d*f1)^2+(x2-d*f2)^2+e*(1-(x1-d*f1))^2;

Dd=diff(D,'d'); dd=solve(Dd); x1=x1-dd*f1; x2=x2-dd*f2;

end

end

if(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))<=0.001)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=0.001)

a(k+1)

b(k+1)

k

f0(k+1)

break;

else

m(k+1)=c*m(k);

end

end

2、用內點法求下列問題的最優解

內點牛頓法

clc

m=zeros(1,50);a=zeros(1,50);b=zeros(1,50);f0=zeros(1,50);

syms x1 x2 e;

m(1)=1;c=0.2;a(1)=2;b(1)=-3;

f=x1^2+x2^2-e*(1/(2*x1+x2-2)+1/(1-x1)); f0(1)=15;

fx1=diff(f,'x1');fx2=diff(f,'x2');fx1x1=diff(fx1,'x1');fx1x2=diff(fx1,'x2');fx2x1=diff(fx2,'x1');fx2x2=diff(fx2,'x2');

for k=1:100

x1=a(k);x2=b(k);e=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))<=0.002)

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]';

x1=X(1,1);x2=X(2,1);

end

end

if(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))<=0.001)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=0.001)

a(k+1)

b(k+1)

k

f0(k+1)

break;

else

m(k+1)=c*m(k);

end

end