鱼C论坛

 找回密码
 立即注册
查看: 15459|回复: 3

[技术交流] 【原创】求解线性规划题目的利器:pymprog库

 关闭 [复制链接]
发表于 2017-9-25 23:04:07 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能^_^

您需要 登录 才可以下载或查看,没有账号?立即注册

x
今天来介绍一个求解线性规划题目的利器:pymprog库

先从一个简单的例子开始,我们来认识一下pymprog。

  1. 已知约束:
  2. x  <=  3  # mountain bike limit
  3. y  <=  4  # racer limit
  4. x + y  <=  5  # frame limit
  5. x >=0, y >=0     # non-negative
  6. 求最大化:
  7. maximize  15 x + 10 y         # profit
复制代码


这就是典型的线性规划的题目,看看用pymprog是怎么求解的。
  1. from pymprog import *
  2. begin('bike production')
  3. x, y = var('x, y') # variables
  4. maximize(15 * x + 10 * y, 'profit')
  5. x <= 3 # mountain bike limit
  6. y <= 4 # racer production limit
  7. x + y <= 5 # metal finishing limit
  8. solve()
复制代码

输出:
GLPK Simplex Optimizer, v4.60
1 row, 2 columns, 2 non-zeros
*     0: obj =  -0.000000000e+00 inf =   0.000e+00 (2)
*     2: obj =   6.500000000e+01 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
结果就是6.500000000e+01 = 65,是不是很简单?

pymprog其实是PythonMathProgram的缩写,它提供了一个现成的模块用来求解多种线性规划问题。

下面让我们来详细认识pymprog库。

1. 安装
  1. pip install pymprog
复制代码

是最简便也是推荐的安装方法,当然你也可以自己去pypi上下载最新的版本手动安装。
注意:pymprog依赖swiglpk库,且版本>=1.3.3。

2. 变量
pymprog中的变量用var(...)来表示,举例
  1. >>> from pymprog import *
  2. >>> begin('test')
  3. model('test') is the default model.
  4. >>> X = var('X')
  5. >>> X # by default, it is non-negative and continuous
  6. X >= 0 continuous
  7. >>> x, y = var('x, y') # many names -> many vars
  8. >>> x, y
  9. (x >= 0 continuous, y >= 0 continuous)
  10. >>> z = var('z', 3)
  11. >>> z # an array of 3 variables
  12. [z[0] >= 0 continuous, z[1] >= 0 continuous, z[2] >= 0 continuous]
  13. >>> z[2] # access the third variable
  14. z[2] >= 0 continuous
  15. >>> v = var('v', kind=bool) # 0/1 variable
  16. >>> v
  17. 0 <= v <= 1 binary
  18. >>> w = var('w', bounds=(0,5)) # specify the bounds
  19. >>> w
  20. 0 <= w <= 5 continuous
  21. >>> colors = ('red', 'green', 'blue') # index set
  22. >>> clr = var('color', colors, bool) # using an index set
  23. >>> clr # a dictionary with keys from the index set
  24. {'blue': 0 <= color['blue'] <= 1 binary, 'green': 0 <= color['green'] <= 1 binary,
  25. 'red': 0 <= color['red'] <= 1 binary}
  26. >>> clr['green']
  27. 0 <= color['green'] <= 1 binary
复制代码

其中,kind关键词可以设置多种类型:
float (default): continuous
int: integer
bool: binary, side-effect: reset bounds to (0,1)
bounds关键词可以指定边界: a pair of numbers, for the lower and upper bounds. If None is used, it means unbounded. The default value is (0, None), so the lower bound is 0, upper bound is none.

3. 组合
pymprog中内置了一种把几个集合(列表或元祖)进行组合的方法:iprod函数。
  1. >>> from pymprog import *
  2. >>> A = (1, 2, 3)
  3. >>> B = (6, 7, 8)
  4. >>> C = iprod(A,B)
  5. >>> for c in C: print(c)
  6. ...
  7. (1, 6) (1, 7) (1, 8) (2, 6) (2, 7) (2, 8) (3, 6) (3, 7) (3, 8)
