数学建模学习---整数规划1


前言

例如:
—这边学习整数规划的定义,分类,当然也有大家最关心的各种代码解法。

一、整数规划是什么?

1.1 定义

  1. 规划中的变量(部分或全部)限制为整数时,称为整数规划。若在线性规划模型中,变量限制为整数,则称为整数线性规划**。目前所流行的求解整数规划的方法,往往只适用于整数线性规划。目前还没有一种方法能有效地求解一切整数规划。

1.2 整数规划的分类

  1. 变量全限制为整数时,称纯(完全)整数规划
  2. 变量部分限制为整数的,称混合整数规划。
  3. 0-1规划,所有决策变量只能取0或1两个整数

1.3 整数规划特点

在这里插入图片描述
其中(2)可理解为x1或X2等的限制区间为(0,1),所以不出现解即无可行解。
其中2是因为简单取整后的值可能不满足约束条件。

二、求解方法

1.求解方法分类

(i)分枝定界法—可求纯或混合整数线性规划。
(ii)割平面法—可求纯或混合整数线性规划。
(iii)隐枚举法—求解“0-1”整数规划:
①过滤隐枚举法;
②分枝隐枚举法。
(iv)匈牙利法—解决指派问题(“0-1”规划特殊情形)。
(v)蒙特卡洛法—求解各种类型规划。
下面将简要介绍常用的几种求解整数规划的方法。

2.1分枝定界法

1.教材: 对有约束条件的最优化问题(其可行解为有限数)的所有可行解空间恰当地进行系统搜索,这就是分枝与定界内容。通常,把全部可行解空间反复地分割为越来越小的子集,称为分枝;并且对每个子集内的解集计算一个目标下界(对于最小值问题),这称为定界。在每次分枝后,凡是界限超出已知可行解集目标值的那些子集不再进一步分枝,样,许多子集可不予考虑,这称剪枝。这就是分枝定界法的主要思路。(看不懂没关系,接下来通俗的讲讲)
2. 通俗:当你用线性规划的方法求出的X中不管是X1还是X2只要有出现不是整数的,那你就选一个不是整数的X,往上取一个最接近的整数或者往下取一个最接近的整数作为一个新的限制条件(例如X1=2.7,往下取一个整数的获得一个上限的限制条件:X1≤2,往上取获得一个下限的限制条件:X2≥3)一直循环,一直增加限制条件即边界,直到算出来X全是整数为止,不就是所谓的分枝定界法嘛。

2.2分枝定界法的matlab代码

  1. 使用intprog()函数进行求解,intprog()在matlab工具箱中不存在需要自己创建
  2. 下面是创建代码:
    首先创建一个branchboun()函数,具体的我就不讲了直接用就行了
 function [newx,newfval,status,newbound] = branchbound(f,A,B,I,x,fval,bound,Aeq,Beq,lb,ub,e)

% 分支定界法求解整数规划
% f,A,B,Aeq,Beq,lb,ub与线性规划相同
% I为整数限制变量的向量
% x为初始解,fval为初始值

options = optimset('display','off');
[x0,fval0,status0]=linprog(f,A,B,Aeq,Beq,lb,ub,[],options);

%递归中的最终退出条件
%无解或者解比现有上界大则返回原解
if status0 <= 0 || fval0 >= bound
    newx = x;
    newfval = fval;
    newbound = bound;
    status = status0;
    return;
end

%是否为整数解,如果是整数解则返回
intindex = find(abs(x0(I) - round(x0(I))) > e);
if isempty(intindex) %判断是否为空值
    newx(I) = round(x0(I));
    newfval = fval0;
    newbound = fval0;
    status = 1;
    return;
end

