一道MATLAB题，求助啊

题主并没有把问题描述清楚。

其实这是一个单一设施选址问题，其中的ai对应的是平樱桐州面上点的坐标，wi为各点的权重。所谓【f是一个特定函数】说的很含糊，其实f是所选点与各已知点距离的加权和，而迭代的目标则是让f达到最小值。

这是一个无约束优化问题，可用fminunc直接求解：

% 常数定义

a1 = [2 0] a2 = [4 10] a3 = [8 3] a4 = [3 5] a5 = [5 7]

w1 = 6 w2 = 2 w3 = 7 w4 = 5 w5 = 4

e = 1e-6

% 为便于统一处理，使用数组

A = [a1 a2 a3 a4 a5]

W = [w1 w2 w3 w4 w5]

% 初值

x0 = rand(2,1) * 10

% 定义求范数和目标函数的匿名函数

n = @ (x) arrayfun( @ (i) norm( A(:,i) - x), 1:length(W) )

f = @ (x) sum( W .* n(x) ) 

x=fminunc(f,x0,optimset('LargeScale','off', 'display', 'off'))

如果用题主贴出的来式子编程序，代码如下：

% 常数定义

a1 = [2 0] a2 = [4 10] a3 = [8 3] a4 = [3 5] a5 = [5 7]

w1 = 6 w2 = 2 w3 = 7 w4 = 5 w5 = 4

e = 1e-6

% 为便于统一处理，使用数组

A = [a1 a2 a3 a4 a5]

W = [w1 w2 w3 w4 w5]

% 初值

x0 = rand(2,1) * 10

% 定义求范数和目标函数的匿名函数

n = @ (x) arrayfun( @ (i) norm( A(:,i) - x), 1:length(W) )

f = @ (x) sum( W .* n(x) )

% 迭代

while true

    N = n(x0)

    x = sum( [W W] .* A ./ [N N], 2 ) / sum( W ./ N, 2 )

    

    % 若满足精度要求，退出循环

    if abs( f(x) - f(x0) ) <= f(x0) * e

        break

    end

    x0 = x

end

x

其实，迭代的终止条件除题主所给的目标函数相对误差（TolFun）之外，常见的还有变量本身的允许误差（TolX）以及迭代次数（MaxIter）等。

下面给出一个增加了绘图示意的版本。等值线为目标函数，各已知点根据脊蔽权重以不同大小表示，动画示意轮颂了从不同初值经迭代到达最优点的过程。

% 常数定义

a1 = [2 0] a2 = [4 10] a3 = [8 3] a4 = [3 5] a5 = [5 7]

w1 = 6 w2 = 2 w3 = 7 w4 = 5 w5 = 4

e = 1e-6

% 为便于统一处理，使用数组

A = [a1 a2 a3 a4 a5]

W = [w1 w2 w3 w4 w5]

% 初值（随机生成）

x0 = rand(2,1) * 10

% 定义求范数和目标函数的匿名函数

n = @ (x) arrayfun( @ (i) norm( A(:,i) - x), 1:length(W) )

f = @ (x) sum( W .* n(x) )

% 绘图

clf

ezcontour(@(x,y)arrayfun(@(x,y)f([xy]),x,y),[-3 12 -1 11])

hold on

for i = 1 : length(W)

    plot(A(1,i), A(2,i), 'o', 'MarkerSize', W(i)*2, 'MarkerFace', 'c' )

    text(A(1,i)+0.4, A(2,i), sprintf('(%g, %g)', A(:,i)))

end

h2 = plot(x0(1), x0(2), ':.', 'color', [1 1 1]*0.8)

h1 = plot(x0(1), x0(2), 'r*')

h3 = title('')

axis equal

colorbar

X = x0

% 迭代

while true

    N = n(x0)

    x = sum( [W W] .* A ./ [N N], 2 ) / sum( W ./ N, 2 )

    

    % 更新绘图

    X = [X x]

    set(h2, 'XData', X(1,:), 'YData', X(2,:))

    set(h1, 'XData', x(1), 'YData', x(2))

    set(h3, 'str', sprintf('Current point: (%.4f, %.4f), f(x) = %.6f', x, f(x)) )

    drawnow

    pause(0.5)

    

    % 若满足精度要求，退出循环，否则以新的值继续迭代

    if abs( f(x) - f(x0) ) <= f(x0) * e

        set(h1, 'Color', 'g', 'Marker', 'p', 'MarkerFace', 'g')

        break

    end

    x0 = x

end

参考：

Weiszfeld 算法

1、最短蠢伍路问题

两个指定顶点之间的最短路径敬圆。

例如，给出了一个连接若干个城镇的铁路网络，在这个网络的两个指定城镇间，找一条最短铁路线。

