Github仓库地址：https://github.com/dekrt/Reports

数学建模上机实验

一、问题重述

​ 某医院每天各时间段内需要的值班护士数如表 1 所示：

$$ \[\begin{aligned} &\begin{array}{|c|c|} \hline \text { 时间区段 } & \text { 护士数量 } \\ \hline 6: 00 \sim 10: 00 & 18 \\ \hline 10: 00 \sim 14: 00 & 20 \\ \hline 14: 00 \sim 18: 00 & 19 \\ \hline 18: 00 \sim 22: 00 & 17 \\ \hline 22: 00 \sim 6: 00 \text { (次日) } & 12 \\ \hline \end{array} \end{aligned}\]

$$ ​ 该医院护士上班分五个班次，每班8小时，具体上班时间为第一班2:00~10:00，第二班6:00~14:00，第三班10:00~18:00，第四班14:00~22:00，第五班18:00~2:00（次日）。每名护士每周上5个班，并被安排在不同的日子，由一名总护士长负责护士的值班安排。值班方案要做到在人员或经济上比较节省，又做到尽可能合情合理。下面是一些正在考虑中的值班方案：

​ 方案1：每名护士连续上班5天，休息2天，并从上班第一天起按从上第一班到第五班顺序安排。

​ 方案2：考虑到方案1中每名护士在周末（周六、周日）两天内休息安排不均匀，于是规定每名护士在周六、周日两天内安排一天、且只安排一天休息，再在周一至周五期间安排4个班，同样上班的5天内分别顺序安排5个不同班次。

​ 在对方案1、2建立线性规划模型并求解后发现，方案2虽然在安排周末休息上比较合理，但所需值班人员要比方案1有较多增加，经济上不太合算，于是又提出了第3方案。

​ 方案3：在方案2的基础上，动员一部分护士放弃周末休息，即每周在周一至周五间由总护士长给安排三天值班，加周六周日共上五个班，同样五个班分别安排不同班次。作为奖励，规定放弃周末休息的护士，其工资和奖金总额比其他护士增加a%。

​ 根据上述方案，帮助总护士长分析研究：

​ （1）对方案1、2建立使值班护士人数为最少的线性规划模型并求解。

​ （2）对方案3，同样建立使值班护士人数为最少的线性规划模型并求解，然后回答a的值为多大时，第3方案较第2方案更经济。

二、问题分析

​ 从该医院各时间段护士值班表可看出：五个时间段所需护士人数分别为 18,20,19,17,12。每个护士每周值五个班且安排在不同日子。为了使人员和经济较节省，并且安排合情合理。现制定以下三种方案，通过建立线性规划模型并求解，分析各种方案的最佳安排方式。

​ 方案一：要求每名护士连续工作 5 天，接着休息两天，从第一天起，按第一班到第五班的顺序依次安排工作。例如从周三起上班，即周三上第一班，周四上第二班....建立线性规划模型并求解使值班护 士人数最少。

​ 方案二：由于方案 1 在周末两天内休息安排不均匀，于是重新规定每名护士在周六周日两天内至少安排一天休息，再通过从第一天起，按第一班到第五班的顺序依次安排工作。同样的建立线性规划模 型求解使值班护士人数最少。

​ 方案三：通过方案 1，方案 2 的比较，由于在经济上不太合算，于是又提出了方案 3。首先动员部 分护士放弃周末休息，即每周在周一至周五期间工作三天，加周六周日共五个班，同样五个班在不同班次。其他护士每周六周日均休息一天，规定周末不休息的护士其工资和奖金总额比其他护士增加 a%作为奖励。建立线性规划模型求解后，得到 a 的值为多大时，方案三较方案二更经济。

三、模型假设

  1. 每位护士在一个工作日内只能被安排一个班次。
  2. 所有护士都有相同的工作效率，即每位护士在任何班次中提供的工作量是相同的。
  3. 每个班次都必须有足够的护士来满足表格中的需求。
  4. 护士的工资是根据工作的小时数来计算的，而不考虑具体的班次时间。
  5. 对于方案3中提到的奖励，假设增加的工资和奖金总额是足够吸引一部分护士放弃周末休息的。

