devbox
Gausst松鼠会
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

数学建模之matlab中线性规划_matlab线性规划

作者：Gausst松鼠会 | 2024-04-03 06:59:41

赞

踩

matlab线性规划

目录

一、线性规划的标准形式

二、整数规划

二、整数规划之分支定界

1.概念

2、代码实现

三、整数规划之割平面法

1、基本思想

 2、代码实现

四、整数规划之匈牙利算法（0-1）

1、适用情况

①0-1变量的使用

② 互斥问题

 ③固定费用问题

④指派问题

 2、指派问题中匈牙利法

①步骤

②举例

 3、代码实现

总结

一、线性规划的标准形式

求的是max，所以f中要加符号。f中是目标函数的系数

a表示的是不等式约束的系数（有不等号的，并且不等号为<）

b表示不等号后面的数

aeq表示线性约束的系数（有等号的系数）

beq表示等号后面的数

lb和ub表示最小值和最大值

（这里的zeros（3，1）表示创建3*1的数组，值都是0）

二、整数规划

1、概念

 

 2、一般形式

 

 

二、整数规划之分支定界

1.概念

 

 举例：

2、代码实现

matlab中代码：首先定义两个函数，然后再写一个测试脚本

  1. function [x,fval,status] = intprog(f,A,B,I,Aeq,Beq,lb,ub,e)
  2. %整数规划求解函数 intprog()
  3. %     其中 f为目标函数向量
  4. %     A和B为不等式约束 Aeq与Beq为等式约束
  5. %     I为整数约束
  6. %     lb与ub分别为变量下界与上界
  7. %     x为最优解，fval为最优值
  8. %例子:
  9. %        maximize 20 x1 + 10 x2 
  10. %        S.T.
  11. % 	             5 x1 + 4 x2 <=24
  12. %                2 x1 + 5 x2 <=13
  13. %                   x1, x2 >=0 
  14. %                   x1, x2是整数
  15. % f=[-20, -10];
  16. % A=[ 5  4; 2 5];
  17. % B=[24; 13];
  18. % lb=[0 0];
  19. % ub=[inf inf];
  20. % I=[1,2];
  21. % e=0.000001;
  22. % [x v s]= IP(f,A,B,I,[],[],lb,ub,,e)
  23. % x = 4     1  v = -90.0000   s = 1
  24.  
  25. % 控制输入参数
  26. if nargin < 9, e = 0.00001;
  27.    if nargin < 8, ub = []; 
  28.       if nargin < 7, lb = []; 
  29.          if nargin < 6, Beq = []; 
  30.             if nargin < 5, Aeq = [];
  31.                if nargin < 4, I = [1:length(f)];
  32.                end, end, end, end, end, end
  33.   
  34. %求解整数规划对应的线性规划,判断是否有解
  35. options = optimset('display','off');
  36. [x0,fval0,exitflag] = linprog(f,A,B,Aeq,Beq,lb,ub,[],options);
  37. if exitflag < 0
  38.     disp('没有合适整数解');
  39.     x = x0;
  40.     fval = fval0;
  41.     status = exitflag;
  42.     return;
  43. else
  44.     %采用分支定界法求解
  45.     bound = inf;
  46.     [x,fval,status] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);
  47. end
  48.  
  1. function [newx,newfval,status,newbound] = branchbound(f,A,B,I,x,fval,bound,Aeq,Beq,lb,ub,e)
  2.  
  3. % 分支定界法求解整数规划
  4. % f,A,B,Aeq,Beq,lb,ub与线性规划相同
  5. % I为整数限制变量的向量
  6. % x为初始解，fval为初始值
  7.  
  8. options = optimset('display','off');
  9. [x0,fval0,status0]=linprog(f,A,B,Aeq,Beq,lb,ub,[],options);
  10.  
  11. %递归中的最终退出条件
  12. %无解或者解比现有上界大则返回原解
  13. if status0 <= 0 || fval0 >= bound
  14.     newx = x;
  15.     newfval = fval;
  16.     newbound = bound;
  17.     status = status0;
  18.     return;
  19. end
  20.  
  21. %是否为整数解,如果是整数解则返回
  22. intindex = find(abs(x0(I) - round(x0(I))) > e);
  23. if isempty(intindex) %判断是否为空值
  24.     newx(I) = round(x0(I));
  25.     newfval = fval0;
  26.     newbound = fval0;
  27.     status = 1;
  28.     return;
  29. end
  30.  
  31. %当有非整可行解时，则进行分支求解
  32. %此时必定会有整数解或空解
  33. %找到第一个不满足整数要求的变量
  34. n = I(intindex(1));
  35. addA = zeros(1,length(f));
  36. addA(n) = 1;
  37. %构造第一个分支 x<=floor(x(n))
  38. A = [A;addA];
  39. B = [B,floor(x(n))];%向下取整
  40. [x1,fval1,status1,bound1] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);
  41. A(end,:) = [];
  42. B(:,end) = [];
  43. %解得第一个分支，若为更优解则替换，若不是则保持原状
  44.  
  45. status = status1;
  46. if status1 > 0 && bound1 < bound
  47.     newx = x1;
  48.     newfval = fval1;
  49.     bound = fval1;
  50.     newbound = bound1;
  51. else
  52.     newx = x0;
  53.     newfval = fval0;
  54.     newbound = bound;
  55. end
  56.  
  57. %构造第二分支
  58. A = [A;-addA];
  59. B = [B,-ceil(x(n))];%向上取整
  60. [x2,fval2,status2,bound2] = branchbound(f,A,B,I,x0,fval0,bound,Aeq,Beq,lb,ub,e);    
  61. A(end,:) = [];
  62. B(:,end) = [];
  63.  
  64. %解得第二分支，并与第一分支做比较，如果更优则替换
  65. if status2 > 0 && bound2 < bound
  66.     status = status2;
  67.     newx = x2;
  68.     newfval = fval2;
  69.     newbound = bound2;
  70. end

 测试：

  1. f=[-20 -10];
  2. A=[5 4;2 5];
  3. B=[24 13];
  4. lb=[0 0];
  5.  
  6. [x,fval,status]= intprog(f,A,B,[1,2],[],[],lb)

