用Powell法优化设计程序与一维搜索黄金分割法组合.doc

用Powell法优化设计程序与一维搜索黄金分割法组合
编程求解函数

的极小点x，初始点x0=[1,1]T，迭代精度ε=0.001。
解：已知，初始点x0=[1,1]T，迭代精度ε=0.001。
在该优化设计过程中，黄金分割搜索法作为POWELL算法主程序中的一部分。在POWELL算法运行过程中会多次调用黄金分割搜索算法程序。这样可以缩短优化设计计算时间。
1. MATLAB源程序代码
1.1 关于α的目标函数源代码
function m = y(x1,x2,d1,d2,alpha)%建立关于α的一元函数y(α)
m = (x1+alpha*d1)^2+2*(x2+alpha*d2)^2-4*(x1+alpha*d1)-2*(x1+alpha*d1)*(x2+alpha*d2);
1.2 一元函数最小值区间函数源代码
function [a,b] = section(x1,x2,d1,d2)%采用外推法求解一元函数的最小值区间
x11 = x1;x22 = x2;d11 = d1;d22 = d2;h0 = 1;h = h0;
alpha1 = 0;
y1 = y(x11,x22,d11,d22,alpha1);
alpha2 = h;
y2 = y (x11,x22,d11,d22,alpha2);
t = 0;
if y2>y1
h = -h;alpha3 = alpha1;y3 = y1;t = 1;
end
while(1)
if t == 1
alpha1 = alpha2;y1 = y2;
alpha2 = alpha3;y2 = y3;
else t = 1;
end
alpha3 = alpha2+h;y3 = y(x11,x22,d11,d22,alpha3);
if y3<y2
h = 2*h;
else
break;
end
end
if alpha1>alpha3
tem = alpha1;alpha1 = alpha3;alpha3 = tem;
a = alpha1;b = alpha3;
else a=alpha1;b = alpha3;
end
1.3黄金分割搜索法函数源代码
function alpha = ALPHA(x1,x2,d1,d2,A,B)%利用黄金分割法求解关于α的函数y(α)的极小点α*
x11 = x1;x22 = x2;d11 = d1;d22 = d2;
a = A;b = B;
ep = 0.001;r = 0.618;
alpha1 = b-r*(b-a);
y1 = y(x11,x22,d11,d22,alpha1);
alpha2 = a+r*(b-a);
y2 = y(x11,x22,d11,d22,alpha2);
while(1)
if y1>=y2
a = alpha1;alpha1 = alpha2;
y1 = y2;
alpha2 = a+r*(b-a);
y2 = y(x11,x22,d11,d22,alpha2);
else
b = alpha2;alpha2 = alpha1;
y2 = y1;
alpha1 = b-r*(b-a);
y1 = y(x11,x22,d11,d22,alpha1);
end
if abs(b-a)<ep && abs(y2-y1)<ep
break;
end
end
alpha = 0.5*(a+b);
1.4 POWELL算法程序源代码
function [z,fmin] = powell(f) %根据无约束优化方法POWELL算法，在二维空间中求解f(x1,x2)的最小值点，求解结果返回变量坐标(x1,x2)和极小值fmin
k = 0;n = 2;
x = [0;0;];
ff(1) = f(x(1),x(2));
ep = 0.001;
d = [1;0;0;1];
while(1)
x00 = [x(1);x(2);];
for i = 1:n
[a(i),b(i)]=section(x(2*i-1),x(2*i),d(2*i-1),d(2*i));
alpha(i)=ALPHA(x(2*i-1),x(2*i),d(2*i-1),d(2*i),a(i),b(i));
x(2*i+1)=x(2*i-1)+alpha(i)*d(2*i-1);
x(2*i+2)=x(2*i)+alpha(i)*d(2*i);
ff(i+1)=f(x(2*i+1),x(2*i+2));
end
for i=1:n
Delta(i)=ff(i)-ff(i+1);
end
delta=max(Delta);
for i = 1:n
if delta == Delta(i)
m = i;
break;
end
end
d(2*n+1) = x(2*n+1)-x(1);
d(2*n+2)=x(2*n+2)-x(2);
x(2*n+3)=2*x(2*n+1)-x(1);
x(2*n+4)=2*x(2*n+2)-x(2);
ff(n+2)=f(x(2*n+3),x(2*n+4));
f0=ff(1);f2=ff(n+1);f3=ff(n+2);
k=k+1;
R(k,:)=[k,x',d',ff];
if f3<f0 && (f0-2*f2+f3)*(f0-f2-delta)^2<0.5*delta*(f0-f3)^2
[a(n+1),b(n+1)]=section(x(2*n+1),x(2*n+2),d(2*n+1),d(2*n+2));
alpha(n+1)=ALPHA(x(2*n+1),x(2*n+2),d(2*n+1),d(2*n+2),a(n+1),b(n+1));
x(1)=x(2*n+1)+alpha(n+1)*d(2*n+1);
x(2)=x(2*n+2)+alpha(n+1)*d(2*n+2);
for i=m:n
d(2*i-1)=d(2*i+1);
d(2*i)=d(2*i+2);
end
else
if f2<f3
x(1)=x(2*n+1);
x(2)=x(2*n+2);
else
x(1)=x(2*n+3);
x(2)=x(2*n+4);
end
end
RR(k,:)=alpha;
ff(1)=f(x(1),x(2));
if (((x(2*n+1)-x00(1))^2+(x(2*n+2)-x00(2))^2)^(1/2))<ep
break;
end
end
z = [x(1);x(2)];
fmin = f(x(1),x(2));
2. 运行程序计算
>> f = @(x1,x2) x1^2+2*x2^2-4*x1-2*x1*x2
f =
@(x1,x2)x1^2+2*x2^2-4*x1-2*x1*x2
>> [z,fmin] = powell(f)
z =
3.9998
1.9998
fmin =
-8.0000
>>
经过计算可知，极小值点为(3.9998,1.9998)，极小值为-8。
3. 验算
通过MATLAB软件内置fminsearch函数进行验算。
3.1 绘图
>> x = linspace(-10,10,50);y = linspace(-10,10,50);
[x,y] = meshgrid(x,y);
z = x.^2+2*y.^2-4*x-2*x.*y;
subplot(1,2,1);%在一幅框中显示两张图，显示第一张图
cs = contour(x,y,z);clabel(cs); %绘制等高线图
xlabel('x1');ylabel('x2');
title('(a)Contour plot');grid;
>> subplot(1,2,2); %在一幅框中显示两张图，显示第二张图
cs = surfc(x,y,z); %绘制空间曲面图
zmin = floor(min(z));
zmax = ceil(max(z));
xlabel('x1');ylabel('x2');zlabel('f(x1,x2)');
title('(b) Mesh plot');
>>

通过观察可以发现极小值点在(4,2)附近。
3.2 fminsearch函数验算
>> f = @(x) x(1)^2+2*x(2)^2-4*x(1)-2*x(1)*x(2)
f =
@(x)x(1)^2+2*x(2)^2-4*x(1)-2*x(1)*x(2)
>> [x,fval] = fminsearch(f,[1,1])
x =
4.0000 2.0000
fval =
-8.0000
>>
通过fminsearch函数验算，确认极小值点为(4,2)，极小值为-8。
4. 结论
函数的极小值点为(4,2)，极小值为-8。
5. 程序流程图
黄金分割法程序框图
开 始
给定a、b、c

是
否

是

结 束
否
POWELL算法程序框图
开 始

是
否

判别条件是否满足
是

否

是
否

是

结 束
否