四、模型的建立与求解

4.1 方案一模型的建立与求解

方案一通过顺序安排各个班次的方法来排班的。因此，由题意得出每天各个班次上 班的护士人数如下表所示：（xi 表示从周 i 开始工作的护工人数）由每个护士上 5 个班次可以看出，每组护士上 5 个班次，而此方案有 35 个班次，共需要 7 组护士。

班次/星期 周一 周二 周三 周四 周五 周六 周日
2:00~10:00 X1 X2 X3 X4 X5 X6 X7
6:00~14:00 X7 X1 X2 X3 X4 X5 X6
10:00~18:00 X6 X7 X1 X2 X3 X4 X5
14:00~22:00 X5 X6 X7 X1 X2 X3 X4
18:00~2:00（次日） X4 X5 X6 X7 X1 X2 X3

为此，我们可以得到如下的约束关系：

班次/星期 一 二 三 四 五 六 日
6：00-10：00 X1+x7>=18 X2+x1>=18 x3+x2>=18 x4+x3>=18 X5+x4>=18 X6+x5>=18 X7+x6>=18
10：00-14：00 X1+x6>=20 X1+x7>=20 X1+x2>=20 X3+x2>=20 X4+x3>=20 x5+x4>=20 X6+x5>=20
14：00-18：00 X6+x5>=19 X7+x6>=19 X1+x7>=19 X2+x1>=19 X3+x2>=19 X4+x3>=19 X5+x4>=19
18：00-22：00 X5+x4>=17 X6+x5>=17 X7+x6>=17 X1+x7>=17 X2+x1>=17 X3+x2>=17 X4+x3>=17
22：00-6：00 X4>=12
X2>=12		 X5>=12
X3>=12		 X6>=12
X4>=12		 X7>=12
X5>=12		 X1>=12
X6>=12		 X2>=12
X7>=12		 X3>=12
X1>=12

为了满足题意，我们使用线性规划模型建模如下： \[ \begin{aligned} & \min f(x)=x_1+x_2+x_3+x_4+x_5+x_6+x_7 \\ & x_1+x_7 \geq 20 \\ & x_7+x_6 \geq 20 \\ & x_6+x_5 \geq 20 \\ & x_5+x_4 \geq 20 \\ & x_4+x_3 \geq 20 \\ & x_3+x_2 \geq 20 \\ & x_2+x_1 \geq 20 \\ & x_i \geq 12(i=1,2, \cdots, 7) \end{aligned} \] 使用matlab编写程序求解（源程序见附录），得到结果如下：

1
2
3
4
5
6
7
8
9
10
Solution found:
Objective function value: 84
Variables:
    12
    12
    12
    12
    12
    12
    12

结果表明，最优解为星期一上第一班的班组的人数为12人，星期二上第一班的班组的人数为12人，星期三上第一班的班组的人数为12人，星期四上第一班的班组的人数为12人，星期五上第一班的班组的人数为12人，星期六上第一班的班组的人数为12人，星期日上第一班的班组的人数为12人。总人数84人。

4.2 方案二模型的建立与求解

因为每名护士在周六、周日两天里必须工作一天，安排休息一天。周一到周五连续安排4个班，所以可以先安排周末的护士值班情况：周六、周末两天共10个班次，用Xj (j=1,2,3,…10)表示周六周末两天10个班次的护士人数，其中 X1-X5分别代表周六第1个到第5个班次的护士人数，X6-X10分别代表周日从第1个到第5个班次的护士人数。其值班安排表如下：

星期班次 一 二 三 四 五 六 日
2：00——10：00 X₁₀ X₅+X₉ X₄+X₈ X₃+X₇ X₂ X₁ X₆
6：00——14：00 X₆ X₁+X₁₀ X₅+X₉ X₄+X₈ X₃ X₂ X₇
10：00——18：00 X₇ X₂+X₆ X₁+X₁₀ X₅+X₉ X₄ X₃ X₈
14：00——22：00 X₈ X₃+X₇ X₂+X₆ X₁+X₁₀ X₅ X₄ X₉
18：00——2：00 X₉ X₄+X₈ X₃+X₇ X₂+X₆ X₁ X₅ X₁₀

