用Python求解线性规划问题

线性规划简介及数学模型表示线性规划简介一个典型的线性规划问题线性规划模型的三要素线性规划模型的数学表示图解法和单纯形法图解法单纯形法使用python求解简单线性规划模型编程思路求解案例例1：使用scipy求解例2：包含非线性项的求解从整数规划到0-1规划整数规划模型0-1规划模型案例：投资的收益和风险问题描述与分析建立与简化模型

线性规划简介及数学模型表示

线性规划简介

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

而随着计算机的发展，线性规划的方法被应用于广泛的领域，已成为数学建模里最为经典，最为常用的模型之一。线性规划模型可用于求解利润最大，成本最小，路径最短等最优化问题。

一个典型的线性规划问题

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

问题分析

这个问题是一个十分典型的线性规划问题，首先对问题提取出关键信息：

决策：生产几台甲、乙机床 优化目标：总利润最大 约束：生产机床的使用时间有限

将上诉三个要素写成数学表达式，就是一个典型的线性规划模型：

规划问题的分类

• 线性规划： 在一组线性约束条件的限制下，求一线性目标函数最大或最小的问题；
• 整数规划：当约束条件加强，要求所有的自变量必须是整数时，成为整数规划（特别地，自变量只能为0或1时称为0-1规划）；
• 非线性规划：无论是约束条件还是目标函数出现非线性项，那么规划问题就变成了非线性规划；
• 多目标规划：在一组约束条件的限制下，求多个目标函数最大或最小的问题；
• 动态规划：将优化目标函数分多阶段，利用阶段间的关系逐一进行求解的方法；

应用举例：旅行商问题、车辆路径规划问题、运输问题、最短路问题、最大流问题、中国邮递员问题

线性规划模型的三要素

线性规划模型主要包括三个部分：决策变量、目标函数、约束条件

决策变量 决策变量是指问题中可以改变的量，例如生产多少货物，选择哪条路径等；线性规划的目标就是找到最优的决策变量。 在线性规划中决策变量包括实数变量，整数变量，0-1变量等。

目标函数

线性规划模型的数学表示

图解法和单纯形法

图解法

对于较为简单且只有两个决策变量的线性规划问题可以使用图解法。

单纯形法

对于决策变量比较多的线性规划模型，图解法不再适用。 单纯形法是1947 年G. B. Dantzig提出的一种十分有效的求解方法，极大地推广了线性规划的应用，直到今日也在一些线性规划的求解器中使用。

从图解法的例子中，我们可以看出，约束条件所围成的区域为一个凸多边形，当决策变量多于两个时，约束条件围成的区域为一个凸多面体，称之为可行域。其中每一个面（称之为超平面）即代表一个约束条件。

可以证明：线性规划的最优解一定在可行域的边界上

单纯形法的思路就是在可行域的一个顶点处找到一个初始可行解，判断该解是不是最优，若不是，则迭代到下一个顶点处进行重复判断。因为最优解的搜索范围从整个可行域缩小到了可行域的有限个顶点，算法的效率得到了极大的提升。

具体的找初始可行解的方法，判断解是否最优的条件，如何进行迭代这里不做详细展开，有兴趣可以查阅相关资料

此外，求解线性规划的方法还有椭球法、卡玛卡算法、内点法等。 其中内点法因为求解效率更高，在决策变量多，约束多的情况下能取得更好的效果，目前主流线性规划求解器都是使用的内点法。

使用python求解简单线性规划模型

编程思路

1. 选择适当的决策变量

在解决实际问题时，把问题归结成一个线性规划数学模型是很重要的一步，但往往也是困难的一步，模型建立得是否恰当，直接影响到求解。而选适当的决策变量，是我们建立有效模型的关键之一。

2.将求解目标简化为求一个目标函数的最大/最小值

能把要求解的问题简化为一个最值问题是能否使用线性规划模型的关键，如果这一点不能达到，之后的工作都有没有意义的。

3. 根据实际要求写出约束条件（正负性，资源约束等）

线性规划的约束条件针对不同的问题有不同的形式，总结来说有以下三种：等式约束、不等式约束、符号约束

