|
马上注册,结交更多好友,享用更多功能^_^
您需要 登录 才可以下载或查看,没有账号?立即注册
x
今天来介绍一个求解线性规划题目的利器:pymprog库
先从一个简单的例子开始,我们来认识一下pymprog。
- 已知约束:
- x <= 3 # mountain bike limit
- y <= 4 # racer limit
- x + y <= 5 # frame limit
- x >=0, y >=0 # non-negative
- 求最大化:
- maximize 15 x + 10 y # profit
复制代码
这就是典型的线性规划的题目,看看用pymprog是怎么求解的。
- from pymprog import *
- begin('bike production')
- x, y = var('x, y') # variables
- maximize(15 * x + 10 * y, 'profit')
- x <= 3 # mountain bike limit
- y <= 4 # racer production limit
- x + y <= 5 # metal finishing limit
- 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. 安装
是最简便也是推荐的安装方法,当然你也可以自己去pypi上下载最新的版本手动安装。
注意:pymprog依赖swiglpk库,且版本>=1.3.3。
2. 变量
pymprog中的变量用var(...)来表示,举例
- >>> from pymprog import *
- >>> begin('test')
- model('test') is the default model.
- >>> X = var('X')
- >>> X # by default, it is non-negative and continuous
- X >= 0 continuous
- >>> x, y = var('x, y') # many names -> many vars
- >>> x, y
- (x >= 0 continuous, y >= 0 continuous)
- >>> z = var('z', 3)
- >>> z # an array of 3 variables
- [z[0] >= 0 continuous, z[1] >= 0 continuous, z[2] >= 0 continuous]
- >>> z[2] # access the third variable
- z[2] >= 0 continuous
- >>> v = var('v', kind=bool) # 0/1 variable
- >>> v
- 0 <= v <= 1 binary
- >>> w = var('w', bounds=(0,5)) # specify the bounds
- >>> w
- 0 <= w <= 5 continuous
- >>> colors = ('red', 'green', 'blue') # index set
- >>> clr = var('color', colors, bool) # using an index set
- >>> clr # a dictionary with keys from the index set
- {'blue': 0 <= color['blue'] <= 1 binary, 'green': 0 <= color['green'] <= 1 binary,
- 'red': 0 <= color['red'] <= 1 binary}
- >>> clr['green']
- 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函数。
- >>> from pymprog import *
- >>> A = (1, 2, 3)
- >>> B = (6, 7, 8)
- >>> C = iprod(A,B)
- >>> for c in C: print(c)
- ...
- (1, 6) (1, 7) (1, 8) (2, 6) (2, 7) (2, 8) (3, 6) (3, 7) (3, 8)
复制代码
对于生成的C,我们可以直接设定变量:
- >>> begin('test')
- >>> x = var('x', C)
- >>> type(x)
- <type 'dict'>
- >>> x[2,7].name
- 'x[2,7]'
复制代码
这是非常高效,也是非常有用的生成组合的方法,在后面会经常用到,所以必须掌握。
4. 参数
pymprog的参数设置用par(...)来表示,举例:
参数可以用列表形式
- >>> from pymprog import *
- >>> k = par('k', [2, 3, 4])
- >>> type(k)
- <type 'list'>
- >>> k[1]
- (k[1]:3)
复制代码
也可以用字典形式
- >>> p = par('P', {'east':0, 'west':2, 'south':3, 'north':1})
- >>> type(p)
- <type 'dict'>
- >>> p['east']
- (P['east']:0)
复制代码
更复杂的,可以两者结合
- >>> r = [{(3,4):3, (1,2):4}, 5]
- >>> R = par('R', r)
- >>> R
- [{(1, 2): (R[0][1,2]:4), (3, 4): (R[0][3,4]:3)}, (R[1]:5)]
- >>> r[0][3,4]
- 3
- >>> R[0][3,4]
- (R[0][3,4]:3)
- >>> r[1]
- 5
- >>> R[1]
- (R[1]:5)
复制代码
既然是参数,那么参数的值也是可以变化的,对参数赋值
看一个例子:
- from pymprog import *
- begin('trader')
- x = var('x', 3)
- c = par('c', [100, 300, 50])
- b = par('b', [93000, 101, 201])
- maximize(sum(c[i]*x[i] for i in range(3)), 'Profit')
- 300*x[0] + 1200*x[1] + 120*x[2] <= b[0]
- 0.5*x[0] + x[1] + 0.5*x[2] <= b[1]
- solve()
- sensitivity()
- b[1].value += 1
- solve()
- end()
复制代码
输出:
- GLPK Simplex Optimizer, v4.60
- 2 rows, 3 columns, 6 non-zeros
- * 0: obj = -0.000000000e+00 inf = 0.000e+00 (3)
- * 2: obj = 2.560000000e+04 inf = 0.000e+00 (0)
- OPTIMAL LP SOLUTION FOUND
- PyMathProg 1.0 Sensitivity Report Created: 2016/12/10 Sat 15:01PM
- ================================================================================
- Variable Activity Dual.Value Obj.Coef Range.From Range.Till
- --------------------------------------------------------------------------------
- *x[0] 94 0 100 87.5 150
- *x[1] 54 0 300 200 366.667
- x[2] 0 -20 50 -inf 70
- ================================================================================
- ================================================================================
- Constraint Activity Dual.Value Lower.Bnd Upper.Bnd RangeLower RangeUpper
- --------------------------------------------------------------------------------
- R1 93000 0.166667 -inf 93000 60600 121200
- R2 101 100 -inf 101 77.5 155
- ================================================================================
- GLPK Simplex Optimizer, v4.60
- 2 rows, 3 columns, 6 non-zeros
- * 2: obj = 2.570000000e+04 inf = 0.000e+00 (0)
- OPTIMAL LP SOLUTION FOUND
复制代码
注:sensitivity函数是非常有用的,它会输出详细的报告,帮助线性规划问题的求解分析。
5. 约束
pymprog的约束用st(...)表示,举例
- begin('illustration')
- x = var('x', 3)
- y = var('y')
- c = [6,4,3]
- st( sum(p*q for p,q in zip(c,x)) <= y )
- st( x[0] + x[2] >= x[1] )
- st( x[0] + x[2] <= 1 )
复制代码
其实,这里的st是可以省略的,这样写也同样成立。
- begin('illustration')
- x = var('x', 3)
- y = var('y')
- c = [6,4,3]
- sum(p*q for p,q in zip(c,x)) <= y
- x[0] + x[2] >= x[1]
- x[0] + x[2] <= 1
复制代码
再看个例子,
- from pymprog import *
- fruits = ('apple', 'pear', 'orange')
- nutrit = ('fat', 'fibre', 'vitamin')
- ntrmin = ( 10, 50, 30 ) # min nutrition intake
- #nutrition ingredients of each fruit
- ingred = ('apple': (3,4,5), 'pear': (2,4,1),
- 'orange': (2,3,4))
- diet = var('diet', fruits, int)
- for n, ntr in enumerate(nutrit):
- sum( diet[frt] * ingred[frt][n]
- for frt in fruits) >= ntrmin[n]
复制代码
双重比较也是可以接受的,例如
- >>> begin('model')
- >>> x, y = var('x, y')
- >>> 3 <= 4 * x + 5 * y <= 6
复制代码
约束还可以这样写(这语法在纯python环境中是不接受的)
- >>> begin('model')
- model('model') is the default model.
- >>> x, y = var('x, y')
- >>> R = 3 * x + y >= 5
- >>> R
- R1: 3 * x + y >= 5
- >>> R.name
- 'R1'
- >>> R.name = 'Sugar Level'
- >>> R
- Sugar Level: 3 * x + y >= 5
复制代码
6. 求解
上面讲了那么多条件设定,最终就是为了求解题目。
最简单的,直接solve(),程序会自动调用最匹配的求解器进行求解。
当然,我们也可以设定我们想要使用的求解器。
- >>> solver('interior', msg_lev=glpk.GLP_MSG_OFF)
- {'msg_lev': 0}
复制代码
或者
- >>> solver(int, br_tech=glpk.GLP_BR_PCH)
复制代码
关于求解器部分,可以用help命令自行查找,普通求解的话,可以直接用默认,无需设置。
7. 模块
pymprog可以直接建模用begin(...),也可以设置模型m=model(...),以后可以直接调用m来执行模型。
对于模型中的元素(参数或约束)的删除,可以用x.delete(),例如
- from pymprog import *
- begin('trader')
- x = var('x', 3)
- c = par('c', [100, 300, 50])
- b = par('b', [93000, 101, 201])
- maximize(sum(c[i]*x[i] for i in range(3)), 'Profit')
- 300*x[0] + 1200*x[1] + 120*x[2] <= b[0]
- 0.5*x[0] + x[1] + 0.5*x[2] <= b[1]
- r = x[0] + x[1] + x[2] <= b[2]
- solve()
- sensitivity()
- r.delete()
- # deleting a basic varriable destroys the basis
- x[1].delete()
- # restore the standard basis
- std_basis()
- solve()
- sensitivity()
- end()
复制代码
模型的保存,可以用save(...)表示。例如
- from pymprog import *
- begin('save')
- x = var('x', 3)
- x[2].kind = int
- c = par('c', [100, 300, 50])
- b = par('b', [93000, 101, 201])
- maximize(sum(c[i]*x[i] for i in range(3)), 'Profit')
- 300*x[0] + 1200*x[1] + 120*x[2] <= b[0]
- 0.5*x[0] + x[1] + 0.5*x[2] <= b[1]
- r = x[0] + x[1] + x[2] <= b[2]
- solve()
- save(mps='_save.mps', sol='_save.sol',
- clp='_save.clp', glp='_save.glp',
- sen='_save.sen', ipt='_save.ipt',
- mip='_save.mip')
- end()
复制代码
错误分析可以调用KKT()函数,例如
- from pymprog import *
- c = (10, 6, 4)
- A = [ ( 1, 1, 1),
- ( 9, 4, 5),
- ( 2, 2, 6) ]
- b = (10, 60, 30)
- begin('basic') # begin modelling
- verbose(True) # be verbose
- x = var('x', 3) #create 3 variables
- maximize(sum(c[i]*x[i] for i in range(3)))
- for i in range(3):
- sum(A[i][j]*x[j] for j in range(3)) <= b[i]
- solve() # solve the model
- print(KKT())
- end() #Good habit: do away with the model
复制代码
输出:
- Max : 10 * x[0] + 6 * x[1] + 4 * x[2]
- R1: x[0] + x[1] + x[2] <= 10
- R2: 9 * x[0] + 4 * x[1] + 5 * x[2] <= 60
- R3: 2 * x[0] + 2 * x[1] + 6 * x[2] <= 30
- GLPK Simplex Optimizer, v4.60
- 3 rows, 3 columns, 9 non-zeros
- * 0: obj = -0.000000000e+00 inf = 0.000e+00 (3)
- * 2: obj = 7.600000000e+01 inf = 0.000e+00 (0)
- OPTIMAL LP SOLUTION FOUND
- Karush-Kuhn-Tucker optimality conditions:
- =========================================
- Solver used for this solution: simplex
- 1) Error for Primal Equality Constraints:
- ----------------------------------------
- Largest absolute error: 0.000000 (row id: 0)
- Largest relative error: 0.000000 (row id: 0)
- 2) Error for Primal Inequality Constraints:
- -------------------------------------------
- Largest absolute error: 0.000000 (row id: 0)
- Largest relative error: 0.000000 (row id: 0)
- 1) Error for Dual Equality Constraints:
- ----------------------------------------
- Largest absolute error: 0.000000 (var id: 0)
- Largest relative error: 0.000000 (var id: 0)
- 2) Error for Dual Inequality Constraints:
- -------------------------------------------
- Largest absolute error: 0.000000 (var id: 0)
- Largest relative error: 0.000000 (var id: 0)
- __del__ is deleting problem: basic
复制代码 |
|