因此，满足如下的约束关系：

星期 一 二 三 四 五 六 日
6：00-10：00 x10+x6>=18 X1+x5+x9+x10>=18 X4+x8+x5+x9>=18 X3+x7+x4+x8>=18 X2+x3>=18 X1+x2>=18 X6+x7>=18
10：00-14：00 X6+x7>=20 X1+x2+x6+10>=20 X1+x5+x9+x10>=20 X4+x8+x5+x9>=20 X3+x4>=20 X2+x3>=20 X7+x8>=20
14：00-18：00 X7+x8>=19 X2+x6+x3+x7>=19 X1+x2+x6+x10>=19 X1+x5+x9+x10>=19 X4+x5>=19 X3+x4>=19 X8+x9>=19
18：00-22：00 X8+x9>=17 x3+x7+x4+x8>=17 X2+x6+x3+x7>=17 X1+x3+x6+x10>=17 X5+x1>=17 X4+x5>=17 X9+x10>=17
22：00-6：00 X9>=12
X5+x9>=12		 X4+x8>=12 X3+x7>=12 X3+x6>=12
X2>=12		 X1>=12 X5>=12
X6>=12		 X10>=12

因此，我们可以建立下面的数学规划模型： \[ \begin{aligned} & \min z=x_1+x_2+x_3+x_4+x_5+x_6+x_7+x_8+x_9+x_{10} \\ & x_6+x_{10} \geqslant 18 \\ & x_6+x_7 \geqslant 20 \\ & x_3+x_4 \geqslant 20 \\ & x_4+x_5 \geqslant 19 \\ & x_1+x_5 \geqslant 17 \\ & x_1+x_2 \geqslant 18 \\ & x_2+x_3 \geqslant 20 \\ & x_7+x_8 \geqslant 20 \\ & x_8+x_9 \geqslant 19 \\ & x_9+x_{10} \geqslant 17 \\ & x_1+x_2+x_6+x_{10} \geqslant 20 \\ & x_1+x_5+x_9+x_{10} \geqslant 20 \\ & x_4+x_5+x_8+x_9 \geqslant 20 \\ & x_4+x_8 \geqslant 12 \\ & x_3+x_7 \geqslant 12 \\ & x_1 \geqslant 12 \\ & x_2 \geqslant 12 \\ & x_5 \geqslant 12 \\ & x_6 \geqslant 12 \\ & x_9 \geqslant 12 \\ & x_{10} \geqslant 12 \\ & x_j \geqslant 0 ; j=1,2,3 \ldots 10 \end{aligned} \]

使用matlab求解，结果如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
Solution found:
Objective function value: 112
Variables:
	12
    12
     8
    12
    12
    12
    13
     7
    12
    12

4.3 方案三模型的建立与求解

分析方案3的突破口主要有以下几点：1、一部分护士周末两天都上班，另外一部分护士周末只上一天。2、连续上班5天，休息2天。3、同样5个班分别安排在不同的班次。因此，先安排周末的值班，设：X1- X5周末两天都上班。X6-X15周末只上一天。 对方案3进行分析，以表格的形式将方案3的护士值班安排表示如下表所示：

星期班次 一 二 三 四 五 六 七
2:00-10:00 x4+x15 x3+x14+x10 x2+x13+x9 x12+x8 X7 x1+x6 x5+x11
6:00-14:00 x5+x11 x4+x15+x6 x3+x14+x10 x13+x9 X8 x2+x7 x1+x12
10:00-18:00 x1+x12 x5+x11+x7 x4+x15+x6 x14+x10 X9 x3+x8 x2+x13
14:00-22:00 x2+x13 x1+x12+x8 x5+x11+x7 x15+x6 X10 x4+x9 x3+x14
18:00-2:00 x3+x14 x2+x13+x9 x1+x12+x8 x11+x7 x6 x5+x10 x4+x15

满足下面的约束条件：

