输入:待排列的全部圆的半径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] = circle_arrangement(R)
% R：待排列的圆的半径，大小为 n×1 的向量
% Vmin：最小圆排列距离
% Path：最优排列方案

n = length(R);
Vmin = inf; % 初始距离为无穷大
Path = []; % 初始最优排列方案为空

% 计算代价函数的值
function [val, bound] = cost_func(X, k)
    if k == n % 叶子节点
        val = sum(X) - 1; % 当前价值为已选圆的数量减1
        bound = val;
    else % 非叶子节点
        val = sum(X(1:k)); % 当前价值为已选圆的数量
        remain_R = R(~X(1:k)); % 剩余圆的半径
        [~, idx] = sort(remain_R, 'ascend'); % 按半径从小到大排序
        i = 1;
        cap = 0;
        while i <= length(remain_R) && cap + remain_R(idx(i)) <= 1 - sum(X(1:k))
            cap = cap + remain_R(idx(i));
            i = i + 1;
        end
        if i <= length(remain_R)
            bound = val + cap/remain_R(idx(i));
        else
            bound = val + cap/remain_R(end);
        end
    end
end

ANode = zeros(n+1, n); % 待处理节点矩阵
ANode(1,:) = 1:n; % 第一行是圆的编号
ANode(n+1,:) = ones(1,n); % 最后一行是已确定部分数量
[AVal, ABound] = cost_func(ANode(1,:), ANode(n+1,:)); % 计算初始节点的价值和上界

while ~isempty(ANode)
    [~, loc] = min(ABound); % 找到价值上界最小的节点
    X = ANode(1:n, loc); % 当前节点对应的排列方案
    k = ANode(n+1, loc); % 已确定部分的数量
    if ABound(loc) >= Vmin % 剪枝
        ANode(:, loc) = [];
        AVal(loc) = [];
        ABound(loc) = [];
    else
        for i = setdiff(1:n, X(1:k)) % 枚举下一个圆的编号
            next_X = [X; i]; % 加入下一个圆
            next_k = k + 1;
            [next_val, next_bound] = cost_func(next_X, next_k); % 计算价值和上界
            if next_k < n && next_bound < Vmin % 非叶子节点且上界小于当前最优解
                next_ANode = [next_X; next_k]; % 新节点
                ANode = [ANode, next_ANode]; % 加入新节点
                AVal = [AVal, next_val]; % 记录新节点的价值
                ABound = [ABound, next_bound]; % 记录新节点的上界
            elseif next_k == n % 叶子节点
                next_val = next_val + sum(next_X);
                if next_val < Vmin % 更新最优解
                    Vmin = next_val;
                    Path = next_X;
                end
            end
        end
        ANode(:, loc) = [];
        AVal(loc) = [];
        ABound(loc) = [];
    end
end

end

注：该代码仅供参考，实际使用时需要根据具体情况进行修改。