三、整数规划之割平面法

1、基本思想

 

实例：

引入松弛变量 x7 x8 x9 x10

 2、代码实现

 定义一个函数DIvidePlane

  1. function  [intx,intf] = DividePlane(A,c,b,baseVector)
  2. %功能：用割平面法求解整数规划
  3. %调用格式：[intx,intf]=DividePlane(A,c,b,baseVector)
  4. %其中, A：约束矩阵；
  5. %      c：目标函数系数向量；
  6. %      b：约束右端向量；
  7. %      baseVector：初始基向量；
  8. %      intx：目标函数取最值时的自变量值；
  9. %      intf：目标函数的最值；
  10. sz = size(A);
  11. nVia = sz(2);%获取有多少决策变量
  12. n = sz(1);%获取有多少约束条件
  13. xx = 1:nVia;
  14.  
  15. if length(baseVector) ~= n
  16.     disp('基变量的个数要与约束矩阵的行数相等！');
  17.     mx = NaN;
  18.     mf = NaN;
  19.     return;
  20. end
  21.  
  22. M = 0;
  23. sigma = -[transpose(c) zeros(1,(nVia-length(c)))];
  24. xb = b;
  25.  
  26. %首先用单纯形法求出最优解
  27. while 1   
  28.     [maxs,ind] = max(sigma);
  29. %--------------------用单纯形法求最优解--------------------------------------
  30.     if maxs <= 0   %当检验数均小于0时，求得最优解。      
  31.         vr = find(c~=0 ,1,'last');
  32.         for l=1:vr
  33.             ele = find(baseVector == l,1);
  34.             if(isempty(ele))
  35.                 mx(l) = 0;
  36.             else
  37.                 mx(l)=xb(ele);
  38.             end
  39.         end
  40.         if max(abs(round(mx) - mx))<1.0e-7  %判断最优解是否为整数解，如果是整数解。
  41.             intx = mx;
  42.             intf = mx*c;
  43.             return;
  44.         else  %如果最优解不是整数解时，构建切割方程
  45.             sz = size(A);
  46.             sr = sz(1);
  47.             sc = sz(2);
  48.             [max_x, index_x] = max(abs(round(mx) - mx));
  49.             [isB, num] = find(index_x == baseVector);
  50.             fi = xb(num) - floor(xb(num));
  51.             for i=1:(index_x-1)
  52.                 Atmp(1,i) = A(num,i) - floor(A(num,i));
  53.             end
  54.             for i=(index_x+1):sc
  55.                 Atmp(1,i) = A(num,i) - floor(A(num,i));
  56.             end
  57.             
  58.             Atmp(1,index_x) = 0; %构建对偶单纯形法的初始表格
  59.             A = [A zeros(sr,1);-Atmp(1,:) 1];
  60.             xb = [xb;-fi];
  61.             baseVector = [baseVector sc+1];
  62.             sigma = [sigma 0];
  63.          
  64.             %-------------------对偶单纯形法的迭代过程----------------------
  65.             while 1
  66.                 %----------------------------------------------------------
  67.                 if xb >= 0    %判断如果右端向量均大于0，求得最优解
  68.                     if max(abs(round(xb) - xb))<1.0e-7   %如果用对偶单纯形法求得了整数解，则返回最优整数解
  69.                         vr = find(c~=0 ,1,'last');
  70.                         for l=1:vr
  71.                             ele = find(baseVector == l,1);
  72.                             if(isempty(ele))
  73.                                 mx_1(l) = 0;
  74.                             else
  75.                                 mx_1(l)=xb(ele);
  76.                             end
  77.                         end
  78.                         intx = mx_1;
  79.                         intf = mx_1*c;
  80.                         return;
  81.                     else   %如果对偶单纯形法求得的最优解不是整数解，继续添加切割方程
  82.                         sz = size(A);
  83.                         sr = sz(1);
  84.                         sc = sz(2);
  85.                         [max_x, index_x] = max(abs(round(mx_1) - mx_1));
  86.                         [isB, num] = find(index_x == baseVector);
  87.                         fi = xb(num) - floor(xb(num));
  88.                         for i=1:(index_x-1)
  89.                             Atmp(1,i) = A(num,i) - floor(A(num,i));
  90.                         end
  91.                         for i=(index_x+1):sc
  92.                             Atmp(1,i) = A(num,i) - floor(A(num,i));
  93.                         end
  94.                         Atmp(1,index_x) = 0;  %下一次对偶单纯形迭代的初始表格
  95.                         A = [A zeros(sr,1);-Atmp(1,:) 1];
  96.                         xb = [xb;-fi];
  97.                         baseVector = [baseVector sc+1];
  98.                         sigma = [sigma 0];
  99.                         continue;
  100.                     end
  101.                 else   %如果右端向量不全大于0，则进行对偶单纯形法的换基变量过程
  102.                     minb_1 = inf;
  103.                     chagB_1 = inf;
  104.                     sA = size(A);
  105.                     [br,idb] = min(xb);
  106.                     for j=1:sA(2)
  107.                         if A(idb,j)<0
  108.                             bm = sigma(j)/A(idb,j);
  109.                             if bm<minb_1
  110.                                 minb_1 = bm;
  111.                                 chagB_1 = j;
  112.                             end
  113.                         end
  114.                     end
  115.                     sigma = sigma -A(idb,:)*minb_1;
  116.                     xb(idb) = xb(idb)/A(idb,chagB_1);
  117.                     A(idb,:) = A(idb,:)/A(idb,chagB_1);
  118.                     for i =1:sA(1)
  119.                         if i ~= idb
  120.                             xb(i) = xb(i)-A(i,chagB_1)*xb(idb);
  121.                             A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:);
  122.                         end
  123.                     end
  124.                     baseVector(idb) = chagB_1;
  125.                 end
  126.               %------------------------------------------------------------
  127.             end 
  128.             %--------------------对偶单纯形法的迭代过程---------------------    
  129.         end     
  130.     else     %如果检验数有不小于0的，则进行单纯形算法的迭代过程
  131.         minb = inf;
  132.         chagB = inf;
  133.         for j=1:n
  134.             if A(j,ind)>0
  135.                 bz = xb(j)/A(j,ind);
  136.                 if bz<minb
  137.                     minb = bz;
  138.                     chagB = j;
  139.                 end
  140.             end
  141.         end
  142.         sigma = sigma -A(chagB,:)*maxs/A(chagB,ind);
  143.         xb(chagB) = xb(chagB)/A(chagB,ind);
  144.         A(chagB,:) = A(chagB,:)/A(chagB,ind);
  145.         for i =1:n
  146.             if i ~= chagB
  147.                 xb(i) = xb(i)-A(i,ind)*xb(chagB);
  148.                 A(i,:) = A(i,:) - A(i,ind)*A(chagB,:);
  149.             end
  150.         end
  151.         baseVector(chagB) = ind;
  152.     end
  153.     M = M + 1;
  154.     if (M == 1000000)
  155.         disp('找不到最优解！');
  156.         mx = NaN; 
  157.         minf = NaN;
  158.         return;
  159.     end
  160. end

