More Advanced Examples

Here some more advanced examples are presented with brief explanations. There are dozens more examples in the models folder of the project if you downloaded the zipped source files from its sourceforge project page. The source code of all the examples here are placed in the docs folder of this project. Some of these examples were originally written in GNU MathProg by Andrew Makhorin (thanks to Andrew), and were adopted into PyMathProg.

  1. Re-optimize the assignment
  2. Interaction with queens
  3. The traveling sales man(TSP)
  4. Zebra: using string as index
  5. Magic numbers
  6. Sudoku: want some AI?
  7. Super Sudoku
  8. Machine Scheduling
  9. Scheduling for Revenue

Re-optimize the assignment

In this example, we first solve the assignment problem, then we change one of the cost parameter, and re-optimize. See how easily this is done in PyMathProg.

 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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
# Assignment Problem
# Written in pymprog by Yingjie Lan <ylan@umd.edu>

# The assignment problem is one of the fundamental combinatorial
#   optimization problems.

#   In its most general form, the problem is as follows:

#   There are a number of agents and a number of tasks. Any agent can be
#   assigned to perform any task, incurring some cost that may vary
#   depending on the agent-task assignment. It is required to perform all
#   tasks by assigning exactly one agent to each task in such a way that
#   the total cost of the assignment is minimized.
#   (From Wikipedia, the free encyclopedia.)

#problem data
m = 8 # agents
M = range(m) #set of agents
n = 8 # tasks
N = range(n) #set of tasks
c = [ #cost
(13,21,20,12,8,26,22,11),
(12,36,25,41,40,11,4,8),
(35,32,13,36,26,21,13,37),
(34,54,7,8,12,22,11,40),
(21,6,45,18,24,34,12,48),
(42,19,39,15,14,16,28,46),
(16,34,38,3,34,40,22,24),
(26,20,5,17,45,31,37,43)]

from pymprog import *

begin("assign")
#verbose(True) # for model output
A = iprod(M, N) # Descartan product 
x = var('x', A) # assignment decisions
# use parameters for automatic model update
c = par('c', c) # when their values change
minimize(sum(c[i][j]*x[i,j] for i,j in A))
# each agent is assigned to at most one task
for k in M: sum(x[k,j] for j in N)<=1 
# each task must be assigned to somebody
for k in N: sum(x[i,k] for i in M)==1 

def report():
    print("Total Cost = %g"%vobj())
    assign = [(i,j) for i in M for j in N 
                if x[i,j].primal>0.5]
    for i,j in assign:
        print("Agent %d gets Task %d"%(i, j))
    return assign

solve()
assign = report()
i,j = assign[0]
# model will be updated for the value change
c[i][j].value += 10 
print("Set cost c%r to higher value %r"%
          ([i,j],c[i][j].value))

solve()
report()
end()

Below is the output, note how the optimal assignment changed after increasing the cost, and how the second time of solving the problem seems much easier.

GLPK Simplex Optimizer, v4.60
16 rows, 64 columns, 128 non-zeros
      0: obj =   0.000000000e+00 inf =   8.000e+00 (8)
      8: obj =   1.750000000e+02 inf =   0.000e+00 (0)
*    24: obj =   7.600000000e+01 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Total Cost = 76
Agent 0 gets Task 0
Agent 1 gets Task 7
Agent 2 gets Task 6
Agent 3 gets Task 4
Agent 4 gets Task 1
Agent 5 gets Task 5
Agent 6 gets Task 3
Agent 7 gets Task 2
Set cost c[0, 0] to higher value 23
GLPK Simplex Optimizer, v4.60
16 rows, 64 columns, 128 non-zeros
*    25: obj =   7.800000000e+01 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Total Cost = 78
Agent 0 gets Task 7
Agent 1 gets Task 0
Agent 2 gets Task 6
Agent 3 gets Task 4
Agent 4 gets Task 1
Agent 5 gets Task 5
Agent 6 gets Task 3
Agent 7 gets Task 2

Interaction with queens

