数学建模——规划问题

线性规划

 运筹学对于线性规划问题直接使用图解法，单纯形法利用求解。在python中可以直接使用scipy.optimize模块的linprog函数求解。
 linprog函数的调用方式：

		scipy.optimize.linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None,
		 method='highs', callback=None, options=None, x0=None, integrality=None)

常用参数解释:
(1)  c:价格向量
(2)  A_ub:不等式约束技术系数矩阵
(3)  b_ub:不等式约束资源向量
(4)  A_eq:等式约束技术系数矩阵
(5)  b_eq:等输约束资源向量
(6)  bounds:表示决策变量的上下界的元组。这里值得注意的是每个决策变量的上下界默认为(0,None),None在linprog函数中表示正/负无穷。如果有变量不是这个范围都要记得修改bounds参数
将数学模型化为linprog支持的形式

示例
 求解的难点就是把上面的各个参数表示出来

# 导入linprog函数和numpy包
import numpy as np
from scipy.optimize import linprog
#注意，对于shape==(n,)的数组，可以理解为列向量，也可以理解为行向量
c=np.array([-1,2,3])
A_ub=np.array([[-2,1,1],[3,-1,-2]])
A_eq=np.array([[4,-2,-1]])
b_ub=np.array([9,-4])
b_eq=np.array([-6])
low=[-10,0,None];high=[None]*3
bounds=tuple(zip(low,high))
res=linprog(c,A_ub,b_ub,A_eq,b_eq,bounds)
print("求解状态:",res.status)
print("最优解:",res.x)
print("最优值:",-res.fun)

整数规划

一般的整数线性规划问题

 什么是整数规划：限制全部或部分决策变量的取值为整数的线性规划问题。
 不比线性规划有单纯形法杀穿一切，没有一种算法可以有效解决一切整数规划问题。但是对于使用python直接用cvxpy无脑求解就行了。
 cvxpy是一种用于凸优化问题的Python嵌入式建模语言。它会自动将问题转换为标准形式，调用解算器，并解开结果。
 关于cvxpy的一些基础操作：

1.用Variables创建决策变量
语法：class cvxpy.expressions.variable.Variable(shape: int | Iterable[int, ...] = (), name: str | None = None, var_id: int | None = None, **kwargs: Any)
常用参数：
（1）shape:决定创建的决策变量的形状，决策变量可以是标量，向量，矩阵。比如我需要5个决策变量：x1,x2,
x3,x4,x5。我就可以这样创建决策变量：x=cp.Varible(5) ,x是一个向量（x1,x2,x3,x4,x5）
（2）integer：设定决策变量为整数。x=cp.Varible(5,integer=True)
2.用Problem生成问题
语法：class cvxpy.Problem(objective: Union[Minimize, Maximize], constraints: Optional[List[Constraint]] = None)
就是把目标（objective）和约束条件（constraints）传入就可以了
3.用Minimize或者Maximize函数生成objective
语法：class cvxpy.Maximize(expr)  其中expr必须是标量
	  class cvxpy.Minimize(expr)
4.用solve求解问题
语法：solve(*args, **kwargs)
常用参数：
（1）solver：求解器，不同的问题可能需要选择不同的求解器下面给出常用求解器以及适用范围

示例

import cvxpy as cp
import numpy as np
c=np.array([40,90])
x=cp.Variable(2,integer=True)
cons1=[np.array([9,7])@x<=56]
cons2=[np.array([7,20])@x>=70]
cons3=[x>=0]
cons=cons1+cons2+cons3
obj=cp.Minimize(c@x)
problem=cp.Problem(obj,cons)
problem.solve(solver="GLPK_MI")
print("求解状态:",problem.status)
print("最优解:",x.value)
print('最优值:',problem.value)

0-1整数规划

 0-1整数规划是指在整数规划的基础上限制决策变量只能取0-1，求解方式与一般的整数规划大同小异，只不过给决策变量加一个限制：既要大于等于0又要小于等于1.于是逼迫决策变量为整数时只能取0或1。
 指派问题的标准形式经常表现为0-1规划模型，标准指派问题的特点是：一个人只可以干一项工作，一项工作只能被一个人干。基于这种特点可以想到指派问题的解是一个置换矩阵（各行各列和为1），因为这样才能保证上述特征。