举例：

 
>> c = [-1;-1]; % 不要加松弛变量
A = [-1 1 1 0;3 1 0 1]; % 加上松弛变量
b = [1;4];
[x fval] = DividePlane(A,c,b,[3 4]); % 松弛变量 3 4

>> x

x =

    1.0000    1.0000

>> maxz=-fval%求得最大值，加负号

maxz =

     2

四、整数规划之匈牙利算法（0-1）

1、适用情况

①0-1变量的使用

② 互斥问题

 

 ③固定费用问题

④指派问题

 

 

 

 2、指派问题中匈牙利法

①步骤

 

 

 

②举例

 

 

 

 3、代码实现

例子1

  1. c = [3 8 2 10 3;8 7 2 9 7;6 4 2 7 5;8 4 2 3 5;9 10 6 9 10] 
  2. c = c(:); % 矩阵转换为向量
  3. a = zeros(10,25);
  4. for i = 1:5
  5. a(i,(i-1)*5+1:5*i) = 1;
  6. a(5+i,i:5:25) = 1;
  7. end % 循环将指派问题转换为线性规划问题
  8. b= ones(10,1); % 10个约束(5*2)
  9. [x y] = linprog(c,[],[],a,b,zeros(25,1),ones(25,1));
  10.  
  11. X = reshape(x,5,5)
  12. opt = y
  13.

 

 例子2：

 %% 指派问题（选择队员去进行游泳接力比赛）
