数学建模十大算法03—线性规划、整数规划、非线性规划、多目标规划

一、线性规划（Linear Programming，LP）

1.1 引例

在人们的生产实践中，经常会遇到如何利用现有资源来安排生产，以取得最大经济效益的问题。 此类问题构成了运筹学的一个重要分支一数学规划，而线性规划（Linear Programming, LP） 则是数学规划的一个重要分支。

简而言之，线性规划就是在有限的条件下，获得最大的收益。

【例】 某机床厂生产甲、乙两种机床，每台销售后的利润分别为4000元与3000元。生产甲机床需用A、B机器加工，加工时间分别为每台2h和1h；生产乙机床需用A、B、C三种机器加工，加工时间为每台各1h。若每天可用于加工的机器时数分别为A机器10h、B机器8h和C机器7h，问该厂应生产甲、乙机床各几台，才能使总利润最大?

上述问题的数学模型为：（其中，$x_1$ 和 $x_2$ 是决策变量）

总之，线性规划问题是在一组线性约束条件的限制下，求一线性目标函数最大或最小的问题。 “线性”意味着所有变量都是一次方。
在解决实际问题时，把问题归结成一个线性规划数学模型是很重要的一步，往往也是很困难的一步，模型建立得是否恰当，直接影响到求解。而选择适当的决策变量，是建立有效模型的关键之一。

1.2 线性规划问题的解

一般线性规划问题的（数学）标准型为：

1.3 Matlab标准形式

线性规划的目标函数可以是求最大值，也可以是求最小值，约束条件的不等号可以是
小于等号也可以是大于等号。为了避免这种形式多样性带来的不便，Matlab中规定线性规划的标准形式为：

Matlab中求解线性规划的命令为：

[x, fval] = linprog(f, A, b)  % A和b对应线性不等式约束
[x, fval] = linprog(f, A, b, Aeq, beq)  % Aeq和beq对应线性等式约束
[x, fval] = linprog(f, A, b, Aeq, beq, lb, ub)  % lb和ub分别对应决策变量的下界向量和上界向量
1
2
3

其中，$x$ 返回决策向量的取值；$fval$ 返回目标函数的最优值。

如果线性规划问题的解为max形式，可以通过以下转换，再应用Matlab提供的方法：

【例】 求解下列线性规划问题。

将其转换为Matlab标准形式：

Matlab代码实现为：

f = [-40;-30]  % 目标函数中变量的系数矩阵
a = [1,1;-1,0;0,-1;240,120]  % 小于等于的约束条件中变量系数矩阵
b = [6;-1;-1;1200]  % 小于等于的约束条件中常数项矩阵

[x,y] = linprog(f, a, b)  % 从等式约束开始后都没有，可以都不写
y = -y  % 注意要还原成求最大值！！！！
1
2
3
4
5
6

求解结果为：

【例】 求解下列线性规划问题。

代码如下：

clc;clear
c=[2;3;1];
a=[1,4,2;3,2,0];
b=[8;6];
[x,y]=linprog(c,-a,-b,[],[],zeros(3,1))  %这里没有等式约束，对应的矩阵为空矩阵
1
2
3
4
5

求解结果为：

1.4 投资的收益和风险（模型建立与分析）

1.4.1 符号规定和基本假设

符号规定：

基本假设：

注意，这里将M假设为1，可以方便后面的计算。（建模比赛中可以适当进行假设）

1.4.2 模型的分析与建立

（1）总体风险用所投资的 $s_i$ 中最大的一个风险来衡量，即

（2）购买 $s_i（i=1,2,...,n)$ 所付交易费是一个分段函数，即

而题目中所给的定值 $u_i$ （单位：元）相对总投资 $M$ 很少， $p_i{u}_i$ 更小，这样购买 $s_i$ 的净收益可以简化为 $(r_i-p_i)x_i$ 。

（3）要使净收益尽可能大，总体风险尽可能小，这是一个多目标规划模型 。

目标函数为：