星期 一 二 三 四 五 六 日
6：00-10：00 X4+x15+x5+x11>=18 X3+x14+x10+x4+x15+x6>=18 x2+x13+x9+X3+x14+x10>=18 X12+x8+x13+x9>=18 X7+x8>=18 X1+x2+x6+x7>=18 X5+x11+x1+x12>=18
10：00-14：00 X1+x5+x11+x12>=20 X4+x15+x6+x5+x11+x7>=20 X3+x14+x10+x4+x15+x6>=20 X13+x9+x14+x10>=20 X8+x9>=20 X2+x7+x3+x8>=20 X1+x12+x2+x13>=20
14：00-18：00 X1+x12+x2+x13>=19 X5+x11+x7+x1+x12+x8>=19 X4+x15+x6+x5+x11+x7>=19 X14+x10+x15+x6>=19 X9+x10>=19 X3+x8+x4+x9>=19 X2+x13++x3+x14>=19
18：00-22：00 X2+x13+x3+x14>=17 x1+x12+x8+x2+x13+x9>=17 X5+x11+x7+x1+x12+x8>=17 X15+x6+x11+x7>=17 X10+x6>=17 x4+x9+x5+x10>=17 X3+x14+x4+x15>=17
22：00-6：00 X3+x14>=12
X3+x14+x10>=12		 X2+x13+x9>=12 X1+x12+x8>=12
X12+x8>=12		 X11+x7>=12
X7>=12		 X6>=12
X1+x6>=12		 X5+x10>=12
X5+x11>=12		 X4+x15>=12

因此，我们可以建立下面的线性规划模型： \[ \begin{aligned} & \min z=x_1+x_2+x_3+x_4+x_5+x_6+x_7+x_8+x_9+x_{10}+x_{11}+x_{12}=x_{13}+x_{14}+x_{15} \\ & x_4+x_5+x_{11}+x_{15} \geqslant 18 \\ & x_1+x_5+x_{11}+x_{12} \geqslant 20 \\ & x_4+x_5+x_6+x_7+x_{11}+x_{15} \geqslant 20 \\ & x_1+x_5+x_7+x_8+x_{11}+x_{12} \geqslant 19 \\ & x_3+x_4+x_6+x_{10}+x_{14}+x_{15} \geqslant 20 \\ & x_9+x_{10}+x_{13}+x_{14} \geqslant 20 \\ & x_6+x_{10}+x_{14}+x_{15} \geqslant 19 \\ & x_6+x_7+x_{11}+x_{15} \geqslant 17 \\ & x_7+x_8 \geqslant 18 \\ & x_8+x_9 \geqslant 20 \\ & x_9+x_{10} \geqslant 19 \\ & x_6+x_{10} \geqslant 17 \\ & x_1+x_2+x_6+x_7 \geqslant 18 \\ & x_2+x_3+x_7+x_8 \geqslant 20 \\ & x_1+x_2+x_{12}+x_{13} \geqslant 20 \\ & x_2+x_3+x_{13}+x_{14} \geqslant 19 \\ & x_3+x_4+x_{14}+x_{15} \geqslant 17 \\ & x_8+x_{12} \geqslant 12 \\ & x_5+x_{11} \geqslant 12 \\ & x_3+x_{14} \geqslant 12 \\ & x_2+x_9+x_{13} \geqslant 12 \\ & x_5+x_{10} \geqslant 12 \\ & x_4+x_{15} \geqslant 12 \\ & x_6 \geqslant 12 \\ & x_7 \geqslant 12 \\ & x_j \geqslant 0, j=1,2 \ldots 15 \end{aligned} \] 使用matlab编程求解，结果如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
Solution found:
Objective function value: 84
Variables:105
	 7
     7
     0
    12
     7
    12
    12
     6
    14
     5
     5
     6
     0
    12
     0

五、结果分析

5.1 方案一的结果

