A Dive-in Tutorial

A dive-in example

Suppose you have an easy linear program. A 3x3 matrix is given as the matrix for constraints, along with 2 vectors of length 3, one for the objective function and the other for the right-hand-sides. You need to solve the problem and give a sensitivity report. This is very easy with PyMathProg:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
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("###>Objective value: %f"%vobj())
sensitivity() # sensitivity report
end() #Good habit: do away with the model

The code is quite straight-forward. Note there are three basic building blocks: it starts with a data block, followed by a model block, and finished with a report block.

  1. Data block (lines 2-6):

    The matrix and the vectors are defined. Of course, you can also read data from some external database, as Python can do that easily.

  2. Model block (lines 7-12):

    It all starts by a function begin(name) on line 7. It creates a new model instance with the given name for later modelling steps to build upon. It is also possible to handle the model directly, we will come to that later.

    Line 8 turns on the verbosity option, which enables PyMathProg to provide feedbacks on each building step.

    Line 9 defines the three variables and organize them in a list: you simply provide the name for the group, and the number of variables to create, which is 3 in this case. By default, these variables are continuous and non-negative.

    Line 10 defines the objective: maximize the summation of the terms b[i]*x[i], where i goes from 0 to 3.

    Lines 11-12 define the constraints with a for loop.

    That’s it for modelling. Now the code moves on to

  3. Report block (lines 13-16):

    Line 13 issues a call solve() to solve the model.

    Line 14 prints the objective value.

    Line 15 produces the sensitivity report.

    Finally, line 16 calls the function end() to do away with the model.

And here is the output (with the verbosity turned on):

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
###>Objective value: 76.000000

PyMathProg 1.0 Sensitivity Report Created: 2016/12/08 Thu 12:55PM
================================================================================
Variable            Activity   Dual.Value     Obj.Coef   Range.From   Range.Till
--------------------------------------------------------------------------------
*x[0]                      4            0           10            6         13.5
*x[1]                      6            0            6      4.44444           10
 x[2]                      0         -2.8            4         -inf          6.8
================================================================================
Note: rows marked with a * list a basic variable.

================================================================================
Constraint       Activity Dual.Value  Lower.Bnd  Upper.Bnd RangeLower RangeUpper
--------------------------------------------------------------------------------
 R1                    10        2.8       -inf         10    6.66667         15
 R2                    60        0.8       -inf         60         40         90
*R3                    20          0       -inf         30         20         20
================================================================================
Note: normally, RangeLower is the min for the binding bound, and RangeUpper
gives the max value. However, when neither bounds are binding, the row is
marked with a *, and RangeLower is the max for Lower.Bnd(whose min is -inf),
and RangeUpper is the min for Upper.Bnd(whose max value is inf). Then the
columns of RangeLower, RangeUpper and Activity all have identical values.

__del__ is deleting problem: basic

A parallel example

In the previous example, we use global functions to create and manipulate a default model – we did not directly handle it. PyMathProg also provides a way of direct handling of the model, as demonstrated by the code below that does exactly the same thing as before:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
from pymprog import model 
c = (10, 6, 4)
A = [ ( 1, 1, 1),     
      ( 9, 4, 5),   
      ( 2, 2, 6) ]   
b = (10, 60, 30)
p = model('basic')  
p.verbose(True)
x = p.var('x', 3) 
p.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] 
p.solve() # solve the model
print("###>Objective value: %f"%p.vobj())
p.sensitivity() # sensitivity report
p.end()

Clearly, there is a line to line correspondence to the previous code. In the first line, this code only imports model from pymprog. Then on line 7, instead of call begin(.), this code uses p=model(.) to create a new model instance. From then on, each function call is converted into a method invocation on the created model instance.

Zero-sum Two-player Game

This example is employed to show the use of the ‘primal’ and ‘dual’ values. A zero-sum two-player game is a game between two players, where the gain of one player is the loss of the other, hence their pay-offs always sums up to zero. The value a[i,j] is the pay-off for player one when player one plays strategy i and player two plays strategy j. Here is an LP formulation to find the equilibrium mixed strategies and the value of the game.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
from pymprog import *

#####Solve this 2-player zero-sum game:
##
##     Gain for player 1 
##    (Loss for player 2)
##   
##            ||  Player  2
##   Player 1 ||  B1     B2
##      A1    ||  5      9
##      A2    ||  8      6
##
##############################

begin('game')
# gain of player 1, a free variable
v = var('game_value', bounds=(None,None))
# mixed strategy of player 2
p = var('p', 2) 
# probability sums to 1
sum(p) == 1
# player 2 chooses p to minimize v
minimize(v) 
# player 1 chooses the better value 
r1 =  v >= 5*p[0] + 9*p[1] 
r2 =  v >= 8*p[0] + 6*p[1]
solve()
print('Game value: %g'% v.primal)
print("Mixed Strategy for player 1:")
print("A1: %g, A2: %g"%(r1.dual, r2.dual))
print("Mixed Strategy for player 2:")
print("B1: %g, B2: %g"%(p[0].primal, p[1].primal))
end()

In this block of code, two variables r1 and r2 are employed to save the constraints for the sake of reporting. Note that in this model, the primal value of the variables gives the probability for player 2’s mixed strategy, and the dual value of the constraints r1 and r2 gives the mixed strategy of player 1. And the output is as follows:

GLPK Simplex Optimizer, v4.60
3 rows, 3 columns, 8 non-zeros
    0: obj =   0.000000000e+00 inf =   1.000e+00 (1)
    3: obj =   7.000000000e+00 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Game value: 7
Mixed Strategy for player 1:
A1: 0.333333, A2: 0.666667
Mixed Strategy for player 2:
B1: 0.5, B2: 0.5

Hope this dive-in tutorial has served to give you the basic idea of how PyMathProg works. Before working on more advanced examples, please read on!