懲罰函式法(內點法、外點法)求解約束優化問題最優值 matlab
阿新 • • 發佈:2018-12-02
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
轉:https://blog.csdn.net/soaringlee_fighting/article/details/72885138