根据前文所述，方案一的结果如下： \[ \begin{array}{|l|l|l|l|l|l|l|l|} \hline & \text { 星期一 } & \text { 星期二 } & \text { 星期三 } & \text { 星期四 } & \text { 星期五 } & \text { 星期六 } & \text { 星期天 } \\ \hline \text { 第一班 } & 12 & 12 & 12 & 12 & 12 & 12 & 12 \\ \hline \text { 第二班 } & 12 & 12 & 12 & 12 & 12 & 12 & 12 \\ \hline \text { 第三班 } & 12 & 12 & 12 & 12 & 12 & 12 & 12 \\ \hline \text { 第四班 } & 12 & 12 & 12 & 12 & 12 & 12 & 12 \\ \hline \text { 第五班 } & 12 & 12 & 12 & 12 & 12 & 12 & 12 \\ \hline \end{array} \]

5.2 方案二的结果

根据前文所述，方案二的结果如下： \[ \begin{array}{|l|c|c|c|c|c|c|c|} \hline & \text { 星期一 } & \text { 星期二 } & \text { 星期三 } & \text { 星期四 } & \text { 星期五 } & \text { 星期六 } & \text { 星期天 } \\ \hline \text { 第一班 } & 12 & 12+12 & 12+7 & 8+13 & 12 & 12 & 12 \\ \hline \text { 第二班 } & 12 & 12+12 & 12+12 & 12+7 & 8 & 12 & 13 \\ \hline \text { 第三班 } & 13 & 12+12 & 12+12 & 12+12 & 12 & 8 & 7 \\ \hline \text { 第四班 } & 7 & 8+13 & 12+12 & 12+12 & 12 & 12 & 12 \\ \hline \text { 第五班 } & 12 & 12+7 & 8+13 & 12+12 & 12 & 12 & 12 \\ \hline \end{array} \]

5.3 方案三的结果

\[ \begin{array}{|c|c|c|c|c|c|c|c|} \hline & \text { 星期一 } & \text { 星期二 } & \text { 星期三 } & \text { 星期四 } & \text { 星期五 } & \text { 星期六 } & \text { 星期天 } \\ \hline \text { 第一班 } & 12+0 & 0+12+5 & 7+0+14 & 6+6 & 12 & 7+12 & 7+5 \\ \hline \text { 第二班 } & 7+5 & 12+0+12 & 0+12+5 & 0+14 & 6 & 7+12 & 7+6 \\ \hline \text { 第三班 } & 7+6 & 7+5+12 & 12+0+12 & 12+5 & 14 & 0+6 & 7+0 \\ \hline \text { 第四班 } & 7+0 & 7+6+6 & 7+5+12 & 0+12 & 5 & 12+14 & 0+12 \\ \hline \text { 第五班 } & 0+12 & 7+0+14 & 7+6+6 & 5+12 & 12 & 7+5 & 12+0 \\ \hline \end{array} \]

作为奖励，规定放弃周末休息的护士，其工资和奖金总额比其他 护士增加 a%，我们有： \[ 33 \times (1 + a\%) + 72 \leq 112 \] 解得： \[ a \leq 21.2121\% \] 结果表明：放弃周末休息的护士，其工资和奖金总额比其他护士增加的值小于等于21.21%时，方案三较方案二更经济。

六、附录

6.1 问题一的matlab代码：

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
% Objective function coefficients
f = [1 1 1 1 1 1 1];

% Inequality constraints
A = [
    -1  0  0  0  0  0 -1;
    -1 -1  0  0  0  0  0;
     0 -1 -1  0  0  0  0;
     0  0 -1 -1  0  0  0;
     0  0  0 -1 -1  0  0;
     0  0  0  0 -1 -1  0;
     0  0  0  0  0 -1 -1;
];
b = [-20; -20; -20; -20; -20; -20; -20];

% Bounds
lb = [12; 12; 12; 12; 12; 12; 12];
ub = [];

% Integer constraints (all variables are integers)
intcon = 1:7;

% Options (optional)
options = optimoptions('intlinprog','Display','off');

% Solve the problem
[x,fval,exitflag,output] = intlinprog(f,intcon,A,b,[],[],lb,ub,options);

% Display the results
if exitflag == 1
    disp('Solution found:');
    disp(['Objective function value: ', num2str(fval)]);
    disp('Variables:');
    disp(x);
