# 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.

- Re-optimize the assignment
- Interaction with queens
- The traveling sales man(TSP)
- Zebra: using string as index
- Magic numbers
- Sudoku: want some AI?
- Super Sudoku
- Machine Scheduling
- 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]))
``` |