This example will randomly pick a board size and then randomly place some queens on the board. The rest of the queens are then placed by the LP program to see if all queens can be placed without attaching each other. If that’s impossible, then the random positions are said to be bad. Note that in that LP, no objective function is needed.

 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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
from pymprog import model, iprod, glpk

# The Queens Problem is to place as many queens as possible on the nxn
#   chess board in a way that they do not fight
#   each other. This problem is probably as old as the chess game itself,
#   and thus its origin is not known, but it is known that Gauss studied
#   this problem.

def queens(n): # n: size of the chess board
   p = model('queens')
   iboard = iprod(range(n), range(n)) #create indices
   x = p.var('x', iboard, bool) #create variables
   sum(x[t] for t in iboard) == n
   for i in range(n): # row-wise
       sum(x[i,j] for j in range(n)) <= 1
   for j in range(n): # column-wise
       sum(x[i,j] for i in range(n)) <= 1
   for k in range(2-n, n-1): # diagonal '\' wise
       sum(x[i,j] for i,j in iboard if i-j == k) <= 1 
   for k in range(1, n+n-2): # anti-diagonal '/' wise
       sum(x[i,j] for i,j in iboard if i+j == k) <= 1 
   return p,x

import random
n = random.randint(6, 11)
print("Board size: %i X %i"%(n,n))
def randpair():
    m = random.randint(0, n*n-1)
    return m%n, m//n
def randpos(k):
    while True:
       pos = [randpair() for i in range(k)]
       if len(set(i for i,j in pos))<k: continue
       if len(set(j for i,j in pos))<k: continue
       if len(set(i+j for i,j in pos))<k: continue
       if len(set(i-j for i,j in pos))<k: continue
       return pos

p, x = queens(n)

def try_out(pos):
   for r,c in pos:
     x[r,c] == 1
   p.solve()
   for r,c in pos:
     x[r,c].reset(0,1)
   return p.status(int) != glpk.GLP_OPT

p.solver(float, msg_lev=glpk.GLP_MSG_OFF)
p.solver(int, msg_lev=glpk.GLP_MSG_OFF)

bads = 0
k = 2
print("randomly put %i queens."%k)
for ii in range(100):
    pos = randpos(k)
    #print(pos)
    bads += try_out(pos)
print('found %i bad positions out of 100'%bads)

The traveling sales man(TSP)

Really, the too well known problem.

  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
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
# TSP, Traveling Salesman Problem
# Model in pymprog by Yingjie Lan <ylan@pku.edu.cn>

# The Traveling Salesman Problem (TSP) is a famous problem.
#   Let a directed graph G = (V, E) be given, where V = {1, ..., n} is
#   a set of nodes, E <subset of> V x V is a set of arcs. Let also each arc
#   e = (i,j) be assigned a number c[i,j], which is the length of the
#   arc e. The problem is to find a closed path of minimal length going
#   through each node of G exactly once. 
from __future__ import print_function