else
    disp('No feasible solution found.');
end

6.2 方案二的matlab代码：

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
A = [
    0, 0, 0, 0, 0, 1, 0, 0, 0, 1;
    0, 0, 0, 0, 0, 1, 1, 0, 0, 0;
    0, 0, 1, 1, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 1, 1, 0, 0, 0, 0, 0;
    1, 0, 0, 0, 1, 0, 0, 0, 0, 0;
    1, 1, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 1, 1, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 1, 1, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 1, 1, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 1, 1;
    1, 1, 0, 0, 0, 1, 0, 0, 0, 1;
    1, 0, 0, 0, 1, 0, 0, 0, 1, 1;
    0, 0, 0, 1, 1, 0, 0, 1, 1, 0;
    0, 0, 0, 1, 0, 0, 0, 1, 0, 0;
    0, 0, 1, 0, 0, 0, 1, 0, 0, 0;
];
b = [18, 20, 20, 19, 17, 18, 20, 20, 19, 17, 20, 20, 20, 12, 12]';
lb = [12, 12, 0, 0, 12, 12, 0, 0, 12, 12]';
f = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]';

% Solve the problem
[x,fval] = linprog(f,-A,-b,[],[],lb);

disp('Solution found:');
disp(['Objective function value: ', num2str(fval)]);
disp('Variables:');
disp(x)

8.3 方案三的matlab代码

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
% 定义约束矩阵 A 和向量 b
A = [
    -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1;
    -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, -1, 0, 0, 0;
    0, 0, 0, -1, -1, -1, -1, 0, 0, 0, -1, 0, 0, 0, -1;
    -1, 0, 0, 0, -1, 0, -1, -1, 0, 0, -1, -1, 0, 0, 0;
    0, 0, -1, -1, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, -1;
    0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 0, 0, -1, -1, 0;
    0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, -1;
    0, 0, 0, 0, 0, -1, -1, 0, 0, 0, -1, 0, 0, 0, -1;
    0, 0, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0;
    -1, -1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0;
    0, -1, -1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0;
    -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 0, 0;
    0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 0;
    0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1;
    0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0;
    0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0;
    0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0;
    0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0;
    0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1;
    0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0;
];

b = -[
    18, 20, 20, 19, 20, 20, 19, 17, 18, 20, 19, 17, 18, 20, 20,
    19, 17, 12, 12, 12, 12, 12, 12, 12, 12
]';

% 定义目标函数系数向量 f
f = ones(15, 1);

% 定义变量的下界 lb
lb = zeros(15, 1); % 因为所有的 x_j >= 0
lb(6) = 12; % x_6 >= 12
lb(7) = 12; % x_7 >= 12

% Options (optional)
options = optimoptions('intlinprog','Display','off');

% Solve the problem
[x,fval,exitflag,output] = intlinprog(f,intcon,A,b,[],[],lb,ub,options);

disp('Solution found:');
disp(['Objective function value: ', num2str(fval)]);
disp('Variables:');
disp(x)

2021级最优化作业题

这里选择问题E,结合具体案例对最速下降法、牛顿法、 DFP算法三个算法进行编程实现。

一、案例

求函数\(f(x) = 2x_1^2+x_2^2\)的极小点。设初始点为\(x_0=(1, 1)\)，\(\epsilon = \cfrac{1}{10}\)。

二、最速下降法

2.1 原理

最速下降法（又称梯度下降法）是一种优化算法，它通过在每一步沿着目标函数梯度的负方向更新参数来寻找函数的极小点。对于函数 $ f() = 2x_1^2 + x_2^2 $，梯度向量可以通过对每个变量求偏导来计算：

\[ \nabla f(\mathbf{x}) = \left( \frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2} \right)^T \]

对于 $ f() = 2x_1^2 + x_2^2 $，偏导数是：

\[ \frac{\partial f}{\partial x_1} = 4x_1, \quad \frac{\partial f}{\partial x_2} = 2x_2 \]

因此，梯度向量为 $ f() = (4x_1, 2x_2)^T $。