复制代码

对于生成的C,我们可以直接设定变量:
  1. >>> begin('test')
  2. >>> x = var('x', C)
  3. >>> type(x)
  4. <type 'dict'>
  5. >>> x[2,7].name
  6. 'x[2,7]'
复制代码

这是非常高效,也是非常有用的生成组合的方法,在后面会经常用到,所以必须掌握。

4. 参数
pymprog的参数设置用par(...)来表示,举例:
参数可以用列表形式
  1. >>> from pymprog import *
  2. >>> k = par('k', [2, 3, 4])
  3. >>> type(k)
  4. <type 'list'>
  5. >>> k[1]
  6. (k[1]:3)
复制代码

也可以用字典形式
  1. >>> p = par('P', {'east':0, 'west':2, 'south':3, 'north':1})
  2. >>> type(p)
  3. <type 'dict'>
  4. >>> p['east']
  5. (P['east']:0)
复制代码

更复杂的,可以两者结合
  1. >>> r = [{(3,4):3, (1,2):4}, 5]
  2. >>> R = par('R', r)
  3. >>> R
  4. [{(1, 2): (R[0][1,2]:4), (3, 4): (R[0][3,4]:3)}, (R[1]:5)]
  5. >>> r[0][3,4]
  6. 3
  7. >>> R[0][3,4]
  8. (R[0][3,4]:3)
  9. >>> r[1]
  10. 5
  11. >>> R[1]
  12. (R[1]:5)
复制代码


既然是参数,那么参数的值也是可以变化的,对参数赋值
  1. >>> p.value = new_value
复制代码

看一个例子:
  1. from pymprog import *
  2. begin('trader')
  3. x = var('x', 3)
  4. c = par('c', [100, 300, 50])
  5. b = par('b', [93000, 101, 201])
  6. maximize(sum(c[i]*x[i] for i in range(3)), 'Profit')

  7. 300*x[0] + 1200*x[1] + 120*x[2] <= b[0]
  8. 0.5*x[0] +      x[1] + 0.5*x[2] <= b[1]

  9. solve()
  10. sensitivity()

  11. b[1].value += 1
  12. solve()

  13. end()
复制代码

输出:
  1. GLPK Simplex Optimizer, v4.60
  2. 2 rows, 3 columns, 6 non-zeros
  3. *     0: obj =  -0.000000000e+00 inf =   0.000e+00 (3)
  4. *     2: obj =   2.560000000e+04 inf =   0.000e+00 (0)
  5. OPTIMAL LP SOLUTION FOUND

  6. PyMathProg 1.0 Sensitivity Report Created: 2016/12/10 Sat 15:01PM
  7. ================================================================================
  8. Variable            Activity   Dual.Value     Obj.Coef   Range.From   Range.Till
  9. --------------------------------------------------------------------------------
  10. *x[0]                     94            0          100         87.5          150
  11. *x[1]                     54            0          300          200      366.667
  12. x[2]                      0          -20           50         -inf           70
  13. ================================================================================
  14. ================================================================================
  15. Constraint       Activity Dual.Value  Lower.Bnd  Upper.Bnd RangeLower RangeUpper
  16. --------------------------------------------------------------------------------
  17. R1                 93000   0.166667       -inf      93000      60600     121200
  18. R2                   101        100       -inf        101       77.5        155
  19. ================================================================================
  20. GLPK Simplex Optimizer, v4.60
  21. 2 rows, 3 columns, 6 non-zeros
  22. *     2: obj =   2.570000000e+04 inf =   0.000e+00 (0)
  23. OPTIMAL LP SOLUTION FOUND
复制代码

注:sensitivity函数是非常有用的,它会输出详细的报告,帮助线性规划问题的求解分析。

