Question

一道MATLAB题，求助啊

怎么用MATLAB写出这个算法???
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
Weiszfeld 算法如图所示

如果出现下图二所示的情况，则运算终止（f是一个特定函数,不等式右边是极小量10^(-6)）





Answer 1

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

其实这是一个单一设施选址问题，其中的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([x;y]),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 算法

Answer 2

QQ号给我，我知道