在最速下降法中，我们需要迭代以下步骤来更新 x ：

  1. 计算当前点$ $的梯度 $ f() $。
  2. 确定步长$ $，这通常通过线搜索来完成，以满足一定的准则，例如Wolfe条件。
  3. 更新$ \(为\) - f() $。
  4. 重复这个过程，直到$ |f()| < \(，在这里\) = $。

我们可以开始实现这个过程，设初始点为$ _0 = (1, 1)^T $。我们需要一个合适的方法来选择步长 $ \(，一个简单的选择是使用固定的步长，例如\) = 0.1 $，或者实现一个简单的回溯线搜索。在这里，我们将使用一个较小的固定步长并观察收敛情况。

2.2 解答

使用python编程如下：

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
import numpy as np

# 定义函数f(x)
def f(x):
    return 2*x[0]**2 + x[1]**2

# 定义梯度计算函数
def grad_f(x):
    return np.array([4*x[0], 2*x[1]])

# 最速下降法
def gradient_descent(initial_x, epsilon, step_size):
    x = initial_x
    history = [x] # 用于存储迭代过程中的x值
    while np.linalg.norm(grad_f(x)) >= epsilon:
        # 更新x
        x = x - step_size * grad_f(x)
        history.append(x)
    return x, history

# 初始点
initial_x = np.array([1, 1])
# 收敛阈值
epsilon = 1/10
# 步长
step_size = 0.1

# 执行最速下降法
optimal_x, history = gradient_descent(initial_x, epsilon, step_size)

print("极小点:", optimal_x)
print("迭代次数:", len(history))

运行上述程序，得到如下输出：

1
2
极小点: [0.00078364 0.04398047]
迭代次数: 15

使用最速下降法，我们找到了函数$ f(x) = 2x_1^2 + x_2^2 \(的极小点大约在\) (0.0008, 0.0440) \(处。这个过程经历了 15 次迭代才满足梯度的模长小于\) = $的收敛条件。

由于这个函数是一个凸函数，并且有一个明显的全局最小值在 (0, 0) ，我们可以推断该方法找到的极小点非常接近真实的全局最小值，即使是在使用了固定步长的简化情况下。

三、牛顿法

3.1 原理

牛顿法是一种寻找函数零点或极值点的迭代算法。对于寻找函数极小点的问题，牛顿法迭代公式如下： \[ \mathbf{x}_{n+1} = \mathbf{x}_n - H^{-1}(\mathbf{x}_n) \nabla f(\mathbf{x}_n) \]

其中，$ _n $ 是当前迭代值，$ H(_n) $ 是函数 $ f $ 在 $ _n $ 处的Hessian矩阵，$ f(_n) $ 是函数 $ f $ 的梯度。

对于二次函数 $ f(x) = 2x_1^2 + x_2^2 $，Hessian矩阵是一个常数，因为二次项的二阶导数是常数。对于这个函数，Hessian矩阵是：

\[ H(\mathbf{x}) = \begin{bmatrix} \cfrac{\partial^2 f}{\partial x_1^2} & \cfrac{\partial^2 f}{\partial x_1 \partial x_2} \\ \cfrac{\partial^2 f}{\partial x_2 \partial x_1} & \cfrac{\partial^2 f}{\partial x_2^2} \end{bmatrix} = \begin{bmatrix} 4 & 0 \\ 0 & 2 \end{bmatrix} \]

由于Hessian矩阵是对角矩阵且不依赖于 $ $，所以其逆矩阵也是容易计算的。在这个特殊情况下，牛顿法会在一步之内直接到达极小点，因为我们是在用二次函数的牛顿法，它直接找到了函数的极小值。

3.2 解答

使用python编程如下：

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

# 定义梯度计算函数
def grad_f(x):
    return np.array([4*x[0], 2*x[1]])

# 牛顿法求解极小点
def newton_method(initial_x):
    # Hessian矩阵的逆矩阵，对于给定函数是常数
    H_inv = np.linalg.inv(np.array([[4, 0], [0, 2]]))
    # 进行一次牛顿迭代
    x_newton = initial_x - H_inv.dot(grad_f(initial_x))
    return x_newton

