数学建模之线性规划(含MATLAB代码)_matlab线性规划代码

数学建模之线性规划

1. 线性规划

1.1 matlab中的标准形式

---

>> c = [-2;-3;5]; % 等同于 c = [-2,-3,5]'; 求max取相反
>> A = [-2,5,-1;1,3,1]; % 不等式
>> b = [-10;12]; % 不等式
>> Aeq = [1,1,1]; % 等式
>> beq = 7; % 等式
>> [x,fval] = linprog(c,A,b,Aeq,beq,zeros(3,1)); % 得结果


>> x

x =

    6.4286
    0.5714
         0
         
>> z = -fval

z =

   14.5714
   
% 说明：在x1 x2 x3 = 6.4286 0.5714 0 的情况下，maxz = 14.5714
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

1.2 可转换为线性规划问题

例题：

解： 变量替换，化作线性规划模型：

c =

     1     2     3     4

>> c = [c,c]'

c =

     1
     2
     3
     4
     1
     2
     3
     4

>> Aeq = [1,-1,-1,1;1,-1,1,-3;1,-1,-2,3]

Aeq =

     1    -1    -1     1
     1    -1     1    -3
     1    -1    -2     3

>> Aeq = [Aeq,-Aeq]

Aeq =

     1    -1    -1     1    -1     1     1    -1
     1    -1     1    -3    -1     1    -1     3
     1    -1    -2     3    -1     1     2    -3
     
>> beq = [0;1;-1/2]

beq =

         0
    1.0000
   -0.5000
   
>> [uv,fval] = linprog(c,[],[],Aeq,beq,zeros(8,1));

Optimal solution found.

>> x = uv(1:4) - uv(5:end)

x =

    0.2500
         0
         0
   -0.2500

>> minz = fval

minz =

    1.2500

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

2. 整数规划

整数规划的整数仅仅针对于决策变量而言，与最优解是否为整数无关。

2.1 分支定界算法

通常，把全部可行解空间反复地分割为越来越小的子
集，称为分枝；并且对每个子集内的解集计算一个目标下界（对于最小值问题），这称
为定界。在每次分枝后，凡是界限超出已知可行解集目标值的那些子集不再进一步分枝，这样，许多子集可不予考虑，这称剪枝。这就是分枝定界法的主要思路。

分枝定界法可用于解纯整数或混合的整数规划问题。

2.1.1 分支定界举例

>> c = [-3;-2];
>> A = [2,3;2,1];
>> b = [14;9];
>> [x,fval] = linprog(c,A,b,[],[],zeros(2,1));

Optimal solution found.

>> x

x =

    3.2500
    2.5000

>> z = - fval

z =

   14.7500

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

可以看到，x并不是整数，下面做分支定界：

>> [x,fval] = linprog(c,A,b,[],[],zeros(2,1),[3,inf]);

Optimal solution found.

>> x

x =

    3.0000
    2.6667

>> z = -fval

z =

   14.3333

>>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

继续做分支定界：

​

>> [x,fval] = linprog(c,A,b,[],[],[4,0]);

Optimal solution found.

>> x

x =

    4.0000
    1.0000

>> z = -fval

z =

   14.0000
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

14.3333比14.0000大，所以在x1≤3下再做分支定界：

>> [x,fval] = linprog(c,A,b,[],[],zeros(2,1),[3,2]);

Optimal solution found.

>> x

x =

     3
     2

>> z = -fval

z =

    13

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

继续分支定界：

>> [x,fval] = linprog(c,A,b,[],[],[0,3],[3,inf]);

Optimal solution found.

>> x

x =

    2.5000
    3.0000

>> z = -fval

z =

   13.5000

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

所以整数规划的最优解为14。

2.1.2 matlab代码实现

branchbound.m

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

intprog.m

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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

---

分支定界流程：

2.1.3 intlinprog函数求解整数规划

>> c = [-40;-90];
>> A = [9,7;7,20];
>> b = [56;70];
>> [x,fval] = intlinprog(c,[1 2],A,b,[],[],zeros(2,1));
>> x

x =

    4.0000
    2.0000

>> -fval

ans =

   340
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

x = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub)：函数相比linprog多了一个参数intcon，用于标定整数变量的位置：x1、x2为整数，即intcon = [1 2]。

2.2 割平面算法

解：两个不等式，这里引入两个松弛变量x3和x4；

2.2.1 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160

2.2.2 割平面算法应用

>> 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

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

2.3 匈牙利算法(0-1规划)

0-1规划问题也可以看做区间在[0,1]的整数规划，下面利用intlinprog函数进行计算：