以各城镇为图G的顶点，两城镇间的直通铁路为图G相应两顶点间的边，得图G。对G的每一边e，赋以一个实数)(ew—直通铁路的长度，称为e的权，得到赋权图G。G的子图的权是指子图的各边的权和。问题就是求赋权图G中指定的两个顶点00,vu间的具最小权的轨。这条轨叫做00,vu间的最短路，它的权叫做00,vu间的距离，亦记作),(00vud。

求最短路已有成熟的算法：迪克斯特拉（Dijkstra）算法，其基本思想是按距0u从近到远为顺序，依次求得0u到G的各顶点的最短路和距离，直至0v（或直至G的所有顶点），算法结束。为避免重复并保留每一步的计算信息，采用了标号算法。下面是该算法。

(i) 令0)(0ul，对0uv，令)(vl，}{00uS，0i。

(ii) 对每个iSv（iiSVS\），用

)}

()(),({minuvwulvli

Su

代替)(vl。计算)}({minvli

Sv，把达到这个最小值的一个顶点记为1iu，令

《计量地理学》（徐建华，高等教育出版社，2005）配套实习指导

140

}{11iiiuSS。

(iii). 若1||Vi，停止；若1||Vi，亮档塌用1i代替i，转(ii)。

算法结束时，从0u到各顶点v的距离由v的最后一次的标号)(vl给出。在v进入iS之前的标号)(vl叫T标号，v进入iS时的标号)(vl叫P标号。算法就是不断修改各项点的T标号，直至获得P标号。若在算法运行过程中，将每一顶点获得P标号所由来的边在图上标明，则算法结束时，0u至各项点的最短路也在图上标示出来了。

2、选址问题－以中位点选址为例

中位点选址问题的质量判据为：使最佳选址为止所在的定点到网络图中其他顶点的最短路径距离的总和（或者以各个顶点的载荷加权求和）达到最小。

function [ROUTES,PL,Tau]=ACASP(G,Tau,K,M,S,E,Alpha,Beta,Rho,Q)

%% ---------------------------------------------------------------

% ACASP.m

% 蚁群算法动态寻路算法

% ChengAihua,PLA Information Engineering University,ZhengZhou,China

% Email:aihuacheng@gmail.com

% All rights reserved

%% ---------------------------------------------------------------

% 输入参数列表

% G地形图为01矩阵，如果为1表示障碍物

% Tau 初始信息素矩阵（认为前面的觅食活动中有残留的信息素）

% K迭代次数（指蚂蚁出动多少波）

% M蚂蚁个数（每一波蚂蚁有多少个）

% S起始点（最短路径的起始点）

% E终止点（最短路径的目的点）

% Alpha表征信息素重要程度的参数

% Beta 表征启发式因子重要程度的参数

% Rho 信息素蒸发系数

% Q信息素增加强度系数

%

% 输出参数列表

% ROUTES 每一代的每一只蚂蚁的爬行路线

% PL 每一代的每一只蚂蚁的爬行路线长度

% Tau 输出动态修正过的信息素

%% --------------------变量初始化----------------------------------

%load

D=G2D(G)

N=size(D,1)%N表示问题的规模（象素个数）

MM=size(G,1)

a=1%小方格象素的边长

Ex=a*(mod(E,MM)-0.5)%终止点横坐闭绝标

if Ex==-0.5

Ex=MM-0.5

end

Ey=a*(MM+0.5-ceil(E/MM))%终止点纵坐标

Eta=zeros(1,N)%启发式信息，取为至目标点的直线距离的倒数

%下面构造启发式信息矩阵

for i=1:N

if ix==-0.5

ix=MM-0.5

end

iy=a*(MM+0.5-ceil(i/MM))

if i~=E

Eta(1,i)=1/((ix-Ex)^2+(iy-Ey)^2)^0.5

else

Eta(1,i)=100

end

end

ROUTES=cell(K,M)%用细胞结构存储每如态液一代的每一只蚂蚁渣物的爬行路线

PL=zeros(K,M)%用矩阵存储每一代的每一只蚂蚁的爬行路线长度

%% -----------启动K轮蚂蚁觅食活动，每轮派出M只蚂蚁--------------------

for k=1:K

disp(k)

for m=1:M

%% 第一步：状态初始化

W=S%当前节点初始化为起始点

Path=S%爬行路线初始化

PLkm=0%爬行路线长度初始化

TABUkm=ones(1,N)%禁忌表初始化

TABUkm(S)=0%已经在初始点了，因此要排除

DD=D%邻接矩阵初始化

%% 第二步：下一步可以前往的节点

DW=DD(W,:)

DW1=find(DW

for j=1:length(DW1)

if TABUkm(DW1(j))==0

DW(j)=inf

end

end

LJD=find(DW

Len_LJD=length(LJD)%可选节点的个数

%% 觅食停止条件：蚂蚁未遇到食物或者陷入死胡同

while W~=E&&Len_LJD>=1

%% 第三步：转轮赌法选择下一步怎么走

PP=zeros(1,Len_LJD)

for i=1:Len_LJD

PP(i)=(Tau(W,LJD(i))^Alpha)*(Eta(LJD(i))^Beta)

end

PP=PP/(sum(PP))%建立概率分布

Pcum=cumsum(PP)

Select=find(Pcum>=rand)

%% 第四步：状态更新和记录

Path=[Path,to_visit]%路径增加

PLkm=PLkm+DD(W,to_visit)%路径长度增加

W=to_visit%蚂蚁移到下一个节点

for kk=1:N

if TABUkm(kk)==0

DD(W,kk)=inf

DD(kk,W)=inf

end

end

TABUkm(W)=0%已访问过的节点从禁忌表中删除

for j=1:length(DW1)

if TABUkm(DW1(j))==0

DW(j)=inf

end

end

LJD=find(DW

Len_LJD=length(LJD)%可选节点的个数

end

%% 第五步：记下每一代每一只蚂蚁的觅食路线和路线长度

ROUTES{k,m}=Path

if Path(end)==E

PL(k,m)=PLkm

else

PL(k,m)=inf

end

end

%% 第六步：更新信息素

Delta_Tau=zeros(N,N)%更新量初始化

for m=1:M

if PL(k,m)ROUT=ROUTES{k,m}

TS=length(ROUT)-1%跳数

PL_km=PL(k,m)

for s=1:TS

x=ROUT(s)

Delta_Tau(x,y)=Delta_Tau(x,y)+Q/PL_km

Delta_Tau(y,x)=Delta_Tau(y,x)+Q/PL_km

end

end

end

Tau=(1-Rho).*Tau+Delta_Tau%信息素挥发一部分，新增加一部分

end

%% ---------------------------绘图--------------------------------

plotif=1%是否绘图的控制参数

if plotif==1

%绘收敛曲线

meanPL=zeros(1,K)

minPL=zeros(1,K)

for i=1:K

PLK=PL(i,:)

Nonzero=find(PLK

PLKPLK=PLK(Nonzero)

meanPL(i)=mean(PLKPLK)

minPL(i)=min(PLKPLK)

end

figure(1)

plot(minPL)

hold on

plot(meanPL)

grid on

title('收敛曲线（平均路径长度和最小路径长度）')

xlabel('迭代次数')

ylabel('路径长度')

%绘爬行图

figure(2)

axis([0,MM,0,MM])

for i=1:MM

for j=1:MM

if G(i,j)==1

x1=j-1y1=MM-i

x2=jy2=MM-i

x3=jy3=MM-i+1

x4=j-1y4=MM-i+1

fill([x1,x2,x3,x4],[y1,y2,y3,y4],[0.2,0.2,0.2])

hold on

else

x1=j-1y1=MM-i

x2=jy2=MM-i

x3=jy3=MM-i+1

x4=j-1y4=MM-i+1

fill([x1,x2,x3,x4],[y1,y2,y3,y4],[1,1,1])

hold on

end

end

end

hold on

ROUT=ROUTES{K,M}

LENROUT=length(ROUT)

Rx=ROUT

Ry=ROUT

for ii=1:LENROUT

Rx(ii)=a*(mod(ROUT(ii),MM)-0.5)

if Rx(ii)==-0.5

Rx(ii)=MM-0.5

end

Ry(ii)=a*(MM+0.5-ceil(ROUT(ii)/MM))

end

plot(Rx,Ry)

end

plotif2=1%绘各代蚂蚁爬行图

if plotif2==1

figure(3)

axis([0,MM,0,MM])

for i=1:MM

for j=1:MM

if G(i,j)==1

x1=j-1y1=MM-i

x2=jy2=MM-i

x3=jy3=MM-i+1

x4=j-1y4=MM-i+1

fill([x1,x2,x3,x4],[y1,y2,y3,y4],[0.2,0.2,0.2])

hold on

else

x1=j-1y1=MM-i

x2=jy2=MM-i

x3=jy3=MM-i+1

x4=j-1y4=MM-i+1

fill([x1,x2,x3,x4],[y1,y2,y3,y4],[1,1,1])

hold on

end

end

end

for k=1:K

PLK=PL(k,:)

minPLK=min(PLK)

pos=find(PLK==minPLK)

m=pos(1)

ROUT=ROUTES{k,m}

LENROUT=length(ROUT)

Rx=ROUT

Ry=ROUT

for ii=1:LENROUT

Rx(ii)=a*(mod(ROUT(ii),MM)-0.5)

if Rx(ii)==-0.5

Rx(ii)=MM-0.5

end

Ry(ii)=a*(MM+0.5-ceil(ROUT(ii)/MM))

end

plot(Rx,Ry)

hold on

end

end

将上述算法应用于机器人路径规划，优化效果如下图所示

---