用Octave计算整数规划问题（割平面法）
广告
{{v.name}}
割平面法的核心思想：
1. 松弛：去掉整数约束，求解对应的线性规划（称为松弛问题）。
2. 构造割平面约束：若松弛解非整数，从单纯形表的行中提取整数性条件，推导一个新的线性约束。这个约束满足：
- 割去当前的非整数松弛解；
- 不割去任何整数可行解。
3. 迭代求解：将割平面约束加入原问题，形成新的LP，重复步骤1-2，直到松弛解为整数。
代码如下：
function [x_int, z_int, status] = cutting_plane(f, A, b, lb, ub, ctype, sense)
    % 割平面法求解纯整数线性规划
    % 输入：同 glpk 函数参数
    % 输出：x_int-整数最优解, z_int-最优值, status-求解状态(0=最优)

    max_iter = 100;  % 最大迭代次数，防止死循环
    iter = 0;
    param.lpsolver = 1;
    param.itlim = 100;
    vartype = repmat('C', length(f), 1);

    while iter < max_iter
        iter = iter + 1;

        % 步骤1：求解当前松弛LP
        [x_relax, z_relax, status] = glpk(f, A, b, lb, ub, ctype, vartype, sense, param);
        if status ~= 0
            disp('松弛问题无可行解');
            x_int = [];
            z_int = [];
            return;
        end

        % 步骤2：判断是否为整数解（允许微小数值误差）
        if all(abs(x_relax - round(x_relax)) < 1e-6)
            x_int = round(x_relax);
            z_int = z_relax;
            disp(['迭代 ', num2str(iter), ' 次找到整数解']);
            return;
        end

        % 步骤3：调用 glpk 的 extra 输出，获取单纯形表（简化版：直接用非基变量构造割平面）
        % 简化处理：选择第一个非整数变量对应的行构造割平面
        k = find(abs(x_relax - round(x_relax)) > 1e-6, 1);
        [m, n] = size(A);

        % 构造松弛LP的标准型，获取单纯形表行（这里简化，直接用原约束推导，实际需从单纯形表提取）
        % 注：完整实现需解析单纯形表，此处为演示，假设 x_k 对应的约束行已知
        % 简化版割平面：针对变量 x_k，直接添加 x_k ≤ floor(x_relax(k)) 或 x_k >= ceil(x_relax(k))
        % 严格割平面需从单纯形表推导，此处简化便于理解
        f_k = x_relax(k) - floor(x_relax(k));
        if f_k > 0.5
            % 添加 x_k >= ceil(x_relax(k))
            new_row = zeros(1, n);
            new_row(k) = 1;
            A = [A; new_row];
            b = [b; ceil(x_relax(k))];
            ctype = [ctype, 'D'];  % 约束类型为 >=
        else
            % 添加 x_k ≤ floor(x_relax(k))
            new_row = zeros(1, n);
            new_row(k) = 1;
            A = [A; new_row];
            b = [b; floor(x_relax(k))];
            ctype = [ctype, 'U'];  % 约束类型为 ≤
        end
    end

    disp('达到最大迭代次数，未找到整数解');
    x_int = [];
    z_int = [];
    status = -1;
end
例子：
\( \max\ Z = 10x_1 + 6x_2 + 4x_3 \\ \)
\( s.t. \begin{cases} x_1 + x_2 + x_3 \le 100 \\ 10x_1 + 4x_2 + 5x_3 \le 600 \\ 2x_1 + 2x_2 + 6x_3 \le 300 \\ x_1,x_2,x_3 \ge 0且为整数 \end{cases} \)
定义目标函数系数（max Z = 10x1 + 6x2 + 4x3），代码如下：
>> c = [10, 6, 4]';
定义约束系数矩阵 A，代码如下：
>> A = [ 1, 1, 1;
     10, 4, 5;
      2, 2, 6];
定义约束右端项 b，代码如下：
>> b = [100, 600, 300]';
定义变量下界，代码如下：
>> lb = [0, 0, 0]';
变量无上界，代码如下：
>> ub = [];
定义约束类型：3个约束都是 ≤，对应 'UUU'，代码如下：
>> ctype = "UUU";
定义优化方向：最大化模型（-1），代码如下：
>> s = -1;
求解，代码如下：
>> [x_relax, z_relax, status] = ...
   cutting_plane (c, A, b, lb, ub, ctype, s)
迭代 2 次找到整数解
x_relax =

   33
   67
    0

z_relax = 732
status = 0
友链