目录
一、线性规划求解
二、非线性规划问题
三、整数规划(包括0-1规划)
四、最大最小化模型
五、多目标规划模型
注意:代码文件仅供参考,一定不要直接用于自己的数模论文中
国赛对于论文的查重要求非常严格,代码雷同也算作抄袭
如何修改代码避免查重的方法:https://www.bilibili.com/video/av59423231 //清风数学建模
在给定条件下,按照某一衡量指标(目标函数)求计划。
一、线性规划求解
根据题弄A x b等矩阵,列式子用matlab求解
[x fval] = linprog(c, A, b, Aeq, beq, lb)
fval = -fval % 注意这个fval要取负号(原来是求最大值,我们添加负号变成了最小值问题)
例题1:生产决策问题
%% 生产决策问题
format long g %可以将Matlab的计算结果显示为一般的长数字格式(默认会保留四位小数,或使用科学计数法)
% (1) 系数向量
c = zeros(9,1); % 初始化目标函数的系数向量全为0
c(1) = 1.25 -0.25 -300/6000*5; % x1前面的系数是c1
c(2) = 1.25 -0.25 -321/10000*7;
c(3) = -250 / 4000 * 6;
c(4) = -783/7000*4;
c(5) = -200/4000 * 7;
c(6) = -300/6000*10;
c(7) = -321 / 10000 * 9;
c(8) = 2-0.35-250/4000*8;
c(9) = 2.8-0.5-321/10000*12-783/7000*11;
c = -c; % 我们求的是最大值,所以这里需要改变符号
% (2) 不等式约束
A = zeros(5,9);
A(1,1) = 5; A(1,6) = 10;
A(2,2) = 7; A(2,7) = 9; A(2,9) = 12;
A(3,3) = 6; A(3,8) = 8;
A(4,4) = 4; A(4,9) = 11;
A(5,5) = 7;
b = [6000 10000 4000 7000 4000]';
% (3) 等式约束
Aeq = [1 1 -1 -1 -1 0 0 0 0;
0 0 0 0 0 1 1 -1 0];
beq = [0 0]';
%(4)上下界
lb = zeros(9,1);
% 进行求解
[x fval] = linprog(c, A, b, Aeq, beq, lb)
fval = -fval
% fval =
% 1146.56650246305
% 注意,本题应该是一个整数规划的例子,我们在后面的整数规划部分再来重新求解。
intcon = 1:9;
[x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb)
fval = -fval
例题二:投料问题
%% 投料问题
clear,clc
format long g %可以将Matlab的计算结果显示为一般的长数字格式(默认会保留四位小数,或使用科学计数法)
% (1) 系数向量
a=[1.25 8.75 0.5 5.75 3 7.25]; % 工地的横坐标
b=[1.25 0.75 4.75 5 6.5 7.25]; % 工地的纵坐标
x = [5 2]; % 料场的横坐标
y = [1 7]; % 料场的纵坐标
c = []; % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量)
for j =1:2
for i = 1:6
c = [c; sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)]; % 每循环一次就在c的末尾插入新的元素
end
end
% (2) 不等式约束
A =zeros(2,12);
A(1,1:6) = 1;
A(2,7:12) = 1;
b = [20,20]';
% (3) 等式约束
Aeq = zeros(6,12);
for i = 1:6
Aeq(i,i) = 1; Aeq(i,i+6) = 1;
end
% Aeq = [eye(6),eye(6)] % 两个单位矩阵横着拼起来
beq = [3 5 4 7 6 11]'; % 每个工地的日需求量
%(4)上下界
lb = zeros(12,1);
% 进行求解
[x fval] = linprog(c, A, b, Aeq, beq, lb)
x = reshape(x,6,2) % 将x变为6行2列便于观察(reshape函数是按照列的顺序进行转换的,也就是第一列读完,读第二列,即x1对应x_1,1,x2对应x_2,1)
二、非线性规划问题
非线性代表一般有如x^2这种非一元线性。这种题也算是数学规划中最难的一类,一般代码放一个文件夹里,有许多函数需要写。
如果函数出现函数或变量 ‘XXX’ 无法识别的问题:
例题1选址问题(上面的投料问题第二问)
%% 使用蒙特卡罗的方法来找初始值(推荐)
clc,clear;
n=10000000; %生成的随机数组数
x1=unifrnd(-100,100,n,1); % 生成在[-100,100]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(-100,100,n,1); % 生成在[-100,100]之间均匀分布的随机数组成的n行1列的向量构成x2
fmin=+inf; % 初始化函数f的最小值为正无穷(后续只要找到一个比它小的我们就对其更新)
for i=1:n
x = [x1(i), x2(i)]; %构造x向量, 这里千万别写成了:x =[x1, x2]
if ((x(1)-1)^2-x(2)<=0) & (-2*x(1)+3*x(2)-6 <= 0) % 判断是否满足条件
result = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ; % 如果满足条件就计算函数值
if result < fmin % 如果这个函数值小于我们之前计算出来的最小值
fmin = result; % 那么就更新这个函数值为新的最小值
x0 = x; % 并且将此时的x1 x2更新为初始值
end
end
end
disp('蒙特卡罗选取的初始值为:'); disp(x0)
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)
fval = -fval
function f = fun1(x)
% 注意:这里的f实际上就是目标函数,函数的返回值也是f
% 输入值x实际上就是决策变量,由x1和x2组成的向量
% fun1是函数名称,到时候会被fmincon函数调用, 可以任意取名
% 保存的m文件和函数名称得一致,也要为fun1.m
% max f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2
f = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ;
end
function [c,ceq] = nonlfun1(x)
% 注意:这里的c实际上就是非线性不等式约束,ceq实际上就是非线性等式约束
% 输入值x实际上就是决策变量,由x1和x2组成的一个向量
% 返回值有两个,一个是非线性不等式约束c,一个是非线性等式约束ceq
% nonlfun1是函数名称,到时候会被fmincon函数调用, 可以任意取名,但不能和目标函数fun1重名
% 保存的m文件和函数名称得一致,也要为nonlfun1.m
% -(x1-1)^2 +x2 >= 0
c = [(x(1)-1)^2-x(2)]; % 千万別写成了: (x1-1)^2 -x2
ceq = []; % 不存在非线性等式约束,所以用[]表示
end
飞行管理问题:非常难,不再粘贴代码了,看更新12
三、整数规划(包括0-1规划)
例题像动态规划问题,如背包问题、指派问题(游泳接力)、钢管切割问题。这里就放一个背包问题。
%% 背包问题(货车运送货物的问题)
c = -[540 200 180 350 60 150 280 450 320 120]; % 目标函数的系数矩阵(最大化问题记得加负号)
intcon=[1:10]; % 整数变量的位置(一共10个决策变量,均为0-1整数变量)
A = [6 3 4 5 1 2 3 5 4 2]; b = 30; % 线性不等式约束的系数矩阵和常数项向量(物品的重量不能超过30)
Aeq = []; beq =[]; % 不存在线性等式约束
lb = zeros(10,1); % 约束变量的范围下限
ub = ones(10,1); % 约束变量的范围上限
%最后调用intlinprog()函数
[x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)
fval = -fval
四、最大最小化模型
这类问题就是在最糟糕的条件下选择最优解。
%% 最大最小化模型 : min{max[f1,f2,···,fm]}
x0 = [6, 6]; % 给定初始值
lb = [3, 4]; % 决策变量的下界
ub = [8, 10]; % 决策变量的上界
[x,feval] = fminimax(@Fun,x0,[],[],[],[],lb,ub)
max(feval)
% x =
% 8.0000 8.5000
% feval =
% 13.5000 5.5000 5.5000 12.5000 8.5000 8.5000 5.5000 13.5000 9.5000 0.5000
% 结论:
% 在坐标为(8,8.5)处建立供应中心可以使该点到各需求点的最大距离最小,最小的最大距离为13.5单位。
function f = Fun(x)
a=[1 4 3 5 9 12 6 20 17 8];
b=[2 10 8 18 1 4 5 10 8 9];
% 函数向量
f=zeros(10,1);
for i = 1:10
f(i) = abs(x(1)-a(i))+abs(x(2)-b(i));
end
% f(1) = abs(x(1)-a(1))+abs(x(2)-b(1));
% f(2) = abs(x(1)-a(2))+abs(x(2)-b(2));
% f(3) = abs(x(1)-a(3))+abs(x(2)-b(3));
% f(4) = abs(x(1)-a(4))+abs(x(2)-b(4));
% f(5) = abs(x(1)-a(5))+abs(x(2)-b(5));
% f(6) = abs(x(1)-a(6))+abs(x(2)-b(6));
% f(7) = abs(x(1)-a(7))+abs(x(2)-b(7));
% f(8) = abs(x(1)-a(8))+abs(x(2)-b(8));
% f(9) = abs(x(1)-a(9))+abs(x(2)-b(9));
% f(10) = abs(x(1)-a(10))+abs(x(2)-b(10));
end
五、多目标规划模型
多个目标处理,如企业在保证利润最大化时候要保证污染最小,但比方不能无脑利润最大化,因此需要引入权重(例如层次分析法)
%% 多目标规划问题
w1 = 0.4; w2 = 0.6; % 两个目标函数的权重 x1 = 5 x2 = 2
w1 = 0.5; w2 = 0.5; % 两个目标函数的权重 x1 = 5 x2 = 2
w1 = 0.3; w2 = 0.7; % 两个目标函数的权重 x1 = 1 x2 = 6
c = [w1/30*2+w2/2*0.4 ;w1/30*5+w2/2*0.3]; % 线性规划目标函数的系数
A = [-1 -1]; b = -7; % 不等式约束
lb = [0 0]'; ub = [5 6]'; % 上下界
[x,fval] = linprog(c,A,b,[],[],lb,ub)
f1 = 2*x(1)+5*x(2)
f2 = 0.4*x(1) + 0.3*x(2)
%% 敏感性分析
clear;clc
W1 = 0.1:0.001:0.5; W2 = 1- W1;
n =length(W1);
F1 = zeros(n,1); F2 = zeros(n,1); X1 = zeros(n,1); X2 = zeros(n,1); FVAL = zeros(n,1);
A = [-1 -1]; b = -7; % 不等式约束
lb = [0 0]; ub = [5 6]; % 上下界
for i = 1:n
w1 = W1(i); w2 = W2(i);
c = [w1/30*2+w2/2*0.4 ;w1/30*5+w2/2*0.3]; % 线性规划目标函数的系数
[x,fval] = linprog(c,A,b,[],[],lb,ub);
F1(i) = 2*x(1)+5*x(2);
F2(i) = 0.4*x(1) + 0.3*x(2);
X1(i) = x(1);
X2(i) = x(2);
FVAL(i) = fval;
end
% 「Matlab」“LaTex字符汇总”讲解:https://blog.csdn.net/Robot_Starscream/article/details/89386748
% 在图上可以加上数据游标,按住Alt加鼠标左键可以设置多个数据游标出来。
figure(1)
plot(W1,F1,W1,F2)
xlabel('f_{1}的权重')
ylabel('f_{1}和f_{2}的取值')
legend('f_{1}','f_{2}')
figure(2)
plot(W1,X1,W1,X2)
xlabel('f_{1}的权重')
ylabel('x_{1}和x_{2}的取值')
legend('x_{1}','x_{2}')
figure(3)
plot(W1,FVAL) % 看起来是两个直线组合起来的下半部分
xlabel('f_{1}的权重')
ylabel('综合指标的值')