# 初始点
initial_x = np.array([1, 1])

# 执行牛顿法
optimal_x = newton_method(initial_x)

print("极小点:", optimal_x)

程序输出结果如下：

1
极小点: [0. 0.]

用牛顿法，我们立即找到了函数 $ f(x) = 2x_1^2 + x_2^2 $ 的极小点，位于 $ (0, 0) $。这是因为牛顿法对于二次函数是非常有效的，特别是当Hessian矩阵是常数时，它可以直接在一次迭代中找到极值点。在这种情况下，牛顿法展示了其对于二次问题的解析求解能力。

四、DFP算法

4.1 原理

DFP（Davidon-Fletcher-Powell）算法是一种求解无约束优化问题的迭代算法，属于拟牛顿方法的一种。DFP算法通过不断更新逆Hessian矩阵的估计来寻找函数的极小点，其基本步骤如下：

  1. 选择一个初始点 $ _0 $ 和初始逆Hessian估计 $ H_0 $，通常 $ H_0 $ 可以是单位矩阵。
  2. 计算梯度 $ _0 = f(_0) $。
  3. 对于 $ k = 0, 1, 2, $ 直到收敛：
    • 计算搜索方向 $ _k = -H_k _k $。
    • 用线搜索找到步长 $ _k $，使得 $ f(_k + _k _k) $ 相对 $ _k $ 最小化。
    • 更新 $ _{k+1} = _k + _k _k $。
    • 计算新梯度 $ {k+1} = f({k+1}) $。
    • 如果 $ |_{k+1}| $ 足够小，则停止迭代。
    • 计算 $ k = {k+1} - _k $ 和 $ k = {k+1} - _k $。
    • 更新逆Hessian估计 $ H_{k+1} $。

函数 $ f() = 2x_1^2 + x_2^2 $ 的梯度和Hessian矩阵都是很简单的。我们可以使用固定的步长，因为这是一个凸二次函数，这样可以简化线搜索过程。

4.2 解答

使用python编程如下：

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
import numpy as np

# 定义函数f(x)
def f(x):
    return 2*x[0]**2 + x[1]**2

# 定义梯度计算函数
def grad_f(x):
    return np.array([4*x[0], 2*x[1]])

# DFP算法实现
def dfp_method(f, grad_f, initial_x, epsilon, max_iter=1000):
    x = initial_x
    n = len(x)
    H = np.eye(n)  # 初始化逆Hessian矩阵为单位矩阵
    g = grad_f(x)
    history = [x]  # 用于存储迭代过程中的x值

    for _ in range(max_iter):
        if np.linalg.norm(g) < epsilon:  # 检查收敛性
            break

        # 计算搜索方向
        p = -np.dot(H, g)
        
        # 这里我们使用一个简单的固定步长
        alpha = 0.1
        
        # 更新x值
        x_new = x + alpha * p
        g_new = grad_f(x_new)
        s = x_new - x
        y = g_new - g
        
        # 避免除以0，加入一个小的正数
        sy = np.dot(s, y)
        if sy == 0:
            sy += 1e-10

        # 更新逆Hessian矩阵
        Hs = np.dot(H, s)
        ys = np.dot(y, s)
        Hys = np.dot(H, y)
        sHs = np.dot(s, Hs)

        # DFP更新公式
        H = H + np.outer(s, s) / sy - np.outer(Hs, Hs) / sHs

        # 为下一次迭代做准备
        x = x_new
        g = g_new
        history.append(x)

    return x, history

# 初始点
initial_x = np.array([1, 1])
# 收敛阈值
epsilon = 1/10

# 执行DFP算法
optimal_x, history = dfp_method(f, grad_f, initial_x, epsilon)

print("极小点:", optimal_x)

结果如下：

1
极小点: [0.0014867  0.04557862]

使用DFP算法，我们找到了函数 $ f(x) = 2x_1^2 + x_2^2 $ 的极小点大约在 $ (0.0015, 0.0456) $ 处。这个结果与最速下降法和牛顿法得到的结果是一致的，都非常接近于函数的真实极小点 $ (0, 0) $。