无约束优化选址问题

1、无约束优化选址问题超市选址问题问题的提出：怎么选择超市的地址，使得居民区离超市所在位置距离最近。居民区位置用二维坐标表示，(X,Yi=1,2,.n 。此问题的优化模型为：n22min D = (x- Xi)(y - yj i d实例分析某投资商想在城市居民区新建一个超市，已知其五个居民区的位置坐标如下表:XY11032-41136-242145-51为使超市离居民区距离之和达到最小，该如何选择超市的位置?居民区坐标位置图绘制如下:则此问题的优化模型为2 2 ； 2 2minD= ,(x-io) (y-3)+ ,(x 4) (yi)+”6)2 (y 2)+(xR)2(y4)+(x 5)2(y-1

2、)2函数用MATLAB画出曲面图以及等高线图X,Y=meshgrid(-10:0.1:10);Z=sqrt(X-10).A2+(Y-3).A2)+sqrt(X+4).A2+(Y-11).A2)+sqrt(X-6).A2+(Y+2).A2)+.sqrt(X-2)A2+(Y-14)A2)+sqrt(X+5)A2+(Y-1)A2);surf(X ,Y,Z)shad ing xlabel( ylabel( zlabel( ti tle(in terpx );'Y')；'Z');'surf of surface')ar eurfac?X,Y=meshgri

3、d(-10:0.1:10);Z=sqrt(X-10)A2+(Y-3)A2)+sqrt(X+4)A2+(Y-11)A2)+sqrt(X-6)A2+(Y+2).A2)+ .sqrt(X-2).A2+(Y-14).A2)+sqrt(X+5).A2+(Y-1).A2);con tour(X,Y,Z,100)shad ingin terpxlabel('X');ylabel('Y');title( 'con tour of surface')10cfflitKkjr of surface6我们可以用无约束优化算法中的最速下降法进行求解Grad.mfun c

4、ti onx,val,k=grad(fu n,gfun, xO)%功能：用最速下降法求解无约束问题：min f(x)%输入：xO是初始点，fun, gfun分别是目标函数和梯度%输岀：x, val分别是近似最优点和最优值，k是迭代次数maxk=5000;%最大迭代次数rho=0.5;sigma=0.4;k=0; epsil on=1e-5;while (k<maxk)g=feval(gfu n, x0);%计算梯度d=-g;%计算搜索方向if (norm(d)<epsilon),break ; endm=0; mk=0;while (m<20) %Armijo 搜索if (f

5、eval(fun,xO+rhoAm*d)<feval(fun,xO)+sigma*rhoAm*g'*d) mk=m;break ;endm=m+1;endxO=xO+rhoAmk*d;k=k+1;endx=xO;val=feval(fu n, x0);Fun .mfun cti onf=fun(x)f=sqrt(x(1)-10)A2+(x(2)-3)A2)+sqrt(x(1)+4)A2+(x(2)-11)A2)+sqrt(x(1)-6)A2+(x(2)+2)A2)+.sqrt(x(1)-2)A2+(x(2)-14)A2)+sqrt(x(1)+5)A2+(x(2)-1)A2);Gf

6、un .mfun cti ong=gfu n(x)x =-10)A2+ (x(2)- 3)A2)A(1/2),-1)A2)A(1/2) + (2*x(2) +g= (2*x(1) + 10)/(2*(x(1) + 5)A2 + (x(2) - 1)A2)A(1/2) + (2*x(1)-4)/(2*(x(1) - 2)A2 + (x(2) - 14)A2)A(1/2) + (2*x(1) - 12)/(2*(x(1)-6)A2 + (x(2) + 2)A2)A(1/2) + (2*x(1) + 8)/(2*(x(1) + 4)人2 + (x(2)-11)A2)A(1/2)+ (2*x(1)-

7、20)/(2*(x(1)(2*x(2) - 2)/(2*(x(1) + 5)A2 + (x(2)4)/(2*(x(1)-6)A2+ (x(2)+ 2)A2)A(1/2)+ (x(2)-3)A2)A(1/2)+(2*x(2)11)A2)A(1/2)+ (2*x(2)- 28)/(2*(x(1)+ (2*x(2)- 6)/(2*(x(1)22)/(2*(x(1)+ 4)A2+-2)A2+ (x(2)- 14)A2)A(1/2).'-10)A2(x(2)we.mx0=0 0'x,val,k=grad(运行结果：'fun','gfun' ,x0)1.51

8、465.4837val =41.8050k =4444所以超市最佳的位置是(1.5146 , ,54837 )。距离为41.8050，算法迭代的次数为 次。运行时间画数名称调用总时问目用时闾*总时间圉廡色各帯二目用时间)we10.093 $0.014 $10.079 s0.075 s450X04 s0.004 s1fun890 50.000 s然后再用阻尼牛顿法解决此问题fun cti onf=fun(x)f=sqrt(x(1)-10)A2+(x(2)-3)A2)+sqrt(x(1)+4)A2+(x(2)-11)A2)+sqrt(x(1 )-6)A2+(x(2)+2)A2)+.sqrt(x(1

9、)-2)A2+(x(2)-14)A2)+sqrt(x(1)+5)A2+(x(2)-1)A2);fun cti ong=gfu n(x) g= (2*x(1) + 10)/(2*(x(1) + 5)A2 + (x(2) - 1)人2)人(1/2) + (2*x(1)- 4)/(2*(x(1) - 2)A2 + (x(2) - 14)人2)人(1/2) + (2*x(1) - 12)/(2*(x(1)-6)A2 + (x(2) + 2)A2)A(1/2) + (2*x(1) + 8)/(2*(x(1) + 4)人2 + (x(2)-11)A2)A(1/2)+ (2*x(1)- 20)/(2*(x(

10、1)- 10)A2+ (x(2)- 3)人2)人(1/2),.(2*x(2) - 2)/(2*(x(1) + 5)A2 + (x(2) - 1)人2)人(1/2) + (2*x(2) +4)/(2*(x(1)- 6)A2+ (x(2)+ 2)A2)A(1/2)+ (2*x(2)- 6)/(2*(x(1)- 10)人2+ (x(2) - 3)A2)A(1/2) + (2*x(2) - 22)/(2*(x(1) + 4)人2 + (x(2)-11)A2)A(1/2)+ (2*x(2)- 28)/(2*(x(1)- 2)人2+ (x(2)- 14)人2)人(1/2).'fun cti onh