约束条件为：

（4）模型简化。

模型一： 在实际投资中，投资者承受风险的程度不一样，若给定风险一个界限 a，使最大的一个风险率为a，即 $q_i{x}_i/M\le a(i=1,2,\dots ,n)$，可找到相应的投资方案。这样以来就可以把多目标规划变成一个目标的线性规划。

即固定风险水平，优化收益：

通俗理解，就是把原来多目标优化中第二个 $min$ 目标函数转变为约束条件，只要风险小于 $a$ 都是投资者可接受的风险范围。

模型二： 若投资者希望总盈利至少达到水平 $k$ 以上，在风险最小的情况下寻求相应的投资组合。

即固定盈利水平，极小化风险：


和模型一类似，这里将 $max$ 目标函数转换为了只要总收益大于 $k$ ，投资者都可以接受。

模型三： 投资者在权衡资产风险和预期收益两方面时，希望选择一个令自己满意的投资组合。因此对风险收益分别赋予权重 $s$ 和 $(1-s)$ ，$s$ 称为投资偏好系数。

这里可以简单理解成投资者是更注重哪一个目标，投资者倾向的那个权重设置的更大即可。

上述模型对应的两条思路也是在数学建模中经常会使用到的：①将多目标转换为单目标规划；②将多目标函数赋予相应的权重进行搜索。

1.4.3 模型一的求解

模型一为：

注意：这里将原目标函数 $max$ 转换成了Matlab标准型的 $min$ 形式，就要在所有值前面加负号。
【数值细节解释】 $x_0$ 对应的-0.05代表的是题目中假定的同期银行存款利率 $r_0=0.5$ ，既无交易费又无风险的条件。 $x_1$ 对应的是当投资第一个资产 $s_1$ 时，其对应的 $-(r_i-p_i)x_i$ 值。同理可得其他值。

由于 $\alpha$ 是任意给定的风险度，不同的投资者有不同的风险度。下面从 $\alpha =0$ 开始，以步长 $△\alpha =0.001$ 进行循环搜索。Matlab求解代码如下：

clc,clear
a=0;
hold on
while a<0.05
    c=[-0.05,-0.27,-0.19,-0.185,-0.185];
    A=[zeros(4,1),diag([0.025,0.015,0.055,0.026])];  % 生成对角矩阵
    b=a*ones(4,1);  % 不等式约束右边的值
    Aeq=[1,1.01,1.02,1.045,1.065];
    beq=1;
    LB=zeros(5,1);  % 下界
    [x,Q]=linprog(c,A,b,Aeq,beq,LB);
    Q=-Q;  % 注意求得是max！！！！！！！！！！！！
    plot(a,Q,'*k');
    a=a+0.001;  % 以步长α=0.001进行循环搜索
end
xlabel('a'),ylabel('Q')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

运行结果为：

1.4.4 结果分析

风险 $\alpha$ 与收益 $Q$ 之间的关系如上图所示，我们可以看出：
　　(1) 风险大，收益也大。
　　(2) 当投资越分散时，投资者承担的风险越小，这与题意一致。冒险的投资者会出现集中投资的情况，保守的投资者则尽量分散投资。
　　(3) 在 $a=0.006$ 附近有一个转折点，在这一点左边，风险增加很少时，利润增长很快；在这一点右边，风险增加很大时，利润增长很缓慢。所以对于风险和收益没有特殊偏好的投资者来说，应该选择曲线的转折点作为最优投资组合，大约是 $a=0.6$%，$Q=20$%。
　　
所对应的投资方案我们可以在代码中加入a=0.6%查看具体的 $x$ 和 $Q$ ：

     while a == 0.006
        disp(x);
        disp(Q);
        break;
    end
1
2
3
4
5

得到具体的投资方案为：

二、整数规划

2.1 引例

数学规划中的变量（部分或全部）限制为整数时，称为整数规划。若在线性规划模型中，变量限制为整数，则称为整数线性规划。目前所流行的求解整数规划的方法，往往只适用于整数线性规划。目前还没有一种方法能有效地求解一切整数规划。