n = 16 #number of nodes
V = range(1,n+1) #set of notes
#cost or each arc, format: (start, end):cost
c = {(1,2):509, (1,3):501, (1,4):312, (1,5):1019, (1,6):736, (1,7):656, 
     (1,8): 60, (1,9):1039, (1,10):726, (1,11):2314, (1,12):479, 
     (1,13):448, (1,14):479, (1,15):619, (1,16):150, 
(2,1):509, (2,3):126, (2,4):474, (2,5):1526, (2,6):1226, (2,7):1133, 
     (2,8):532, (2,9):1449, (2,10):1122, (2,11):2789, (2,12):958, 
     (2,13):941, (2,14):978, (2,15):1127, (2,16):542, 
(3,1):501, (3,2):126, (3,4):541, (3,5):1516, (3,6):1184, (3,7):1084, 
     (3,8):536, (3,9):1371, (3,10):1045, (3,11):2728, (3,12):913, 
     (3,13):904, (3,14):946, (3,15):1115, (3,16):499, 
(4,1):312, (4,2):474, (4,3):541, (4,5):1157, (4,6):980, (4,7):919, 
     (4,8):271, (4,9):1333, (4,10):1029, (4,11):2553, (4,12):751, 
     (4,13):704, (4,14):720, (4,15):783, (4,16):455, 
(5,1):1019, (5,2):1526, (5,3):1516, (5,4):1157, (5,6):478, (5,7):583, 
     (5,8):996, (5,9):858, (5,10):855, (5,11):1504, (5,12):677, 
     (5,13):651, (5,14):600, (5,15):401, (5,16):1033, 
(6,1):736, (6,2):1226, (6,3):1184, (6,4):980, (6,5):478, (6,7):115, 
     (6,8):740, (6,9):470, (6,10):379, (6,11):1581, (6,12):271, 
     (6,13):289, (6,14):261, (6,15):308, (6,16):687, 
(7,1):656, (7,2):1133, (7,3):1084, (7,4):919, (7,5):583, (7,6):115, 
     (7,8):667, (7,9):455, (7,10):288, (7,11):1661, (7,12):177, 
     (7,13):216, (7,14):207, (7,15):343, (7,16):592, 
(8,1): 60, (8,2):532, (8,3):536, (8,4):271, (8,5):996, (8,6):740, 
     (8,7):667, (8,9):1066, (8,10):759, (8,11):2320, (8,12):493, 
     (8,13):454, (8,14):479, (8,15):598, (8,16):206, 
(9,1):1039, (9,2):1449, (9,3):1371, (9,4):1333, (9,5):858, (9,6):470, 
     (9,7):455, (9,8):1066, (9,10):328, (9,11):1387, (9,12):591, 
     (9,13):650, (9,14):656, (9,15):776, (9,16):933, 
(10,1):726, (10,2):1122, (10,3):1045, (10,4):1029, (10,5):855, 
     (10,6):379, (10,7):288, (10,8):759, (10,9):328, (10,11):1697, 
     (10,12):333, (10,13):400, (10,14):427, (10,15):622, (10,16):610, 
(11,1):2314, (11,2):2789, (11,3):2728, (11,4):2553, (11,5):1504, 
     (11,6):1581, (11,7):1661, (11,8):2320, (11,9):1387, (11,10):1697, 
     (11,12):1838, (11,13):1868, (11,14):1841, (11,15):1789, (11,16):2248, 
(12,1):479, (12,2):958, (12,3):913, (12,4):751, (12,5):677, (12,6):271, 
     (12,7):177, (12,8):493, (12,9):591, (12,10):333, (12,11):1838, 
     (12,13): 68, (12,14):105, (12,15):336, (12,16):417, 
(13,1):448, (13,2):941, (13,3):904, (13,4):704, (13,5):651, (13,6):289, 
     (13,7):216, (13,8):454, (13,9):650, (13,10):400, (13,11):1868, 
     (13,12): 68, (13,14): 52, (13,15):287, (13,16):406, 
(14,1):479, (14,2):978, (14,3):946, (14,4):720, (14,5):600, (14,6):261, 
     (14,7):207, (14,8):479, (14,9):656, (14,10):427, (14,11):1841, 
     (14,12):105, (14,13): 52, (14,15):237, (14,16):449, 
(15,1):619, (15,2):1127, (15,3):1115, (15,4):783, (15,5):401, (15,6):308, 
     (15,7):343, (15,8):598, (15,9):776, (15,10):622, (15,11):1789, 
     (15,12):336, (15,13):287, (15,14):237, (15,16):636, 
(16,1):150, (16,2):542, (16,3):499, (16,4):455, (16,5):1033, (16,6):687, 
     (16,7):592, (16,8):206, (16,9):933, (16,10):610, (16,11):2248, 
     (16,12):417, (16,13):406, (16,14):449, (16,15):636}
#set of arcs: (i,j) repr an arc from i to j
E = c.keys()

from pymprog import model
p = model("tsp")
x = p.var('x', E, bool) # created over E.
#minize the total travel distance
p.min(sum(c[t]*x[t] for t in E), 'totaldist')
#subject to: leave each city exactly once
for k in V: sum(x[k,j] for j in V if (k,j) in E)==1 
#subject to: enter each city exactly once
for k in V: sum(x[i,k] for i in V if (i,k) in E)==1 