11、e=Hess(x)n=len gth(x);he=zeros( n,n);he= 1/(x(1) + 5)A2 + (x(2) - 1)人2)人(1 /2) - (2*x(1) - 4)A2/(4*(x(1)-2)A2 + (x(2) - 14)A2)A(3/2) - (2*x -12)人2/(4*(&-6)人2 + (x(2)+ 2)A2)A(3/2)- (2*x(1)+ 8)A2/(4*(x(1)+ 4)人2+ (x(2)- 11)人2)人(3/2)-(2*x(1) - 20)A2/(4*(x(1) - 10)A2 + (x(2) - 3)人2)人(3/2) - (2*x(1) +

12、10)A2/(4*(x(1) + 5)A2 + (x(2) - 1)A2)A(3/2) + 1/(x(1) - 6)人2 + (x(2)+ 2)A2)A(1/2)+ 1/(x(1) -10)A2+ (x(2)- 3)A2)A(1/2)+ 1/(x(1)+ 4)人2+ (x(2) - 11)A2)A(1/2) + 1/(x -2)A2 + (x(2) - 14)人2)人(1/2),-(2*x + 10)*(2*x (2) - 2)/(4*(x(1) + 5)A2 + (x(2) - 1)人2)人(3/2)-(2*x -12)*(2*x (2) + 4)/(4*(x(1) - 6)人2 + (x(

13、2) + 2)人2)人(3/2)-(2*x -20)*(2*x (2) - 6)/(4*(x(1) - 10)人2 + (x(2) - 3)人2)人(3/2)-(2*x(1) + 8)*(2*x(2) - 22)/(4*(x(1) + 4)人2 + (x(2) - 11)人2)人(3/2)-(2*x(1)- 4)*(2*x(2)- 28)/(4*(x(1)- 2)人2+ (x(2)- 14)人2)人(3/2);-(2*x(1)+ 10)*(2*x(2)- 2)/(4*(x(1)+ 5)A2+ (x(2)- 1)人2)人(3/2)-(2*x -12)*(2*x (2) + 4)/(4*(x(1)

14、 - 6)人2 + (x(2) + 2)人2)人(3/2)-(2*x -20)*(2*x (2) - 6)/(4*(x(1) - 10)人2 + (x(2) - 3)人2)人(3/2)-(2*x(1) + 8)*(2*x(2) - 22)/(4*(x(1) + 4)A2 + (x(2) - 11)A2)A(3/2)-(2*x(1)- 4)*(2*x(2)- 28)/(4*(x(1)- 2)人2+ (x(2)1/(x(1) + 5)A2 + (x(2) - 1)A2)A(1/2) - (2*x(2) + 4)A2/(4*(x(1) - 6)人2+ (x(2) + 2)A2)A(3/2) - (2

15、*x(2) - 6)人2/(4*©(1) - 10)人2 + (x(2)-3)A2)A(3/2) - (2*x(2) - 22)A2/(4*(x(1) + 4)人2 + (x(2) - 11)人2)人(3/2)-(2*x(2) - 28)A2/(4*(x(1) - 2)人2 + (x(2) - 14)人2)人(3/2) - (2*x(2)-2)A2/(4*(x(1) + 5)A2 + (x(2) - 1)人2)人(3/2) + 1/(x -6)人2 + (x(2)+ 2)A2)A(1/2)+ 1/(x(1) -10)A2+ (x(2)- 3)A2)A(1/2)+ (x(2) - 11

16、)A2)A(1/2) + 1/(x -2)A2 + (x(2) - 14)人2)人(1/2);Su.m-14)A2)A(3/2),+ 1/(x(1)+ 4)A2x0=0 0'x,val,k=damp nm('fun' , 'gfun' , 'Hess' ,x0)运行结果:>> sux =1.51455.4837val =41.80504运行时间分析葩容称调用总时间目用时间*总时间囹糜色条带=目用时间）10.011 $0.002 sdampnm10.009 s0j002 sHess50.00 s0.005 s口 flirtSOj