clear;clc
c = [66.8 75.6 87 58.6 57.2 66 66.4 53 78 67.8 84.6 59.4 70 74.2 69.6 57.2 67.4 71 83.8 62.4]';  % 目标函数的系数矩阵（先列后行的写法）
intcon = [1:20];  % 整数变量的位置(一共20个决策变量，均为0-1整数变量)
% 线性不等式约束的系数矩阵和常数项向量（每个人只能入选四种泳姿之一，一共五个约束）
A = [1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
       0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
       0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1];
% A = zeros(5,20);
% for i = 1:5
%     A(i, (4*i-3): 4*i) = 1;
% end
b = [1;1;1;1;1];
% 线性等式约束的系数矩阵和常数项向量 （每种泳姿有且仅有一人参加，一共四个约束）
Aeq = [1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0;
          0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0;
          0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0;
          0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1];
% Aeq = [eye(4),eye(4),eye(4),eye(4),eye(4)];  % 或者写成 repmat(eye(4),1,5)  
beq = [1;1;1;1];
lb = zeros(20,1);  % 约束变量的范围下限
ub = ones(20,1);  % 约束变量的范围上限
%最后调用intlinprog()函数
[x,fval] = intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)
% reshape(x,4,5)'
%      0     0     0     1    甲自由泳
%      1     0     0     0    乙蝶泳
%      0     1     0     0    丙仰泳
%      0     0     1     0    丁蛙泳
%      0     0     0     0    戊不参加

 （感觉这个例子很实用，就引用过来了，要是有啥侵权，告诉我我删掉，不好意思）

总结

以上为线性规划中算法代码，图片来自数学建模老哥课上ppt，仅为笔记。

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/Gausst松鼠会/article/detail/354686?site
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号