1、无约束优化选址问题无约束优化选址问题超市选址问题 问题的提出：怎么选择超市的地址，使得居民区离超市所在位置距离最近。居民区位置用二维坐标表示，（X,Y）i=1,2,.n 。此问题的优化模型为： min D = 实例分析某投资商想在城市居民区新建一个超市，已知其五个居民区的位置坐标如下表：XY11032-41136-242145-51为使超市离居民区距离之和达到最小，该如何选择超市的位置？居民区坐标位置图绘制如下： 则此问题的优化模型为 minD=+函数用MATLAB画出曲面图以及等高线图X,Y=meshgrid(-10:0.1:10); Z=sqrt(X-10).2+(Y-3).2)+sqrt
2、(X+4).2+(Y-11).2)+sqrt(X-6).2+(Y+2).2)+. sqrt(X-2).2+(Y-14).2)+sqrt(X+5).2+(Y-1).2); surf(X,Y,Z) shading interpxlabel(X);ylabel(Y);zlabel(Z);title(surf of surface)X,Y=meshgrid(-10:0.1:10); Z=sqrt(X-10).2+(Y-3).2)+sqrt(X+4).2+(Y-11).2)+sqrt(X-6).2+(Y+2).2)+. sqrt(X-2).2+(Y-14).2)+sqrt(X+5).2+(Y-1).2)
3、;contour(X,Y,Z,100)shading interpxlabel(X);ylabel(Y);title(contour of surface)我们可以用无约束优化算法中的最速下降法进行求解Grad.mfunction x,val,k=grad(fun,gfun,x0)% 功能: 用最速下降法求解无约束问题: min f(x)%输入: x0是初始点, fun, gfun分别是目标函数和梯度%输出: x, val分别是近似最优点和最优值, k是迭代次数.maxk=5000; %最大迭代次数rho=0.5;sigma=0.4;k=0; epsilon=1e-5;while(kmaxk)
4、 g=feval(gfun,x0); %计算梯度 d=-g; %计算搜索方向 if(norm(d)epsilon), break; end m=0; mk=0; while(m20) %Armijo搜索 if(feval(fun,x0+rhom*d) sux = 1.5145 5.4837val = 41.8050k = 4运行时间分析修正牛顿法function x,val,k=revisenm(fun,gfun,Hess,x0)%功能: 用修正牛顿法求解无约束问题: min f(x)%输入: x0是初始点, fun, gfun, Hess 分别是求% 目标函数,梯度,Hesse 阵的函数%输
5、出: x, val分别是近似最优点和最优值, k是迭代次数.n=length(x0); maxk=150;rho=0.55;sigma=0.4; tau=0.0;k=0; epsilon=1e-5;while(kmaxk) gk=feval(gfun,x0); % 计算梯度 muk=norm(gk)(1+tau); Gk=feval(Hess,x0); % 计算Hesse阵 Ak=Gk+muk*eye(n); dk=-Akgk; %解方程组Gk*dk=-gk, 计算搜索方向 if(norm(gk)epsilon), break; end %检验终止准则 m=0; mk=0; while(m20
6、) %用Armijo搜索求步长 if(feval(fun,x0+rhom*dk) ex = 1.5145 5.4837val = 41.8050k =12运行时间：最速下降法阻尼牛顿法修正牛顿法最优点1.5146 5.48371.51455.48371.51455.4837最优值41.805041.805041.8050迭代次数44412运行时间0.0790.0090.020修正牛顿法是最速下降法和牛顿法的结合，从表格可以看出，阻尼牛顿法的收敛速度比修正牛顿法还要好。问题扩展n个居民区到超市的最短距离问题。其目标函数的坐标值我们用for循环实现对数据的读取。Fun.mfunction f=fu
7、n(x)n=6;A=1 2 3 4 5 6;B=0 1 2 3 4 5;f=0;for i=1:n f=f+sqrt(x(1)-A(i)2+(x(2)-B(i)2); endf;梯度函数Gfun.mfunction g=gfun(x)g= (2*x(1) - 2)/(2*(x(1) - 1)2 + x(2)2)(1/2) + (2*x(1) - 4)/(2*(x(1) - 2)2 + (x(2) - 1)2)(1/2) + (2*x(1) - 6)/(2*(x(1) - 3)2 + (x(2) - 2)2)(1/2) + (2*x(1) - 8)/(2*(x(1) - 4)2 + (x(2)
8、- 3)2)(1/2) + (2*x(1) - 10)/(2*(x(1) - 5)2 + (x(2) - 4)2)(1/2) + (2*x(1) - 12)/(2*(x(1) - 6)2 + (x(2) - 5)2)(1/2),.x(2)/(x(1) - 1)2 + x(2)2)(1/2) + (2*x(2) - 2)/(2*(x(1) - 2)2 + (x(2) - 1)2)(1/2) + (2*x(2) - 4)/(2*(x(1) - 3)2 + (x(2) - 2)2)(1/2) + (2*x(2) - 6)/(2*(x(1) - 4)2 + (x(2) - 3)2)(1/2) + (2*x(2) - 8)/(2*(x(1) - 5)2 + (x(2) - 4)2)(1/2) + (2*x(2) - 10)/(2*(x(1) - 6)2 + (x(2) - 5)2)(1/2).;u.mx0=0 0;x,val,k=grad(fun,gfun,x0)运行结果：x = 3.1967 2.1967val = 12.7279k = 10