整数规划的分类：
如不加特殊说明，则一般指整数线性规划。整数线性规划模型大致分为两类：
1）变量全限制为整数时，称纯（完全）整数规划。
2）变量部分限制为整数时，称混合整数规划。

整数规划特点：
1）原线性规划有最优解，当自变量限制为整数后，其整数规划解出现下述情况。
　　①原线性规划最优解全是整数，则整数规划最优解与线性规划最优解一致。
　　②整数规划无可行解。

③有可行解（当然就存在最优解），但最优解值变差。

2）整数规划最优解不能按照实数最优解简单取整而获得。

求解方法分类：

注：蒙特卡罗算法请看之前的笔记：数学建模十大算法01-蒙特卡洛算法（Monte Carlo）

2.2 “0-1 型”整数规划

0-1型整数规划是整数规划中的特殊情形，它的变量 $x_j$ 仅取值0或1。这时 $x_j$ 称为0-1变量，或称二进制变量。 $x_j$ 仅取值0或1这个条件可由下述约束条件 $$0\le x_j\le 1且x_j为整数$$所代替，是和一般整数规划的约束条件形式一致的。在实际问题中，如果引人0-1变量，就可以把有各种情况需要分别讨论的数学规划问题统一在一个问题中讨论了。下面先介绍引入0-1变量的实际问题。

2.2.1 相互排斥的约束条件

有两种运输方式可供选择，但只能选择一运输方式，或者用车运输，或者用船运输。
用车运输的约束条件为$5x_1+4x_2\le 24$，用船运输的约束条件为$7x_1+3x_2\le 45$。即有两个相互排斥的约束条件$$5x_1+4x_2\le 24或7x_1+3x_2\le 45$$
为了统一在一个问题中，引入0-1变量：

则上述约束条件可改写为：

把相互排斥的约束条件改成普通的约束条件，未必需要引进充分大的正实数，例如相互排斥的约束条件：

如果有 $m$ 个互相排斥的约束条件：

2.2.2 指派问题的数学模型

拟分配 $n$ 人去做 $n$ 项工作，每人做且仅做1项工作，若分配第 $i$ 人去做第 $j$ 项工作，需花费 $c_{ij}$ 单位时间，问应该如何分配工作才能使工人花费的总时间最少？

思路： 只需要给出矩阵 $C=c_{ij}$ （指派问题的系数矩阵）。

引入0-1变量：

上述指派问题的数学模型为：

上述指派问题的可行解可以用一个矩阵表示，其每行、每列均有且只有一个元素为1，其余元素均为0；还可以用 $1,...,n$ 中的一个置换表示。

Matlab中指派问题可用以下语法格式实现：

[x,fval]=intlinprog(f,intcon,a,b,aeq,beq,lb,ub)
1

对应下列数学模型：

【例】 求解下列指派问题，已知指派矩阵为：

这里需要把二维决策变量 $x_{ij}(i,j=1,...,5)$ 变为一维决策变量 $y_k(k=1,...,25)$。可以理解成将每行每列只有一个“1”转换为一维考虑。

Matlab实现代码如下：

clc,clear
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(:);  % 将矩阵c中的每列（！！！）合并成一个长的列向量
a = zeros(10,25);
intcon = 1:25;

for i = 1:5
    a(i,(i-1)*5+1:5*i) = 1;  % 第i行，连续5个值为1----> 
    a(5+i,i:5:25) = 1;  % 第5+i行，每隔5个赋值1---->
end

b = ones(10,1);
lb = zeros(25,1);
ub = ones(25,1);
x = intlinprog(c,intcon,[],[],a,b,lb,ub)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

把5个物品放入一个5×5的方阵中，且每行、每列都必须且只能有一个物品。使用混合整数线性规划的函数intlinprog求解，就要把这5x5的位置看成25个优化变量，每个变量只能取0或1。a和b表示等式约束，10行表示10个等式，其中前5个表示每行有一个物品，后5行表示每列有一个物品。