%当有非整可行解时,则进行分支求解
%此时必定会有整数解或空解
%找到第一个不满足整数要求的变量
n = I(intindex(1));
addA = zeros(1,length(f));
addA(n) = 1;
%构造第一个分支 x<=floor(x(n))
A = [A;addA];
B = [B,floor(x(n))];%向下取整
[x1,fval1,status1,bound1] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);
A(end,:) = [];
B(:,end) = [];
%解得第一个分支,若为更优解则替换,若不是则保持原状

status = status1;
if status1 > 0 && bound1 < bound
    newx = x1;
    newfval = fval1;
    bound = fval1;
    newbound = bound1;
else
    newx = x0;
    newfval = fval0;
    newbound = bound;
end

%构造第二分支
A = [A;-addA];
B = [B,-ceil(x(n))];%向上取整
[x2,fval2,status2,bound2] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);    
A(end,:) = [];
B(:,end) = [];

%解得第二分支,并与第一分支做比较,如果更优则替换
if status2 > 0 && bound2 < bound
    status = status2;
    newx = x2;
    newfval = fval2;
    newbound = bound2;
end

然后创建一个intprog()函数

function [x,fval,status] = intprog(f,A,B,I,Aeq,Beq,lb,ub,e)
%整数规划求解函数 intprog()
%     其中 f为目标函数向量
%     A和B为不等式约束 Aeq与Beq为等式约束
%     I为整数约束
%     lb与ub分别为变量下界与上界
%     x为最优解,fval为最优值
%例子:
%        maximize 20 x1 + 10 x2 
%        S.T.
% 	             5 x1 + 4 x2 <=24
%                2 x1 + 5 x2 <=13
%                   x1, x2 >=0 
%                   x1, x2是整数
% f=[-20, -10];
% A=[ 5  4; 2 5];
% B=[24; 13];
% lb=[0 0];
% ub=[inf inf];
% I=[1,2];
% e=0.000001;
% [x v s]= IP(f,A,B,I,[],[],lb,ub,,e)
% x = 4     1  v = -90.0000   s = 1

% 控制输入参数
if nargin < 9, e = 0.00001;
   if nargin < 8, ub = []; 
      if nargin < 7, lb = []; 
         if nargin < 6, Beq = []; 
            if nargin < 5, Aeq = [];
               if nargin < 4, I = [1:length(f)];
               end, end, end, end, end, end
  
%求解整数规划对应的线性规划,判断是否有解
options = optimset('display','off');
[x0,fval0,exitflag] = linprog(f,A,B,Aeq,Beq,lb,ub,[],options);
if exitflag < 0
    disp('没有合适整数解');
    x = x0;
    fval = fval0;
    status = exitflag;
    return;
else
    %采用分支定界法求解
    bound = inf;
    [x,fval,status] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);
end

最后测试一下

%例子1
% f = [-40 -90];
 %A = [9 7;7 20];
 %B = [56 70];
% lb = [0 0]';
%例子2
 f = [-20 -10];
 A = [5 4;2 5];
 B = [24 13];
 lb = [0 0];


[x,fval,status] = intprog(f,A,B,[1 2],[],[],lb)

得到最终的intprog()函数:

 [x,fval,status] = intprog(f,A,B,I,Aeq,Beq,lb,ub,e)
1. 其中 f为目标函数向量  
2.  A和B为不等式约束 Aeq与Beq为等式约束  
3.   I为整数约束    
4.  lb与ub分别为变量下界与上界  	
5.  5.  x为最优解,fval为最优值
6.  e为误差

以下例题加代码运用:

	例题:maximize 20 x1 + 10 x2 
 S.T.
         5 x1 + 4 x2 <=24
         2 x1 + 5 x2 <=13
          x1, x2 >=0 
          x1, x2是整数

在这里插入图片描述

3.1割平面法

3.1.1定义

  1. 如果线性规划问题无解,则整数规划无解;
  2. 如果线性规划问题的最优解为整数向量,则也是整数规划的最优解;
  3. 如果线性规划问题的解含有非整数分量,则对线性规划问题增加割平面条件:即对线性规划问题增加一个线性约束,将线性规划问题的可行区域割掉一块,使得非整数解恰好在割掉的一块中, 但又没有割掉原问题整数规划的可行解,得到一个新的问题,重 复上述的过程