>> c = [-3 2 -5]';
>> A = [1 2 -1;1 4 1];
>> b = [2;4];
>> c = [-3 2 -5]';
>> A = [1 2 -1;1 4 1;1 1 0;0 4 1];
>> b = [2;4;3;6];
>> [x fval] = intlinprog(c,[1 2 3],A,b,[],[],zeros(3,1),ones(3,1));
>> x

x =

     1
     0
     1

>> maxz = -fval

maxz =

     8

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

---

匈牙利算法解决0-1规划问题：

2.3.1 投资问题

2.3.2 互斥约束条件问题

2.3.3 固定费用问题

2.3.4 指派问题

举例

2.3.5 非标准形式的指派问题

2.3.6 匈牙利算法说明

每一行找出最小值让该行各个元素减去最小值，每一列找出最小值让该列各个元素减去最小值，从只有一个0元素的行开始，将0元素圈圈，划掉该0元素所在列的其他0元素，接着再找只有一个0元素的行，重复操作，当行中0元素没有一个的时候，找只有一个0元素的列，圈圈，划掉该0元素所在行的其他0元素。

最后，若圈出的0元素的个数小于阶数n，进行如下操作：(否则找到解)

对没有圈0元素的行打√，对打√号的行中的所有划掉0元素的所在列打√，对打√的列中含圈0元素的所在行打√，重复上述3步，直到得不出新的打√号的行列为止。

下面，对没有打√的行画横线，对打√的列画纵线，从而得到覆盖所有0元素的最少横线数，应当等于圈0元素的个数，然后进行如下操作，变换当前矩阵来增加0元素，以便找到n个独立0元素：

找到没有被直线覆盖的所有元素中的最小元素值，将打√的各行减去最小元素值，打√的各列加上最小元素值，从而得到一个新的指派矩阵，用这个新的指派矩阵重复上述所有步骤，直到找出n个独立0元素结果为止。

例题

实现步骤

2.3.7 匈牙利算法代码实现

>> 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] 

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

>> c = c(:); % 矩阵转换为向量
>> a = zeros(10,25);
>> for i = 1:5
a(i,(i-1)*5+1:5*i) = 1;
a(5+i,i:5:25) = 1;
end % 循环将指派问题转换为线性规划问题
>> b = ones(10,1); % 10个约束(5*2)
>> [x y] = linprog(c,[],[],a,b,zeros(25,1),ones(25,1));

Optimal solution found.

>> X = reshape(x,5,5)

X =

     0     0     0     0     1
     0     0     1     0     0
     0     1     0     0     0
     0     0     0     1     0
     1     0     0     0     0

>> opt = y

opt =

    21

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

>> [x y] = linprog(-c,[],[],a,b,zeros(25,1),ones(25,1));

Optimal solution found.

>> X = reshape(x,5,5)

X =

     0     1     0     0     0
     0     0     0     1     0
     0     0     1     0     0
     1     0     0     0     0
     0     0     0     0     1

>> opt = -y

opt =

    37

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

>> c = [6 7 11 2;4 5 9 8;3 1 10 4;5 9 8 2];
>> a = zeros(8,16);
>>  for i = 1:4
a(i,(i-1)*4+1:4*i) = 1;
a(4+i,i:4:16) = 1;
end
>> b = ones(8,1);
>>  [x y] = linprog(c,[],[],a,b,zeros(16,1),ones(16,1));

Optimal solution found.

>> X = reshape(x,4,4)

X =

     0     0     0     1
     1     0     0     0
     0     1     0     0
     0     0     1     0

>> opt = y

opt =

    15

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

>> [x y]= linprog(-c,[],[],a,b,zeros(25,1),ones(25,1));

Optimal solution found.

>> X = reshape(x,5,5)

X =

     0     1     0     0     0
     0     0     0     1     0
     0     0     1     0     0
     1     0     0     0     0
     0     0     0     0     1

>> opt = -y

opt =

    37

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

>> c = [6 7 11 2;4 5 9 8;3 1 10 4;5 9 8 2];
>> a = zeros(8,16);
>>  for i = 1:4
a(i,(i-1)*4+1:4*i) = 1;
a(4+i,i:4:16) = 1;
end
>> b = ones(8,1);
>>  [x y] = linprog(c,[],[],a,b,zeros(16,1),ones(16,1));

Optimal solution found.

>> X = reshape(x,4,4)

X =

     0     0     0     1
     1     0     0     0
     0     1     0     0
     0     0     1     0

>> opt = y

opt =

    15

>> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

---

3. 参考资料

• 《数学建模算法与应用》(司守奎)
• 《数学建模算法与应用习题解答》(司守奎)
• >> bilibili数学建模学习教程

---