17、OOI sOJOOI sfung0X01 e0.001 s自用时问昱指固數耗辭的町间.但不包括囲数的子固：t谊福的时闾*目用时间还包掏袍过保中产生的开销。修正牛顿法fun cti onx,val,k=revise nm(fu n,gfu n, Hess,xO)%功能：用修正牛顿法求解无约束问题：min f(x)%输入:x0 是初始点，fun, gfun, Hess分别是求%目标函数，梯度,Hesse阵的函数%输岀：x, val分别是近似最优点和最优值，k是迭代次数.n=len gth(xO); maxk=150;rho=0.55;sigma=0.4; tau=0.0;k=0; epsil on

18、=1e-5;while (k<maxk)gk=feval(gfu n, x0);% 计算梯度muk=no rm(gk)A(1+tau);Gk=feval(Hess,xO);% 计算 Hesse 阵Ak=Gk+muk*eye( n);dk=-Akgk;%解方程组Gk*dk=-gk, 计算搜索方向if (norm(gk)<epsilon),break ; end%检验终止准则m=0; mk=0;while (m<20)%用Armijo 搜索求步长if (feval(fun,x0+rhoAm*dk)<feval(fun,x0)+sigma*rhoAm*gk'*dk)b

19、reak'fun' , 'gfun' , 'Hess' ,x0)mk=m;endm=m+1;endx0=x0+rhoAmk*dk; k=k+1;endx=x0;val=feval(fu n, x);E.mx0=0 0'x,val,k=revise nm(运行结果:>> e1.51455.4837val =41.8050k =12运行时间:国救宕称调用总时问目用时间沪总时间圏蘇色各帯二目用时间）10.057 £0.037 5revisenm10j020 s0.004 sHees130JO12 s0.012 百130.0

20、03 s0.003 s1fun250401 50.001 s1自用时间昱指QU康蛊旳盯间但不包括函数盯匸十廉前时间.目用盯间辻包掏社曲星中产生的开欝.最速下降法阻尼牛顿法修正牛顿法最优点1.51465.48371.51455.48371.51455.4837最优值41.805041.805041.8050迭代次数44412运行时间0.0790.0090.020修正牛顿法是最速下降法和牛顿法的结合，从表格可以看出，阻尼牛顿法的收敛速度比修正牛顿法还要好。问题扩展n个居民区到超市的最短距离问题。其目标函数的坐标值我们用for循环实现对数据的读取。Fun .mfun cti onf=fun(x)n=6;A=1 2 3 4 5 6;B=0 1 2 3 4 5;f=0;for i=1: nf=f+sqrt(x(1)-A(i)f2+(x (2)-B(i)A2);endf；梯度函数Gfun .mfun cti ong=