#some flow constraints to eliminate subtours.
#y: the number of cars carried: city 1 has n cars.
#exactly one car will be distributed to each city.
y=p.var('y', E) 
for t in E: (n-1)*x[t] >= y[t] 
for k in V: (
  sum(y[k,j] for j in V if (k,j) in E) # cars out
  - sum(y[i,k] for i in V if (i,k) in E) # cars in
  ==  (n if k==1 else 0) - 1 )

p.solve(float) #solve as LP only.
print("simplex done: %r"% p.status())
p.solve(int) #solve the IP problem
print(p.vobj())
tour = [t for t in E if x[t].primal>.5]
cat, car = 1, n
print("This is the optimal tour with [cars carried]:")
for k in V: 
   print(cat, end='')
   for i,j in tour: 
      if i==cat: 
         print("[%g]"%y[i,j].primal, end='')
         cat=j
         break
print(cat)
p.end()

Zebra: using string as index

A logic puzzle, see how strings are used as indices here to enhance the readability of the model.

  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
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
# ZEBRA, Who Owns the Zebra? 
# Adapted from GLPK to pymprog by Yingjie Lan <ylan@pku.edu.cn>

########################################################################
#  The Zebra Puzzle is a well-known logic puzzle.
#
#  It is often called "Einstein's Puzzle" or "Einstein's Riddle"
#  because it is said to have been invented by Albert Einstein as a boy,
#  with the common claim that Einstein said "only 2 percent of the
#  world's population can solve". It is also sometimes attributed to
#  Lewis Carroll. However, there is no known evidence for Einstein's or
#  Carroll's authorship.
#
#  There are several versions of this puzzle. The version below is
#  quoted from the first known publication in Life International
#  magazine on December 17, 1962.
#
#   1. There are five houses.
#   2. The Englishman lives in the red house.
#   3. The Spaniard owns the dog.
#   4. Coffee is drunk in the green house.
#   5. The Ukrainian drinks tea.
#   6. The green house is immediately to the right of the ivory house.
#   7. The Old Gold smoker owns snails.
#   8. Kools are smoked in the yellow house.
#   9. Milk is drunk in the middle house.
#  10. The Norwegian lives in the first house.
#  11. The man who smokes Chesterfields lives in the house next to the
#      man with the fox.
#  12. Kools are smoked in the house next to the house where the horse
#      is kept.
#  13. The Lucky Strike smoker drinks orange juice.
#  14. The Japanese smokes Parliaments.
#  15. The Norwegian lives next to the blue house.
#
#  Now, who drinks water? Who owns the zebra?
#
#  In the interest of clarity, it must be added that each of the five
#  houses is painted a different color, and their inhabitants are of
#  different national extractions, own different pets, drink different
#  beverages and smoke different brands of American cigarettes. One
#  other thing: In statement 6, right means your right.
#
#  (From Wikipedia, the free encyclopedia.)
########################################################################

HOUSE=range(1,6)

COLOR = ("blue", "green", "ivory", "red", "yellow")

NATIONALITY = ("Englishman", "Japanese", "Norwegian",
               "Spaniard", "Ukranian")

DRINK = ("coffee", "milk", "orange_juice", "tea", "water")

SMOKE = ("Chesterfield", "Kools", "Lucky_Strike", 
          "Old_Gold", "Parliament")

PET = ("dog", "fox", "horse", "snails", "zebra")

from pymprog import *
#Uncomment to forbid comparison with a variable as an expression.
#pymprog.variable.comparable=False
zeb = model("zebra")

color=zeb.var('color', iprod(HOUSE, COLOR), bool)
for h in HOUSE: sum(color[h,c] for c in COLOR)==1
for c in COLOR: sum(color[h,c] for h in HOUSE)==1

nationality=zeb.var('nation', iprod(HOUSE, NATIONALITY), bool)
for h in HOUSE: sum(nationality[h,n] for n in NATIONALITY)==1
for n in NATIONALITY: sum(nationality[h,n] for h in HOUSE)==1 

