用matlab编程求解powell法问题。

1/目标函数程序清单

function m=f(x1,x2)

m=x1^2+2*x2^2-4*x1-2*x1*x2

2/关于α的目标函数

function m=y(x1,x2,d1,d2,alpha)

m=(x1+alpha*d1)^2+(x2+alpha*d2)^2-(x1+alpha*d1)*(x2+alpha*d2)-10*(x1+alpha*d1)-4*(x2+

alpha*d2)+60

3\

function [a,b]=section(x1,x2,d1,d2)

x11=x1%给出起始点坐标x1，x2和搜索方向d1，d2

x22=x2

d11=d1

d22=d2

h0=1 %初始化

h=h0alpha1=0

y1=y(x11,x22,d11,d22,alpha1)%代入α1求解y1

alpha2=hy2=y(x11,x22,d11,d22,alpha2)

t=0

if y2>y1

h=-halpha3=alpha1y3=y1t=1%如果y2>y1，则改变搜索方向

end

while(1)

if t==1

alpha1=alpha2y1=y2%实现交换

alpha2=alpha3y2=y3

else t=1

end

alpha3=alpha2+hy3=y(x11,x22,d11,d22,alpha3)

if y3<y2

h=2*h%改变搜索步长，将其加倍

else

break

end

end

if alpha1>alpha3

tem=alpha1alpha1=alpha3alpha3=tem%比较大小，保证输出区间为【a，b】

a=alpha1b=alpha3

else a=alpha1b=alpha3

end

4\

function alpha=ALPHA(x1,x2,d1,d2,A,B)

x11=x1x22=x2 %给出起始点坐标x1，x2和搜索方向d1，d2

d11=d1d22=d2

a=Ab=B %获取区间

epsilon=0.000001%初始化，给定进度

r=0.618

alpha1=b-r*(b-a)

y1=y(x11,x22,d11,d22,alpha1)%代入α1求解y1

alpha2=a+r*(b-a)

y2=y(x11,x22,d11,d22,alpha2)

while(1)

if y1>=y2 %根据区间消去法原理缩短搜索空间

a=alpha1alpha1=alpha2

y1=y2

alpha2=a+r*(b-a)

y2=y(x11,x22,d11,d22,alpha2)

else

b=alpha2alpha2=alpha1

y2=y1

alpha1=b-r*(b-a)

y1=y(x11,x22,d11,d22,alpha1)

end

if abs(b-a)<epsilon &

abs(y2-y1)<epsilon %判断是否满足进度要求

break%满足进度要求则退出循环迭代过程

end

end

alpha=0.5*(a+b)%返回值

5\主程序

clear all

clc

k=0n=2

x=[00]ff(1)=f(x(1),x(2)) %初始化

epsilon=0.00001

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))%调用section（）函数求解一元函数有y（α）最小值时的区间

alpha(i)=ALPHA(x(2*i-1),x(2*i),d(2*i-1),d(2*i),a(i),b(i))%调用ALPHA（）函数求解y（α）最小值点α*

x(2*i+1)=x(2*i-1)+alpha(i)*d(2*i-1)%沿di方向搜索

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) %求出反射点搜索方向dn+1

d(2*n+2)=x(2*n+2)-x(2)

x(2*n+3)=2*x(2*n+1)-x(1) %搜索到反射点Xn+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 %根据函数值之差最大值的下标值m，对原方向组进行替换

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)) %计算下一轮迭代过程需要的f0值

if

(((x(2*n+1)-x00(1))^2+(x(2*n+2)-x00(2))^2)^(1/2))<epsilon %判断是否满足精度要求

break %满足进度要求则退出循环迭代过程

end

end

xx=[x(1)x(2)]

fmin=f(x(1),x(2)) %显示最小值及其所对应的坐标

最数下降法解无约束优化的程序

先建立一维搜索的m文件：minWP.m如下

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,minf]=minWP(f,XMAX,c1,c2,alpha,tol)

%一维搜索的Wolfe-Powell法

%作者:龚纯 王正林<<精通 Matlab 最优化计算>>

%f:目标函数

%XMAX:搜索最大值

%c1:可接受系数1

%c2:可接受系数2

%alpha:增大步长倍数

%tol:精度

%x:极小值点

%minf:极小值点处的函数值

format long

if nargin==5

tol=1.0e-6

end

if ~(c1>0)||~(c1<c2)||~(c2<1)||~(XMAX>0)||~(alpha>1)

error('参数不对')

end

var=findsym(f)

df=diff(f)

f0=subs(f,var,0)

df0=subs(df,var,0)

a=0

b=XMAX

if b<inf

t=(a+b)/2

else

t=10

end

while 1

ft=subs(f,var,t)

fu=f0+c1*t*df0

if ft<=fu

dft=subs(df,var,t)

dfl=c2*t*df0

if dft>dfl

x=t

break

else

a=t

if b==inf

t=t*alpha

else

t=(a+b)/2

end

continue

end

else

b=t

t=(a+b)/2

continue

end

end

minf=subs(f,var,t)

format short

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

然后建立最速下降法的m文件如下：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%minFD.m

function [x,minf]=minFD(f,x0,var,tol)

%最速下降法

%作者：龚纯 王正林《精通Matlab最优化计算》

%目标函数：f

%初始点：x0

%自变量向量：var

%精度：tol

%所求的驻点：x

%驻点处的函数值

format long

if nargin==3

tol=1.0e-6

end

gradf=jacobian(f,var) %f的梯度

wucha=1

syms lamda

while wucha>tol

d=-subs(gradf,var,x0)

wucha=norm(d)

y=x0+lamda*d

yf=subs(f,var,y)

%l=minHJ(yf,0,100)

l=minWP(yf,10,0.2,0.6,2)

x1=x0+l*d

x0=x1

end

x=x1

minf=subs(f,var,x)

format short

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

保存后就可调用了：

>>clear

>>syms x1 x2

>>f=1.5*x1^2+0.5*x2^2-x1*x2-2*x1

f =

3/2*x1^2+1/2*x2^2-x1*x2-2*x1

>>[x,mf]=minFD(f,[-2 4],[x1,x2])

x =

1.00001.0000

mf =

-1.0000

---