示例

import cvxpy as cp
import numpy as np
x=cp.Variable((5,5),integer=True)
c=np.array([[ 4,8,7,15,12],
             [ 7,9,17,14,10],
             [ 6,9,12,8,7.],
             [ 6,7,14,6,10.],
             [ 6,9,12,10,6.]])
# 用（*）做元素级乘法报错了，所以改用mutiply。我不太清楚为什么不能用（*）做
obj=cp.Minimize(cp.sum(cp.multiply(c,x)))
cons1=[cp.sum(x,axis=0)==1]
cons2=[cp.sum(x,axis=1)==1]
cons3=[x>=0,x<=1]
cons=cons1+cons2+cons3
prob=cp.Problem(obj,cons)
prob.solve(solver='GLPK_MI')
print(prob.status)
print(prob.value)
print(x.value)

广义指派问题

 广义指派问题和标准指派问题差异在与后者可能是目标函数求极大值，人数与任务量不等，一人可有多个任务，某任务不能有某人完成。要手算的话这种问题需要化为标准形式再用匈牙利算法。但是用python直接当做一般的整数规划问题就行了。
示例

import cvxpy as cp
import numpy as np
x=cp.Variable((2,7),integer=True)
l=np.array([48.7,52.0,61.3,72.0,48.7,52.0,64.0])
w=np.array([2000,3000,1000,500,4000,2000,1000])
a=np.array([8,7,9,6,6,4,8])
# 不太目标为什么cvxpy中(7,)和（2,7）不能广播,在numpy中是可以广播的。改成（1,7）后没问题。
total_load=cp.sum(cp.multiply(l.reshape(1,7),x))
obj=cp.Maximize(total_load)
cons1=[cp.sum(x,axis=0)<=a]
cons2=[cp.multiply(l.reshape(1,7),x)<=1020]
cons3=[cp.multiply(w.reshape(1,7),x)<=40000]
cons4=[cp.sum(cp.multiply(l[4:].reshape(1,3),x[:,4:]))<=302.7]
cons5=[x>=0]
cons=cons1+cons2+cons3+cons4+cons5
prob=cp.Problem(obj,cons)
prob.solve(solver='GLPK_MI')
print(prob.status)
print(prob.value)
print(x.value)

import cvxpy as cp
import numpy as np
# 决策变量X={x_ij | i=1,2 and j=1,2,3,4,5,6,7}, x_ij为整数
X=cp.Variable((2,7),integer=True)
# 厚度系数向量
L=np.array([48.7,52.0,61.3,72.0,48.7,52.0,64.0])
# 重量系数向量
W=np.array([2000,3000,1000,500,4000,2000,1000])
# 件数向量
A=np.array([8,7,9,6,6,4,8])
# 限制条件1:两个平板车上j种货物的装载量之和小于j种货物数量
cons1=[cp.sum(X,axis=0)<=A]
# 限制条件2：一个平板车上各种货物的厚度和小于1020
cons2=[cp.multiply(L.reshape(1,7),X)<=1020]
# 限制条件3：一个平板车上各种货物的重量和小于40000
cons3=[cp.multiply(W.reshape(1,7),X)<=40000]
# 限制条件4：两辆平板车上5,6,7类货箱的厚度和不超过302.7
cons4=[cp.sum(cp.multiply(X[:,4:],L[4:].reshape(1,3)))<=302.7]
# 限制条件5：X>=0
cons5=[X>=0]
# 限制条件
constraints=cons1+cons2+cons3+cons4+cons5
# 目标:两辆平板车装的重量最大
objective=cp.Maximize(cp.sum(cp.multiply(W.reshape(1,7),X)))
# 创建问题
problem=cp.Problem(objective,constraints)
# 求解
problem.solve(solver='GLPK_MI')
print(problem.status)
print(problem.value)
print(X.value)

