来源：数学建模清风学习内容整理
一、概述
（1）什么是数学规划？
- 求目标函数在一定约束条件下的极值问题
数学规划是运筹学的一个分支，用来研究：在给定的条件下（约束条件），如何按照某一衡量指标（目标函数）来寻求计划、管理工作中的最优方案。
（2）数学规划的一般形式
（3）数学规划的分类
① 线性规划（Linear programming）
如果目标函数
f
(
x
)
f(x)
f(x)和约束条件均是决策变量的线性表达式，那么此时的数学规划问题就属于线性规划。
② 非线性规划（nonlinear programming）
当目标函数
f
(
x
)
f(x)
f(x)或者约束条件中有一个是决策变量
x
x
x的非线性表达式，那么此时的数学规划问题就属于非线性规划。
③ 整数规划（integer programming）
整数规划是一类要求变量取整数值的数学规划（线性整数规划和非线性整数规划）。
④ 0-1规划（0-1 programming）
整数规划的特例，整数变量的取值只能为0和1。
二、线性规划问题的求解
（1）Matlab中线性规划的标准型
% 1题代码
[x,fval]=linprog(c,A,b,[],[],lb) %不写ub意味着不存在上界约束
% 2题代码
[x,fval]=linprog(c,A,b,Aeq,beg,lb)
% 3题代码
[x,fval]=linprog(c,A,b,Aeq,beg,lb)
fval=-fval
（2）Matlab求解线性规划的函数
[x,fval]=linprog(c,A,b,Aeq,beq,lb,ub,xo)
① xo表示给定Matlab迭代求解的初始值（一般不用给）；
② c，A，b，Aeq，beq，lb，ub的意义和标准型中的意义一致 ；
③ 若不存在不等式约束，可用 [ ] 替代 A 和 b ；
④ 若不存在等式约束，可用 [ ] 替代 Aeq 和 beq ；
⑤ 若某个
x
i
x_i
xi无下届或上届，则设置 lb(i)=-inf，ub(i)=+inf ；
⑥ 返回的
x
x
x表示最小值处的
x
x
x取值；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
（2）例题二：投料问题
%% 投料问题
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）
% fval =
% 135.281541790676
四、非线性规划问题的求解
（1）Matlab中非线性规划的标准型
（2）Matlab中求解非线性规划的命令
[x,fval]=fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)
（1） 非线性规划中对于初始值
x
0
x_0
x0 的选取非常重要，因为非线性规划的算法求解出来的是一个局部最优解。
（线性规划不存在这个问题）
（2） 如果要求“全局最优解”，有两种思路：
① 给定不同的初始值，在里面找到最优解；
② 先用蒙特卡罗模拟，得到一个蒙特卡罗解，再将这个解作为初始值来求最优解；
（3）“option” 选项可以给定求解的算法，一共有四种：
interior-point（内点法）
sqp（序列二次规划法）
active-set（有效集法）
trust-region-reflective（信赖域反射算法）
（4） 不同的算法有其各自的优缺点和适用情况，我们可以改变求解的算法来看求解的结果是否变好；
（5 ）“ @fun ” 表示目标函数，要编写一个独立的“.m”文件储存目标函数：
funciton f=fun(x) % 注：1.fun可以任意取名；
f= ... % 2.f也可以任意取名，但返回的f和函数内部的f得完全一致；
end % 3.这里的x实际上是表示决策变量的向量，其行列方向取决于初始值x；
% 4.调用函数：fmincon(@fun, ... ) 求解
（6）“@nonlfun”表示非线性部分的约束，同样得编写一个独立的 “.m” 文件储存非线性约束条件：
function[c,ceq]=nonlfun(x) % 1. nonlfun同样可任意取名，不和上面的fun相同就可，保存的.m文件也得是这个名字；
c=[非线性不等式约束1； % 2. c、ceq中可能有多个约束，因此写成列向量的形式；
非线性不等式约束p；] % 3. 若不存在非线性不等式约束，则可令 c=[];
ceq=[非线性等式约束1； % 4. 调用函数：fmincon（...,@nonlfun,options）求解
非线性等式约束q；
end
（7） 注意要把下标改写为括号，例如：
f
=
x
1
2
+
3
x
2
f=x_1^2+3x_2
f=x12+3x2 写成matlab能识别的形式：f=x(1)^2+3*x(2)；
（8） 若不存在某种约束，则可用 “ [ ] ” 替代，若后面全为 “ [ ] ” 且不指定option（使用默认的求解方法），则 “ [ ] ” 也可以省略掉；
对于例题的matlab使用讲解
- 主函数
%% 非线性规划的函数
% [x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)
% x0表示给定的初始值（用行向量或者列向量表示），必须得写
% A b表示线性不等式约束
% Aeq beq 表示线性等式约束
% lb ub 表示上下界约束
% @fun表示目标函数
% @nonlfun表示非线性约束的函数
% option 表示求解非线性规划使用的方法
clear;clc
format long g %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
%% 例题1的求解
% max f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2
% s.t. -(x1-1)^2 +x2 >= 0 ; 2x1-3x2+6 >= 0
x0 = [0 0]; %任意给定一个初始值
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1) % 注意 fun1.m文件和nonlfun1.m文件都必须在当前文件夹目录下
fval = -fval
% 一个值得讨论的地方，能不能把线性不等式约束Ax <= b也写到nonlfun1函数中？
% 先把nonlfun1中的c改为下面这样：
% c = [(x(1)-1)^2-x(2);
% -2*x(1)+3*x(2)-6];
% [x,fval] = fmincon(@fun1,x0,[],[],[],[],[],[],@nonlfun1)
% 结果也是可以计算出来的，但并不推荐这样做~
%% 使用其他算法对例题1求解
% edit fmincon % 查看fmincon的“源代码”
% Matlab2017a默认使用的算法是'interior-point' 内点法
% 使用interior point算法 （内点法）
option = optimoptions('fmincon','Algorithm','interior-point')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)
fval = -fval
% 使用SQP算法 （序列二次规划法）
option = optimoptions('fmincon','Algorithm','sqp')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)
fval = -fval %得到-4.358，远远大于内点法得到的-1,猜想是初始值的影响
% 改变初始值试试
x0 = [1 1]; %任意给定一个初始值
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option) % 最小值为-1，和内点法相同（这说明内点法的适应性要好）
fval = -fval
% 使用active set算法 （有效集法）
option = optimoptions('fmincon','Algorithm','active-set')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)
fval = -fval
% 使用trust region reflective (信赖域反射算法)
option = optimoptions('fmincon','Algorithm','trust-region-reflective')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)
fval = -fval
% this algorithm does not solve problems with the constraints you have specified.
% 这说明这个算法不适用我们这个约束条件，所以以后遇到了不能求解的情况，记得更换其他算法试试！！！
%% 选取初始值得到的结果可能会不满足限定条件，出现了一个Bug 因此选择的初始值很重要
x0 = [40.8, 10.8];
option = optimoptions('fmincon','Algorithm','interior-point')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)
fval = -fval
% https://cn.mathworks.com/help/optim/ug/fmincon.html
%% 生成不同的随机初始值来优化代码，有一定几率会触发上面那个Bug，因此不推荐
n = 10; % 重复n次
Fval = +inf; X = [0,0]; %初始化最优的结果
A = [-2 3]; b = 6;
for i = 1:n
x0 = [rand()*10 , rand()*10]; %用随机数生成一个初始值（随机数的范围自己根据题目条件设置）
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option); % 注意 fun1.m文件和nonlfun1.m文件都必须在当前文件夹目录下
if fval < Fval % 如果找到了更小的值，那么就代替最优的结果
Fval = fval;
X = x;
end
end
Fval = -Fval
X
%% 使用蒙特卡罗的方法来找初始值(推荐）
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
%% 例题二的求解
x0 = [1 1 1]; %任意给定一个初始值
lb = [0 0 0]; % 决策变量的下界
[x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2) % 注意 fun2.m文件和nonfun2.m文件都必须在当前文件夹目录下
% x =
% 0.552167405729277 1.20325915507969 0.947824046150443
% fval =
% 10.6510918606939
%% 使用蒙特卡罗的方法来找初始值(推荐）
clc,clear;
n=1000000; %生成的随机数组数
x1= unifrnd(0,2,n,1); % 生成在[0,2]之间均匀分布的随机数组成的n行1列的向量构成x1
x2 = sqrt(2-x1); % 根据非线性等式约束用x1计算出x2
x3 = sqrt((3-x2)/2); % 根据非线性等式约束用x2计算出x3
fmin=+inf; % 初始化函数f的最小值为正无穷（后续只要找到一个比它小的我们就对其更新）
for i=1:n
x = [x1(i), x2(i), x3(i)]; %构造x向量, 这里千万别写成了：x =[x1, x2, x3]
if (-x(1)^2+x(2)-x(3)^2<=0) & (x(1)+x(2)^2+x(3)^2-20<=0) % 判断是否满足条件
result =sum(x.*x) + 8 ; % 如果满足条件就计算函数值
if result < fmin % 如果这个函数值小于我们之前计算出来的最小值
fmin = result; % 那么就更新这个函数值为新的最小值
x0 = x; % 并且将此时的x1 x2 x3更新为初始值
end
end
end
disp('蒙特卡罗选取的初始值为：'); disp(x0)
lb = [0 0 0]; % 决策变量的下界
[x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2) % 注意 fun2.m文件和nonfun2.m文件都必须在当前文件夹目录下
%% 例题三的求解(蒙特卡罗模拟那一讲的例题)
clear;clc
% 蒙特卡罗模拟得到的最大值为3445.6014
% 最大值处x1 x2 x3的取值为：
% 22.5823101903968 12.5823101903968 12.1265223966757
A = [1 -2 -2; 1 2 2]; b = [0 72];
x0 = [ 22.58 12.58 12.13];
Aeq = [1 -1 0]; beq = 10;
lb = [-inf 10 -inf]; ub = [inf 20 inf];
[x,fval] = fmincon(@fun3,x0,A,b,Aeq,beq,lb,ub,[]) % 注意没有非线性约束，所以这里可以用[]替代，或者干脆不写
fval = -fval
% % 注意：代码文件仅供参考，一定不要直接用于自己的数模论文中
% % 国赛对于论文的查重要求非常严格，代码雷同也算作抄袭
- 被调用的函数：fun1 （例题一的目标函数）
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
- 被调用的函数：fun2 （例题二的目标函数）
function f = fun2(x)
% f = x(1)^2+x(2)^2 +x(3)^2+8 ;
f = sum(x.*x) + 8; % 可别忘了x实际上是一个向量，我们可以使用矩阵的运算符号对其计算
end
- 被调用的函数：fun3 （例题三的目标函数）
function f = fun3(x)
f = -prod(x); % 可别忘了x实际上是一个向量（prod表示连乘符号，用法和sum类似）
end
- 被调用的函数： nonlfun1 （例题一的非线性约束条件）
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
- 被调用的函数：nonlfun2 （例题二的非线性约束条件）
function [c,ceq] = nonlfun2(x)
% 非线性不等式约束
c = [-x(1)^2+x(2)-x(3)^2; % 一定要注意写法的规范，再次强调这里的x是一个向量！不能把x(1)写成x1
x(1)+x(2)^2+x(3)^2-20];
% 非线性等式约束
ceq = [-x(1)-x(2)^2+2;
x(2)+2*x(3)^2-3];
end
五、非线性规划的典型例题
（1）例题一：选址问题
- 主函数
%% 选址问题
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,16); % 注意这里要改成16
A(1,1:6) = 1;
A(2,7:12) = 1;
b = [20,20]';
% (3) 等式约束
Aeq = zeros(6,16); % 注意这里要改成16
for i = 1:6
Aeq(i,i) = 1; Aeq(i,i+6) = 1;
end
beq = [3 5 4 7 6 11]'; % 每个工地的日需求量
%（4）上下界
lb = zeros(16,1);
% lb = [zeros(12,1); -inf*ones(4,1)]; 两个新料场坐标的下界可以设为-inf
% 进行求解
% 注意哦，这里我们只尝试了这一个初始值，大家可以试试其他的初始值，有可能能够找到更好的解。
% 还可以用遗传算法求解
x0 = [3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7]; % 用第一问的结果作为初始值
[x,fval] = fmincon(@fun5,x0,A,b,Aeq,beq,lb) % 注意没有非线性约束，所以这里可以用[]替代，或者干脆不写
reshape(x(1:12),6,2) % 将x的前12个元素变为6行2列便于观察（reshape函数是按照列的顺序进行转换的，也就是第一列读完，读第二列，即x1对应x_1,1，x2对应x_2,1）
% 新坐标（5.74，4.99） （7.25，7.25）
% fval =
% 89.9231692432933
% 第一问的fval =
% 135.281541790676
135.281541790676 - 89.9231692432933 % 45.3583725473827
- 被调用的函数
function f = fun5(xx) % 注意为了避免和下面的x同号，我们把决策变量的向量符号用xx表示（注意xx的长度为16）
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 = [xx(13) xx(15)]; % 新料场的横坐标
y = [xx(14) xx(16)]; % 新料场的纵坐标
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
% 下面我们要求吨千米数，注意c是列向量，我们计算非线性规划时给定的初始值x0是行向量
f = xx(1:12) * c;
end
（2）例题二：飞行管理问题
计算新坐标的示意图
- 主函数
%% 飞行管理问题
format long g
%% （1）画六架飞机的位置
clear;clc
figure(1) % 生成一个图形
box on % 显示为封闭的盒子
% 绘制飞机的初始位置
data = [150 140 243;
85 85 236;
150 155 220.5;
145 50 159;
130 150 230;
0 0 52];
plot(data(:,1),data(:,2),'.r')
axis([0 160,0,160]);% 设置坐标轴刻度范围
hold on;
% 在图上标上注释
for i = 1:6
txt = ['飞机',num2str(i)];
text(data(i,1)+2,data(i,2)+2,txt,'FontSize',8)
end
% 把Matlab做出来的图可以导出，然后再放到PPT中画出飞机飞行方向的箭头
%% 求解非线性规划问题
x0 = [0 0 0 0 0 0]; % 初始值
lb = -pi/6*ones(6,1);
ub = pi/6*ones(6,1);
[x,fval] = fmincon(@fun6,x0,[],[],[],[],lb,ub,@nonlfun6)
x = x * 180 / pi % 将弧度转换为度数
% 定义一：fval = 3.7315°
% 定义二: fval = 6.9547((°)^2)
- 被调用函数：fun6
function f = fun6(delta) % 决策变量delta为六架飞机调整的角度
% f =sum(abs(delta)) * 180 / pi; % 目标函数第一种定义：绝对值的和（将弧度转换为度数）
f = sum(delta .* delta) * (180 / pi)^2; % 目标函数第二种定义：平方和（将弧度转换为度数）
end
- 被调用函数：nonlfun6
function [c,ceq] = nonlfun6(delta) % 决策变量delta为六架飞机调整的角度
x = [150 85 150 145 130 0]; % 飞机初始位置的横坐标
y = [140 85 155 50 150 0]; % 飞机初始位置的纵坐标
theta = [243 236 220.5 159 230 52] * pi / 180; % 飞机初始的飞行方向角
v = 800; % 飞机速度
co = cos(theta + delta); % 包含6个元素的向量
si = sin(theta + delta); % 包含6个元素的向量
% 下面开始计算飞机i和j之间的最短距离（只需要计算矩阵的一半即可）
d = zeros(6); % 初始化飞机两两之间的最短距离矩阵
for i = 2: 6
for j = 1: i-1
% 套用我们推导出来的公式计算飞机i和飞机j相距最近的时间
fenzi = ((y(j)-y(i))*(si(j)-si(i)) +(x(j)-x(i))*(co(j)-co(i))) ; % 分子
fenmu = v * ((co(j)-co(i))^2 + (si(j)-si(i))^2); % 分母
t(i,j) =- fenzi / fenmu;
if t(i,j) <0
d(i, j) = 1000; % 此时最初的位置就是相距最近的点，因为最初的时候所有飞机两两之间的距离就大于8，因此未来绝不会相撞，我们令它们的距离为一个特别大的数
else
d(i, j) = sqrt((x(j)-x(i)+v*t(i,j)*(co(j)-co(i)))^2+(y(j)-y(i)+v*t(i,j)*(si(j)-si(i)))^2);
end
end
end
% 非线性不等式约束
c =ones(15,1)*8.000001 - [d(2,1); d(3,1:2)'; d(4,1:3)'; d(5,1:4)'; d(6,1:5)'];
% 12个非线性不等式约束: “最短距离>8” 等价于 “8 - 最短距离<0”
% 注意： 由于Matlab标准型中取的是小于等于号，因此这里取一个比8略大的数：8.000001-最短距离<=0
ceq = []; % 没有非线性等式约束
end
六、整数规划
（1）Matlab中线性整数规划求解
[x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)
%% 线性整数规划问题
%% 例1
c=[-20,-10]';
intcon=[1,2]; % x1和x2限定为整数
A=[5,4;
2,5];
b=[24;13];
lb=zeros(2,1);
[x,fval]=intlinprog(c,intcon,A,b,[],[],lb)
fval = -fval
%% 例2
c=[18,23,5]';
intcon=3; % x3限定为整数
A=[107,500,0;
72,121,65;
-107,-500,0;
-72,-121,-65];
b=[50000;2250;-500;-2000];
lb=zeros(3,1);
[x,fval]=intlinprog(c,intcon,A,b,[],[],lb)
%% 例3
c=[-3;-2;-1]; intcon=3; % x3限定为整数
A=ones(1,3); b=7;
Aeq=[4 2 1]; beq=12;
lb=zeros(3,1); ub=[+inf;+inf;1]; %x(3)为0-1变量
[x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)
（2）Matlab中线性 0 - 1 规划求解
- 0 -1 规划无非就是把上界定为1，下届定为0的整数规划 （限制 lb、un）
七、整数规划的典型例题
（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
（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 戊不参加
（3）例题三：钢管切割问题
%% 钢管切割问题
%% (1)枚举法找出同一个原材料上所有的切割方法
for i = 0: 2 % 2.9m长的圆钢的数量
for j = 0: 3 % 2.1m长的圆钢的数量
for k = 0:6 % 1m长的圆钢的数量
if 2.9*i+2.1*j+1*k >= 6 && 2.9*i+2.1*j+1*k <= 6.9
disp([i, j, k])
end
end
end
end
% 老的MATLAB版本，会出现浮点数计算的误差
% 只需要将上面的if这一行进行适当的放缩即可。
% if 2.9*i+2.1*j+1*k >= 6-0.0000001 && 2.9*i+2.1*j+1*k <= 6.9+0.0000001
%% (2) 线性整数规划问题的求解
c = ones(7,1); % 目标函数的系数矩阵
intcon=[1:7]; % 整数变量的位置(一共7个决策变量，均为整数变量)
A = -[1 2 0 0 0 0 1;
0 0 3 2 1 0 1;
4 1 0 2 4 6 1]; % 线性不等式约束的系数矩阵
b = -[100 100 100]'; % 线性不等式约束的常数项向量
lb = zeros(7,1); % 约束变量的范围下限
[x,fval]=intlinprog(c,intcon,A,b,[],[],lb)
八、最大化最小化模型
（1）模型的一般形式
[x,feval] = fminimax(@Fun,x0,A,b,Aeq,beq,lb,ub，@nonlfun,option)
（2）典型例题
（3）模型的求解
[x,feval] = fminimax(@Fun,x0,A,b,Aeq,beq,lb,ub，@nonlfun,option)
- 主函数
%% 最大最小化模型 : 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单位。
- 被调用的函数：Fun
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
九、多目标规划模型
（1）求解思路
- 先统一量纲和讲目标函数统一为最大化和最小化；
- 通过加权求和把多目标规划转换为单目标规划
- 求解得到符号约束条件的最优解，可以再通过最优解还原多目标函数
（2）典型例题
%% 多目标规划问题
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('综合指标的值')
% % 注意：代码文件仅供参考，一定不要直接用于自己的数模论文中
% % 国赛对于论文的查重要求非常严格，代码雷同也算作抄袭
十、拓展例题
（1）护士值班问题
%% 护士排班问题
clear
% (1) 系数向量
c = ones(6,1);
% (2) 不等式约束
A =zeros(6,6);
A(1,1) = -1; A(1,6) = -1;
for i = 1:5
A(i+1, i) = -1; A(i+1,i+1) = -1;
end
b = -[60 70 60 50 20 30]';
%（3）上下界
lb = zeros(6,1);
% 注意，这题应该是一个整数规划问题哦
intcon = [1:6];
[x fval] = intlinprog(c,intcon, A, b, [], [], lb)
% fval =
% 150
% 注：本题的x可能有多个解！
（2）货舱装运问题
%% 货机装货问题
clear
format long g %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
% 输入已知数据
w = [10 16 8]'; % 前、中、后仓的重量限制
v = [6800 8700 5300]'; % 前、中、后仓的体积限制
a = [18 15 23 12]'; % 货物的重量
b = [480 650 580 390]'; % 货物单位重量占用的空间
c = [3100 3800 3500 2850]'; % 货物单位重量的利润
% (1) 系数向量
cc = []; % 因为c这个符号被占用了，因此我们用cc表示系数向量
for i = 1:4
for j = 1:3
cc = [cc;c(i)]; % 用循环的方法生成cc
end
end
cc = -cc; % 因为要求最大值，所以这里对系数向量变号
% (2) 不等式约束
A = zeros(10,12); % 10个不等式约束
for i = 1:4
A(i,3*i-2) = 1; A(i,3*i-1) = 1; A(i,3*i) = 1;
end
for i = 5:7
j = i-4;
A(i,j) = 1; A(i,j+3) = 1; A(i,j+6) = 1; A(i,j+9) = 1;
end
for i = 8:10
j = i-7;
A(i,j) = b(1); A(i,j+3) = b(2); A(i,j+6) = b(3); A(i,j+9) = b(4);
end
b = [a;w;v];
% (3) 等式约束
Aeq = zeros(2,12);
for i = 1:4
Aeq(1,3*i-2) = 1/10;
Aeq(1,3*i-1) = -1/16;
end
for i = 1:4
Aeq(2,3*i-2) = 1/10;
Aeq(2,3*i) = -1/8;
end
beq = [0,0]';
%（4）上下界
lb = zeros(12,1);
% 进行求解
[x fval] = linprog(cc, A, b, Aeq, beq, lb)
fval = -fval
% fval =
% 121515.789473684
% 注：本题的x可能有多个解！
% 下面将x重新排列成4*3的矩阵
% 注意这里不能使用reshape函数，reshape函数是按照列的顺序进行转换的，也就是第一列读完，读第二列，即x1对应x_1,1，x2对应x_2,1
% 我们这题之前在定义x1-x12时，x1对应的是x_1,1, x2对应的是x_1,2
% 我们可以这样操作：先变形后转置
x = reshape(x,3,4)'
% 如果不能理解上面的操作，可以看下面这个简单的例题
% y = [1,2,3,4,5,6]; % 想要把y里面的元素按照行的顺序来转换成大小为3*2的矩阵
% yy = reshape(y,2,3) % 先按照列的顺序转换为2*3的矩阵
% yyy = yy' % 对上一步操作得到的结果来进行转置
（3）非线性规划的求解
- 主函数
%% 拓展三：非线性规划的求解
clear;clc
A=[-1 -2 0]; b=-1; % 线性不等式约束
x0 = [1 0 1]; % 初始值
lb = [0 -inf -inf]; % 下界
[x,fval] = fmincon(@fun,x0,A,b,[],[],lb,[],@nonlfun)
fval = -fval
% 先使用蒙特卡罗算法确定初始值，然后再进行计算
N = 1000000;
x1 = unifrnd(0,3,N,1); % x1在0~3之间均匀分布
x2 = unifrnd(-8,7,N,1); % x2在-8~7之间均匀分布
x3 = 2-x1.^2; % x1^2 + x3 = 2, 注意这里要用“.^”，表示每个元素分别平方
fmax=-inf; % 初始化函数f的最大值为负无穷（后续只要找到一个比它大的我们就对其更新）
for i=1:N
x = [x1(i), x2(i), x3(i)]; %构造x向量, 这里千万别写成了：x =[x1, x2, x3]
if (x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10<=0) & (x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50<=0) & (2*x(1)+x(1)^2+2*x(2)+x(3)-40<=0) & (x(1)+2*x(2)>=1) % 判断是否满足条件
result =2*x(1)+3*x(1)^2+3*x(2)+x(2)^2+x(3) ; % 如果满足条件就计算函数值
if result > fmax % 如果这个函数值大于我们之前计算出来的最大值
fmax = result; % 那么就更新这个函数值为新的最大值
x0 = x; % 并且将此时的x1 x2 x3更新为初始值
end
end
end
disp('蒙特卡罗算法选取的初始值为：'); disp(x0)
[x,fval] = fmincon(@fun,x0,A,b,[],[],lb,[],@nonlfun)
fval = -fval
- 被调用函数：fun （目标函数）
function f=fun(x)
f=2*x(1)+3*x(1)^2+3*x(2)+x(2)^2+x(3);
f=-f; % 最大值加负号改为最小值
end
- 被调用函数：nonlfun （非线性约束条件）
function [c,ceq]=nonlfun(x)
% 非线性不等式约束
c=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10
x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50
2*x(1)+x(1)^2+2*x(2)+x(3)-40];
% 非线性等式约束
ceq=x(1)^2+x(3)-2;
end
（4）覆盖问题
%% 第四题：覆盖问题
c = ones(1,6); % 目标函数的系数矩阵
intcon=[1:6]; % 整数变量的位置(一共6个决策变量，均为0-1整数变量)
% 线性不等式约束的系数矩阵和常数项向量
A =- [1 1 1 0 0 0;
0 1 0 1 0 0;
0 0 1 0 1 0;
0 0 0 1 0 1;
0 0 0 0 1 1;
1 0 0 0 0 0;
0 1 0 1 0 1];
b = -ones(7,1);
Aeq = []; beq = []; % 不存在线性等式约束
lb = zeros(1,6); % 约束变量的范围下限
ub = ones(1,6); % 约束变量的范围上限
%最后调用intlinprog()函数
[x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
