python之建模规划篇

一.线性规划

线性规划求解需清晰两部分，目标函数（max,min）和约束条件（s.t.）,求解前应转化为标准形式：

举例：求解下列动态规划问题

1.MatLab求解：

 [x,fval]=linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS)

LB和UB分别为x的上界和下界
===============================================================================================

2. scipy库求解：

from scipy import optimize
import nunmpy as np
# 求解函数
res=optimize.linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS)
#目标函数最小值
print(res.fun)
#最优解
print(res.x)

Python代码如下：

输出的结果：

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

结果这样看：
1.fun: 最值
2.x: 最优解

注意事项：
1.符号要一致，对于这个问题做了以下调整：
2X1-5X2+X3>=10 => -2X1=5X2-X3<=-10
2.最大值，最小值要看清：
默认求最小值，所以求得最大值时要加 “-”

3.pulp库求解：

Python代码如下：

二.整数规划

1.整数规划的模型与线性规划基本相同，只是额外增加了部分变量为整数的约束。
2.整数规划求解的基本框架是分支定界法，首先去除整数约束得到“松弛模型”，使用线性规划的方法求解。
3.若有某个变量不是整数，在松弛模型上分别添加约束:x<floor(A)和x ≥ceil(A),然后再分别求解，这个过程叫做分支。当节点求解结果中所有变量都是整数时，停止分支。这样不断迭代，形成了一颗树。所谓定界，指的是叶子节点产生后，相当于给问题定了一个下界。之后在求解过程中一旦某个节点的目标函数值小于这个下界，那就直接pass，不再进行分支了；每次新产生叶子节点，则更新下界。

1.分支定界代码：

from scipy.optimize import linprog
import numpy as np
import math
import sys
from queue import Queue


class ILP():
    def __init__(self, c, A_ub, b_ub, A_eq, b_eq, bounds):
        # 全局参数
        self.LOWER_BOUND=-sys.maxsize
        self.UPPER_BOUND = sys.maxsize
        self.opt_val = None
        self.opt_x = None
        self.Q = Queue()

        # 这些参数在每轮计算中都不会改变，因为求最大值所以c=-c
        self.c = -c
        self.A_eq = A_eq
        self.b_eq = b_eq
        self.bounds = bounds

        # 首先计算一下初始问题
        r = linprog(-c, A_ub, b_ub, A_eq, b_eq, bounds)

        # 若最初问题线性不可解
        if not r.success:
            raise ValueError('Not a feasible problem!')

        # 将解和约束参数放入队列
        self.Q.put((r, A_ub, b_ub))

    def solve(self):
        while not self.Q.empty():
            # 取出当前问题
            res, A_ub, b_ub = self.Q.get(block=False)

            # 当前最优值小于总下界，则排除此区域
            if -res.fun < self.LOWER_BOUND:
                continue

            # 若结果 x 中全为整数，则尝试更新全局下界、全局最优值和最优解
            if all(list(map(lambda f: f.is_integer(), res.x))):
                if self.LOWER_BOUND < -res.fun:
                    self.LOWER_BOUND = -res.fun

                if self.opt_val is None or self.opt_val < -res.fun:
                    self.opt_val = -res.fun
                    self.opt_x = res.x

                continue

            # 进行分枝
            else:
                # 寻找 x 中第一个不是整数的，取其下标 idx
                idx = 0
                for i, x in enumerate(res.x):
                    if not x.is_integer():
                        break
                    idx += 1

                # 构建新的约束条件（分割
                new_con1 = np.zeros(A_ub.shape[1]) #返回长度为2的一维数组[0,0]
                new_con1[idx] = -1                 #此时new_con1=[-1,0]
                new_con2 = np.zeros(A_ub.shape[1]) #返回长度为2的一维数组[0,0]
                new_con2[idx] = 1                  #此时new_con2=[1,0]

                #添入新的约束条件，此时new_A_ub_1=[[ 9  7][ 7 20][-1  0]],不懂的可以参照numpy的insert函数用法
                new_A_ub1 = np.insert(A_ub, A_ub.shape[0], new_con1, axis=0)
                # 添入新的约束条件，此时new_A_ub_2=[[ 9  7][ 7 20][1  0]]
                new_A_ub2 = np.insert(A_ub, A_ub.shape[0], new_con2, axis=0)

                #此时new_b_ub1=[[56,70],[-5,0]]
                new_b_ub1 = np.insert(
                    b_ub, b_ub.shape[0], -math.ceil(res.x[idx]), axis=0)
                # 此时new_b_ub2=[[56,70],[4,0]]
                new_b_ub2 = np.insert(
                    b_ub, b_ub.shape[0], math.floor(res.x[idx]), axis=0)

                # 将新约束条件加入队列，先加最优值大的那一支
                r1 = linprog(self.c, new_A_ub1, new_b_ub1, self.A_eq,
                             self.b_eq, self.bounds)
                r2 = linprog(self.c, new_A_ub2, new_b_ub2, self.A_eq,
                             self.b_eq, self.bounds)
                if not r1.success and r2.success:
                    self.Q.put((r2, new_A_ub2, new_b_ub2))
                elif not r2.success and r1.success:
                    self.Q.put((r1, new_A_ub1, new_b_ub1))
                elif r1.success and r2.success:
                    if -r1.fun > -r2.fun:
                        self.Q.put((r1, new_A_ub1, new_b_ub1))
                        self.Q.put((r2, new_A_ub2, new_b_ub2))
                    else:
                        self.Q.put((r2, new_A_ub2, new_b_ub2))
                        self.Q.put((r1, new_A_ub1, new_b_ub1))


def test1():
    """ 此测试的真实，最优值为340，最优解为 [4, 2] """
    c = np.array([40, 90])
    A = np.array([[9, 7], [7, 20]])
    b = np.array([56, 70])
    Aeq = None
    beq = None
    bounds = [(0, None), (0, None)]

    solver = ILP(c, A, b, Aeq, beq, bounds)
    solver.solve()

    print("Test 1's result:", solver.opt_val, solver.opt_x)
    print("Test 1's true optimal x: [4, 2]\n")



if __name__ == '__main__':
    test1()

摘选自知乎：https://zhuanlan.zhihu.com/p/386049272

2.非线性规划：

(1).非线性规划可以简单分两种，目标函数为凸函数or非凸函数 (2).凸函数的非线性规划，比如fun = x2+y+xy,有很多常用库完成，比如cvxpy (3).非凸函数的非线性规划（求极值），可以尝试以下方法: 纯数学方法，求导求极值 (3).神经网络、深度学习（反向传播算法中链式求导过程) (4).scipy.optimize.minimize (minimize局部最优)