最终求得的最优指派方案为：用x = reshape(x,[5,5])查看：

【例】 求解如下的混合整数规划问题：

Matlab代码如下：

clc,clear
f=[-3;-2;-1]; intcon=3;  % 整数变量的地址
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=intlinprog(f,intcon,a,b,aeq,beq,lb,ub)
z=dot(f,x)  % 目标函数的最优值为-12
1
2
3
4
5
6
7

求解结果为：

【例】 某公司新购置了某种设备6台，欲分配给下属的4个企业，已知各企业获得这种设备后年创利润如表所列（单位：千万元）。问应如何分配这些设备能使年创总利润最大，最大利润是多少？

思路：用 $j=1,2,3,4$ 分别表示甲乙丙丁四个企业，$c_{ij}$ 表示第 $i(i=1,...,6)$ 台设备分配给第 $j$ 个企业创造的利润，引进0-1变量：

则问题的数学模型为：

Matlab代码如下：（和上述代码不同，这里用的是基于问题的求解方法）

% 非常规指派问题
clc,clear
prob = optimproblem('ObjectiveSense','max');  % 基于问题的求解方法：ObjectiveSense可以是max和min，代表优化最大值还是最小值，默认是min
x = optimvar('x',6,4,'TYPE','integer','LowerBound',0,'UpperBound',1);  % x是变量名，x为6行4列。该函数所属类型为int。下界为0，上界为1.
c = [4 2 3 4;6 4 5 5;7 6 7 6;7 8 8 6;7 9 8 6;7 10 8 6];
M=c.*x;
prob.Objective = sum(sum(M));  % 目标函数，目标函数需要得到一个标量数值，不是矩阵向量
% 约束条件 .con是标签，可以自己命名，多个约束条件时，必须标签不能一样。
prob.Constraints.con2 =sum(x)>=1;  % 每列求和，结果大于等于1
prob.Constraints.con1 =sum(x,2)==1;  % 每行求和，结果等于1
[sol,flav,flag] = solve(prob);  % fval是最优值，sol.x是决策变量的值，当多个决策变量时，可以sol.y，flag在线性规划中不用在意，在非线性规划中注意不能为负值。
xx = sol.x
sum(sum(c.*xx))
1
2
3
4
5
6
7
8
9
10
11
12
13

Python实现代码如下所示：

# 决策变量
n = 8
# nonneg参数，变量是否为非负
x = cp.Variable(n,nonneg = True)

# 约束1
A1 = np.array([[1,1,1,0,0,0,0,0],
              [0,1,0,1,0,0,0,0],
              [0,0,1,0,1,0,0,0],
              [0,0,0,1,0,1,0,0],
              [0,0,0,0,1,1,0,0],
              [1,0,0,0,0,0,0,0],
              [0,1,0,1,0,1,0,0]])
b1 = np.array([1,1,1,1,1,1,1])

# 约束2
A2 = np.ones((8, 8))
for i in range(A2.shape[0]):
    for j in range(A2.shape[1]):
        if i == j:
            pass
        else:
            A2[i, j] = A2[i, j]*0

b2 = np.array([0,0,0,0,0,0,0,0])
# 约束3
A3 = np.ones((8, 8))
for i in range(A3.shape[0]):
    for j in range(A3.shape[1]):
        if i == j:
            pass
        else:
            A3[i, j] = A3[i, j]*0
b3 = np.array([1,1,1,1,1,1,1,1])

# 目标函数
c = np.array([1,1,1,1,1,1,1,1])

# 定义问题，添加约束条件
prob = cp.Problem(cp.Minimize(c.T @ x),
                  [A1 @ x >= b1,
                  A2 @ x >= b2,
                  A3 @ x <= b3])

# 求解
ans = prob.solve(solver=cp.GLOP)
# 输出结果
print("目标函数最小值:", ans)
# 对x向量各元素取整数后再输出
print(x.value)

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