5. 约束
pymprog的约束用st(...)表示,举例
  1. begin('illustration')
  2. x = var('x', 3)
  3. y = var('y')
  4. c = [6,4,3]
  5. st( sum(p*q for p,q in zip(c,x)) <= y )
  6. st( x[0] + x[2] >= x[1] )
  7. st( x[0] + x[2] <= 1 )
复制代码

其实,这里的st是可以省略的,这样写也同样成立。
  1. begin('illustration')
  2. x = var('x', 3)
  3. y = var('y')
  4. c = [6,4,3]
  5. sum(p*q for p,q in zip(c,x)) <= y
  6. x[0] + x[2] >= x[1]
  7. x[0] + x[2] <= 1
复制代码


再看个例子,
  1. from pymprog import *
  2. fruits = ('apple', 'pear', 'orange')
  3. nutrit = ('fat', 'fibre', 'vitamin')
  4. ntrmin = ( 10, 50, 30 ) # min nutrition intake
  5. #nutrition ingredients of each fruit
  6. ingred = ('apple': (3,4,5), 'pear': (2,4,1),
  7.           'orange': (2,3,4))
  8. diet = var('diet', fruits, int)
  9. for n, ntr in enumerate(nutrit):
  10.    sum( diet[frt] * ingred[frt][n]
  11.           for frt in fruits) >= ntrmin[n]
复制代码


双重比较也是可以接受的,例如
  1. >>> begin('model')
  2. >>> x, y = var('x, y')
  3. >>> 3 <= 4 * x + 5 * y <= 6
复制代码


约束还可以这样写(这语法在纯python环境中是不接受的)
  1. >>> begin('model')
  2. model('model') is the default model.
  3. >>> x, y = var('x, y')
  4. >>> R = 3 * x + y >= 5
  5. >>> R
  6. R1: 3 * x + y >= 5
  7. >>> R.name
  8. 'R1'
  9. >>> R.name = 'Sugar Level'
  10. >>> R
  11. Sugar Level: 3 * x + y >= 5
复制代码


6. 求解
上面讲了那么多条件设定,最终就是为了求解题目。
最简单的,直接solve(),程序会自动调用最匹配的求解器进行求解。

当然,我们也可以设定我们想要使用的求解器。
  1. >>> solver('interior', msg_lev=glpk.GLP_MSG_OFF)
  2. {'msg_lev': 0}
复制代码

或者
  1. >>> solver(int, br_tech=glpk.GLP_BR_PCH)
复制代码

关于求解器部分,可以用help命令自行查找,普通求解的话,可以直接用默认,无需设置。

7. 模块
pymprog可以直接建模用begin(...),也可以设置模型m=model(...),以后可以直接调用m来执行模型。
对于模型中的元素(参数或约束)的删除,可以用x.delete(),例如
  1. from pymprog import *
  2. begin('trader')
  3. x = var('x', 3)
  4. c = par('c', [100, 300, 50])
  5. b = par('b', [93000, 101, 201])
  6. maximize(sum(c[i]*x[i] for i in range(3)), 'Profit')

  7. 300*x[0] + 1200*x[1] + 120*x[2] <= b[0]
  8. 0.5*x[0] +      x[1] + 0.5*x[2] <= b[1]
  9. r = x[0] +          x[1] +     x[2] <= b[2]

  10. solve()
  11. sensitivity()

  12. r.delete()
  13. # deleting a basic varriable destroys the basis
  14. x[1].delete()
  15. # restore the standard basis
  16. std_basis()
  17. solve()
  18. sensitivity()

  19. end()
复制代码


