无约束优化选址问题.docx
《无约束优化选址问题.docx》由会员分享，可在线阅读，更多相关《无约束优化选址问题.docx（12页珍藏版）》请在冰豆网上搜索。
无约束优化选址问题
无约束优化选址问题
超市选址问题
问题的提出：
怎么选择超市的地址，使得居民区离超市所在位置距离最近。
居民区位置用二维坐标表示，（X
Y
）i=1,2,..n。
此问题的优化模型为：
minD=
实例分析
某投资商想在城市居民区新建一个超市，已知其五个居民区的位置坐标如下表：
X
Y
1
10
3
2
-4
11
3
6
-2
4
2
14
5
-5
1
为使超市离居民区距离之和达到最小，该如何选择超市的位置？
居民区坐标位置图绘制如下：
则此问题的优化模型为
minD=
+
+
+
+
函数用MATLAB画出曲面图以及等高线图
[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）;
surf（X,Y,Z）
shadinginterp
xlabel（'X'）;
ylabel（'Y'）;
zlabel（'Z'）;
title（'surfofsurface'）
[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）;
contour（X,Y,Z,100）
shadinginterp
xlabel（'X'）;
ylabel（'Y'）;
title（'contourofsurface'）
我们可以用无约束优化算法中的最速下降法进行求解
Grad.m
function[x,val,k]=grad（fun,gfun,x0）
%功能:
用最速下降法求解无约束问题:
minf（x）
%输入:
x0是初始点,fun,gfun分别是目标函数和梯度
%输出:
x,val分别是近似最优点和最优值,k是迭代次数.
maxk=5000;%最大迭代次数
rho=0.5;sigma=0.4;
k=0;epsilon=1e-5;
while（kg=feval（gfun,x0）;%计算梯度
d=-g;%计算搜索方向
if（norm（d）m=0;mk=0;
while（m<20）%Armijo搜索
if（feval（fun,x0+rho^m*d）mk=m;break;
end
m=m+1;
end
x0=x0+rho^mk*d;
k=k+1;
end
x=x0;
val=feval（fun,x0）;
Fun.m
functionf=fun（x）
f=sqrt（（x
（1）-10）^2+（x
（2）-3）^2）+sqrt（（x
（1）+4）^2+（x
（2）-11）^2）+sqrt（（x
（1）-6）^2+（x
（2）+2）^2）+...
sqrt（（x
（1）-2）^2+（x
（2）-14）^2）+sqrt（（x
（1）+5）^2+（x
（2）-1）^2）;
Gfun.m
functiong=gfun（x）
g=[（2*x
（1）+10）/（2*（（x
（1）+5）^2+（x
（2）-1）^2）^（1/2））+（2*x
（1）-4）/（2*（（x
（1）-2）^2+（x
（2）-14）^2）^（1/2））+（2*x
（1）-12）/（2*（（x
（1）-6）^2+（x
（2）+2）^2）^（1/2））+（2*x
（1）+8）/（2*（（x
（1）+4）^2+（x
（2）-11）^2）^（1/2））+（2*x
（1）-20）/（2*（（x
（1）-10）^2+（x
（2）-3）^2）^（1/2））,...
（2*x
（2）-2）/（2*（（x
（1）+5）^2+（x
（2）-1）^2）^（1/2））+（2*x
（2）+4）/（2*（（x
（1）-6）^2+（x
（2）+2）^2）^（1/2））+（2*x
（2）-6）/（2*（（x
（1）-10）^2+（x
（2）-3）^2）^（1/2））+（2*x
（2）-22）/（2*（（x
（1）+4）^2+（x
（2）-11）^2）^（1/2））+（2*x
（2）-28）/（2*（（x
（1）-2）^2+（x
（2）-14）^2）^（1/2））].';
we.m
x0=[00]';
[x,val,k]=grad（'fun','gfun',x0）
运行结果：
x=
1.5146
5.4837
val=
41.8050
k=
44
所以超市最佳的位置是（1.5146，,54837）。
距离为41.8050，算法迭代的次数为44次。
运行时间
然后再用阻尼牛顿法解决此问题
functionf=fun（x）
f=sqrt（（x
（1）-10）^2+（x
（2）-3）^2）+sqrt（（x
（1）+4）^2+（x
（2）-11）^2）+sqrt（（x
（1）-6）^2+（x
（2）+2）^2）+...
sqrt（（x
（1）-2）^2+（x
（2）-14）^2）+sqrt（（x
（1）+5）^2+（x
（2）-1）^2）;
functiong=gfun（x）
g=[（2*x
（1）+10）/（2*（（x
（1）+5）^2+（x
（2）-1）^2）^（1/2））+（2*x
（1）-4）/（2*（（x
（1）-2）^2+（x
（2）-14）^2）^（1/2））+（2*x
（1）-12）/（2*（（x
（1）-6）^2+（x
（2）+2）^2）^（1/2））+（2*x
（1）+8）/（2*（（x
（1）+4）^2+（x
（2）-11）^2）^（1/2））+（2*x
（1）-20）/（2*（（x
（1）-10）^2+（x
（2）-3）^2）^（1/2））,...
（2*x
（2）-2）/（2*（（x
（1）+5）^2+（x
（2）-1）^2）^（1/2））+（2*x
（2）+4）/（2*（（x
（1）-6）^2+（x
（2）+2）^2）^（1/2））+（2*x
（2）-6）/（2*（（x
（1）-10）^2+（x
（2）-3）^2）^（1/2））+（2*x
（2）-22）/（2*（（x
（1）+4）^2+（x
（2）-11）^2）^（1/2））+（2*x
（2）-28）/（2*（（x
（1）-2）^2+（x
（2）-14）^2）^（1/2））].';
functionhe=Hess（x）
n=length（x）;
he=zeros（n,n）;
he=[1/（（x
（1）+5）^2+（x
（2）-1）^2）^（1/2）-（2*x
（1）-4）^2/（4*（（x
（1）-2）^2+（x
（2）-14）^2）^（3/2））-（2*x
（1）-12）^2/（4*（（x
（1）-6）^2+（x
（2）+2）^2）^（3/2））-（2*x
（1）+8）^2/（4*（（x
（1）+4）^2+（x
（2）-11）^2）^（3/2））-（2*x
（1）-20）^2/（4*（（x
（1）-10）^2+（x
（2）-3）^2）^（3/2））-（2*x
（1）+10）^2/（4*（（x
（1）+5）^2+（x
（2）-1）^2）^（3/2））+1/（（x
（1）-6）^2+（x
（2）+2）^2）^（1/2）+1/（（x
（1）-10）^2+（x
（2）-3）^2）^（1/2）+1/（（x
（1）+4）^2+（x
（2）-11）^2）^（1/2）+1/（（x
（1）-2）^2+（x
（2）-14）^2）^（1/2）,-（（2*x
（1）+10）*（2*x
（2）-2））/（4*（（x
（1）+5）^2+（x
（2）-1）^2）^（3/2））-（（2*x
（1）-12）*（2*x
（2）+4））/（4*（（x
（1）-6）^2+（x
（2）+2）^2）^（3/2））-（（2*x
（1）-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）^2+（x
（2）-1）^2）^（3/2））-（（2*x
（1）-12）*（2*x
（2）+4））/（4*（（x
（1）-6）^2+（x
（2）+2）^2）^（3/2））-（（2*x
（1）-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））,1/（（x
（1）+5）^2+（x
（2）-1）^2）^（1/2）-（2*x
（2）+4）^2/（4*（（x
（1）-6）^2+（x
（2）+2）^2）^（3/2））-（2*x
（2）-6）^2/（4*（（x
（1）-10）^2+（x
（2）-3）^2）^（3/2））-（2*x
（2）-22）^2/（4*（（x
（1）+4）^2+（x
（2）-11）^2）^（3/2））-（2*x
（2）-28）^2/（4*（（x
（1）-2）^2+（x
（2）-14）^2）^（3/2））-（2*x
（2）-2）^2/（4*（（x
（1）+5）^2+（x
（2）-1）^2）^（3/2））+1/（（x
（1）-6）^2+（x
（2）+2）^2）^（1/2）+1/（（x
（1）-10）^2+（x
（2）-3）^2）^（1/2）+1/（（x
（1）+4）^2+（x
（2）-11）^2）^（1/2）+1/（（x
（1）-2）^2+（x
（2）-14）^2）^（1/2）];
Su.m
x0=[00]';
[x,val,k]=dampnm（'fun','gfun','Hess',x0）
运行结果：
>>su
x=
1.5145
5.4837
val=
41.8050
k=
4
运行时间分析
修正牛顿法
function[x,val,k]=revisenm（fun,gfun,Hess,x0）
%功能:
用修正牛顿法求解无约束问题:
minf（x）
%输入:
x0是初始点,fun,gfun,Hess分别是求
%目标函数,梯度,Hesse阵的函数
%输出:
x,val分别是近似最优点和最优值,k是迭代次数.
n=length（x0）;maxk=150;
rho=0.55;sigma=0.4;tau=0.0;
k=0;epsilon=1e-5;
while（kgk=feval（gfun,x0）;%计算梯度
muk=norm（gk）^（1+tau）;
Gk=feval（Hess,x0）;%计算Hesse阵
Ak=Gk+muk*eye（n）;
dk=-Ak\gk;%解方程组Gk*dk=-gk,计算搜索方向
if（norm（gk）m=0;mk=0;
while（m<20）%用Armijo搜索求步长
if（feval（fun,x0+rho^m*dk）mk=m;break;
end
m=m+1;
end
x0=x0+rho^mk*dk;
k=k+1;
end
x=x0;
val=feval（fun,x）;
E.m
x0=[00]';
[x,val,k]=revisenm（'fun','gfun','Hess',x0）
运行结果:
>>e
x=
1.5145
5.4837
val=
41.8050
k=
12
运行时间：
最速下降法
阻尼牛顿法
修正牛顿法
最优点
1.5146
5.4837
1.5145
5.4837
1.5145
5.4837
最优值
41.8050
41.8050
41.8050
迭代次数
44
4
12
运行时间
0.079
0.009
0.020
修正牛顿法是最速下降法和牛顿法的结合，从表格可以看出，阻尼牛顿法的收敛速度比修正牛顿法还要好。
问题扩展
n个居民区到超市的最短距离问题。
其目标函数的坐标值我们用for循环实现对数据的读取。
Fun.m
functionf=fun（x）
n=6;
A=[123456];
B=[012345];
f=0;
fori=1:
n
f=f+sqrt（（x
（1）-A（i））^2+（x
（2）-B（i））^2）;
end
f;
梯度函数
Gfun.m
functiong=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）-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.m
x0=[00]';
[x,val,k]=grad（'fun','gfun',x0）
运行结果：
x=
3.1967
2.1967
val=
12.7279
k=
10
升级会员 