drink=zeb.var('drink', iprod(HOUSE, DRINK), bool)
for h in HOUSE: sum(drink[h,d] for d in DRINK)==1 
for d in DRINK: sum(drink[h,d] for h in HOUSE)==1 

smoke=zeb.var('smoke',iprod(HOUSE, SMOKE),bool)
for h in HOUSE: sum(smoke[h,s] for s in SMOKE)==1 
for s in SMOKE: sum(smoke[h,s] for h in HOUSE)==1 

pet=zeb.var('pet', iprod(HOUSE, PET),bool)
for h in HOUSE: sum(pet[h,p] for p in PET)==1 
for p in PET: sum(pet[h,p] for h in HOUSE)==1 

## the Englishman lives in the red house */
for h in HOUSE: nationality[h,"Englishman"]==color[h,"red"] 

## the Spaniard owns the dog */
for h in HOUSE: nationality[h,"Spaniard"]==pet[h,"dog"] 

## coffee is drunk in the green house */
for h in HOUSE: drink[h,"coffee"]==color[h,"green"] 

## the Ukrainian drinks tea */
for h in HOUSE: nationality[h,"Ukranian"]==drink[h,"tea"] 

## the green house is immediately to the right of the ivory house */
for h in HOUSE:
   color[h,"green"] == (color[h-1, "ivory"] if h>1 else 0)
## the Old Gold smoker owns snails */
for  h in HOUSE: smoke[h,"Old_Gold"]==pet[h,"snails"] 

## Kools are smoked in the yellow house */
for h in HOUSE: smoke[h,"Kools"]==color[h,"yellow"] 

## milk is drunk in the middle house */
drink[3,"milk"] == 1

## the Norwegian lives in the first house */
nationality[1,"Norwegian"] == 1

# the man who smokes Chesterfields lives in the house 
#   next to the man with the fox */
for h in HOUSE:(
   1 - smoke[h,"Chesterfield"] +
   (0 if h == 1 else pet[h-1,"fox"]) +
   (0 if h == 5 else pet[h+1,"fox"]) >= 1 )

# Kools are smoked in the house next to the house 
#   where the horse is kept */
for h in HOUSE:(
   1 - smoke[h,"Kools"] +
   (0 if h == 1 else pet[h-1,"horse"]) +
   (0 if h == 5 else pet[h+1,"horse"]) >= 1)

## the Lucky Strike smoker drinks orange juice */
for h in HOUSE:
   smoke[h,"Lucky_Strike"] == drink[h,"orange_juice"]


## the Japanese smokes Parliaments */
for h in HOUSE:
   nationality[h,"Japanese"] == smoke[h,"Parliament"]


#zeb.verb = True
## the Norwegian lives next to the blue house */
for h in HOUSE:(
   1 - nationality[h,"Norwegian"] +
   (0 if h == 1 else color[h-1,"blue"]) +
   (0 if h == 5 else color[h+1,"blue"]) >= 1)

zeb.solve()

def answ(h):
   for c in COLOR: 
      if color[h,c].primal>0: break
   for n in NATIONALITY:
      if nationality[h,n].primal>0: break
   for d in DRINK:
      if drink[h,d].primal>0: break
   for s in SMOKE:
      if smoke[h,s].primal>0: break
   for p in PET:
      if pet[h,p].primal>0: break
   return (h,c,n,d,s,p)

print("HOUSE  COLOR   NATIONALITY  DRINK         SMOKE         PET")
for h in HOUSE:
   print("%5d  %-6s  %-11s  %-12s  %-12s  %-6s"%answ(h))
zeb.end()

Magic numbers

Magic numbers, a very interesting problem.

 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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
"""In recreational mathematics, a magic square of order n is an
   arrangement of n^2 numbers, usually distinct integers, in a square,
   such that n numbers in all rows, all columns, and both diagonals sum
   to the same constant. A normal magic square contains the integers
   from 1 to n^2.

   (From Wikipedia, the free encyclopedia.)

When n=5, we have:

the magic sum must be:  65
solver status: opt
10  2  4 24 25
 9 11 20  6 19
22 14 16  8  5
23 17  7 15  3
 1 21 18 12 13"""