运行结果为：

更多用Python实现的例题可参考：数学建模算法与应用—整数规划（cvxpy包）

指派问题的求解可以使用匈牙利算法、拍卖算法等算法，这里就不讨论了。

三、非线性规划

3.1 引例

如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。简而言之，就是至少有一个变量不是一次方。

一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不像线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。

3.2 Matlab求解

Matlab中非线性规划的的数学模型写成以下形式：

Matlab中的命令是：

[x, fval] = fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon)
1

其中，各项的含义如下表：

注意：这里的fun是单独函数文件里面定义的目标函数。

【例】 求下列非线性规划：

第一步：编写M函数fun1.m定义目标函数：

第二步：编写M函数fun2.m定义非线性约束条件：

第三步：编写主程序文件如下：

求得最优情况下，x和y的取值为：

3.3 数学建模中的应用

非线性规划适用的典型赛题有：
题目中提到“怎么安排/分配”、“尽量多（少）”、“最多（少）”、“利润最大”、“最合理”等词，且变量要非一次方。

• 投资规划： 组合投资、总收益率最大/最大投资方案
• 角度调整： 飞行管理避免相撞；影院最佳视角等；设计三角函数为非线性
• 生产安排： 原材料、设备有限制，总利润最大（目标函数或约束条件含有非线性变量）

★【例】飞行管理问题。
在约10000m高空的某边长160km的正方形区域内，经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞。现假定条件如下：
　　　
　　(1) 不碰撞的标准为任意两架飞机的距离大于8km。
　　(2) 飞机飞行方向角调整的幅度不应超过30°。
　　(3) 所有飞机飞行速度均为800km/h。
　　(4) 进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60km以上。
　　(5) 最多需考虑6架飞机。
　　(6) 不必考虑飞机离开此区域后的状况。
　　
请对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计
算（方向角误差不超过0.01度），要求飞机飞行方向角调整的幅度尽量小。
设该区域4个顶点的坐标为(0,0)、(160,0)、(160,160)、(0,160)。记录数据见表3.1。

注：方向角指飞行方向与 $x$ 轴正向的夹角。

为方便以后的讨论，引进如下记号：

模型一： 目标函数为所有飞机的调整量绝对值之和的最小值。

根据相对运动的观点在考察两架飞机 $i$ 和 $j$ 的飞行时，可以将 $i$ 视为不动，而飞机 $j$ 以相对速度

相对于飞机 $i$ 运动，对上式进行适当的计算，得：

不妨设 ${\theta }_j\ge {\theta }_i$ ，此时相对飞行方向角为，如图所示：


由于两架飞机的初始距离为：

因此只要当相对非新方向角 ${\beta }_{ij}$ 满足：

时，两架飞机就不可能碰撞。

本问题中的优化目标函数可以有不同的形式：如使所有飞机的最大调整量最小，所有飞机的调整量绝对值之和最小等。这里以所有飞机的调整量绝对值之和最小为目标函数，可以得到如下的数学规划模型：

Matlab代码如下所示：

clc,clear
x0=[150 85 150 145 130 0];
y0=[140 85 155 50 150 0];
q=[243 236 220.5 159 230 52];
xy0=[x0; y0];
d0=dist(xy0);   %求矩阵各个列向量之间的距离
d0(find(d0==0))=inf;
a0=asind(8./d0)  %以度为单位的反函数
xy1=x0+i*y0
xy2=exp(i*q*pi/180)
for m=1:6
     for n=1:6
         if n~=m
         b0(m,n)=angle((xy2(n)-xy2(m))/(xy1(m)-xy1(n))); 
         end
     end
end
b0=b0*180/pi;
dlmwrite('txt1.txt',a0,'delimiter', '\t','newline','PC');
dlmwrite('txt1.txt','~','-append');       %往纯文本文件中写LINGO数据的分割符
dlmwrite('txt1.txt',b0,'delimiter', '\t','newline','PC','-append','roffset', 1)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

