数学模型第一次作业思路分享
清明时节雨纷纷,作业让我欲断魂。
第一个作业题 P131 Q4
题目
一家保姆服务公司专门向顾主提供保姆服务。
- 根据估计,下一年的需求是:
季节 | 春季 | 夏季 | 秋季 | 冬季 |
---|---|---|---|---|
人日 | 6000 | 7500 | 5500 | 9000 |
-
公司新招聘的保姆 必须经过5天的培训才能上岗,每个保姆每季度工作(新保姆包括培训)65 天。
-
保姆从该公司而不是从顾主那里得到报酬,每人每月工资800 元。
-
春季开始时公司拥有120 名保姆。在每个季度结束后,将有15%的保姆自动离职。
(1)如果公司不允许解雇保姆,请你为公司制定下一年的招聘计划:
- 哪些季度需求的增加不影响招聘计划?
- 可以增加多少?
(2)如果公司在每个季度结束后允许解雇保姆,请为公司制定下一年的招聘计划。
第一小问
我们先对第一问进行分析,第一问要求我们制定下一年的招聘计划,也就是我们要在不同季度分别招募的保姆人数。
我们不妨先以未知数的形式给出:
季节 | 春季 | 夏季 | 秋季 | 冬季 |
---|---|---|---|---|
招募人数 | x 1 x_1 x1 | x 2 x_2 x2 | x 3 x_3 x3 | x 4 x_4 x4 |
这样的话,每个季度的工作保姆数以及新人数如下表所示:
季节 | 春季 | 夏季 | 秋季 | 冬季 |
---|---|---|---|---|
招募人数(新员工数) | x 1 x_1 x1 | x 2 x_2 x2 | x 3 x_3 x3 | x 4 x_4 x4 |
老员工数量 | y 1 y_1 y1 | y 2 y_2 y2 | y 3 y_3 y3 | y 4 y_4 y4 |
其中:
y
1
=
120
,
y_1=120,
y1=120,
y
i
=
(
y
i
−
1
+
x
i
−
1
)
∗
0.85
,
i
=
2
,
3
,
4
y_i=(y_{i-1}+x_{i-1})*0.85,i=2,3,4
yi=(yi−1+xi−1)∗0.85,i=2,3,4
计算得:
y
1
=
120
;
y_1=120;
y1=120;
y
2
=
120
∗
0.85
+
0.85
∗
x
1
;
y_2=120*0.85+0.85*x_1;
y2=120∗0.85+0.85∗x1;
y
3
=
120
∗
0.8
5
2
+
0.8
5
2
∗
x
1
+
0.85
∗
x
2
y_3=120*0.85^2+0.85^2*x_1+0.85*x_2
y3=120∗0.852+0.852∗x1+0.85∗x2
y
4
=
120
∗
0.8
5
3
+
0.8
5
3
∗
x
1
+
0.8
5
2
∗
x
2
+
0.85
∗
x
3
y_4=120*0.85^3+0.85^3*x_1+0.85^2*x_2+0.85*x_3
y4=120∗0.853+0.853∗x1+0.852∗x2+0.85∗x3
用矩阵表示上式为:
(
y
1
y
2
y
3
y
4
)
=
(
0
0
0
0
0.8
5
1
0
0
0
0.8
5
2
0.8
5
1
0
0
0.8
5
3
0.8
5
2
0.8
5
1
0
)
(
x
1
x
2
x
3
x
4
)
+
120
∗
(
0.8
5
0
0.8
5
1
0.8
5
2
0.8
5
3
)
\left( \begin{array}{c} y_1\\ y_2\\ y_3\\ y_4 \end{array} \right)= \left( \begin{array}{cccc} 0&0&0&0\\ 0.85^1&0&0&0\\ 0.85^2&0.85^1&0&0\\ 0.85^3&0.85^2&0.85^1&0 \end{array} \right) \left( \begin{array}{c} x_1\\ x_2\\ x_3\\ x_4 \end{array} \right)+120* \left( \begin{array}{c} 0.85^0\\ 0.85^1\\ 0.85^2\\ 0.85^3 \end{array} \right)
⎝⎜⎜⎛y1y2y3y4⎠⎟⎟⎞=⎝⎜⎜⎛00.8510.8520.853000.8510.8520000.8510000⎠⎟⎟⎞⎝⎜⎜⎛x1x2x3x4⎠⎟⎟⎞+120∗⎝⎜⎜⎛0.8500.8510.8520.853⎠⎟⎟⎞
而在上述表达式下每个月度的总人日数为:
季节 | 春季 | 夏季 | 秋季 | 冬季 |
---|---|---|---|---|
招募人数(新员工数) | x 1 x_1 x1 | x 2 x_2 x2 | x 3 x_3 x3 | x 4 x_4 x4 |
老员工数量 | y 1 y_1 y1 | y 2 y_2 y2 | y 3 y_3 y3 | y 4 y_4 y4 |
总人日数 | z 1 z_1 z1 | z 2 z_2 z2 | z 3 z_3 z3 | z 4 z_4 z4 |
其中:
z
i
=
60
∗
x
i
+
65
∗
y
i
z_i=60*x_i+65*y_i
zi=60∗xi+65∗yi
(新员工一个季度工作60h,老员工一个季度工作65h).
矩阵表示为:
(
z
1
z
2
z
3
z
4
)
=
65
∗
(
12
/
13
0
0
0
0.8
5
1
12
/
13
0
0
0.8
5
2
0.8
5
1
12
/
13
0
0.8
5
3
0.8
5
2
0.8
5
1
12
/
13
)
(
x
1
x
2
x
3
x
4
)
+
65
∗
120
∗
(
0.8
5
0
0.8
5
1
0.8
5
2
0.8
5
3
)
\left( \begin{array}{c} z_1\\ z_2\\ z_3\\ z_4 \end{array} \right)= 65*\left( \begin{array}{cccc} 12/13&0&0&0\\ 0.85^1&12/13&0&0\\ 0.85^2&0.85^1&12/13&0\\ 0.85^3&0.85^2&0.85^1&12/13 \end{array} \right) \left( \begin{array}{c} x_1\\ x_2\\ x_3\\ x_4 \end{array} \right)+65*120* \left( \begin{array}{c} 0.85^0\\ 0.85^1\\ 0.85^2\\ 0.85^3 \end{array} \right)
⎝⎜⎜⎛z1z2z3z4⎠⎟⎟⎞=65∗⎝⎜⎜⎛12/130.8510.8520.853012/130.8510.8520012/130.85100012/13⎠⎟⎟⎞⎝⎜⎜⎛x1x2x3x4⎠⎟⎟⎞+65∗120∗⎝⎜⎜⎛0.8500.8510.8520.853⎠⎟⎟⎞
到此为止,我们能够表示出每个季度的总人日数,下面就是给定条件寻找最优解了。
条件:
z
1
≥
6000
,
z
2
≥
7500
,
z
3
≥
5500
,
z
4
≥
9000
z_1\ge 6000,z_2\ge 7500,z_3\ge 5500,z_4\ge 9000
z1≥6000,z2≥7500,z3≥5500,z4≥9000
目标:
m
i
n
(
1
,
1
,
1
,
1
)
(
x
1
x
2
x
3
x
4
)
+
(
1
,
1
,
1
,
1
)
(
y
1
y
2
y
3
y
4
)
min{(1,1,1,1)\left( \begin{array}{c} x_1\\ x_2\\ x_3\\ x_4 \end{array} \right)+(1,1,1,1)\left(\begin{array}{c} y_1\\ y_2\\ y_3\\ y_4 \end{array} \right)}
min(1,1,1,1)⎝⎜⎜⎛x1x2x3x4⎠⎟⎟⎞+(1,1,1,1)⎝⎜⎜⎛y1y2y3y4⎠⎟⎟⎞
该目标等价于:
m
i
n
(
1
+
0.85
+
0.8
5
2
+
0.8
5
3
,
1
+
0.85
+
0.8
5
2
,
1
+
0.85
,
1
)
(
x
1
x
2
x
3
x
4
)
min{(1+0.85+0.85^2+0.85^3,1+0.85+0.85^2,1+0.85,1)\left(\begin{array}{c} x_1\\ x_2\\ x_3\\ x_4 \end{array} \right)}
min(1+0.85+0.852+0.853,1+0.85+0.852,1+0.85,1)⎝⎜⎜⎛x1x2x3x4⎠⎟⎟⎞
我们考虑用python的线性规划函数进行求解:
from scipy.optimize import linprog
import numpy as np
if __name__=='__main__':
C = np.array([1+0.85+0.85**2+0.85**3,1+0.85+0.85**2,1+0.85,1])
A = np.array([[-60,0,0,0],[-65*0.85,-60,0,0],[-65*0.85**2,-65*0.85,-60,0],[-65*0.85**3,-65*0.85**2,-65*0.85,-60]])
b = 65*120*np.array([[1],[0.85],[0.85**2],[0.85**3]])-np.array([[6000],[7500],[5500],[9000]])
X1_bounds = [0,None]
X2_bounds = [0,None]
X3_bounds = [0, None]
X4_bounds = [0, None]
res = linprog(C,A,b,bounds=(X1_bounds,X2_bounds,X3_bounds,X4_bounds))
print(res)
得到结果为:
con: array([], dtype=float64)
fun: 96.11572935215932
message: 'Optimization terminated successfully.'
nit: 7
slack: array([ 1.80000001e+03, 1.05913114e-06, 9.36625002e+02, -3.00459578e-06])
status: 0
success: True
x: array([1.94284663e-07, 1.44999998e+01, 3.69380977e-08, 5.88144791e+01])
即是说得到的结果为
X
=
(
0
14.5
0
58.8
)
X=\left(\begin{array}{c} 0\\ 14.5\\ 0\\ 58.8 \end{array} \right)
X=⎝⎜⎜⎛014.5058.8⎠⎟⎟⎞
但是这是整数规划问题,因而我们对附近整数进行计算:
- x = ( 0 , 14 , 0 , 58 ) T x=(0,14,0,58)^T x=(0,14,0,58)T
import numpy as np
if __name__=='__main__':
X=np.array([[0],[14],[0],[58]])
Z=65*np.array([[12/13,0,0,0],[0.85,12/13,0,0],[0.85**2,0.85,12/13,0],[0.85**3,0.85**2,0.85,12/13]])@X+65*120*np.array([[1],[0.85],[0.85**2],[0.85**3]])
print(Z)
结果为:
[[7800. ]
[7470. ]
[6409. ]
[8927.65]]
不满足要求,舍弃。
-
x
=
(
0
,
14
,
0
,
59
)
T
x=(0,14,0,59)^T
x=(0,14,0,59)T
结果为:
[[7800. ]
[7470. ]
[6409. ]
[8987.65]]
不满足要求,舍弃。
3)
x
=
(
0
,
15
,
0
,
58
)
T
x=(0,15,0,58)^T
x=(0,15,0,58)T
结果为:
[[7800. ]
[7530. ]
[6464.25 ]
[8974.6125]]
不满足要求,舍弃。
4)
x
=
(
0
,
15
,
0
,
59
)
T
x=(0,15,0,59)^T
x=(0,15,0,59)T
结果为:
[[7800. ]
[7530. ]
[6464.25 ]
[9034.6125]]
满足要求,保留。
对应的工资开销为:
import numpy as np
if __name__=='__main__':
X=np.array([[0],[15],[0],[59]])
Z=65*np.array([[12/13,0,0,0],[0.85,12/13,0,0],[0.85**2,0.85,12/13,0],[0.85**3,0.85**2,0.85,12/13]])@X+65*120*np.array([[1],[0.85],[0.85**2],[0.85**3]])
print(Z)
Y=np.array([[0,0,0,0],[0.85,0,0,0],[0.85**2,0.85,0,0],[0.85**3,0.85**2,0.85,0]])@X+120*np.array([[1],[0.85],[0.85**2],[0.85**3]])
Sal=3*800*np.array([1,1,1,1])@(X+Y)#一个季度三个月(hhh忘掉了
print(Sal)
结果为:
[[7800. ]
[7530. ]
[6464.25 ]
[9034.6125]]
[1151958.]
因而总工资开销为1,151,958RMB。
KEY:
但是在这种情形下,Y中一些数值是非整数值,这种情况在现实中是没有意义的,因而我们将Y向上取整(预算可以多整,少整可不行啊。)
程序更改为:
Y=np.ceil(np.array([[0,0,0,0],[0.85,0,0,0],[0.85**2,0.85,0,0],[0.85**3,0.85**2,0.85,0]])@X+120*np.array([[1],[0.85],[0.85**2],[0.85**3]]))
算出预算结果为:
[1154400.]
同时,上述计算Z的结果也表示出了那些季节的需求可增加且不影响招募。
第二小问
第二问中我们有了更多可以周旋的余地,我们终于有权利去裁员了哈哈哈(有了当老板的感觉!)。
那么仿照上面的做法,我们假设各个季度结束后我们裁员的人数为:
季节 | 春季 | 夏季 | 秋季 | 冬季 |
---|---|---|---|---|
招募人数(新员工数) | x 1 x_1 x1 | x 2 x_2 x2 | x 3 x_3 x3 | x 4 x_4 x4 |
裁员人数 | r 1 r_1 r1 | r 2 r_2 r2 | r 3 r_3 r3 | |
老员工数 | y 1 y_1 y1 | y 2 y_2 y2 | y 3 y_3 y3 | y 4 y_4 y4 |
其中:
y
1
=
120
,
y_1=120,
y1=120,
y
i
=
(
y
i
−
1
+
x
i
−
1
)
∗
0.85
−
r
i
,
i
=
2
,
3
,
4
y_i=(y_{i-1}+x_{i-1})*0.85-r_i,i=2,3,4
yi=(yi−1+xi−1)∗0.85−ri,i=2,3,4
计算得:
y
1
=
120
;
y_1=120;
y1=120;
y
2
=
120
∗
0.85
+
0.85
∗
x
1
−
r
1
;
y_2=120*0.85+0.85*x_1-r_1;
y2=120∗0.85+0.85∗x1−r1;
y
3
=
120
∗
0.8
5
2
+
0.8
5
2
∗
x
1
+
0.85
∗
x
2
−
r
2
y_3=120*0.85^2+0.85^2*x_1+0.85*x_2-r_2
y3=120∗0.852+0.852∗x1+0.85∗x2−r2
y
4
=
120
∗
0.8
5
3
+
0.8
5
3
∗
x
1
+
0.8
5
2
∗
x
2
+
0.85
∗
x
3
−
r
3
y_4=120*0.85^3+0.85^3*x_1+0.85^2*x_2+0.85*x_3-r_3
y4=120∗0.853+0.853∗x1+0.852∗x2+0.85∗x3−r3
用矩阵表示上式为:
(
y
1
y
2
y
3
y
4
)
=
(
0
0
0
0
0
0
0
0.8
5
1
0
0
0
−
1
0
0
0.8
5
2
0.8
5
1
0
0
0
−
1
0
0.8
5
3
0.8
5
2
0.8
5
1
0
0
0
−
1
)
(
x
1
x
2
x
3
x
4
r
1
r
2
r
3
)
+
120
∗
(
0.8
5
0
0.8
5
1
0.8
5
2
0.8
5
3
)
\left( \begin{array}{c} y_1\\ y_2\\ y_3\\ y_4 \end{array} \right)= \left( \begin{array}{cccc} 0&0&0&0&0&0&0\\ 0.85^1&0&0&0&-1&0&0\\ 0.85^2&0.85^1&0&0&0&-1&0\\ 0.85^3&0.85^2&0.85^1&0&0&0&-1 \end{array} \right) \left( \begin{array}{c} x_1\\ x_2\\ x_3\\ x_4\\ r_1\\ r_2\\ r_3 \end{array} \right)+120* \left( \begin{array}{c} 0.85^0\\ 0.85^1\\ 0.85^2\\ 0.85^3 \end{array} \right)
⎝⎜⎜⎛y1y2y3y4⎠⎟⎟⎞=⎝⎜⎜⎛00.8510.8520.853000.8510.8520000.85100000−10000−10000−1⎠⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛x1x2x3x4r1r2r3⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞+120∗⎝⎜⎜⎛0.8500.8510.8520.853⎠⎟⎟⎞
而总人日计算公式仍为:
Z
=
65
∗
Y
+
60
∗
X
Z=65*Y+60*X
Z=65∗Y+60∗X
条件:
z
1
≥
6000
,
z
2
≥
7500
,
z
3
≥
5500
,
z
4
≥
9000
;
z_1\ge 6000,z_2\ge 7500,z_3\ge 5500,z_4\ge 9000;
z1≥6000,z2≥7500,z3≥5500,z4≥9000;
r
1
≥
0
,
r
2
≥
0
,
r
3
≥
0.
r_1\ge0,r_2\ge0,r_3\ge0.
r1≥0,r2≥0,r3≥0.
目标:
m
i
n
(
1
,
1
,
1
,
1
,
0
,
0
,
0
)
(
x
1
x
2
x
3
x
4
r
1
r
2
r
3
)
+
(
1
,
1
,
1
,
1
)
(
y
1
y
2
y
3
y
4
)
min{(1,1,1,1,0,0,0)\left( \begin{array}{c} x_1\\ x_2\\ x_3\\ x_4 \\ r_1\\ r_2\\ r_3 \end{array} \right)+(1,1,1,1)\left(\begin{array}{c} y_1\\ y_2\\ y_3\\ y_4 \end{array} \right)}
min(1,1,1,1,0,0,0)⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛x1x2x3x4r1r2r3⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞+(1,1,1,1)⎝⎜⎜⎛y1y2y3y4⎠⎟⎟⎞
我们代入程度:
from scipy.optimize import linprog
import numpy as np
if __name__=='__main__':
#X=np.array([[0],[0],[0],[0],[0],[0],[0]])
A=np.array([[0,0,0,0,0,0,0],[0.85,0,0,0,-1,0,0],[0.85**2,0.85,0,0,0,-1,0],[0.85**3,0.85**2,0.85,0,0,0,-1]])
b=np.array([[1],[0.85],[0.85**2],[0.85**3]])
#Y=A@X+120*b
B=np.array([[1,0,0,0,0,0,0],[0,1,0,0,0,0,0],[0,0,1,0,0,0,0],[0,0,0,1,0,0,0]])
#Z=60*B@X+65*Y
C=np.array([1,1,1,1,0,0,0])+np.array([1,1,1,1])@A
#条件
X1_bounds = [0, None]
X2_bounds = [0, None]
X3_bounds = [0, None]
X4_bounds = [0, None]
r1_bounds = [0, None]
r2_bounds = [0, None]
r3_bounds = [0, None]
res = linprog(C, -60*B-65*A, 65*120*b-np.array([[6000],[7500],[5500],[9000]]), bounds=(X1_bounds, X2_bounds, X3_bounds, X4_bounds, r1_bounds, r2_bounds, r3_bounds))
print(res)
得到:
con: array([], dtype=float64)
fun: 81.70611377752874
message: 'Optimization terminated successfully.'
nit: 6
slack: array([ 1.80000000e+03, -4.55738700e-08, 3.51953986e-08, -2.96003236e-07])
status: 0
success: True
x: array([2.91460980e-10, 1.45000000e+01, 6.95358956e-09, 5.88144792e+01,
6.65316200e-09, 1.44096154e+01, 1.77970082e-09])
即是说:
x
1
=
0
,
x
2
=
14.5
,
x
3
=
0
,
x
4
=
58.8
x_1=0,x_2=14.5,x_3=0,x_4=58.8
x1=0,x2=14.5,x3=0,x4=58.8
r
1
=
0
,
r
2
=
14.4
,
r
3
=
0
r_1=0,r_2=14.4,r_3=0
r1=0,r2=14.4,r3=0
经过程序验证:
import numpy as np
if __name__=='__main__':
X=np.array([[0],[15],[0],[59],[0],[14],[0]])
A=np.array([[0,0,0,0,0,0,0],[0.85,0,0,0,-1,0,0],[0.85**2,0.85,0,0,0,-1,0],[0.85**3,0.85**2,0.85,0,0,0,-1]])
b=np.array([[1],[0.85],[0.85**2],[0.85**3]])
Y=A@X+120*b
B=np.array([[1,0,0,0,0,0,0],[0,1,0,0,0,0,0],[0,0,1,0,0,0,0],[0,0,0,1,0,0,0]])
Z=60*B@X+65*Y
C=np.array([1,1,1,1,0,0,0])+np.array([1,1,1,1])@A
print(Z)
当:
x
1
=
0
,
x
2
=
15
,
x
3
=
0
,
x
4
=
59
x_1=0,x_2=15,x_3=0,x_4=59
x1=0,x2=15,x3=0,x4=59
r
1
=
0
,
r
2
=
14
,
r
3
=
0
r_1=0,r_2=14,r_3=0
r1=0,r2=14,r3=0
满足目标条件。
此时的总工资花销为:1,118,358RMB
[[7800. ]
[7530. ]
[5554.25 ]
[9034.6125]]
[1118358.]
同样的,我们考虑取整问题:
Y=np.ceil(A@X+120*b)
得到预算结果为:
[1120800.]
至此,第四题分析结束。
第二个作业题 P131 Q5
在甲乙双方的一场战争中,一部分甲方部队被乙方部队包围长达 4 个月。由于乙方封锁了所有水陆交通通道,被包围的甲方部队只能依靠空中交通维持供给。
运送4个月的供给分别需要 2 次,3 次,3 次,4次飞行,每次飞行编队由 50 架飞机组成(每架飞机需要3名飞行员),可以运送10万吨物资。
每架飞机每个月只能飞行一次,每名飞行员每个月也只能飞行一次。在执行完运输任务后的返回途中有 20%的飞机会被乙方部队击落,相应的飞行员也因此牺牲或失踪。
在第1个月开始时,甲方拥有110 架飞机和 330 名熟练的飞行员。在每个月开始时,甲方可以招聘新飞行员和购买新飞机。新飞机必须经过一个月的检查后才可以投人使用。新飞行员必须在熟练飞行员的指导下经过一个月的训练才能投入飞行。
每名熟练飞行员可以作为教练每个月指导 20 名飞行员(包括他自己在内)进行训练。每名飞行员在完成— 个月的飞行任务后,必须有—个月的带薪假期,假期结束后才能再投人飞行。
- 已知各项费用(单位略去)如下表所示,请你为甲方安排一个飞行计划:
第一个月 | 第二个月 | 第三个月 | 第四个月 | |
---|---|---|---|---|
新飞机价格 | 200.0 | 195.0 | 190.0 | 185.0 |
闲置的熟练飞行员报酬 | 7.0 | 6.9 | 6.8 | 6.7 |
教练和新飞行员报酬(包括培训费用) | 10.0 | 9.9 | 9.8 | 9.7 |
执行飞行任务的熟练飞行员报酬 | 9.0 | 8.9 | 9.8 | 9.7 |
休假期间的熟练飞行员报酬 | 5.0 | 4.9 | 4.8 | 4.7 |
- 如果每名熟练飞行员可以作为教练每个月指导不超过 20 名飞行员(包括他自己在内)进行训练,模型和结果有哪些改变?
第一小问
我们首先来解决第一小问,首先假设如下变量:
月份 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
购入飞机数量(月初) | f 1 f_1 f1 | f 2 f_2 f2 | f 3 f_3 f3 | f 4 f_4 f4 |
招募飞行员数量(月初) | p 1 p_1 p1 | p 2 p_2 p2 | p 3 p_3 p3 | p 4 p_4 p4 |
总的熟练飞行员数量 | T 1 T_1 T1 | T 2 T_2 T2 | T 3 T_3 T3 | T 4 T_4 T4 |
总可飞行飞机数目 | F 1 F_1 F1 | F 2 F_2 F2 | F 3 F_3 F3 | F 4 F_4 F4 |
执行飞行任务的飞机数量(月中) | k 1 k_1 k1 | k 2 k_2 k2 | k 3 k_3 k3 | k 4 k_4 k4 |
执行飞行任务的熟练飞行员数量(月中) | t 1 t_1 t1 | t 2 t_2 t2 | t 3 t_3 t3 | t 4 t_4 t4 |
教练熟练飞行员数量(月中) | i 1 i_1 i1 | i 2 i_2 i2 | i 3 i_3 i3 | i 4 i_4 i4 |
闲置熟练飞行员数量(月中) | b 1 b_1 b1 | b 2 b_2 b2 | b 3 b_3 b3 | b 4 b_4 b4 |
休假熟练飞行员数量 | c 1 c_1 c1 | c 2 c_2 c2 | c 3 c_3 c3 | c 4 c_4 c4 |
被击落飞机数(月末) | L 1 L_1 L1 | L 2 L_2 L2 | L 3 L_3 L3 | L 4 L_4 L4 |
牺牲飞行员数目(月末) | l 1 l_1 l1 | l 2 l_2 l2 | l 3 l_3 l3 | l 4 l_4 l4 |
那么其中存在一些关系:
- 总的熟练飞行员人数如何变化?:
T 1 = 330 , T i = p i − 1 + T i − 1 − l i − 1 , i = 2 , 3 , 4 T_1=330,T_i=p_{i-1}+T_{i-1}-l_{i-1},i=2,3,4 T1=330,Ti=pi−1+Ti−1−li−1,i=2,3,4 - 总的可飞行飞机数量如何变化?:
F 1 = 110 , F i = f i − 1 + F i − 1 − L i − 1 , i = 2 , 3 , 4 F_1=110,F_i=f_{i-1}+F_{i-1}-L_{i-1},i=2,3,4 F1=110,Fi=fi−1+Fi−1−Li−1,i=2,3,4 - 一个教官带20个新生: 20 ∗ i i = p i 20*i_i=p_i 20∗ii=pi
- 休假的飞行员是上个月飞行但未失踪的飞行员:
c 1 = 0 , c i = t i − 1 − l i − 1 , i = 2 , 3 , 4 c_1=0,c_i=t_{i-1}-l_{i-1},i=2,3,4 c1=0,ci=ti−1−li−1,i=2,3,4 - 飞机和人数始终是3:1关系,被击沉飞机占比20%:
3 ∗ k i = t i , 3 ∗ L i = l i , 5 ∗ L i = k i 3*k_i=t_i, 3*L_i=l_i,5*L_i=k_i 3∗ki=ti,3∗Li=li,5∗Li=ki - 总的飞行员数量以及总的飞机数量受限制:
T i = t i + i i + b i + c i , F i ≥ k i , i = 1 , 2 , 3 , 4 T_i=t_i+i_i+b_i+c_i,F_i\ge k_i,i=1,2,3,4 Ti=ti+ii+bi+ci,Fi≥ki,i=1,2,3,4 - 每个月份的飞行要求:
k 1 = 100 , k 2 = 150 , k 3 = 150 , k 4 = 200 k_1=100,k_2=150,k_3=150,k_4=200 k1=100,k2=150,k3=150,k4=200
关系实在是太复杂了啊啊啊!!
所以我们来把有关系的几个变量放一块:
首先是我们的输入变量:
月份 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
购入飞机数量(月初) | f 1 f_1 f1 | f 2 f_2 f2 | f 3 f_3 f3 | f 4 f_4 f4 |
招募飞行员数量(月初) | p 1 p_1 p1 | p 2 p_2 p2 | p 3 p_3 p3 | p 4 p_4 p4 |
其次是秩为1的一组关系:
月份 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
执行飞行任务的飞机数量(月中) | k 1 k_1 k1 | k 2 k_2 k2 | k 3 k_3 k3 | k 4 k_4 k4 |
执行飞行任务的熟练飞行员数量(月中) | t 1 t_1 t1 | t 2 t_2 t2 | t 3 t_3 t3 | t 4 t_4 t4 |
被击落飞机数(月末) | L 1 L_1 L1 | L 2 L_2 L2 | L 3 L_3 L3 | L 4 L_4 L4 |
牺牲飞行员数目(月末) | l 1 l_1 l1 | l 2 l_2 l2 | l 3 l_3 l3 | l 4 l_4 l4 |
再其次是飞机数量变化关系:
月份 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
购入飞机数量(月初) | f 1 f_1 f1 | f 2 f_2 f2 | f 3 f_3 f3 | f 4 f_4 f4 |
总可飞行飞机数目 | F 1 F_1 F1 | F 2 F_2 F2 | F 3 F_3 F3 | F 4 F_4 F4 |
被击落飞机数(月末) | L 1 L_1 L1 | L 2 L_2 L2 | L 3 L_3 L3 | L 4 L_4 L4 |
之后是熟练飞行员的数量关系:
月份 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
招募飞行员数量(月初) | p 1 p_1 p1 | p 2 p_2 p2 | p 3 p_3 p3 | p 4 p_4 p4 |
总的熟练飞行员数量 | T 1 T_1 T1 | T 2 T_2 T2 | T 3 T_3 T3 | T 4 T_4 T4 |
牺牲飞行员数目(月末) | l 1 l_1 l1 | l 2 l_2 l2 | l 3 l_3 l3 | l 4 l_4 l4 |
之后是与不飞行的熟练飞行员数量计算有关的量:
月份 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
休假熟练飞行员数量 | c 1 c_1 c1 | c 2 c_2 c2 | c 3 c_3 c3 | c 4 c_4 c4 |
执行飞行任务的熟练飞行员数量(月中) | t 1 t_1 t1 | t 2 t_2 t2 | t 3 t_3 t3 | t 4 t_4 t4 |
教练熟练飞行员数量(月中) | i 1 i_1 i1 | i 2 i_2 i2 | i 3 i_3 i3 | i 4 i_4 i4 |
招募飞行员数量(月初) | p 1 p_1 p1 | p 2 p_2 p2 | p 3 p_3 p3 | p 4 p_4 p4 |
下面我们来用数学式子表示其中的关系:
∘ 1. O u t P u t = ( 200 , 195 , 190 , 185 ) ∗ f + ( 7.0 , 6.9 , 6.8 , 6.7 ) ∗ b + 21 ∗ ( 10.0 , 9.9 , 9.8 , 9.7 ) ∗ i + ( 9.0 , 8.9 , 9.8 , 9.7 ) ∗ t + ( 5.0 , 4.9 , 4.8 , 4.7 ) ∗ c \circ{1}.OutPut=(200,195,190,185)*f+(7.0,6.9,6.8,6.7)*b+21*(10.0,9.9,9.8,9.7)*i+(9.0,8.9,9.8,9.7)*t+(5.0,4.9,4.8,4.7)*c ∘1.OutPut=(200,195,190,185)∗f+(7.0,6.9,6.8,6.7)∗b+21∗(10.0,9.9,9.8,9.7)∗i+(9.0,8.9,9.8,9.7)∗t+(5.0,4.9,4.8,4.7)∗c
∘ 2. b = T − i − c − t \circ{2}.b=T-i-c-t ∘2.b=T−i−c−t
∘ 2 1 . b = ( − 1 / 20 0 0 0 1 − 1 / 20 0 0 1 1 − 1 / 20 0 1 1 1 − 1 / 20 ) ∗ p + ( 30 60 90 − 150 ) \circ{2_1}.b=\left(\begin{array}{cc} -1/20&0&0&0\\ 1&-1/20&0&0\\ 1&1&-1/20&0\\ 1&1&1&-1/20\\ \end{array} \right)*p+\left(\begin{array}{c} 30\\ 60\\ 90\\ -150 \end{array} \right) ∘21.b=⎝⎜⎜⎛−1/201110−1/201100−1/201000−1/20⎠⎟⎟⎞∗p+⎝⎜⎜⎛306090−150⎠⎟⎟⎞
∘ 3. T = ( 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ) ∗ ( p + T − l ) + ( 330 0 0 0 ) \circ{3}.T=\left(\begin{array}{cc} 0&0&0&0\\ 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ \end{array} \right)*(p+T-l)+\left(\begin{array}{cc} 330\\ 0\\ 0\\ 0\\ \end{array} \right) ∘3.T=⎝⎜⎜⎛0100001000010000⎠⎟⎟⎞∗(p+T−l)+⎝⎜⎜⎛330000⎠⎟⎟⎞
∘ 3 1 . T = ( 1 0 0 0 − 1 1 0 0 0 − 1 1 0 0 0 − 1 1 ) − 1 ∗ ( 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ) ∗ ( p − l ) + ( 1 0 0 0 − 1 1 0 0 0 − 1 1 0 0 0 − 1 1 ) − 1 ∗ ( 330 0 0 0 ) \circ{3_1}. T=\left(\begin{array}{cc} 1&0&0&0\\ -1&1&0&0\\ 0&-1&1&0\\ 0&0&-1&1\\ \end{array} \right)^{-1}*\left(\begin{array}{cc} 0&0&0&0\\ 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ \end{array} \right)*(p-l)+\left(\begin{array}{cc} 1&0&0&0\\ -1&1&0&0\\ 0&-1&1&0\\ 0&0&-1&1\\ \end{array} \right)^{-1}*\left(\begin{array}{cc} 330\\ 0\\ 0\\ 0\\ \end{array} \right) ∘31.T=⎝⎜⎜⎛1−10001−10001−10001⎠⎟⎟⎞−1∗⎝⎜⎜⎛0100001000010000⎠⎟⎟⎞∗(p−l)+⎝⎜⎜⎛1−10001−10001−10001⎠⎟⎟⎞−1∗⎝⎜⎜⎛330000⎠⎟⎟⎞
∘ 3 2 . T = ( 0 0 0 0 1 0 0 0 1 1 0 0 1 1 1 0 ) ∗ p + ( 330 270 180 90 ) \circ{3_2}.T=\left(\begin{array}{cc} 0&0&0&0\\ 1&0&0&0\\ 1&1&0&0\\ 1&1&1&0\\ \end{array} \right)*p+\left(\begin{array}{c} 330\\ 270\\ 180\\ 90 \end{array} \right) ∘32.T=⎝⎜⎜⎛0111001100010000⎠⎟⎟⎞∗p+⎝⎜⎜⎛33027018090⎠⎟⎟⎞
∘ 4. i = 1 / 20 ∗ p \circ{4}.i=1/20*p ∘4.i=1/20∗p
∘ 5. c = ( 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ) ∗ ( t − l ) \circ{5}.c=\left(\begin{array}{cc} 0&0&0&0\\ 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ \end{array} \right)*(t-l) ∘5.c=⎝⎜⎜⎛0100001000010000⎠⎟⎟⎞∗(t−l)
∘ 6. t = 3 ∗ k \circ{6}.t=3*k ∘6.t=3∗k
∘ 7. l = 3 / 5 ∗ k \circ{7}.l=3/5*k ∘7.l=3/5∗k
∘ 1 1 . O u t P u t = ( 200 , 195 , 190 , 185 ) ∗ f + ( 7.0 , 6.9 , 6.8 , 6.7 ) ∗ ( ( − 1 / 20 0 0 0 1 − 1 / 20 0 0 1 1 − 1 / 20 0 1 1 1 − 1 / 20 ) ∗ p + ( 330 330 330 330 ) − ( 3 0 0 0 3 3 0 0 3 / 5 3 3 0 3 / 5 3 / 5 3 3 ) ∗ k ) + 21 / 20 ∗ ( 10.0 , 9.9 , 9.8 , 9.7 ) ∗ p + ( 9.0 , 8.9 , 9.8 , 9.7 ) ∗ 3 ∗ k + ( 5.0 , 4.9 , 4.8 , 4.7 ) ∗ ( 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ) ∗ 12 / 5 ∗ k \circ{1_1}.OutPut=(200,195,190,185)*f+(7.0,6.9,6.8,6.7)*(\left(\begin{array}{cc} -1/20&0&0&0\\ 1&-1/20&0&0\\ 1&1&-1/20&0\\ 1&1&1&-1/20\\ \end{array} \right)*p+\left(\begin{array}{cc} 330\\ 330\\ 330\\ 330\\ \end{array} \right)-\left(\begin{array}{cc} 3&0&0&0\\ 3&3&0&0\\ 3/5&3&3&0\\ 3/5&3/5&3&3\\ \end{array} \right)*k)+21/20*(10.0,9.9,9.8,9.7)*p+(9.0,8.9,9.8,9.7)*3*k+(5.0,4.9,4.8,4.7)*\left(\begin{array}{cc} 0&0&0&0\\ 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ \end{array} \right)*12/5*k ∘11.OutPut=(200,195,190,185)∗f+(7.0,6.9,6.8,6.7)∗(⎝⎜⎜⎛−1/201110−1/201100−1/201000−1/20⎠⎟⎟⎞∗p+⎝⎜⎜⎛330330330330⎠⎟⎟⎞−⎝⎜⎜⎛333/53/50333/500330003⎠⎟⎟⎞∗k)+21/20∗(10.0,9.9,9.8,9.7)∗p+(9.0,8.9,9.8,9.7)∗3∗k+(5.0,4.9,4.8,4.7)∗⎝⎜⎜⎛0100001000010000⎠⎟⎟⎞∗12/5∗k
∘ 8. k = ( 100 , 150 , 150 , 200 ) T \circ{8}.k=(100,150,150,200)^T ∘8.k=(100,150,150,200)T
∘ 9. F = ( 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ) ∗ ( f + F − L ) + ( 110 0 0 0 ) \circ{9}.F=\left(\begin{array}{cc} 0&0&0&0\\ 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ \end{array} \right)*(f+F-L)+\left(\begin{array}{cc} 110\\ 0\\ 0\\ 0\\ \end{array} \right) ∘9.F=⎝⎜⎜⎛0100001000010000⎠⎟⎟⎞∗(f+F−L)+⎝⎜⎜⎛110000⎠⎟⎟⎞
∘ 9 1 . F = ( 0 0 0 0 1 0 0 0 1 1 0 0 1 1 1 0 ) ∗ f + ( 110 90 60 30 ) \circ{9_1}.F=\left(\begin{array}{cc} 0&0&0&0\\ 1&0&0&0\\ 1&1&0&0\\ 1&1&1&0\\ \end{array} \right)*f+\left(\begin{array}{cc} 110\\ 90\\ 60\\ 30\\ \end{array} \right) ∘91.F=⎝⎜⎜⎛0111001100010000⎠⎟⎟⎞∗f+⎝⎜⎜⎛110906030⎠⎟⎟⎞
目标:(最小化)
∘ 1 2 . O u t P u t = ( 200 , 195 , 190 , 185 ) ∗ f + ( 30.55 , 23.55 , 16.65 , 9.85 ) ∗ p + 9030 \circ{1_2}.OutPut=(200,195,190,185)*f+(30.55,23.55,16.65,9.85)*p+9030 ∘12.OutPut=(200,195,190,185)∗f+(30.55,23.55,16.65,9.85)∗p+9030
条件:
b ≥ 0 , F ≥ k b\ge0,F\ge k b≥0,F≥k
from scipy.optimize import linprog
import numpy as np
if __name__ == '__main__':
C = np.array([200,195,190,185,30.55,23.55,16.65,9.85])
A = np.array([[0,0,0,0,0,0,0,0], [-1,0,0,0,0,0,0,0], [-1,-1,0,0,0,0,0,0],
[-1,-1,-1,0,0,0,0,0],[0,0,0,0,1/20,0,0,0], [0,0,0,0,-1,1/20,0,0], [0,0,0,0,-1,-1,1/20,0],
[0,0,0,0,-1,-1,-1,1/20]])
b =np.array([[10], [-60], [-90], [-170],[30], [60], [90], [-150]])
f1_bounds = [0, None]
f2_bounds = [0, None]
f3_bounds = [0, None]
f4_bounds = [0, None]
p1_bounds = [0, None]
p2_bounds = [0, None]
p3_bounds = [0, None]
p4_bounds = [0, None]
res = linprog(C, A, b, bounds=(f1_bounds, f2_bounds, f3_bounds, f4_bounds,p1_bounds, p2_bounds, p3_bounds, p4_bounds))
print(res)
得到结果如下:
con: array([], dtype=float64)
fun: 35547.49572128034
message: 'Optimization terminated successfully.'
nit: 6
slack: array([ 1.00000000e+01, 9.79953937e-06, 1.25172376e-06, -2.05973778e-05,
2.99999999e+01, 6.00000012e+01, 8.25000038e+01, -2.75783846e-05])
status: 0
success: True
x: array([6.00000098e+01, 2.99999915e+01, 7.99999782e+01, 0.00000000e+00,
1.22194725e-06, 1.05037096e-06, 1.49999970e+02, 1.35485963e-06])
月份 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
购入飞机数量(月初) | 60 | 30 | 80 | 0 |
招募飞行员数量(月初) | 0 | 0 | 150 | 0 |
总的熟练飞行员数量 | 330 | T 2 T_2 T2 | T 3 T_3 T3 | T 4 T_4 T4 |
总可飞行飞机数目 | 110 | F 2 F_2 F2 | F 3 F_3 F3 | F 4 F_4 F4 |
执行飞行任务的飞机数量(月中) | 100 | 150 | 150 | 200 |
执行飞行任务的熟练飞行员数量(月中) | 300 | 450 | 450 | 600 |
教练熟练飞行员数量(月中) | 0 | 0 | 8 | 0 |
闲置熟练飞行员数量(月中) | b 1 b_1 b1 | b 2 b_2 b2 | b 3 b_3 b3 | b 4 b_4 b4 |
休假熟练飞行员数量 | 0 | 240 | c 3 c_3 c3 | c 4 c_4 c4 |
被击落飞机数(月末) | 20 | 30 | 30 | 40 |
牺牲飞行员数目(月末) | 60 | 90 | 90 | 120 |
题目:《数学模型》(姜启源,第三版)
文稿:锦帆远航
代码:锦帆远航
图片:百度图片(岁城琉心·鬼蜮咒原)
您的三连支持是航航坚持写下去的不懈动力!