Scipy是Python中一个强大且广泛使用的科学计算库,内置了许多优化算法,可以方便地解决线性规划、非线性规划和整数规划等常见的规划问题。

线性规划

线性规划即目标函数、约束条件均为线性的规划问题。使用Scipy中的linprog函数可方便求解:

import scipy.optimize

c = [2, 3, -5]#目标函数系数 2x1+3x2-5x3
A_ub = [[-2, 5, -1],
        [1, 3, 1]]##约束矩阵系数
b_ub = [-10, 12]##约束条件右值
A_eq = [[1, 1, 1]]###约束矩阵等式系数
b_eq = [7.12]###约束等式右值
bounds = [(0, 3.5), (0, 3.6), (0, 100.1)]#变量区间

res = scipy.optimize.linprog(c=c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
print(res)
print(res.get('x'))

A_ub和b_ub分别指定约束条件矩阵和右值。以上建模写成了矩阵形式，方便读入参数建模。求解结果如下：

       message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
        success: True
         status: 0
            fun: -15.440000000000003
              x: [ 2.880e+00  0.000e+00  4.240e+00]
            nit: 0
          lower:  residual: [ 2.880e+00  0.000e+00  4.240e+00]
                 marginals: [ 0.000e+00  5.000e+01  0.000e+00]
          upper:  residual: [ 6.200e-01  3.600e+00  9.586e+01]
                 marginals: [ 0.000e+00  0.000e+00  0.000e+00]
          eqlin:  residual: [ 0.000e+00]
                 marginals: [-1.200e+01]
        ineqlin:  residual: [ 0.000e+00  4.880e+00]
                 marginals: [-7.000e+00 -0.000e+00]
 mip_node_count: 0
 mip_dual_bound: 0.0
        mip_gap: 0.0
[2.88 0.   4.24]

如果是对矩阵不熟悉，这里转化成常规约束形式：

#目标函数（默认为求最小值，加负号可表示求最大值）
2x1+3x2-5x3
#约束函数
-2x1+5x2-x3<=-10
x1+3x2+x3<=12
x1+x2+x3=7.12
#变量区间
x1,x2,x3>0

非线性规划

非线性规划即目标函数和/或约束条件含有非线性项。Scipy通过minimize函数提供了多种非线性优化算法求解:

import scipy.optimize
import numpy as np

# 'eq'表示表达式等于0，'ineq'表示表达式大于等于0
f = lambda x: x[0]**2 +x[3]**2 + x[0]*x[1]*x[2]-100
K1 = lambda x: x[3]*x[2]+x[1]*x[2]*x[4]-x[0]*x[1]*x[2]*x[3]
K2 = lambda x: x[4]**3-(x[1]-5)**2-x[4]*x[2]- x[2]
K3 = lambda x: x[0]-x[1]+x[2]+x[3]-x[4]-100
constraints = [{'type': 'ineq', 'fun': lambda x: K1(x)},##大于等于0
               {'type': 'ineq', 'fun': lambda x: K2(x)},##大于等于0
               {'type': 'eq', 'fun': lambda x: K3(x)}##等于0   
              ]
bounds = [(1, 100), (1, 100), (1, 100), (1, 100), (1, 100)]

#设置初始猜测值
x0 = np.array([2, 2, 3, 0, 0])

res = scipy.optimize.minimize(f, x0, method='SLSQP', bounds = bounds, constraints=constraints)
print(res)
print("目标值",f(res.x))
print('第一个约束左值',K1(res.x))
print('第二个约束左值',K2(res.x))
print('第三个约束左值',K3(res.x))

minimize函数可以指定目标函数f、初值x0、约束条件constraints，这里选择SLSQP算法求解。 求解结果如下：

 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 112.00974926638202
       x: [ 1.000e+00  1.000e+00  1.000e+02  1.054e+01  1.054e+01]
     nit: 14
     jac: [ 1.020e+02  1.000e+02  1.000e+00  2.107e+01  0.000e+00]
    nfev: 72
    njev: 12
目标值 112.00974926638202
第一个约束左值 1053.6116421944012
第二个约束左值 6.110099093348254e-10
第三个约束左值 0.0

整数规划

整数规划只需要把求解函数改成milp。在设置变量和约束，按照milp函数内部设定要求即可。这边给出代码和注释供大家学习。

import numpy as np
from scipy.optimize import LinearConstraint
from scipy.optimize import milp
c = np.array([1, -3])#x1-x2
A = np.array([[-1, 1], [3, 2], [2, 3]])#约束矩阵
b_u = np.array([1, 12, 12])##右边界
b_l = np.full_like(b_u, 0)##变量均大于0
constraints = LinearConstraint(A, b_l, b_u)
integrality = np.ones_like(c)
res = milp(c=c, constraints=constraints, integrality=integrality)
print(res)
print(res.x)

输出结果：

        message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
        success: True
         status: 0
            fun: -5.0
              x: [ 1.000e+00  2.000e+00]
 mip_node_count: 1
 mip_dual_bound: -5.0

如果不需要限制变量和目标值为整数，可以把设置求解整数规划的条件拿掉。

from scipy.optimize import LinearConstraint
from scipy.optimize import milp
c = np.array([1, -3])#x1-x2
A = np.array([[-1, 1], [3, 2], [2, 3]])#约束矩阵
b_u = np.array([1, 12, 12])##右边界
b_l = np.full_like(b_u, 0)##变量均大于0
constraints = LinearConstraint(A, b_l, b_u)
res = milp(c=c, constraints=constraints,)
print(res)
print(res.x)

输出结果：

 message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
        success: True
         status: 0
            fun: -6.599999999999998
              x: [ 1.800e+00  2.800e+00]
 mip_node_count: None
 mip_dual_bound: None
        mip_gap: None
[1.8 2.8]

以上示例展示了通过Scipy库内置的优化函数可以方便高效地解决常见线性、非线性和整数规划问题。另外，Scipy在求解这些问题其实提供了多种算法可以选择，这里没有详细介绍，感兴趣可以去官网查看。

参考资料：scipy官方文档 https://scipy.org/