Step1: 导入相关库

import numpy as np 
from scipy import optimize as op 

Step2: 定义决策变量

# 给出变量取值范围
x1=(0,None)  
x2=(0,None)
x3=(0,None)

Step3: 将原问题化为标准形式

注意：编程时默认为最小化目标函数，因此这里改为 ；第二个约束为大于等于约束，这里化为小于等于约束；

Step4: 定义目标函数系数和约束条件系数

c=np.array([-2,-3,5])   # 目标函数系数,3x1列向量

A_ub=np.array([[-2,5,-1],[1,3,1]]) # 不等式约束系数A，2x3维矩阵
B_ub=np.array([-10,12])  # 等式约束系数b, 2x1维列向量
A_eq=np.array([[1,1,1]])  # 等式约束系数Aeq，3x1维列向量
B_eq=np.array([7])   # 等式约束系数beq，1x1数值

Step5: 求解

res=op.linprog(c,A_ub,B_ub,A_eq,B_eq,bounds=(x1,x2,x3)) #调用函数进行求解
res

     con: array([0.])
     fun: -14.571428571428571
 message: 'Optimization terminated successfully.'
     nit: 3
   slack: array([0.        , 3.85714286])
  status: 0
 success: True
       x: array([6.42857143, 0.57142857, 0.        ])

求解案例

例1：使用scipy求解

#导入相关库
import numpy as np
from scipy import optimize as op

#定义决策变量范围
x1=(0,None)
x2=(0,None)
x3=(0,None)

#定义目标函数系数
c=np.array([2,3,1])

#定义约束条件系数
A_ub=np.array([[-1,-4,-2],[-3,-2,0]])
B_ub=np.array([-8,-6])

#求解
res=op.linprog(c,A_ub,B_ub,bounds=(x1,x2,x3))
res

     con: array([], dtype=float64)
     fun: 7.0
 message: 'Optimization terminated successfully.'
     nit: 2
   slack: array([0., 0.])
  status: 0
 success: True
       x: array([0.8, 1.8, 0. ])

例2：包含非线性项的求解

Step1:导入相关库

import numpy as np
from scipy.optimize import minimize

Step2:使用函数的形式表示目标和约束

# 定义目标函数
def objective(x):
    return x[0] ** 2 + x[1]**2 + x[2]**2 +8

# 定义约束条件
def constraint1(x):
    return x[0] ** 2 - x[1] + x[2]**2  # 不等约束

def constraint2(x):
    return -(x[0] + x[1]**2 + x[2]**2-20)  # 不等约束

def constraint3(x):
    return -x[0] - x[1]**2 + 2        # 等式约束

def constraint4(x):
    return x[1] + 2*x[2]**2 -3           # 等式约束

Step3:定义约束条件

con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
con3 = {'type': 'eq', 'fun': constraint3}
con4 = {'type': 'eq', 'fun': constraint4}

# 4个约束条件
cons = ([con1, con2, con3,con4])  

# 决策变量的符号约束
b = (0.0, None) #即决策变量的取值范围为大于等于0
bnds = (b, b ,b) 

注意：每一个约束为一个字典，其中 type 表示约束类型：ineq为大于等于，eq为等于；fun 表示约束函数表达式，即step2中的自定义函数。

Step4:求解

x0=np.array([0, 0, 0]) #定义初始值
solution = minimize(objective, x0, method='SLSQP', \
                    bounds=bnds, constraints=cons)

注意：minimize为最小化目标函数，且约束条件中默认为大于等于约束。

Step5：打印求解结果

x = solution.x

print('目标值: ' + str(objective(x)))
print('最优解为')
print('x1 = ' + str(round(x[0],2)))
print('x2 = ' + str(round(x[1],2)))
print('x3 = ' + str(round(x[2],2)))
solution

目标值: 10.651091840572583
最优解为
x1 = 0.55
x2 = 1.2
x3 = 0.95





     fun: 10.651091840572583
     jac: array([1.10433471, 2.40651834, 1.89564812])
 message: 'Optimization terminated successfully.'
    nfev: 86
     nit: 15
    njev: 15
  status: 0
 success: True
       x: array([0.55216734, 1.20325918, 0.94782404])