模型的保存,可以用save(...)表示。例如
  1. from pymprog import *
  2. begin('save')
  3. x = var('x', 3)
  4. x[2].kind = int
  5. c = par('c', [100, 300, 50])
  6. b = par('b', [93000, 101, 201])
  7. maximize(sum(c[i]*x[i] for i in range(3)), 'Profit')

  8. 300*x[0] + 1200*x[1] + 120*x[2] <= b[0]
  9. 0.5*x[0] +      x[1] + 0.5*x[2] <= b[1]
  10. r = x[0] +          x[1] +     x[2] <= b[2]

  11. solve()

  12. save(mps='_save.mps', sol='_save.sol',
  13.      clp='_save.clp', glp='_save.glp',
  14.      sen='_save.sen', ipt='_save.ipt',
  15.      mip='_save.mip')

  16. end()
复制代码


错误分析可以调用KKT()函数,例如
  1. from pymprog import *
  2. c = (10, 6, 4)
  3. A = [ ( 1, 1, 1),     
  4.       ( 9, 4, 5),   
  5.       ( 2, 2, 6) ]   
  6. b = (10, 60, 30)
  7. begin('basic') # begin modelling
  8. verbose(True)  # be verbose
  9. x = var('x', 3) #create 3 variables
  10. maximize(sum(c[i]*x[i] for i in range(3)))
  11. for i in range(3):
  12.   sum(A[i][j]*x[j] for j in range(3)) <= b[i]
  13. solve() # solve the model
  14. print(KKT())
  15. end() #Good habit: do away with the model
复制代码

输出:
  1. Max : 10 * x[0] + 6 * x[1] + 4 * x[2]
  2. R1: x[0] + x[1] + x[2] <= 10
  3. R2: 9 * x[0] + 4 * x[1] + 5 * x[2] <= 60
  4. R3: 2 * x[0] + 2 * x[1] + 6 * x[2] <= 30
  5. GLPK Simplex Optimizer, v4.60
  6. 3 rows, 3 columns, 9 non-zeros
  7. *     0: obj =  -0.000000000e+00 inf =   0.000e+00 (3)
  8. *     2: obj =   7.600000000e+01 inf =   0.000e+00 (0)
  9. OPTIMAL LP SOLUTION FOUND
  10. Karush-Kuhn-Tucker optimality conditions:
  11. =========================================
  12. Solver used for this solution: simplex

  13. 1) Error for Primal Equality Constraints:
  14. ----------------------------------------
  15. Largest absolute error: 0.000000 (row id: 0)
  16. Largest relative error: 0.000000 (row id: 0)

  17. 2) Error for Primal Inequality Constraints:
  18. -------------------------------------------
  19. Largest absolute error: 0.000000 (row id: 0)
  20. Largest relative error: 0.000000 (row id: 0)

  21. 1) Error for Dual Equality Constraints:
  22. ----------------------------------------
  23. Largest absolute error: 0.000000 (var id: 0)
  24. Largest relative error: 0.000000 (var id: 0)

  25. 2) Error for Dual Inequality Constraints:
  26. -------------------------------------------
  27. Largest absolute error: 0.000000 (var id: 0)
  28. Largest relative error: 0.000000 (var id: 0)

  29. __del__ is deleting problem: basic
复制代码
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复

使用道具 举报

 楼主| 发表于 2017-9-25 23:26:28 | 显示全部楼层