3.1.2具体流程

在这里插入图片描述

3.2割平面法matlab解法

function  [intx,intf] = DividePlane(A,c,b,baseVector)
%功能:用割平面法求解整数规划
%调用格式:[intx,intf]=DividePlane(A,c,b,baseVector)
%其中,A:约束矩阵;
%      c:目标函数系数向量;
%      b:约束右端向量;
%      baseVector:初始基向量;
%      intx:目标函数取最小值时的自变量值;
%      intf:目标函数的最小值;
sz = size(A);
nVia = sz(2);%获取有多少决策变量
n = sz(1);%获取有多少约束条件
xx = 1:nVia;

if length(baseVector) ~= n
    disp('基变量的个数要与约束矩阵的行数相等!');
    mx = NaN;
    mf = NaN;
    return;
end
 
M = 0;
sigma = -[transpose(c) zeros(1,(nVia-length(c)))];
xb = b;
 
%首先用单纯形法求出最优解
while 1   
    [maxs,ind] = max(sigma);
%--------------------用单纯形法求最优解--------------------------------------
    if maxs <= 0   %当检验数均小于0时,求得最优解。      
        vr = find(c~=0 ,1,'last');
        for l=1:vr
            ele = find(baseVector == l,1);
            if(isempty(ele))
                mx(l) = 0;
            else
                mx(l)=xb(ele);
            end
        end
        if max(abs(round(mx) - mx))<1.0e-7  %判断最优解是否为整数解,如果是整数解。
            intx = mx;
            intf = mx*c;
            return;
        else  %如果最优解不是整数解时,构建切割方程
            sz = size(A);
            sr = sz(1);
            sc = sz(2);
            [max_x, index_x] = max(abs(round(mx) - mx));
            [isB, num] = find(index_x == baseVector);
            fi = xb(num) - floor(xb(num));
            for i=1:(index_x-1)
                Atmp(1,i) = A(num,i) - floor(A(num,i));
            end
            for i=(index_x+1):sc
                Atmp(1,i) = A(num,i) - floor(A(num,i));
            end
            
            Atmp(1,index_x) = 0; %构建对偶单纯形法的初始表格
            A = [A zeros(sr,1);-Atmp(1,:) 1];
            xb = [xb;-fi];
            baseVector = [baseVector sc+1];
            sigma = [sigma 0];
         
            %-------------------对偶单纯形法的迭代过程----------------------
            while 1
                %----------------------------------------------------------
                if xb >= 0    %判断如果右端向量均大于0,求得最优解
                    if max(abs(round(xb) - xb))<1.0e-7   %如果用对偶单纯形法求得了整数解,则返回最优整数解
                        vr = find(c~=0 ,1,'last');
                        for l=1:vr
                            ele = find(baseVector == l,1);
                            if(isempty(ele))
                                mx_1(l) = 0;
                            else
                                mx_1(l)=xb(ele);
                            end
                        end
                        intx = mx_1;
                        intf = mx_1*c;
                        return;
                    else   %如果对偶单纯形法求得的最优解不是整数解,继续添加切割方程
                        sz = size(A);
                        sr = sz(1);
                        sc = sz(2);
                        [max_x, index_x] = max(abs(round(mx_1) - mx_1));
                        [isB, num] = find(index_x == baseVector);
                        fi = xb(num) - floor(xb(num));
                        for i=1:(index_x-1)
                            Atmp(1,i) = A(num,i) - floor(A(num,i));
                        end
                        for i=(index_x+1):sc
                            Atmp(1,i) = A(num,i) - floor(A(num,i));
                        end
                        Atmp(1,index_x) = 0;  %下一次对偶单纯形迭代的初始表格
                        A = [A zeros(sr,1);-Atmp(1,:) 1];
                        xb = [xb;-fi];
                        baseVector = [baseVector sc+1];
                        sigma = [sigma 0];
                        continue;
                    end
                else   %如果右端向量不全大于0,则进行对偶单纯形法的换基变量过程
                    minb_1 = inf;
                    chagB_1 = inf;
                    sA = size(A);
                    [br,idb] = min(xb);
                    for j=1:sA(2)
                        if A(idb,j)<0
                            bm = sigma(j)/A(idb,j);
                            if bm<minb_1
                                minb_1 = bm;
                                chagB_1 = j;
                            end
                        end
                    end
                    sigma = sigma -A(idb,:)*minb_1;
                    xb(idb) = xb(idb)/A(idb,chagB_1);
                    A(idb,:) = A(idb,:)/A(idb,chagB_1);
                    for i =1:sA(1)
                        if i ~= idb
                            xb(i) = xb(i)-A(i,chagB_1)*xb(idb);
                            A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:);
                        end
                    end
                    baseVector(idb) = chagB_1;
                end
              %------------------------------------------------------------
            end 
            %--------------------对偶单纯形法的迭代过程---------------------    
        end     
    else     %如果检验数有不小于0的,则进行单纯形算法的迭代过程
        minb = inf;
        chagB = inf;
        for j=1:n
            if A(j,ind)>0
                bz = xb(j)/A(j,ind);
                if bz<minb
                    minb = bz;
                    chagB = j;
                end
            end
        end
        sigma = sigma -A(chagB,:)*maxs/A(chagB,ind);
        xb(chagB) = xb(chagB)/A(chagB,ind);
        A(chagB,:) = A(chagB,:)/A(chagB,ind);
        for i =1:n
            if i ~= chagB
                xb(i) = xb(i)-A(i,ind)*xb(chagB);
                A(i,:) = A(i,:) - A(i,ind)*A(chagB,:);
            end
        end
        baseVector(chagB) = ind;
    end
    M = M + 1;
    if (M == 1000000)
        disp('找不到最优解!');
        mx = NaN; 
        minf = NaN;
        return;
    end
