文章目录
数学建模之线性规划
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 可转换为线性规划问题
例题：
解： 变量替换，化作线性规划模型：
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
>>
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
>>
可以看到，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
>>
继续做分支定界：
>> [x,fval] = linprog(c,A,b,[],[],[4,0]);
Optimal solution found.
>> x
x =
4.0000
1.0000
>> z = -fval
z =
14.0000
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
>>
继续分支定界：
>> [x,fval] = linprog(c,A,b,[],[],[0,3],[3,inf]);
Optimal solution found.
>> x
x =
2.5000
3.0000
>> z = -fval
z =
13.5000
>>
所以整数规划的最优解为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
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
分支定界流程：
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
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
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
>>
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
>>
匈牙利算法解决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
>>
>> [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
>>
>> 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
>>
>> [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
>>
>> 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
>>
3. 参考资料
《数学建模算法与应用》(司守奎)《数学建模算法与应用习题解答》(司守奎)- >> bilibili数学建模学习教程
779
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