从整数规划到0-1规划

整数规划模型

规划中的变量（部分或全部）限制为整数时，称为整数规划。若在线性规划模型中，变量限制为整数，则称为整数线性规划。

当决策变量均为整数时，称纯整数规划；

当决策变量中部分为整数，部分为实数时，称混合整数规划；

将第一节中的线性规划图解法的例子添加整数约束，则可行域变为了多边形内的整点，如下图所示：

可以看出，可行域变成了离散的点，这也使得整数规划问题比线性规划问题要更难求解，但现实中的许多决策变量都只能取整数，因此混合整数规划问题也成为了了研究最多的线性规划问题。

注意：整数规划最优解不能按照实数最优解简单取整而获得

整数规划的两个常用求解方法：分支定界算法、割平面法

分枝定界法

step1不考虑整数约束的情况下求解得到最优解 （一般不是整数）；

step2以该解的上下整数界限建立新的约束，将原整数规划问题变为两个问题（分枝）；

step3分别对两个子问题求解（不考虑整数约束），若解刚好为整数解则结束；若不为整数解则继续进行分枝；

step4以最开始的目标函数值作为上界，子问题求解中得到的任一整数解为下界（定界），对子问题进行剪枝，减小问题规模；

step5重复以上步骤直到得到最优解

割平面法

step1不考虑整数约束的情况下求解得到最优解 （一般不是整数）；

step2通过该解做一个割平面（二维情况下为一条直线），缩小可行域；

step3在缩小后的可行域中求最优解（不考虑整数约束）

step4重复步骤2和步骤3，直到最优解满足整数约束

0-1规划模型

当整数规划问题中的整数型决策变量限制为只能取0或1时，称为0-1整数规划，简称为0-1规划。

因为0-1规划问题的解空间比一般的整数规划问题较少，求解起来较为容易，且所有的整数规划问题都可以化为0-1规划问题，所以在建立混合整数规划模型求解实际问题时，应尽量使用0-1决策变量进行建模。

例如：有十个工厂可供决策时，可以使用10个0-1变量，当取值为0时时代表不使用这个工厂，取值为1时使用该工厂。

0-1规划的常用求解方法：分支定界算法、割平面法、隐枚举法

案例：投资的收益和风险

问题描述与分析

符号规定

模型假设

建立与简化模型

对于一个多目标优化模型，常用的考虑方式为先固定其中一个目标，再优化另一个目标。

使用python scipy库求解

#导入相关库
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as op

#定义a的取值
a = 0
profit_list = [] #记录最大收益
a_list = [] #记录a的取值




while a<0.05:
    #定义决策变量取值范围
    x1=(0,None)

    #定义目标函数系数
    c=np.array([-0.05,-0.27,-0.19,-0.185,-0.185]) 
    #定义不等式约束条件左边系数
    A = np.hstack((np.zeros((4,1)),np.diag([0.025,0.015,0.055,0.026])))
    #定义不等式约束条件右边系数
    b=a*np.ones((4,1));
    #定义等式约束条件左边系数
    Aeq=np.array([[1,1.01,1.02,1.045,1.065]])
    #定义等式约束条件右边系数
    beq=np.array([1]);
    #求解
    res=op.linprog(c,A,b,Aeq,beq,bounds=(x1,x1,x1,x1,x1))
    profit = -res.fun
    profit_list.append(profit)
    a_list.append(a)
    a = a+0.001

#绘制风险偏好a与最大收益的曲线图    
plt.figure(figsize=(10,7))
plt.plot(a_list,profit_list)
plt.xlabel('a');plt.ylabel('Profit')

Text(0, 0.5, 'Profit')

从上图中可以看出

1.风险越大，收益也就越大；

2，当投资越分散时，投资者承担的风险越小，这与题意一致。即： 冒险的投资者会出现集中投资的情况，保守的投资者则尽量分散投资。

模型求解的其他思路

在上面的例子中，我们使用固定风险水平来最大化收益的方法来将多目标转化为单目标，也可考虑其他思路：