end

其中basevector为基向量,主要由约束矩阵的秩决定,例如4*5矩阵它的秩为4,则basevector=[2,3,4,5]

  • 5
    点赞
  • 42
    收藏
    觉得还不错? 一键收藏
  • 4
    评论
数学建模动态规划代码是指包括线性规划、整数规划、动态规划等在内的各种数学建模资源的代码。这些代码可以用于解决一系列问题,如最短路线、库存管理、资源分配、设备更新、排序、装载等。动态规划方法在经济管理、生产调度、工程技术和最优控制等领域得到了广泛的应用。对于使用动态规划方法求解问题,我们可以按照以下步骤建立数学模型:首先,将问题分解为若干个阶段;其次,定义阶段状态和决策变量;然后,建立状态转移方程和目标函数;最后,通过计算和迭代求解最优解。需要注意的是,动态规划方法适用于求解以时间划分阶段的动态过程的优化问题,但也可以人为地引入时间因素,将与时间无关的静态规划问题视为多阶段决策过程来求解。在学习和应用动态规划方法时,需要正确理解基本概念和方法,同时具备丰富的想象力和创造性的技巧,以建立适合具体问题的模型并求解最优解。<span class="em">1</span><span class="em">2</span><span class="em">3</span> #### 引用[.reference_title] - *1* [数学建模各种代码](https://download.csdn.net/download/xiaohen123456/11139413)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v93^chatsearchT3_2"}}] [.reference_item style="max-width: 50%"] - *2* *3* [【数学建模数学建模学习4---动态规划(例题+matlab代码实现)](https://blog.csdn.net/qq_44528283/article/details/118975065)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v93^chatsearchT3_2"}}] [.reference_item style="max-width: 50%"] [ .reference_list ]

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值