# square order 
n = 4
# set of numbers
N = range(1,n*n+1)
#the magic sum:
s = sum(t for t in N)//n
print("the magic sum must be: %r "% s)

from pymprog import  model, iprod
S = iprod(range(1,n+1), range(1,n+1))
T = iprod(range(1,n+1), range(1,n+1), N)
p = model('magic')
x = p.var('x', T, bool)
#x[i,j,k] = 1 means that cell (i,j) contains integer k

#each cell must be assigned exactly one integer
for i,j in S: sum(x[i,j,k] for k in N)==1 

#each integer must be assigned exactly to one cell
for k in N: sum(x[i,j,k] for i,j in S)==1 

#the sum in each row must be the magic sum 
for i in range(1, n+1):
  sum(k*x[i,j,k] for j in range(1,n+1) for k in N)==s

#the sum in each column must be the magic sum 
for j in range(1, n+1):
  sum(k*x[i,j,k] for i in range(1,n+1) for k in N)==s

#the sum in the diagonal must be the magic sum
sum(k*x[i,i,k] for i in range(1,n+1) for k in N)==s

#the sum in the co-diagonal must be the magic sum
sum(k*x[i,n+1-i,k] for i in range(1,n+1) for k in N)==s

p.solve()
print("solver status: %s"% p.status())
for i in range(1,n+1):
   print(' '.join("%2g"%sum(x[i,j,k].primal*k for k in N) 
             for j in range(1,n+1)))
p.end()

Sudoku: want some AI?

Do you play Sudoku? This program might help you crack a Sudoku in seconds ;).

 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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
# SUDOKU, Number Placement Puzzle 
from __future__ import print_function
"""Sudoku, also known as Number Place, is a logic-based placement
   puzzle. The aim of the canonical puzzle is to enter a numerical
   digit from 1 through 9 in each cell of a 9x9 grid made up of 3x3
   subgrids (called "regions"), starting with various digits given in
   some cells (the "givens"). Each row, column, and region must contain
   only one instance of each numeral.

   Example:

   +-------+-------+-------+
   | 5 3 . | . 7 . | . . . |
   | 6 . . | 1 9 5 | . . . |
   | . 9 8 | . . . | . 6 . |
   +-------+-------+-------+
   | 8 . . | . 6 . | . . 3 |
   | 4 . . | 8 . 3 | . . 1 |
   | 7 . . | . 2 . | . . 6 |
   +-------+-------+-------+
   | . 6 . | . . . | 2 8 . |
   | . . . | 4 1 9 | . . 5 |
   | . . . | . 8 . | . 7 9 |
   +-------+-------+-------+

   (From Wikipedia, the free encyclopedia.)"""

# the "givens" 
g =(
(5,3,0,0,7,0,0,0,0),
(6,0,0,1,9,5,0,0,0),
(0,9,8,0,0,0,0,6,0),
(8,0,0,0,6,0,0,0,3),
(4,0,0,8,0,3,0,0,1),
(7,0,0,0,2,0,0,0,6),
(0,6,0,0,0,0,2,8,0),
(0,0,0,4,1,9,0,0,5),
(0,0,0,0,8,0,0,7,9))

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

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

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

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

#cells in \-diagonal
#p.st([sum(x[i,i,k] for i in I)==1 for k in K])

#cells in /-diagonal
#p.st([sum(x[i,10-i,k] for i in I)==1 for k in K])

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

p.solve()

for i in I:
   if i in range(1,10,3):
      print(" +-------+-------+-------+")
   print('', end=' ')
   for j in J:
      if j in range(1,10,3): print("|", end=' ')
      print("%g"%sum(x[i,j,k].primal*k for k in K), end=' ')
      if j==9: print("|")
   if i == 9:
      print(" +-------+-------+-------+")
p.end()

Super Sudoku

This example will provide a sample Super Sudoku: in addition to satisfying all the requirements of Sudoku, Super Sudoku also requires that the elements in each diagonal must be distinct, below is a sample Super Sudoku obtained by the code provided afterwards:

