赞
踩
线性规划求解需清晰两部分,目标函数(max,min)和约束条件(s.t.),求解前应转化为标准形式:
举例:求解下列动态规划问题
[x,fval]=linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS)
LB和UB分别为x的上界和下界
===============================================================================================
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.最大值,最小值要看清:
默认求最小值,所以求得最大值时要加 “-”
Python代码如下:
1.整数规划的模型与线性规划基本相同,只是额外增加了部分变量为整数的约束。
2.整数规划求解的基本框架是分支定界法,首先去除整数约束得到“松弛模型”,使用线性规划的方法求解。
3.若有某个变量不是整数,在松弛模型上分别添加约束:x<floor(A)和x ≥ceil(A),然后再分别求解,这个过程叫做分支。当节点求解结果中所有变量都是整数时,停止分支。这样不断迭代,形成了一颗树。所谓定界,指的是叶子节点产生后,相当于给问题定了一个下界。之后在求解过程中一旦某个节点的目标函数值小于这个下界,那就直接pass,不再进行分支了;每次新产生叶子节点,则更新下界。
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
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。