来看看如何用线性规划的方法来求解“数独”问题。
“数独”其实就是一个线性规划问题,约束:当所有横行、竖行和九宫格都满足1~9时,即求得解。
  1. g = [
  2. [0,0,0,0,0,0,0,0,0],
  3. [0,0,0,0,0,3,0,8,5],
  4. [0,0,1,0,2,0,0,0,0],
  5. [0,0,0,5,0,7,0,0,0],
  6. [0,0,4,0,0,0,1,0,0],
  7. [0,9,0,0,0,0,0,0,0],
  8. [5,0,0,0,0,0,0,7,3],
  9. [0,0,2,0,1,0,0,0,0],
  10. [0,0,0,0,4,0,0,0,9]]

  11. from pymprog import model, iprod
  12. p = model("sudoku")
  13. I = range(1,10)
  14. J = range(1,10)
  15. K = range(1,10)
  16. T = iprod(I,J,K) #create Indice tuples
  17. x = p.var('x', T, bool)
  18. #x[i,j,k] = 1 means cell [i,j] is assigned number k
  19. #assign pre-defined numbers using the "givens"
  20. [x[i,j,k] == (1 if g[i-1][j-1] == k else 0)
  21.        for (i,j,k) in T if g[i-1][j-1] > 0 ]

  22. #each cell must be assigned exactly one number
  23. [sum(x[i,j,k] for k in K)==1 for i in I for j in J]

  24. #cells in the same row must be assigned distinct numbers
  25. [sum(x[i,j,k] for j in J)==1 for i in I for k in K]

  26. #cells in the same column must be assigned distinct numbers
  27. [sum(x[i,j,k] for i in I)==1 for j in J for k in K]

  28. #cells in the same region must be assigned distinct numbers
  29. [sum(x[i,j,k] for i in range(r,r+3) for j in range(c, c+3))==1
  30.     for r in range(1,10,3) for c in range(1,10,3) for k in K]
  31.    
  32. #there is no need for an objective function here

  33. p.solve()

  34. for i in I:
  35.    if i in range(1,10,3):
  36.       print(" +-------+-------+-------+")
  37.    print('', end=' ')
  38.    for j in J:
  39.       if j in range(1,10,3): print("|", end=' ')
  40.       print("%g"%sum(x[i,j,k].primal*k for k in K), end=' ')
  41.       if j==9: print("|")
  42.    if i == 9:
  43.       print(" +-------+-------+-------+")
  44. p.end()
复制代码

几乎就是秒出的。
  1. +-------+-------+-------+
  2. | 9 8 7 | 6 5 4 | 3 2 1 |
  3. | 2 4 6 | 1 7 3 | 9 8 5 |
  4. | 3 5 1 | 9 2 8 | 7 4 6 |
  5. +-------+-------+-------+
  6. | 1 2 8 | 5 3 7 | 6 9 4 |
  7. | 6 3 4 | 8 9 2 | 1 5 7 |
  8. | 7 9 5 | 4 6 1 | 8 3 2 |
  9. +-------+-------+-------+
  10. | 5 1 9 | 2 8 6 | 4 7 3 |
  11. | 4 7 2 | 3 1 9 | 5 6 8 |
  12. | 8 6 3 | 7 4 5 | 2 1 9 |
  13. +-------+-------+-------+
  14. [Finished in 0.3s]
复制代码

同样的,八皇后问题也是线性规划的一种。
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2017-11-24 10:30:49 | 显示全部楼层
楼主,写得真棒。有个问题想请教下,如果求解大数据量的混合整型线性规划,我看到控制台上会每隔一分钟就打印信息,说哪个求解器找到了可行解,但是最优解却要很久才能找到,目前情况好像是只有找到最优解才会输出结果。我现在就想获取到可行解,然后输出可行解,怎么做到呢?或者有没有限时求解的方法?谢谢
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

 楼主| 发表于 2017-11-24 10:49:35 | 显示全部楼层
闪电侠没注册 发表于 2017-11-24 10:30
楼主,写得真棒。有个问题想请教下,如果求解大数据量的混合整型线性规划,我看到控制台上会每隔一分钟就打 ...

对于这个库的了解好像目前还不支持多重解的情况,只有得到唯一解才会输出答案。
也有可能是我没有全部了解所有功能,你可以自己对照帮助文档自己深入学习。
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|手机版|Archiver|鱼C工作室 ( 粤ICP备18085999号-1 | 粤公网安备 44051102000585号)

GMT+8, 2024-3-29 06:47

Powered by Discuz! X3.4

© 2001-2023 Discuz! Team.

快速回复 返回顶部 返回列表