求得 ${\alpha }_{ij}^0$ 的值如表所示：

求得 ${\beta }_{ij}^0$ 的值如表所示：

该题还有其他思路，可以参考：【数学建模】通过调整飞行角度使飞机顺利飞行（Matlab）

四、多目标规划

4.1 引例

简单理解就是：既要……，又要……。
【例】某工厂生产产品Ⅰ和产品Ⅱ，有关数据下，若只追求最大化利润，得到模型：

但现在有3个目标：①尽量使产品Ⅰ的产量不超过产品Ⅱ的产量；②尽可能充分利用所有设备；③尽可能使利润不少于56w。

绝对约束： 模型自带的约束条件，必须满足，否则是不可行解。
目标约束： 模型中对不等式右端追求的值允许有偏差。

上述三个目标使我们要尽可能实现的目标。目标1是“不超过”，也就是尽量“≤”；目标2是“充分利用”，也就是尽量“=”；目标3是“不少于”，也就是尽量“>”。

需要衡量每个目标的完成情况，并主观上区分三个目标的重要性，使得整体的完成情况尽量好。

衡量每个目标的完成情况：正负偏差变量，绝对约束和目标约束，优先因子。


加上与目标值差值，减去超过目标值的部分。即多退少补，把目标函数变成了等式约束条件。

优先因子： 主观上确定优先因子 $P_k$ ，优先满足那个目标。

回到最原始的例子，得到多目标规划模型：

其中，第二个等式代表第一个目标：尽量使产品Ⅰ的产量不超过Ⅱ的产量。第三个等式表示第二个目标：尽可能充分利用所有设备。第四个等式表示第三个目标：尽可能使利润不少于56万。注意这里的优化函数中P体现了重要性。

4.2 Matlab求解

fgoalattain函数或者序贯算法或者用Lingo求解。

语法格式如下：

[X,FVAL] = fgoalattain(fun,x0,goal,a,b,Aeq,Beq,LB,UB,nonlcon)；
1

1. X 为最终解 ， FVAL为最终解对应的函数值。 注意：求最大值时，结果FVAL需要取反
2. fun是定义的决策函数，通常通过M文件或者匿名函数进行定义。 注意：当所求为最大值时，系数需要取反。
3. goal为欲达到的目标，通常通过linprog函数先计算得到每个决策函数目标值
4. a 为约束条件中不等式组的系数矩阵 ，a的列数等于f的列数。 注意：当不等号为 > 或 ≥ 时，矩阵需要取反
5. b 为约束条件中不等式组右边的值。 注意：当不等号为 > 或 ≥ 时，矩阵需要取反
6. Aeq 为约束条件中等式组的系数矩阵 ，Aeq的列数等于f的列数
7. Beq 为约束条件中等式组右边的值
8. LB、UB是解的范围
9. nonlcon 为定义的向量函数

【例】 求解多目标线性规划问题。

Matlab求解代码如下：

a = [-1,-1,0,0;0,0,-1,-1;3,0,2,0;0,3,0,2];  % 不等式约束左边x系数
b = [-30;-30;120;48];  % 不等式约束右边值
c1 = [-100,-90,-80,-70];  % 第一个目标函数的系数
c2 = [0,3,0,2];  % 第二个目标函数的系数
fun = @(x)[c1;c2]*x;  % 目标函数
[x1,g1] = linprog(c1,a,b,[],[],zeros(4,1));  % 求解第一个目标
[x2,g2] = linprog(c2,a,b,[],[],zeros(4,1));  % 求解第二个目标

g3 = [g1,g2];  % ★★
[x,fval] = fgoalattain(fun,rand(4,1),g3,abs(g3),a,b,[],[],zeros(4,1))  % abs(g3)为权重项
1
2
3
4
5
6
7
8
9
10

运行结果为：

参考资料

[1] 数学建模——整数规划笔记
[2] 数学建模算法与应用—整数规划（cvxpy包）
[3] 数学建模3.9—多目标规划