+-------+-------+-------+
| 2 1 3 | 4 6 5 | 8 7 9 |
| 4 5 6 | 7 8 9 | 2 1 3 |
| 7 8 9 | 2 1 3 | 4 5 6 |
+-------+-------+-------+
| 1 2 4 | 3 5 6 | 7 9 8 |
| 3 6 8 | 9 7 2 | 5 4 1 |
| 9 7 5 | 8 4 1 | 3 6 2 |
+-------+-------+-------+
| 8 4 2 | 1 9 7 | 6 3 5 |
| 6 3 1 | 5 2 4 | 9 8 7 |
| 5 9 7 | 6 3 8 | 1 2 4 |
+-------+-------+-------+

 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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
from __future__ import print_function
"""Sudoku, also known as Number Place, is a logic-based placement
   puzzle. The aim of the canonical puzzle is to enter a numerical
   digit from 1 through 9 in each cell of a 9x9 grid made up of 3x3
   subgrids (called "regions"), starting with various digits given in
   some cells (the "givens"). Each row, column, and region must contain
   only one instance of each numeral.

   (From Wikipedia, the free encyclopedia.)

   This example will provide a sample Super Sudoku:
   in addition to satisfying all the requirements 
   of Sudoku, Super Sudoku also requires that the
   elements in each diagonal must be distinct."""

from pymprog import *        # Import the module
begin("ssudoku")
I = range(1,10)
J = range(1,10)
K = range(1,10)
T = iprod(I,J,K) #create Indice tuples
#x[i,j,k] = 1 means cell [i,j] is assigned number k 
x = var('x', T, bool) #binary vars
#each cell must be assigned exactly one number
[sum(x[i,j,k] for k in K)==1 for i in I for j in J]
#cells in the same row must be assigned distinct numbers
[sum(x[i,j,k] for j in J)==1 for i in I for k in K]
#cells in the same column must be assigned distinct numbers
[sum(x[i,j,k] for i in I)==1 for j in J for k in K]
#cells in the same region must be assigned distinct numbers
[sum(x[i,j,k] for i in range(r,r+3) for j in range(c, c+3))==1
 for r in range(1,10,3) for c in range(1,10,3) for k in K]

#cells in \-diagonal
[sum(x[i,i,k] for i in I)==1 for k in K]

#cells in /-diagonal
[sum(x[i,10-i,k] for i in I)==1 for k in K]

#there is no need for an objective function here

solve()

for i in I:
   if i in range(1,10,3):
      print(" +-------+-------+-------+")
   print('', end=' ')
   for j in J:
      if j in range(1,10,3): print("|", end=' ')
      print("%g"%sum(x[i,j,k].primal*k for k in K), end=' ')
      if j==9: print("|")
   if i == 9:
      print(" +-------+-------+-------+")

end()

Machine Scheduling

Schedule jobs on a single machine, given the duration, earliest start time, and the latest finish time of each job. Jobs can not be interrupted, and the machine can only handle one job at one time.

 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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
n = 3
N = range(n)
M = [(i,j) for i in N for j in N if i<j]

D = (3,4,2) #duration of each job
L = (0,2,0) #earliest start
U = (9,7,8) #latest finish

from pymprog import *

begin("job-scheduling")
x = var('x',  N) #start time
#MD[i,j] = (D[i]+D[j])/2.0
#T[i] = x[i] + D[i]
#y[i,j]<= |T[i]-x[j]-MD[i,j]|
#y[i,j] < MD[i,j] <==> overlap betw jobs i,j
y = var('y',  M ) 
#w[i,j]: the 'OR' for |T[i]-x[j]-MD[i,j]|
w = var('w', M, kind=bool)
# z[i,j] >= MD[i,j] - y[i,j]
z = var('z',  M )

minimize( sum(z[i,j] for i,j in M) )

for i,j in M:
   ((D[i]+D[j])/2.0 - (x[i]+D[i] - x[j]) + 
       (U[i]-L[j]) * w[i,j] >= y[i,j])

   ((x[i]+D[i] - x[j]) - (D[i]+D[j])/2.0 +
       (U[j]-L[i])*(1-w[i,j]) >= y[i,j])

   (D[i]+D[j])/2.0 - y[i,j] <= z[i,j] 

