%powell·求方程
f=10*(x1+x2-5)^4+(x1-x2+x3)^2
+(x2+x3)^6
的最优解
function MyPowell()
syms x1 x2 x3 a;
f=10*(x1+x2-5)^4+(x1-x2+x3)^2
+(x2+x3)^6;
error=10^(-3);
D=eye(3);
x0=[0 0 0]';
for k=1:1:10^6
MaxLength=0;x00=x0;m=0;
if k==1,s=D;end
for i=1:3
x=x0+a*s(:,i);
ff=subs(f,{x1,x2,x3},{x(1),x(2),x(3)});
t=Divide(ff,a);
%调用了进退法分割区间
aa=OneDemensionslSearch(ff,a,t);
%调用了0.618法进行一维搜索
xx=x0+aa*s(:,i);
fx0=subs(f,{x1,x2,x3},{x0(1),x0(2),x0(3)});
fxx=subs(f,{x1,x2,x3},{xx(1),xx(2),xx(3)});
length=fx0-fxx;
if length>MaxLength,MaxLength=length;m=m+1;end
x0=xx;
end
ss=x0-x00;
ReflectX=2*x0-x00;
f1=subs(f,{x1,x2,x3},{x00(1),x00(2),x00(3)});
f2=subs(f,{x1,x2,x3},{x0(1),x0(2),x0(3)});
f3=subs(f,{x1,x2,x3},{ReflectX(1),ReflectX(2),ReflectX(3)});
if f3
x=x0+a*ss;
ff=subs(f,{x1,x2,x3},{x(1),x(2),x(3)});
t=Divide(ff,a);
aa=OneDemensionslSearch(ff,a,t);
x0=x0+aa*ss;
for j=m:(3-1),s(:,j)=s(:,j+1);end
s(:,3)=ss;
else
if f2>f3,
x0=ReflectX;end
end
if norm(x00-x0)
k
x0
end
opx=x0;
val=subs(f,{x1,x2,x3},{opx(1),opx(2),opx(3)});
disp('最优点:');opx'
disp('最优化值:');val
disp('迭代次数:');k
附录:
0.618法
function output=OneDemensionslSearch(f,x,s,r)
if nargin<4,r=1e-6;end
a=s(1);b=s(2);
a1=a+0.382*(b-a);fa1=subs(f,{x},{a1});
a2=a+0.618*(b-a);fa2=subs(f,{x},{a2});
while abs((b-a)/b)>r
&&
abs((fa2-fa1)/fa2)>r
if fa1
b=a2;a2=a1;fa2=fa1;a1=a+0.382*(b-a);fa1=subs(f,{x},{a1});
else
a=a1;a1=a2;fa1=fa2;a2=a+0.618*(b-a);fa2=subs(f,{x},{a2});
end
end
op=(a+b)/2;
%fop=subs(f,{x},{op});
output=op;
进退法
%对任意一个一维函数函数进行区间分割,使其出现“高—低—高”的型式
function output=Divide(f,x,m,n)
if nargin<4,n=1e-6;end
if nargin<3,m=0;end
step=n;
t0=m;ft0=subs(f,{x},{t0});
t1=t0+step;ft1=subs(f,{x},{t1});
if ft0>=ft1
t2=t1+step;ft2=subs(f,{x},{t2});
while ft1>ft2
t0=t1;
%ft0=ft1;
t1=t2;ft1=ft2;
step=2*step;t2=t1+step;ft2=subs(f,{x},{t2});
end
else
step=-step;
t=t0;t0=t1;t1=t;ft=ft0;
%ft0=ft1;
ft1=ft;
t2=t1+step;ft2=subs(f,{x},{t2});
while ft1>ft2
t0=t1;
%ft0=ft1;
t1=t2;ft1=ft2;
step=2*step;t2=t1+step;ft2=subs(f,{x},{t2});
end
end
output=[t0,t2];