非线性规划

 手算非线性规划问题比较麻烦，用python来算是比较简单的。scipy.optimize模块的minimize函数或cvxpy包配合上支持求解非线性规划的solver都可以来求解非线性规划问题。
 非线性规划是指有非线性约束条件或目标函数的数学规划，是运筹学的一个重要分支。非线性规划研究一个n元实函数在一组等式或不等式的约束条件下的极值问题，且目标函数和约束条件至少有一个是未知量的非线性函数。
 使用minimize函数求解的基础知识（minimize貌似没法解既有不等式约束又有等式约束的非线性规划问题）

1.语法：scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)
2.常用参数：
（1）fun:要被极小化的可调用函数，它含有一个参数x，这个参数应该被视为一个shape==(n,)的向量
（2）x0:最优解的初始估计值,它也是一个向量且x0.shape==x.shape
（3）bounds：决策变量的上下界的元组
（4）constraints:这个constraints是一个字典，比较常用
的两个key是'type'和'fun'。前者用于表明约束条件的类型：
可选‘eq’或‘ineq’。后者是用来定义约束的可调用函数对象。
这里有一个易错点：Equality constraint means that the constraint function result is to be zero whereas inequality means that it is to be non-negative. Note that COBYLA only supports inequality constraints.
意思就是说如果你的约束类型是等式约束，那么这个等式约束右侧应该为0.如果是不等式约束，那么不等式应该是大于0的，比如AX<=b，应该变形为b-AX>=0.于是f(x)=b-AX.

示例

from scipy.optimize import minimize
import numpy as np
# 这道题的obj是没法用矩阵运算表示出来的，所以用用这种方式定义函数。
# 一维序列如list和tuple都可以当做行向量或列向量。于是x还是一个向量。
def obj(x):
    x1,x2,x3=x
    return (2+x1)/(1+x2)-3*x1+4*x3
# 设初始猜测的最优解看着范围猜，别太离谱
x0=[0.6]*3
low=[0.1]*3;high=[0.9]*3
bounds=tuple(zip(low,high))
# 非线性规划可以没有约束条件
res=minimize(obj,x0,bounds=bounds)
res

from scipy.optimize import minimize
# 目标化为一般形式：min z'=-...,最优解不变，最优值互为相反数
# 这题目标函数也可以写为 lambda x:np.array([-1,-1,-3,-4,-2])@x+np.array([8,2,3,1,2])@x
def objective(x):
    return np.array([-1,-1,-3,-4,-2])@x**2+np.array([8,2,3,1,2])@x
# Equality constraint means that the constraint function result is to be zero whereas inequality means that it is to be non-negative.
# 所以定义不等式约束的函数fun>=0
A=[[1,1,1,1,1],
   [1,2,2,1,6],
   [2,1,6,0,0],
   [0,0,1,1,5]]
b=np.array([400,800,200,200])
def fun(x):
    return b-A@x
constraints={'type':'ineq','fun':fun}
# 决策变量的上下界元组
low=[0]*5
high=[99]*5
bounds=tuple(zip(low,high))
x0=np.ones(5)*50
res=minimize(objective,x0,bounds=bounds,constraints=constraints)
print('求解状态：',res.success)
print("最优解：",res.x)
print("最优值：",-res.fun)

二次规划

 二次规划指的是目标为决策变量的二次函数而约束全部是线性的非线性规划问题。这种问题可以用cvxpy包加上支持二次规划的solver进行求解。
示例

import cvxpy as cp
import numpy as np
c1=np.array([1.5,1,0.85]);c2=np.array([3,-8.2,-1.95])
x=cp.Variable(3) # x.shape==(3,)
objective=cp.Minimize(c1@x**2+c2@x)
cons1=[np.array([1,0,1])@x<=2]
cons2=[np.array([-1,2,0])@x<=2]
cons3=[np.array([0,1,2])@x<=3]
cons4=[np.ones(3)@x==3]
cons=cons1+cons2+cons3+cons4
problem=cp.Problem(objective,cons)
problem.solve(solver="ECOS")
print(problem.status)
print(problem.value)
print(x.value)