#set bounds on x
for i in N:
   L[i] <= x[i] <= U[i] - D[i] 

#another way to enforce no overlapping:
#
# x[i] >= T[j] or x[j] >= T[i]
#
# Which can be formulated as:
#
# x[i] + (U[j]-L[i])*w[i,j]>= x[j]+D[j]
# x[j] + (U[i]-L[j])*(1-w[i,j]) >= x[i]+D[i]

solve()

print("status: %r"% status())
print( "overlap: %r"% vobj())
print( "schedule:")

for i in N:
   start = x[i].primal
   print("job %i: %r, %r"%(i, start, start+D[i]))

The output schedule is as follows (your run may not be the same, as the data is randomly generated – you might even get infeasible problem instances):

Duration [7, 4, 1, 7]
Earliest [13, 17, 19, 12]
Latest [22, 41, 28, 38]
status: opt
schedule:
job 0: 13.0 20.0
job 1: 20.0 24.0
job 2: 27.0 28.0
job 3: 31.0 38.0

Scheduling for Revenue

The same scheduling problem as above, but add a new spin to the problem: when it is not feasible to accept all jobs, try accept those that will maximize the total revenue.

 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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
from __future__ import print_function
n = 20
N = range(n)

from random import Random

rand = Random()

D = [rand.randint(1,10) for i in N] # processing time
R = [rand.randint(1,10) for i in N] # revenue
tt = sum(D); tt -= tt//5
L = [rand.randint(1,tt) for i in N] # arrival time
U = [L[i]+D[i]+rand.randint(10,30) for i in N] # duedate

#overlapping jobs
M = [(i,j) for i in N for j in N 
   if i<j and L[i]<U[j] and L[j]<U[i]] 

#i must preceed j and they could bump into each other
P = [(i,j) for i in N for j in N if i!=j
   and L[i]<U[j] and L[j]<U[i]
   #lastest start < earliest finish
   and U[i]-D[i]<L[j]+D[j]]


from pymprog import *

begin("job-scheduling")
x = var('x', N ) #start time
#MD[i,j] = (D[i]+D[j])/2.0
#T[i] = x[i] + D[i]
#y[i,j]<= |T[i]-x[j]-MD[i,j]|
#y[i,j] < MD[i,j] <==> overlap betw jobs i,j
y = var('y', M) 
#w[i,j]: the 'OR' for |T[i]-x[j]-MD[i,j]|
w = var('w', M, bool)
# z[i,j] >= MD[i,j] - y[i,j]
z = var('z', M )
#u[i] = 1 iff job i is scheduled.
u = var('u', N, bool)

maximize(sum(R[i]*u[i] for i in N), 'revenue')

sum(z[i,j] for i,j in M) == 0 #single
#w[i,j]=0  ==> 
#x[j] + D[j]/2 >= x[i]+D[i]/2 + y[i,j]
#That is: i preceed j when w[i,j]=0.
for i,j in M:
   ((D[i]+D[j])/2.0 - (x[i]+D[i]- x[j]) + 
        (U[i]-L[j]) * w[i,j] >= y[i,j])

   ((x[i]+D[i]- x[j]) - (D[i]+D[j])/2.0 +
        (U[j]-L[i])*(1-w[i,j]) >= y[i,j])

   ((D[i]+D[j])/2.0 - y[i,j] <= z[i,j] + 
       (2-u[i]-u[j])*(D[i]+D[j])/2.0)

   if (i,j) in P or (j,i) in P:
       w[i,j] == (0 if (i,j) in P else 1)

#set bounds on x
for i in N:
   L[i] <= x[i] <= U[i] - D[i] 


solver(int, 
    #this branching option often helps 
    br_tech=glpk.GLP_BR_PCH, 
)
solve()

print("status:", status())
print("revenue:", vobj(), 'max:', sum(R))
print("schedule:")
for i in N:
   start = x[i].primal
   used = u[i].primal
   print("job %i:"%i, start, start+D[i],
          "Accept" if used else "Reject",
          R[i]/float(D[i]))