输入:待排列的全部圆的半径R集合，圆的总数n 输出:最小圆排列距离Vmin，最优排列方案Path 算法: 1. 令Vmin为一个理论上绝对足够长的距离, Path=[] 2. 令ANode=zeros(n+1,n), ANode(1,:)=1:n, ANode(n+1,:)=1 (最后一行记录X中已确定部分数量) 3. 根据ANode(:,1)的信息计算出AVal和ABound(利用代价函数)对应的值 4. 当ANode非空时循环 5. --令AVal中最小值下标为loc, 则令X=ANode(1:n,loc), k=ANode(n+1,loc) 6. --如果ABound(loc)小于Vmin 7. ----对X(k+1)={1,2,...,n}-{X(1),X(2),...,X(k)}进行循环 8. ------如果k+1小于n, 即X(1:(k+1))为非叶结点 9. --------令lb为由X的前k+1个部分信息利用代价函数计算出的ABound的值 10. --------若lb小于Vmin，则把[X;k+1]加入ANode，把X对应的当前价值加入AVal，把lb加入Abound 11. ------否则 12. --------计算X的当前价值，并确定是否更新Vmin和Path 13. --从ANode中删去第loc列，从AVal中删去第loc个元素，从ABound中删去第loc个元素,写出其matlab程序

以下是对应的Matlab程序：

function [Vmin, Path] = min_circle_arrange(R)
% 待排列的全部圆的半径R集合，圆的总数n
n = length(R);
% 最小圆排列距离Vmin，最优排列方案Path
Vmin = inf;
Path = [];
% 令ANode=zeros(n+1,n), ANode(1,:)=1:n, ANode(n+1,:)=1
ANode = zeros(n+1, n);
ANode(1, :) = 1:n;
ANode(n+1, :) = 1;
% 根据ANode(:,1)的信息计算出AVal和ABound(利用代价函数)对应的值
AVal = zeros(1, n);
ABound = zeros(1, n);
for j = 1:n
    ANode_j = ANode(:, j);
    X = ANode_j(1:n);
    k = ANode_j(n+1);
    [AVal(j), ABound(j)] = calculate_val_bound(X, k, R);
end
% 当ANode非空时循环
while ~isempty(ANode)
    % 令AVal中最小值下标为loc, 则令X=ANode(1:n,loc), k=ANode(n+1,loc)
    [min_val, loc] = min(AVal);
    X = ANode(1:n, loc);
    k = ANode(n+1, loc);
    % 如果ABound(loc)小于Vmin
    if ABound(loc) < Vmin
        % 对X(k+1)={1,2,...,n}-{X(1),X(2),...,X(k)}进行循环
        for i = setdiff(1:n, X(1:k))
            % 如果k+1小于n, 即X(1:(k+1))为非叶结点
            if k+1 < n
                % 令lb为由X的前k+1个部分信息利用代价函数计算出的ABound的值
                lb = calculate_bound(X(1:k+1), R);
                % 若lb小于Vmin，则把[X;k+1]加入ANode，把X对应的当前价值加入AVal，把lb加入Abound
                if lb < Vmin
                    ANode(:, end+1) = [X;i;k+1];
                    AVal(end+1) = calculate_val(X, i, R);
                    ABound(end+1) = lb;
                end
            % 否则
            else
                % 计算X的当前价值，并确定是否更新Vmin和Path
                val = calculate_val(X, i, R);
                if val < Vmin
                    Vmin = val;
                    Path = [X;i];
                end
            end
        end
    end
    % 从ANode中删去第loc列，从AVal中删去第loc个元素，从ABound中删去第loc个元素
    ANode(:, loc) = [];
    AVal(loc) = [];
    ABound(loc) = [];
end
end

% 计算X中第k+1个圆放置在位置i时的价值
function val = calculate_val(X, i, R)
% 已放置的圆的圆心坐标
C = zeros(length(X), 2);
% 已放置的圆的半径
r = zeros(length(X), 1);
for j = 1:length(X)
    C(j, :) = [R(X(j)) + r(j), 0];
    for k = 1:j-1
        d = norm(C(j,:) - C(k,:));
        r(j) = max(r(j), r(k) + R(X(j)) - d);
    end
end
% 加上第k+1个圆
C = [C; [R(i) + r(k+1), 0]];
r = [r; R(i)];
% 计算圆排列的距离
d = 0;
for j = 1:length(X)+1
    for k = j+1:length(X)+1
        d = d + max(0, norm(C(j,:) - C(k,:)) - r(j) - r(k));
    end
end
val = d;
end

% 计算X的前k+1个部分信息利用代价函数计算出的ABound的值
function lb = calculate_bound(X, R)
% 已放置的圆的圆心坐标
C = zeros(length(X), 2);
% 已放置的圆的半径
r = zeros(length(X), 1);
for j = 1:length(X)
    C(j, :) = [R(X(j)) + r(j), 0];
    for k = 1:j-1
        d = norm(C(j,:) - C(k,:));
        r(j) = max(r(j), r(k) + R(X(j)) - d);
    end
end
% 计算圆排列的距离的下界
lb = 0;
for j = 1:length(X)
    for k = j+1:length(X)
        lb = lb + max(0, norm(C(j,:) - C(k,:)) - r(j) - r(k));
    end
end
% 加上剩余未放置的圆的最小距离
lb = lb + calculate_min_distance(setdiff(1:length(R), X), R);
end

% 计算剩余未放置的圆的最小距离
function d = calculate_min_distance(I, R)
if isempty(I)
    d = 0;
    return;
end
if length(I) == 1
    d = 2*R(I);
    return;
end
d = inf;
for i = I
    J = setdiff(I, i);
    d_i = R(i) + calculate_min_distance(J, R);
    if d_i < d
        d = d_i;
    end
end
end