小白自学数学建模丨M001 - 线性规划（LP）

M001 - 线性规划 (LP)

线性规划的实例与定义

例1.1 某机床厂生产甲、乙两种机床，每台销售后的利润分别为4千元和3千元。生产甲机床需要使用A,B机器加工，加工时间分别为每台2小时和1小时;生产乙机床需要A,B,C三种机器加工，加工时间为每台各一小时。若每天可用于加工的机器时数分别为A机器10小时、B机器8小时和C机器7小时，问该厂应生产甲、乙机床各几台，才能使总利润最大？
上述问题的数学模型: 设该厂生产\(x_1\)台甲机床和\(x_2\)乙机床时总利润z最大，则\(x_1, x_2\)应满足:

\[（1.1） \max \ \ z=4{x_1}+3{x_2} \]

\[（1.2） \begin{cases} 2{x_1}+{x_2} \le 10\\ {x_1}+{x_2} \le 8\\ {x_2} \le 7\\ {x_1},{x_2} \\ \end{cases} \]

变量\(x_1,x_2\)称为决策变量，(1.1)式被称为问题的目标函数，(1.2)中的几个不等式是问题的约束条件，记为s.t.(即subject to)。

线性规划问题的解的概念

\[\max_x c^Tx \]

\[s.t. \begin{cases} Ax \le b\\ Aeq \cdot x = beq\\ lb \le x \le ub \end{cases} \]

其中c和x为n维列向量
A、Aeq为适当维数的矩阵
b、beq为适当维数的列向量

一般线性规划问题的(数学)标准型为

\[（1.3） \max \ \ z=\sum_{j=0}^{n} c_j x_j \]

\[（1.4）s.t. \begin{cases} \displaystyle\sum_{j=1}^{n} a_{ij} x_j = b_j \quad i=1, 2, \cdots , m\\ x_j \ge 0 \quad j = 1, 2, \cdots, n \end{cases} \]

可行解: 满足约束条件(1.4)的解\(x=[x_1, L, x_n]^T\)，称为线性规划问题的可行解，而使目标函数(1.3)达到最大值的可行解叫最优解。
可行域: 所有可行解构成的集合称为问题的可行域，记为R。

线性规划的Matlab标准形式及软件求解

Matlab规定线性规划的标准形式为

\[\max_x c^Tx \]

\[s.t. \begin{cases} Ax \le b\\ Aeq \cdot x = beq\\ lb \le x \le ub \end{cases} \]

其中c、x、b、beq、lb、ub为列向量
f称为价值向量
b称为资源向量
A、Aeq为矩阵

Matlab中求解线性规划的命令如下

[x,fval] = linprog(c,A,b)
[x,fval] = linprog(c,A,b,Aeq,beq)
[x,fval] = linprog(c,A,b,Aeq,beq,lb,ub)

其中
x返回的是决策向量的取值，
fval返回的是目标函数的最优解，
c为价值向量，
A，b对应的是线性不等式约束，
Aeq，beq对应的是线性等式约束，
lb和ub分别对应的是决策向量的下界向量和上界向量

例1.2 求解下列线性规划问题

\[\max \ \ z=2x_1 + 3x_2 - 5x_3 \]

\[s.t. \begin{cases} x_1 + x_2 + x_3 = 7\\ 2x_1 - 5x_2 + x_3 \ge 10\\ x_1 + 3x_2 + x_3 \le 12\\ x_1,x_2,x_3 \ge 0 \end{cases} \]

使用Matlab的解法

f = [-2; -3; 5];
a = [-2, 5, -1; 1, 3, 1];
b = [-10; 12];
Aeq = [1, 1, 1];
beq = 7; 
[x, y] = linprog(f, a, b, Aeq, beq, zeros(3 , 1));
x, y = -y

使用Python的解法

from scipy import optimize
import numpy as np
f = np.array([-2, -3, 5])
a = np.array([[-2, 5, -1], [1, 3, 1]])
b = np.array([-10, 12])
Aeq = np.array([[1, 1, 1]])
beq = np.array(7)
res = optimize.linprog(f, a, b, Aeq, beq)
print(res)

Python结果如下:

     con: array([1.80714554e-09])
     fun: -14.571428565645032
 message: 'Optimization terminated successfully.'
     nit: 5
   slack: array([-2.24602559e-10,  3.85714286e+00])
  status: 0
 success: True
       x: array([6.42857143e+00, 5.71428571e-01, 2.35900788e-10])

例1.3 求解下列线性规划问题

\[\max \ \ z=2x_1 + 3x_2 + x_3 \]

\[s.t. \begin{cases} x_1 + 2x_2 + 4x_3 = 101\\ x_1 + 4x_2 + 2x_3 \ge 8\\ 3x_1 + 2x_2 \ge 6\\ x_1,x_2,x_3 \ge 0 \end{cases} \]

使用Python的解法一

from scipy import optimize
import numpy as np
f = np.array([-2, -3, -1])
a = np.array([[-2, -4, -2], [-3, -2, 0]])
b = np.array([-8, -6])
Aeq = np.array([[1, 2, 4]])
beq = np.array(101)
res = optimize.linprog(f, a, b, Aeq, beq)
print(res)

解法一结果如下所示 :

     con: array([1.46057431e-07])
     fun: -201.99999956122394
 message: 'Optimization terminated successfully.'
     nit: 7
   slack: array([193.99999969, 296.99999902])
  status: 0
 success: True
       x: array([1.01000000e+02, 1.24887119e-07, 3.11058324e-09])

使用Python的解法二

import pulp as pp
z = [2, 3, 1]
a = [[1, 4, 2], [3, 2, 0]]
b = [8, 6]
Aeq = [[1, 2, 4]]
beq = [101]
#确定最大化还是最小化问题
m = pp.LpProblem(sense = pp.LpMaximize)
#定义三个变量放到列表中如x1,x2,x3
x = [pp.LpVariable(f'x{i}', lowBound = 0) for i in[1, 2, 3]]
#定义一个目标函数并将目标函数加入求解的问题中
m += pp.lpDot(z, x)#LpDot用于计算点积
#设置比较条件
for i in range(len(a)):
	m += (pp.lpDot(a[i], x) >= b[i])
#设置等式条件
for i in range(len(Aeq)):
	m += (pp.lpDot(Aeq[i], x) == beq[i])
#求解问题
m.solve()
#输出结果
print(f'优化结果: {pp.value(m.objective)}')
print(f'参数取值: {[pp.value(var) for var in x]}')

解法二结果如下所示 :

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - E:\PyChram Workspace\Learning\venv\lib\site-packages\pulp\apis\..\solverdir\cbc\win\64\cbc.exe  
C:\Users\Lucifer\AppData\Local\Temp\···pulp.mps max timeMode elapsed branch printingOptions all solution 
C:\Users\Lucifer\AppData\Local\Temp\···pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 8 COLUMNS
At line 20 RHS
At line 24 BOUNDS
At line 25 ENDATA
Problem MODEL has 3 rows, 3 columns and 8 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Presolve 3 (0) rows, 3 (0) columns and 8 (0) elements
0  Obj -0 Primal inf 29.25 (3) Dual inf 5.9999997 (3)
0  Obj -0 Primal inf 29.25 (3) Dual inf 5.1666667e+10 (3)
2  Obj 202
Optimal - objective value 202
Optimal objective 202 - 2 iterations time 0.002
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

优化结果: 202.0
参数取值: [101.0, 0.0, 0.0]