用 Scipy 解决优化的线性规划问题

Question

我有以下线性编。我想使用Scipy 解决的问题。

Maximize: x0 * c + x1 * d
Such that:
          x0 * a + b * x1 >= 0
          x0 + y0 = 1
          x0, x1 belong [0,1]

我试过这个：

from scipy.optimize import linprog

c = [c, d]
A = [[-a, -b], [1, 1]]
b = [0, 1]
x0_bounds = (0, 1)
x1_bounds = (0, 1)
res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds])

a, b, c, d = 1.0, -1.0, 1.0, 5.0 的答案应该是 x0 = 0.5, x1 = 0.5，但我知道了：

     con: array([], dtype=float64)
     fun: 1.4471213546122847e-12
 message: 'Optimization terminated successfully.'
     nit: 5
   slack: array([-2.45517088e-13,  1.00000000e+00])
  status: 0
 success: True
       x: array([3.65893190e-14, 2.82106407e-13])

有什么想法吗？

Answer 1

scipy 中的 linprog 有时不方便，因为：

1. 它总是解决一个最小化问题，所以如果你想最大化一个目标函数，你必须做一个解决方法like in this solution 将其转换为一个最小化问题

2. 具有 >= 的方程需要乘以 -1 才能变为

3. 一起创建约束，比如A_ub A_eq，它们是分离的矩阵，所以单独创建

看看docs，他们也有一个很好的例子

linprog 正在实现一个求解器：

因此，如果您需要 A_lb x >= b_lb 形式的不等式，则需要将其转换为形式（表达式）

第一个解决您最初问题的简单解决方案

from scipy.optimize import linprog
import numpy as np
####defining the minimization problem
##### those are your inputs  following scipy notation
a,b,c,d = 10.0,-3.0,10.0,4.666666666666666
c0, c1= a, b # this is your a and b  10x0 - 3x1 is the obj func to minimize
a_ub0, a_ub1= a, b # this is your a and b in "x0*a + b*x1"  
b0_ub = 0 # this is your constraint in "x0*a + b*x1 >= 0" 
a_eq0, a_eq1= 1, 1 # this is your "1" in "1*x0 + 1*y0"
b0_eq = 1 # this is your 1 constraint in "x0 + x1 = 1" 
x0_bounds = (0, 1) # those are the bounds
x1_bounds = (0, 1)

C = [c0, c1] # you're minimizing so no need to multiply by -1 C
A_ub = [[-a_ub0, -a_ub1]] # but still need to invert the signal here
A_eq = [[a_eq0, a_eq1]]
b_ub = [b0_ub]
b_eq = [b0_eq]
res = linprog(C, A_ub=A_ub, b_ub=b_ub,
                 A_eq=A_eq, b_eq=b_eq,
                 bounds=[x0_bounds, x1_bounds])
minimization_objfunc_output = res['fun']
print(f'min value is {res["fun"]} with solutions x: {res["x"]}')


####defining your original problem
##### those are your inputs  following scipy notation
a,b,c,d = 10.0,-3.0,10.0,4.666666666666666
c0, c1= c, d # this is your c and d 
a_ub0, a_ub1= a, b # this is your a and b in "x0*a + b*x1" 
b0_ub = 0 # this is your constraint in "x0*a + b*x1 >= 0" 
a_eq0, a_eq1= 1, 1 # this is your "1" in "1*x0 + 1*y0"
b0_eq = 1 # this is your 1 constraint in "x0 + x1 = 1" 
a_eq2, a_eq3= a, b # this is your "a" and b in "a*x0 + b*y0"
b1_eq = minimization_objfunc_output # this is your minimization constraint in "a*x0 + b* x1 = result of minimization" 
x0_bounds = (0, 1) # those are the bounds
x1_bounds = (0, 1)

C = [-c0, -c1]
A_ub = [[-a_ub0, -a_ub1]]
A_eq = [[a_eq0, a_eq1],[a_eq2, a_eq3]]
b_ub = [b0_ub]
b_eq = [b0_eq, b1_eq]
res = linprog(C, A_ub=A_ub, b_ub=b_ub,
                 A_eq=A_eq, b_eq=b_eq,
                 bounds=[x0_bounds, x1_bounds])
