运筹学——使用matlab实现整数规划割平面法

题目

代码

A = [-1 3 1 0; 7 1 0 1]; 
b = [6 35]'; 
c = [7 9 0 0];
[xstar,fxstar,iter] = Gomory(A,b,c)

function [xstar,fxstar,iter] = Gomory(A,b,c)
format rat
%UNTITLED 此处显示有关此函数的摘要
iter=0;%初始化迭代次数
while true
    [m,n]=size(A);%A矩阵大小
    if min(b)>=0
        [x_opt,fx_opt,CA,Cb] = simplex(A,b,c);
    else
        [x_opt,fx_opt,CA,Cb] = DSimplex_eye(A,b,c);
    end
    %判断是否已经解出了整数最优解
    flag_zhengshu=1;%flag_zhengshu初始化为1，只要有一个不是整数,flag_zhengshu赋值0
    for pos_x = 1:m
        if abs(round(x_opt(pos_x))-x_opt(pos_x))>=1e-3%判断整数条件
            flag_zhengshu=0;
            break;
        end
    end
    if flag_zhengshu==1 %如果解全是整数，满足条件，循环结束
        xstar=x_opt;
        fxstar=fx_opt;
        break;
    end
    iter=iter+1;
    %否则增加约束条件的行
    %找出b中和整数相差最大的数
    %循环遍历
    cha=0;
    row=0;
    for r=1:m
        t=abs(floor(x_opt(r))-x_opt(r));
        if t>cha
            cha=t;
            row=r;%标记当前最大差值的位置
        end
    end
    n=n+1;
    m=m+1;
    iter=iter+1;
    %更新矩阵系数，在原基础上增加一行一列，第（m,n）=1
    tmp_A=zeros(m,n);
    tmp_b=zeros(m,1);
    tmp_c=zeros(1,n);
    for i=1:m-1
        for j=1:n-1
            tmp_A(i,j)=CA(i,j);
        end
        tmp_b(i,1)=Cb(i,1);
    end
    tmp_b(m,1)=floor(Cb(row,1))-Cb(row,1);
    tmp_A(m,n)=1;
    for i =1:n-1
        tmp_c(1,i)=c(i);
    end
    %加上约束条件
    for i=1:n-1
        if tmp_A(row,i)==0
            tmp_A(m,i)=0;
        else
            tmp_A(m,i)=floor(tmp_A(row,i))-tmp_A(row,i);
        end
    end
    A=tmp_A;
    b=tmp_b;
    c=tmp_c;
end         
end

function [x_opt,fx_opt,A,b] = simplex(A,b,c)
% 单纯形法求解标准形线性规划问题: max cx s.t. Ax=b x>=0
% x_opt,fx_opt,iter(最优解、最优函数值、迭代次数）
format rat %元素使用分数表示
[m,n] = size(A);              %m约束条件个数, n决策变量数
ind_B =has_ones(A);
ind_N = setdiff(1:n, ind_B);  %非基变量的索引,原理是返回在1:n中出现而不在ind_B即基变量索引中出现的元素，并从小到大排序
% 循环求解
while true
    x0 = zeros(n,1);
    x0(ind_B) = b;               %初始基可行解，令基变量的位置为方程组右边系数，非基变量取值为0
    cB = c(ind_B);               %计算cB，即目标函数在基变量处对应的系数
    Sigma = zeros(1,n);          %Sigma为检验数向量
    Sigma(ind_N) = c(ind_N) - cB*A(:,ind_N);   %计算检验数（非基变量），因为基变量对应的初始检验数一定为0
    [~, k] = max(Sigma);         %选出最大检验数, 确定进基变量索引k；~表示忽略第一个参数（即最大值），k是索引
    if ~any(Sigma > 0)           %所有检验数都小于0，此基可行解为最优解, any表示存在某个检验数>0        
        x_opt = x0;
        fx_opt = c * x_opt;      %算出最优解
        return
    end
    if all(A(:,k) <= 0)          %表示检验数这一列每个数都<=0，有无界解
        x_opt = [];
        break
    end
    Theta = b ./ A(:,k);         %计算θ（点除，即矩阵中对应元素相除，得到一个新的矩阵）
    Theta(Theta<=0) = 10000;
    [~, q] = min(Theta);         %选出最小θ
    el = ind_B(q);               %确定出基变量在系数矩阵中的列索引el, 主元为A(q,k)
    % 换基
    ind_B(ind_B == el) = k;      %新的基变量索引
    ind_N = setdiff(1:n, ind_B); %非基变量索引
    % 更新A和b
    A(:,ind_N) = A(:,ind_B) \ A(:,ind_N); %基矩阵的逆乘以非基矩阵
    b = A(:,ind_B) \ b;  %基矩阵的逆乘以b
    A(:,ind_B) = eye(m,m);  %基矩阵更新为单位矩阵
end
end

function [x_opt,fx_opt,A,b] = DSimplex_eye(A,b,c)
% 对偶单纯形法求解标准形线性规划问题: max cx s.t. Ax = b x>=0
% 输入参数: c为目标函数系数, A为约束方程组系数矩阵, b为约束方程组常数项
% 输出参数: x_opt为最求解，fx_opt为最优函数值，iter为迭代次数
format rat                    %元素使用分数表示
[m,n] = size(A);              %m约束条件个数, n决策变量数
ind_B = has_ones(A);
ind_N = setdiff(1:n, ind_B);  %非基变量的索引,原理是返回在1:n中出现而不在ind_B即基变量索引中出现的元素，并从小到大排序
while true
    x0 = zeros(n,1);
    x0(ind_B) = b;                %初始基可行解
    cB = c(ind_B);                %计算cB
    if ~any(b < 0)                %此基可行解为最优解, any存在某个<0        
        x_opt = x0;
        fx_opt = c*x_opt;
        return
    end
    index=find(b<0);
    for i = 1:numel(index)
        if all(A(index(i),:)>=0)
            x_opt=[];
            fx_opt = [];         %此时原问题无可行解，对偶问题存在无界解
            return
        end
    end
    Sigma = zeros(1,n);
    Sigma(ind_N) = c(ind_N) - cB*A(:,ind_N);   %计算检验数
    [~,q] = min(b);               %选出最小的b
    r = ind_B(q);                 %确定出基变量索引r
    Theta = Sigma ./ A(q,:);      %计算θ
    Theta(Theta<=0) = 10000;
    [~,s] = min(Theta);           %确定进基变量索引s, 主元为A(q,s)
    % 换基
    ind_B(ind_B == r) = s;        %新的基变量索引
    ind_N = setdiff(1:n, ind_B);  %非基变量索引
    % 更新A和b
    A(:,ind_N) = A(:,ind_B) \ A(:,ind_N);
    b = A(:,ind_B) \ b;
    A(:,ind_B) = eye(m,m);
end
end

function [index]=has_ones(a)  %判断矩阵是否含有单位阵
    [m,n]=size(a);
    b=min(m,n);
    one=eye(b);
    index=[];
    result=zeros(m,1);
    for i =1:n
        for j=1:b
            if a(:,i)==one(:,j)
                index=[index,i];
                break;
            end
        end
    end
    for k=index
        result=result+a(:,k);
    end
    if result ~= ones(m,1)
       index=[];
    end
end

结果

>> Gomoryhomework

xstar =

       4       
       3       
       1       
       4       
       0       
       0       


fxstar =

      55       


iter =

       4       

作者创作不易，麻烦各位看官点赞收藏转发哈~