print(f'max value is {res["fun"]} with solutions x: {res["x"]}')

输出：

min value is 5.210942788380635e-12 with solutions x: [0.23076923 0.76923077]
max value is -5.897435897450446 with solutions x: [0.23076923 0.76923077]

第二个更通用的解决方案可以将任何问题定义适应 LinProg 格式

from scipy.optimize import linprog
import numpy as np
##### Ideal generic inputs
def convertInputs_toLinProg(objective, 
                            c, 
                            A_GtE, b_GtE,
                            A_LtE, b_LtE,
                            A_Eq, b_Eq,
                            list_of_bounds):
    A_GtE = -1 * A_GtE
    b_GtE = -1 * b_GtE
    if not (A_GtE is None) and not (A_LtE is None): 
        A_ub = np.vstack([A_GtE,A_GtE])
    elif not (A_GtE is None) and (A_LtE is None):
        A_ub = A_GtE
    elif (A_GtE is None) and not (A_LtE is None):
        A_ub = A_GtE

    if objective == "maximize":
        c = -1*c
    return {
            "c":c, 
            "A_ub": A_ub,
            "A_eq": A_Eq,
            "b_ub": b_GtE,
            "b_eq": b_Eq,
            "bounds": list_of_bounds 
           }

a,b,c,d = 10.0,-3.0,10.0,4.666666666666666
####defining the minimization problem here  
##### those are your inputs  following scipy notation
# objective function:
#    10*x0 - 3*x1 
c = [a,     b]
c = np.array(c)
#                           10*x0 -3*x1 >= 0
A_greaterEq, b_greaterEq =[[a,   b]],   [0]
A_greaterEq, b_greaterEq =np.array(A_greaterEq), np.array(b_greaterEq)
#                      1*x0 + 1*x1 = 1
A_equals, b_equals = [[1,     1]],   [1]
A_equals, b_equals = np.array(A_equals), np.array(b_equals) 
# bounds are exactly how you define them
x0_bounds, x1_bounds = (0, 1), (0, 1)
inputs = convertInputs_toLinProg(objective="minimize",
                                 c=c, 
                                 A_GtE= A_greaterEq, b_GtE= b_greaterEq,
                                 A_LtE= None,       b_LtE= None,
                                 A_Eq=  A_equals,   b_Eq= b_equals,
                                 list_of_bounds=[x0_bounds, x1_bounds])
res = linprog(**inputs)
minimization_objfunc_output = res['fun']
print(f'min value is {res["fun"]} with solutions x: {res["x"]}')
##### generic definition of your problem
####defining your original problem
##### those are your inputs  following scipy notation
# objective function:
#   1*x0 + 5*x1 
a,b,c,d = 10.0,-3.0,10.0,4.666666666666666
c = [c,d]
c = np.array(c)
#                           1*x0 -1*x1 >= 0
A_greaterEq, b_greaterEq =[[a,   b]],   [0]
A_greaterEq, b_greaterEq =np.array(A_greaterEq), np.array(b_greaterEq)
#                      1*x0 + 1*x1 = 1
A_equals, b_equals = [[1,     1],[a,     b]],   [1,minimization_objfunc_output]
A_equals, b_equals = np.array(A_equals), np.array(b_equals) 
# bounds are exactly how you define them
x0_bounds, x1_bounds = (0, 1), (0, 1)

inputs = convertInputs_toLinProg(objective="maximize",
                                 c=c, 
                                 A_GtE= A_greaterEq, b_GtE= b_greaterEq,
                                 A_LtE= None,       b_LtE= None,
                                 A_Eq=  A_equals,   b_Eq= b_equals,
                                 list_of_bounds=[x0_bounds, x1_bounds])

res = linprog(**inputs)
print(f'max value is {res["fun"]} with solutions x: {res["x"]}')

输出 2

min value is 5.210942788380635e-12 with solutions x: [0.23076923 0.76923077]
max value is 5.897435897450446 with solutions x: [0.23076923 0.76923077]